123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249 |
- /*
- Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
- Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program.
- Use, modification and distribution is subject to the Boost Software License,
- Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
- http://www.boost.org/LICENSE_1_0.txt)
- This file is converted from GeographicLib, https://geographiclib.sourceforge.io
- GeographicLib is originally written by Charles Karney.
- Author: Charles Karney (2008-2017)
- Last updated version of GeographicLib: 1.49
- Original copyright notice:
- Copyright (c) Charles Karney (2009-2015) <charles@karney.com> and
- licensed under the MIT/X11 License. For more information, see
- https://geographiclib.sourceforge.io
- Compute the series expansions for the ellipsoidal geodesic problem.
- References:
- Charles F. F. Karney,
- Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
- https://doi.org/10.1007/s00190-012-0578-z
- Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
- The code below contains minor modifications to conform with
- Boost Geometry style guidelines.
- To run the code, start Maxima and enter
- load("geod.mac")$
- */
- taylordepth:5$
- ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$
- jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
- ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
- computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
- sintegrand:sqrt(1+k2*sin(sigma)^2),
- sintegrandexp:ataylor(
- (1-eps)*subst([k2=4*eps/(1-eps)^2],sintegrand),
- eps,maxpow),
- s:trigreduce(integrate(sintegrandexp,sigma)),
- s:s-subst(sigma=0,s),
- A1:expand(subst(sigma=2*%pi,s)/(2*%pi)),
- tau1:ataylor(s/A1,eps,maxpow),
- for i:1 thru maxpow do C1[i]:coeff(tau1,sin(2*i*sigma)),
- if expand(tau1-sigma-sum(C1[i]*sin(2*i*sigma),i,1,maxpow)) # 0
- then error("left over terms in B1"),
- A1:A1/(1-eps),
- 'done)$
- codeA1(maxpow):=block([tab2:" ",tab3:" "],
- print("// The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
- static inline CT evaluate_series_A1(CT eps) {
- CT eps2 = math::sqr(eps);
- CT t;
- switch (SeriesOrder/2) {"),
- for n:0 thru entier(maxpow/2) do block([
- q:horner(ataylor(subst([eps=sqrt(eps2)],A1*(1-eps)-1),eps2,n)),
- linel:1200],
- print(concat(tab2,"case ",string(n),":")),
- print(concat(tab3,"t = ",string(q),";")),
- print(concat(tab3,"break;"))),
- print(" }
- return (t + eps) / (1 - eps);
- }"),
- 'done)$
- computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
- sintegrand:1/sqrt(1+k2*sin(sigma)^2),
- sintegrandexp:ataylor(
- (1+eps)*subst([k2=4*eps/(1-eps)^2],sintegrand),
- eps,maxpow),
- s:trigreduce(integrate(sintegrandexp,sigma)),
- s:s-subst(sigma=0,s),
- A2:expand(subst(sigma=2*%pi,s)/(2*%pi)),
- tau1:ataylor(s/A2,eps,maxpow),
- for i:1 thru maxpow do C2[i]:coeff(tau1,sin(2*i*sigma)),
- if expand(tau1-sigma-sum(C2[i]*sin(2*i*sigma),i,1,maxpow)) # 0
- then error("left over terms in B2"),
- A2:A2/(1+eps),
- 'done)$
- codeA2(maxpow):=block([tab2:" ",tab3:" "],
- print("// The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
- CT evaluate_series_A2(CT const& eps)
- {
- CT const eps2 = math::sqr(eps);
- CT t;
- switch (SeriesOrder/2) {"),
- for n:0 thru entier(maxpow/2) do block([
- q:horner(ataylor(subst([eps=sqrt(eps2)],A2*(1+eps)-1),eps2,n)),
- linel:1200],
- print(concat(tab2,"case ",string(n),":")),
- print(concat(tab3,"t = ",string(q),";")),
- print(concat(tab3,"break;"))),
- print(" }
- return (t - eps) / (1 + eps);
- }"),
- 'done)$
- computeI3(maxpow):=block([int,intexp,dlam,eta,del,eps,nu,f,z,n],
- maxpow:maxpow-1,
- int:subst([k2=4*eps/(1-eps)^2],
- (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma)^2))),
- int:subst([f=2*n/(1+n)],int),
- intexp:jtaylor(int,n,eps,maxpow),
- dlam:trigreduce(integrate(intexp,sigma)),
- dlam:dlam-subst(sigma=0,dlam),
- A3:expand(subst(sigma=2*%pi,dlam)/(2*%pi)),
- eta:jtaylor(dlam/A3,n,eps,maxpow),
- A3:jtaylor(A3,n,eps,maxpow),
- for i:1 thru maxpow do C3[i]:coeff(eta,sin(2*i*sigma)),
- if expand(eta-sigma-sum(C3[i]*sin(2*i*sigma),i,1,maxpow)) # 0
- then error("left over terms in B3"),
- 'done)$
- codeA3(maxpow):=block([tab2:" ",tab3:" "],
- print("// The scale factor A3 = mean value of (d/dsigma)I3
- static inline void evaluate_series_A3(CT const& n, CT c[])
- {
- switch (SeriesOrder) {"),
- for nn:0 thru maxpow do block(
- [q:if nn=0 then 0 else
- jtaylor(subst([n=n],A3),n,eps,nn-1),
- linel:1200],
- print(concat(tab2,"case ",string(nn),":")),
- for i : 0 thru nn-1 do
- print(concat(tab3,"c[",i,"] = ",
- string(horner(coeff(q,eps,i))),";")),
- print(concat(tab3,"break;"))),
- print(" }
- }"),
- 'done)$
- codeC1(maxpow):=block([tab2:" ",tab3:" "],
- print("// The coefficients C1[l] in the Fourier expansion of B1
- static inline evaluate_coeffs_C1(CT eps, CT c[]) {
- CT eps2 = math::sqr(eps);
- CT d = eps;
- switch (SeriesOrder) {"),
- for n:0 thru maxpow do (
- print(concat(tab2,"case ",string(n),":")),
- for m:1 thru n do block([q:d*horner(
- subst([eps=sqrt(eps2)],ataylor(C1[m],eps,n)/eps^m)),
- linel:1200],
- if m>1 then print(concat(tab3,"d *= eps;")),
- print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
- print(concat(tab3,"break;"))),
- print(" }
- }"),
- 'done)$
- revertI1(maxpow):=block([tau,eps,tauacc:1,sigacc:0],
- for n:1 thru maxpow do (
- tauacc:trigreduce(ataylor(
- -sum(C1[j]*sin(2*j*tau),j,1,maxpow-n+1)*tauacc/n,
- eps,maxpow)),
- sigacc:sigacc+expand(diff(tauacc,tau,n-1))),
- for i:1 thru maxpow do C1p[i]:coeff(sigacc,sin(2*i*tau)),
- if expand(sigacc-sum(C1p[i]*sin(2*i*tau),i,1,maxpow)) # 0
- then error("left over terms in B1p"),
- 'done)$
- codeC1p(maxpow):=block([tab2:" ",tab3:" "],
- print("// The coefficients C1p[l] in the Fourier expansion of B1p
- static inline evaluate_coeffs_C1p(CT eps, CT c[])
- {
- CT const eps2 = math::sqr(eps);
- CT d = eps;
- switch (SeriesOrder) {"),
- for n:0 thru maxpow do (
- print(concat(tab2,"case ",string(n),":")),
- for m:1 thru n do block([q:d*horner(
- subst([eps=sqrt(eps2)],ataylor(C1p[m],eps,n)/eps^m)),
- linel:1200],
- if m>1 then print(concat(tab3,"d *= eps;")),
- print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
- print(concat(tab3,"break;"))),
- print(" }
- }"),
- 'done)$
- codeC2(maxpow):=block([tab2:" ",tab3:" "],
- print("// The coefficients C2[l] in the Fourier expansion of B2
- static inline void evaluate_coeffs_C2(CT const& eps, CT c[])
- {
- CT const eps2 = math::sqr(eps);
- CT d = eps;
- switch (SeriesOrder) {"),
- for n:0 thru maxpow do (
- print(concat(tab2,"case ",string(n),":")),
- for m:1 thru n do block([q:d*horner(
- subst([eps=sqrt(eps2)],ataylor(C2[m],eps,n)/eps^m)),
- linel:1200],
- if m>1 then print(concat(tab3,"d *= eps;")),
- print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
- print(concat(tab3,"break;"))),
- print(" }
- }"),
- 'done)$
- codeC3(maxpow):=block([tab2:" ",tab3:" "],
- print("// The coefficients C3[l] in the Fourier expansion of B3
- static inline void evaluate_coeffs_C3(CT const& n, CT c[])
- {
- const CT n2 = math::sqr(n);
- switch (SeriesOrder) {"),
- for nn:0 thru maxpow do block([c],
- print(concat(tab2,"case ",string(nn),":")),
- c:0,
- for m:1 thru nn-1 do block(
- [q:if nn = 0 then 0 else
- jtaylor(subst([n=n],C3[m]),_n,eps,nn-1),
- linel:1200],
- for j:m thru nn-1 do (
- print(concat(tab3,"c[",c,"] = ",
- string(horner(coeff(q,eps,j))),";")),
- c:c+1)
- ),
- print(concat(tab3,"break;"))),
- print(" }
- }"),
- 'done)$
- maxpow:8$
- computeI1(maxpow)$
- computeI2(maxpow)$
- computeI3(maxpow)$
- revertI1(maxpow)$
- codeA1(maxpow)$
- codeA2(maxpow)$
- codeA3(maxpow)$
- codeC1(maxpow)$
- codeC2(maxpow)$
- codeC3(maxpow)$
- codeC1p(maxpow)$
|