1/* 2Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan. 3 4Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program. 5 6Use, modification and distribution is subject to the Boost Software License, 7Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 8http://www.boost.org/LICENSE_1_0.txt) 9 10This file is converted from GeographicLib, https://geographiclib.sourceforge.io 11GeographicLib is originally written by Charles Karney. 12 13Author: Charles Karney (2008-2017) 14 15Last updated version of GeographicLib: 1.49 16 17Original copyright notice: 18 19Copyright (c) Charles Karney (2009-2015) <charles@karney.com> and 20licensed under the MIT/X11 License. For more information, see 21https://geographiclib.sourceforge.io 22 23Compute the series expansions for the ellipsoidal geodesic problem. 24 25References: 26 27 Charles F. F. Karney, 28 Algorithms for geodesics, J. Geodesy 87, 43-55 (2013), 29 https://doi.org/10.1007/s00190-012-0578-z 30 Addenda: https://geographiclib.sourceforge.io/geod-addenda.html 31 32The code below contains minor modifications to conform with 33Boost Geometry style guidelines. 34 35To run the code, start Maxima and enter 36 37 load("geod.mac")$ 38*/ 39 40taylordepth:5$ 41ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$ 42jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1], 43ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$ 44 45computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps], 46 sintegrand:sqrt(1+k2*sin(sigma)^2), 47 sintegrandexp:ataylor( 48 (1-eps)*subst([k2=4*eps/(1-eps)^2],sintegrand), 49 eps,maxpow), 50 s:trigreduce(integrate(sintegrandexp,sigma)), 51 s:s-subst(sigma=0,s), 52 A1:expand(subst(sigma=2*%pi,s)/(2*%pi)), 53 tau1:ataylor(s/A1,eps,maxpow), 54 for i:1 thru maxpow do C1[i]:coeff(tau1,sin(2*i*sigma)), 55 if expand(tau1-sigma-sum(C1[i]*sin(2*i*sigma),i,1,maxpow)) # 0 56 then error("left over terms in B1"), 57 A1:A1/(1-eps), 58 'done)$ 59 60codeA1(maxpow):=block([tab2:" ",tab3:" "], 61print("// The scale factor A1-1 = mean value of (d/dsigma)I1 - 1 62static inline CT evaluate_series_A1(CT eps) { 63 CT eps2 = math::sqr(eps); 64 CT t; 65 switch (SeriesOrder/2) {"), 66 for n:0 thru entier(maxpow/2) do block([ 67 q:horner(ataylor(subst([eps=sqrt(eps2)],A1*(1-eps)-1),eps2,n)), 68 linel:1200], 69 print(concat(tab2,"case ",string(n),":")), 70 print(concat(tab3,"t = ",string(q),";")), 71 print(concat(tab3,"break;"))), 72 print(" } 73 return (t + eps) / (1 - eps); 74}"), 75'done)$ 76 77computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps], 78 sintegrand:1/sqrt(1+k2*sin(sigma)^2), 79 sintegrandexp:ataylor( 80 (1+eps)*subst([k2=4*eps/(1-eps)^2],sintegrand), 81 eps,maxpow), 82 s:trigreduce(integrate(sintegrandexp,sigma)), 83 s:s-subst(sigma=0,s), 84 A2:expand(subst(sigma=2*%pi,s)/(2*%pi)), 85 tau1:ataylor(s/A2,eps,maxpow), 86 for i:1 thru maxpow do C2[i]:coeff(tau1,sin(2*i*sigma)), 87 if expand(tau1-sigma-sum(C2[i]*sin(2*i*sigma),i,1,maxpow)) # 0 88 then error("left over terms in B2"), 89 A2:A2/(1+eps), 90 'done)$ 91 92codeA2(maxpow):=block([tab2:" ",tab3:" "], 93print("// The scale factor A2-1 = mean value of (d/dsigma)I2 - 1 94CT evaluate_series_A2(CT const& eps) 95{ 96 CT const eps2 = math::sqr(eps); 97 CT t; 98 switch (SeriesOrder/2) {"), 99 for n:0 thru entier(maxpow/2) do block([ 100 q:horner(ataylor(subst([eps=sqrt(eps2)],A2*(1+eps)-1),eps2,n)), 101 linel:1200], 102 print(concat(tab2,"case ",string(n),":")), 103 print(concat(tab3,"t = ",string(q),";")), 104 print(concat(tab3,"break;"))), 105 print(" } 106 return (t - eps) / (1 + eps); 107}"), 108'done)$ 109 110computeI3(maxpow):=block([int,intexp,dlam,eta,del,eps,nu,f,z,n], 111 maxpow:maxpow-1, 112 int:subst([k2=4*eps/(1-eps)^2], 113 (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma)^2))), 114 int:subst([f=2*n/(1+n)],int), 115 intexp:jtaylor(int,n,eps,maxpow), 116 dlam:trigreduce(integrate(intexp,sigma)), 117 dlam:dlam-subst(sigma=0,dlam), 118 A3:expand(subst(sigma=2*%pi,dlam)/(2*%pi)), 119 eta:jtaylor(dlam/A3,n,eps,maxpow), 120 A3:jtaylor(A3,n,eps,maxpow), 121 for i:1 thru maxpow do C3[i]:coeff(eta,sin(2*i*sigma)), 122 if expand(eta-sigma-sum(C3[i]*sin(2*i*sigma),i,1,maxpow)) # 0 123 then error("left over terms in B3"), 124 'done)$ 125 126codeA3(maxpow):=block([tab2:" ",tab3:" "], 127print("// The scale factor A3 = mean value of (d/dsigma)I3 128static inline void evaluate_series_A3(CT const& n, CT c[]) 129{ 130 switch (SeriesOrder) {"), 131 for nn:0 thru maxpow do block( 132 [q:if nn=0 then 0 else 133 jtaylor(subst([n=n],A3),n,eps,nn-1), 134 linel:1200], 135 print(concat(tab2,"case ",string(nn),":")), 136 for i : 0 thru nn-1 do 137 print(concat(tab3,"c[",i,"] = ", 138 string(horner(coeff(q,eps,i))),";")), 139 print(concat(tab3,"break;"))), 140 print(" } 141}"), 142'done)$ 143 144codeC1(maxpow):=block([tab2:" ",tab3:" "], 145 print("// The coefficients C1[l] in the Fourier expansion of B1 146static inline evaluate_coeffs_C1(CT eps, CT c[]) { 147 CT eps2 = math::sqr(eps); 148 CT d = eps; 149 switch (SeriesOrder) {"), 150 for n:0 thru maxpow do ( 151 print(concat(tab2,"case ",string(n),":")), 152 for m:1 thru n do block([q:d*horner( 153 subst([eps=sqrt(eps2)],ataylor(C1[m],eps,n)/eps^m)), 154 linel:1200], 155 if m>1 then print(concat(tab3,"d *= eps;")), 156 print(concat(tab3,"c[",string(m),"] = ",string(q),";"))), 157 print(concat(tab3,"break;"))), 158 print(" } 159}"), 160'done)$ 161 162revertI1(maxpow):=block([tau,eps,tauacc:1,sigacc:0], 163 for n:1 thru maxpow do ( 164 tauacc:trigreduce(ataylor( 165 -sum(C1[j]*sin(2*j*tau),j,1,maxpow-n+1)*tauacc/n, 166 eps,maxpow)), 167 sigacc:sigacc+expand(diff(tauacc,tau,n-1))), 168 for i:1 thru maxpow do C1p[i]:coeff(sigacc,sin(2*i*tau)), 169 if expand(sigacc-sum(C1p[i]*sin(2*i*tau),i,1,maxpow)) # 0 170 then error("left over terms in B1p"), 171 'done)$ 172 173codeC1p(maxpow):=block([tab2:" ",tab3:" "], 174 print("// The coefficients C1p[l] in the Fourier expansion of B1p 175static inline evaluate_coeffs_C1p(CT eps, CT c[]) 176{ 177 CT const eps2 = math::sqr(eps); 178 CT d = eps; 179 switch (SeriesOrder) {"), 180 for n:0 thru maxpow do ( 181 print(concat(tab2,"case ",string(n),":")), 182 for m:1 thru n do block([q:d*horner( 183 subst([eps=sqrt(eps2)],ataylor(C1p[m],eps,n)/eps^m)), 184 linel:1200], 185 if m>1 then print(concat(tab3,"d *= eps;")), 186 print(concat(tab3,"c[",string(m),"] = ",string(q),";"))), 187 print(concat(tab3,"break;"))), 188 print(" } 189}"), 190'done)$ 191 192codeC2(maxpow):=block([tab2:" ",tab3:" "], 193print("// The coefficients C2[l] in the Fourier expansion of B2 194static inline void evaluate_coeffs_C2(CT const& eps, CT c[]) 195{ 196 CT const eps2 = math::sqr(eps); 197 CT d = eps; 198 switch (SeriesOrder) {"), 199 for n:0 thru maxpow do ( 200 print(concat(tab2,"case ",string(n),":")), 201 for m:1 thru n do block([q:d*horner( 202 subst([eps=sqrt(eps2)],ataylor(C2[m],eps,n)/eps^m)), 203 linel:1200], 204 if m>1 then print(concat(tab3,"d *= eps;")), 205 print(concat(tab3,"c[",string(m),"] = ",string(q),";"))), 206 print(concat(tab3,"break;"))), 207print(" } 208}"), 209'done)$ 210 211codeC3(maxpow):=block([tab2:" ",tab3:" "], 212print("// The coefficients C3[l] in the Fourier expansion of B3 213static inline void evaluate_coeffs_C3(CT const& n, CT c[]) 214{ 215 const CT n2 = math::sqr(n); 216 switch (SeriesOrder) {"), 217 for nn:0 thru maxpow do block([c], 218 print(concat(tab2,"case ",string(nn),":")), 219 c:0, 220 for m:1 thru nn-1 do block( 221 [q:if nn = 0 then 0 else 222 jtaylor(subst([n=n],C3[m]),_n,eps,nn-1), 223 linel:1200], 224 for j:m thru nn-1 do ( 225 print(concat(tab3,"c[",c,"] = ", 226 string(horner(coeff(q,eps,j))),";")), 227 c:c+1) 228 ), 229 print(concat(tab3,"break;"))), 230 print(" } 231}"), 232'done)$ 233 234maxpow:8$ 235 236computeI1(maxpow)$ 237computeI2(maxpow)$ 238computeI3(maxpow)$ 239 240revertI1(maxpow)$ 241codeA1(maxpow)$ 242codeA2(maxpow)$ 243codeA3(maxpow)$ 244 245codeC1(maxpow)$ 246codeC2(maxpow)$ 247codeC3(maxpow)$ 248 249codeC1p(maxpow)$ 250