/* Compute the series expansions for the ellipsoidal geodesic problem. Copyright (c) Charles Karney (2009-2015) and licensed under the MIT/X11 License. For more information, see https://geographiclib.sourceforge.io/ 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 There are 4 sections in this file (1) Functions to compute the expansions (2) Functions to print C++ code (3) Functions to display the results (4) Calls to the above. Edit the section at the end, to modify what is done. As distributed this code computes the 8th order series. This takes about 10 secs. If you want to compute accurate geodesics using geodesic.mac, then you need alse to uncomment the last line of this file so that the series get saved as geodNN.lsp. To run the code, start Maxima and enter load("geod.mac")$ */ /* EXPANSIONS FOR INTEGRALS */ 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)))))$ /* Express I1 = integrate( sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma ) as a series A1 * ( sigma + sum(C1[l] * sin(2*l*sigma), l, 1, maxpow) ) valid for k2 small. It is convenient to write k2 = 4 * eps / (1 - eps)^2 and to expand (1 - eps) * I1 retaining terms up to order eps^maxpow in A1 and C1[l]. This leads to a series where half the terms drop out. */ computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps], sintegrand:sqrt(1+k2*sin(sigma)^2), /* Multiplicative factor 1/(1-eps) */ 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)$ /* Write tau1 = sigma + sum(C1[l] * sin(2*l*sigma), l, 1, maxpow) and revert this to obtain sigma = tau1 + sum(C1p[l] * sin(2*tau1), l, 1, maxpow) retaining terms up to order eps^maxpow in tp[l]. Write tau = sigma + B1(sigma) sigma = tau + B1p(tau) B1(sigma) = sum(C1[l] * sin(2*l*sigma), l, 1, inf) B1p(tau) = sum(C1p[l] * sin(2*tau), l, 1, inf) Then the Lagrange Inversion Theorem J. L. Lagrange, Nouvelle methode pour resoudre les equations litterales par le moyen des series, Mem. de l'Acad. Roy. des Sciences de Berlin 24, 251-326 (1768, publ. 1770), Sec. 16, https://books.google.com/books?id=YywPAAAAIAAJ&pg=PA25 gives B1p(tau) = sum( (-1)^n/n! * diff( B1(tau)^n, tau, n-1 ), n, 1, inf) Call this after computeI1(maxpow)$ */ 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)$ /* Express I2 = integrate( 1/sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma ) as a series A2 * ( sigma + sum(C2[l] * sin(2*l*sigma), l, 1, maxpow) ) valid for k2 small. It is convenient to write k2 = 4 * eps / (1 - eps)^2 and to expand 1/(1 - eps) * I2 retaining terms up to order eps^maxpow in A2 and C2[l]. This leads to a series where half the terms drop out. */ computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps], sintegrand:1/sqrt(1+k2*sin(sigma)^2), /* Multiplicative factor 1/(1+eps) */ 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)$ /* Express I3 = integrate( (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma1)^2)), sigma1, 0, sigma ) as a series A3 * ( sigma + sum(C3[l] * sin(2*l*sigma), l, 1, maxpow-1) ) valid for f and k2 small. It is convenient to write k2 = 4 * eps / (1 - eps)^2 and f = 2*n/(1+n) and expand in eps and n. This procedure leads to a series where the coefficients of eps^j are terminating series in n. */ 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)$ /* Express I4 = -integrate( (t(ep2) - t(k2*sin(sigma1)^2)) / (ep2 - k2*sin(sigma1)^2) * sin(sigma1)/2, sigma1, pi/2, sigma ) where t(x) = sqrt(1+1/x)*asinh(sqrt(x)) + x as a series sum(C4[l] * cos((2*l+1)*sigma), l, 0, maxpow-1) ) valid for ep2 and k2 small. It is convenient to write k2 = 4 * eps / (1 - eps)^2 and ep2 = 4 * n / (1 - n)^2 and to expand in eps and n. This procedure leads to a series which converges even for very eccentric ellipsoids. */ computeI4(maxpow):=block([int,t,intexp,area, x,ep2,k2], maxpow:maxpow-1, t : sqrt(1+1/x) * asinh(sqrt(x)) + x, int:-(tf(ep2) - tf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2) * sin(sigma)/2, int:subst([tf(ep2)=subst([x=ep2],t), tf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],t)], int), int:subst([abs(sin(sigma))=sin(sigma)],int), int:subst([k2=4*eps/(1-eps)^2,ep2=4*n/(1-n)^2],int), intexp:jtaylor(int,n,eps,maxpow), area:trigreduce(integrate(intexp,sigma)), area:expand(area-subst(sigma=%pi/2,area)), for i:0 thru maxpow do C4[i]:coeff(area,cos((2*i+1)*sigma)), if expand(area-sum(C4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0 then error("left over terms in I4"), 'done)$ /* Call all of the above */ computeall():=( computeI1(maxpow), revertI1(maxpow), computeI2(maxpow), computeI3(maxpow), computeI4(maxpow))$ /* FORMAT FOR C++ */ /* If nA1, nC1, nC1p, nA2, nA3, nC3 are compile-time constants indicating the required order, the compiler will include only the needed code. GEOGRAPHICLIB_STATIC_ASSERT is a macro to cause a compile-time error if the assertion is false. */ codeA1(maxpow):=block([tab2:" ",tab3:" "], print(" // The scale factor A1-1 = mean value of (d/dsigma)I1 - 1 Math::real Geodesic::A1m1f(real eps) { real eps2 = Math::sq(eps), t; switch (nA1_/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(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nA1_ >= ",string(0), " && nA1_ <= ",string(maxpow),", \"Bad value of nA1_\");")), print(concat(tab3,"t = 0;")), print(" } return (t + eps) / (1 - eps); }"), 'done)$ codeC1(maxpow):=block([tab2:" ",tab3:" "], print(" // The coefficients C1[l] in the Fourier expansion of B1 void Geodesic::C1f(real eps, real c[]) { real eps2 = Math::sq(eps), d = eps; switch (nC1_) {"), 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(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nC1_ >= ",string(0), " && nC1_ <= ",string(maxpow),", \"Bad value of nC1_\");")), print(" } }"), 'done)$ codeC1p(maxpow):=block([tab2:" ",tab3:" "], print(" // The coefficients C1p[l] in the Fourier expansion of B1p void Geodesic::C1pf(real eps, real c[]) { real eps2 = Math::sq(eps), d = eps; switch (nC1p_) {"), 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(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nC1p_ >= ",string(0), " && nC1p_ <= ",string(maxpow),", \"Bad value of nC1p_\");")), print(" } }"), 'done)$ codeA2(maxpow):=block([tab2:" ",tab3:" "], print(" // The scale factor A2-1 = mean value of (d/dsigma)I2 - 1 Math::real Geodesic::A2m1f(real eps) { real eps2 = Math::sq(eps), t; switch (nA2_/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(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nA2_ >= ",string(0), " && nA2_ <= ",string(maxpow),", \"Bad value of nA2_\");")), print(concat(tab3,"t = 0;")), print(" } return (t - eps) / (1 + eps); }"), 'done)$ codeC2(maxpow):=block([tab2:" ",tab3:" "], print(" // The coefficients C2[l] in the Fourier expansion of B2 void Geodesic::C2f(real eps, real c[]) { real eps2 = Math::sq(eps), d = eps; switch (nC2_) {"), 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(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nC2_ >= ",string(0), " && nC2_ <= ",string(maxpow),", \"Bad value of nC2_\");")), print(" } }"), 'done)$ codeA3(maxpow):=block([tab2:" ",tab3:" "], print(" // The scale factor A3 = mean value of (d/dsigma)I3 void Geodesic::A3coeff() { switch (nA3_) {"), 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,"_A3x[",i,"] = ", string(horner(coeff(q,eps,i))),";")), print(concat(tab3,"break;"))), print(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nA3_ >= ",string(0), " && nA3_ <= ",string(maxpow),", \"Bad value of nA3_\");")), print(" } }"), 'done)$ codeC3(maxpow):=block([tab2:" ",tab3:" "], print(" // The coefficients C3[l] in the Fourier expansion of B3 void Geodesic::C3coeff() { switch (nC3_) {"), 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,"_C3x[",c,"] = ", string(horner(coeff(q,eps,j))),";")), c:c+1) ), print(concat(tab3,"break;"))), print(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nC3_ >= ",string(0), " && nC3_ <= ",string(maxpow),", \"Bad value of nC3_\");")), print(" } }"), 'done)$ codeC4(maxpow):=block([tab2:" ",tab3:" "], print(" // The coefficients C4[l] in the Fourier expansion of I4 void Geodesic::C4coeff() { switch (nC4_) {"), for nn:0 thru maxpow do block([c], print(concat(tab2,"case ",string(nn),":")), c:0, for m:0 thru nn-1 do block( [q:jtaylor(subst([n=_n],C4[m]),_n,eps,nn-1), linel:1200], for j:m thru nn-1 do ( print(concat(tab3,"_C4x[",c,"] = ", string(horner(coeff(q,eps,j))),";")), c:c+1) ), print(concat(tab3,"break;"))), print(concat(tab2,"default:")), print(concat(tab3,"GEOGRAPHICLIB_STATIC_ASSERT(nC4_ >= ",string(0), " && nC4_ <= ",string(maxpow),", \"Bad value of nC4_\");")), print(" } }"), 'done)$ printcode():=( print(""), print(concat(" // Generated by Maxima on ",timedate())), print(""), codeA1(maxpow), print(""), codeC1(maxpow), print(""), codeC1p(maxpow), print(""), codeA2(maxpow), print(""), codeC2(maxpow), print(""), codeA3(maxpow), print(""), codeC3(maxpow), print(""), codeC4(maxpow))$ /* FORMAT FOR DISPLAY */ dispA1(ord):=block( [tt:ataylor(A1*(1-eps),eps,ord),ttt,linel:1200], for j:2 step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat(if j = 2 then "A1 = (1 " else " ", if ttt>0 then "+ " else "- ", string(abs(ttt)), " * ", string(eps^j), if j=ord or j = ord-1 then ") / (1 - eps);" else ""))))$ dispC1(ord):=for i:1 thru ord do block([tt:ataylor(C1[i],eps,ord),ttt,linel:1200], print(), for j:i step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat( if j = i then concat("C1[",string(i),"] = ") else " ", if ttt>0 then "+ " else "- ", string(abs(ttt)), " * ", string(eps^j), if j=ord or j=ord-1 then ";" else ""))))$ dispC1p(ord):=for i:1 thru ord do block([tt:ataylor(C1p[i],eps,ord),ttt,linel:1200], print(), for j:i step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat( if j = i then concat("C1p[",string(i),"] = ") else " ", if ttt>0 then "+ " else "- ", string(abs(ttt)), " * ", string(eps^j), if j=ord or j=ord-1 then ";" else ""))))$ dispA2(ord):=block( [tt:ataylor(A2*(1+eps),eps,ord),ttt,linel:1200], for j:2 step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat(if j = 2 then "A2 = (1 " else " ", if ttt>0 then "+ " else "- ", string(abs(ttt)), " * ", string(eps^j), if j=ord or j = ord-1 then ") / (1 + eps);" else ""))))$ dispC2(ord):=for i:1 thru ord do block([tt:ataylor(C2[i],eps,ord),ttt,linel:1200], print(), for j:i step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat( if j = i then concat("C2[",string(i),"] = ") else " ", if ttt>0 then "+ " else "- ", string(abs(ttt)), " * ", string(eps^j), if j=ord or j=ord-1 then ";" else ""))))$ dispA3(ord):=(ord:ord-1,block( [tt:jtaylor(A3,n,eps,ord),ttt,t4,linel:1200,s], for j:1 thru ord do (ttt:expand(coeff(tt,eps,j)), if ttt # 0 then block([a:taylor(ttt+n^(ord+1),n,0,ord+1),paren,s], paren : is(length(a) > 2), s:if j=1 then "A3 = 1" else " ", if subst([n=1],part(a,1)) > 0 then s:concat(s," + ") else (s:concat(s," - "), a:-a), if paren then s:concat(s,"("), for k:1 thru length(a)-1 do block([term:part(a,k),nn], nn:subst([n=1],term), term:term/nn, if nn > 0 and k > 1 then s:concat(s," + ") else if nn < 0 then (s:concat(s," - "), nn:-nn), if lopow(term,n) = 0 then s:concat(s,string(nn)) else ( if nn # 1 then s:concat(s,string(nn)," * "), s:concat(s,string(term)) )), if paren then s:concat(s,")"), s:concat(s," * ", string(eps^j)), print(concat(s,if j = ord then ";" else ""))))))$ dispC3(ord):=(ord:ord-1,for i:1 thru ord do block([tt:jtaylor(C3[i],eps,n,ord), ttt,t4,linel:1200], for j:i thru ord do ( ttt:coeff(tt,eps,j), if ttt # 0 then block([a:taylor(ttt+n^(ord+1),n,0,ord+1),paren,s], paren : is(length(a) > 2), s:if j = i then concat("C3[",i,"] = ") else " ", if subst([n=1],part(a,1)) > 0 then s:concat(s,"+ ") else (s:concat(s,"- "), a:-a), if paren then s:concat(s,"("), for k:1 thru length(a)-1 do block([term:part(a,k),nn], nn:subst([n=1],term), term:term/nn, if nn > 0 and k > 1 then s:concat(s," + ") else if nn < 0 then (s:concat(s," - "), nn:-nn), if lopow(term,n) = 0 then s:concat(s,string(nn)) else ( if nn # 1 then s:concat(s,string(nn)," * "), s:concat(s,string(term)) )), if paren then s:concat(s,")"), s:concat(s," * ", string(eps^j)), print(concat(s,if j = ord then ";" else ""))))))$ dispC4(ord):=(ord:ord-1,for i:0 thru ord do block([tt:jtaylor(C4[i],n,eps,ord), ttt,t4,linel:1200], for j:i thru ord do ( ttt:coeff(tt,eps,j), if ttt # 0 then block([a:taylor(ttt+n^(ord+1),n,0,ord+1),paren,s], paren : is(length(a) > 2), s:if j = i then concat("C4[",i,"] = ") else " ", if subst([n=1],part(a,1)) > 0 then s:concat(s,"+ ") else (s:concat(s,"- "), a:-a), if paren then s:concat(s,"("), for k:1 thru length(a)-1 do block([term:part(a,k),nn], nn:subst([n=1],term), term:term/nn, if nn > 0 and k > 1 then s:concat(s," + ") else if nn < 0 then (s:concat(s," - "), nn:-nn), if lopow(term,n) = 0 then s:concat(s,string(nn)) else ( if nn # 1 then s:concat(s,string(nn)," * "), s:concat(s,string(term)) )), if paren then s:concat(s,")"), if j>0 then s:concat(s," * ", string(eps^j)), print(concat(s,if j = ord then ";" else ""))))))$ dispseries():=( print(""), print(concat("// Generated by Maxima on ",timedate())), print(""), dispA1(maxpow), print(""), dispC1(maxpow), print(""), dispC1p(maxpow), print(""), dispA2(maxpow), print(""), dispC2(maxpow), print(""), dispA3(maxpow), print(""), dispC3(maxpow), print(""), dispC4(maxpow), print(""))$ /* CALL THE FUNCTIONS */ /* Timings for computeall(n) n time(s) 8 8 10 19 12 43 20 571 30 6671 (111m) */ maxpow:8$ computeall()$ /* printcode()$ */ dispseries()$ load("polyprint.mac")$ printgeod():= block([macro:if simplenum then "GEOGRAPHICLIB_GEODESIC_ORDER" else "GEOGRAPHICLIB_GEODESICEXACT_ORDER"], value1('(A1*(1-eps)-1),'eps,2,0), array1('C1,'eps,2,0), array1('C1p,'eps,2,0), value1('(A2*(1+eps)-1),'eps,2,0), array1('C2,eps,2,0), value2('A3,'n,'eps,1), array2('C3,'n,'eps,1), array2('C4,'n,'eps,1))$ printgeod()$ /* Save the values needed for geodesic.mac This is commented out here to avoid accidentally overwriting files in a user's directory. */ /* (file:concat("geod",maxpow,".lsp"), save(file, values, arrays))$ */