610 lines
19 KiB
Plaintext
610 lines
19 KiB
Plaintext
/*
|
|
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
|
|
|
|
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),
|
|
int:subst([abs(eps-1)=1-eps,abs(n-1)=1-n],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. */
|
|
|
|
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,"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,"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,"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,"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,"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,"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,"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,"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)
|
|
|
|
For computeI4(n)
|
|
n time(s)
|
|
8 1.5
|
|
16 29
|
|
24 215
|
|
30 702 = 11.7 m
|
|
32 1040 = 17.3 m
|
|
40 3851
|
|
48 11521
|
|
64 69960 = 19.4 h
|
|
|
|
*/
|
|
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))$ */
|