forked from boostorg/geometry
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'feature/geodesic_direct' into feature/karney_inverse
- Loading branch information
Showing
2 changed files
with
246 additions
and
285 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,232 @@ | ||
/* | ||
Compute the series expansions for the ellipsoidal geodesic problem. | ||
|
||
Copyright (c) Charles Karney (2009-2015) <charles@karney.com> 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 | ||
|
||
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)$ |
Oops, something went wrong.