139 lines
4.7 KiB
Plaintext
139 lines
4.7 KiB
Plaintext
/*
|
|
Compute the series expansion for the great ellipse area.
|
|
|
|
Copyright (c) Charles Karney (2014) <charles@karney.com> and licensed
|
|
under the MIT/X11 License. For more information, see
|
|
https://geographiclib.sourceforge.io/
|
|
|
|
Area of great ellipse quad
|
|
|
|
area from equator to phi
|
|
dA = dlambda*b^2*sin(phi)/2*
|
|
(1/(1-e^2*sin(phi)^2) + atanh(e*sin(phi))/(e*sin(phi)))
|
|
Total area = 2*pi*(a^2 + b^2*atanh(e)/e)
|
|
|
|
convert to beta using
|
|
sin(phi)^2 = sin(beta)^2/(1-e^2*cos(beta)^2)
|
|
|
|
dA = dlambda*sin(beta)/2*
|
|
( a^2*sqrt(1-e^2*cos(beta)^2)
|
|
+ b^2*atanh(e*sin(beta)/sqrt(1-e^2*cos(beta)^2))
|
|
/(e*sin(beta))) )
|
|
|
|
subst for great ellipse
|
|
sin(beta) = cos(gamma0)*sin(sigma)
|
|
dlambda = dsigma * sin(gamma0) / (1 - cos(gamma0)^2*sin(sigma)^2)
|
|
(1-e^2*cos(beta)^2) = (1-e^2)*(1+k^2*sin(sigma)^2)
|
|
k^2 = ep^2*cos(gamma0)^2
|
|
e*sin(beta) = sqrt(1-e^2)*k*sin(sigma)
|
|
|
|
dA = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*sqrt(1+ep^2)*
|
|
( (1+k^2*sin(sigma)^2)
|
|
+ atanh(k*sin(sigma)/sqrt(1+k^2*sin(sigma)^2)) / (k*sin(sigma))) )
|
|
/ (ep^2 - k^2*sin(sigma)^2)
|
|
|
|
Spherical case radius c, c^2 = a^2/2 + b^2/2*atanh(e)/e
|
|
dA0 = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*
|
|
( a^2 + b^2*atanh(e)/e)
|
|
/ (1 - cos(gamma0)^2*sin(sigma)^2)
|
|
= dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*sqrt(1+ep^2)
|
|
( sqrt(1+ep^2) + atanh(ep/sqrt(1+ep^2))/ep )
|
|
/ (ep^2 - k^2*sin(sigma)^2)
|
|
|
|
dA-dA0 = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*sqrt(1+ep^2)*
|
|
-( ( sqrt(1+ep^2)
|
|
+ atanh(ep/sqrt(1+ep^2))/ep ) -
|
|
( sqrt(1+k^2*sin(sigma)^2)
|
|
+ atanh(k*sin(sigma)/sqrt(1+k^2*sin(sigma)^2)) / (k*sin(sigma)) ) ) /
|
|
/ (ep^2 - k^2*sin(sigma)^2)
|
|
|
|
atanh(y/sqrt(1+y^2)) = asinh(y)
|
|
|
|
dA-dA0 = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*
|
|
- ( ( sqrt(1+ep^2) + asinh(ep)/ep ) -
|
|
( sqrt(1+k^2*sin(sigma)^2)
|
|
+ asinh(k*sin(sigma)) / (k*sin(sigma)) ) ) /
|
|
/ (ep^2 - k^2*sin(sigma)^2)
|
|
|
|
r(x) = sqrt(1+x) + asinh(sqrt(x))/sqrt(x)
|
|
dA-dA0 = e^2*a^2*dsigma*sin(gamma0)*cos(gamma0)*
|
|
- ( r(ep^2) - r(k^2*sin(sigma)^2)) /
|
|
( ep^2 - k^2*sin(sigma)^2 ) *
|
|
sin(sigma)/2*sqrt(1+ep^2)*
|
|
|
|
subst
|
|
ep^2=4*n/(1-n)^2 -- second eccentricity in terms of third flattening
|
|
ellipse semi axes = [a, a * sqrt(1-e^2*cos(gamma0)^2)]
|
|
third flattening for ellipsoid
|
|
eps = (1 - sqrt(1-e^2*cos(gamma0)^2)) / (1 + sqrt(1-e^2*cos(gamma0)^2))
|
|
e^2*cos(gamma0)^2 = 4*eps/(1+eps)^2 -- 1st ecc in terms of 3rd flattening
|
|
k2=((1+n)/(1-n))^2 * 4*eps/(1+eps)^2
|
|
|
|
Taylor expand in n and eps, integrate, trigreduce
|
|
*/
|
|
|
|
taylordepth:5$
|
|
jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
|
|
ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
|
|
|
|
/* Great ellipse via r */
|
|
computeG4(maxpow):=block([int,r,intexp,area, x,ep2,k2],
|
|
maxpow:maxpow-1,
|
|
r : sqrt(1+x) + asinh(sqrt(x))/sqrt(x),
|
|
int:-(rf(ep2) - rf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2)
|
|
* sin(sigma)/2*sqrt(1+ep2),
|
|
int:subst([rf(ep2)=subst([x=ep2],r),
|
|
rf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],r)],
|
|
int),
|
|
int:subst([abs(sin(sigma))=sin(sigma)],int),
|
|
int:subst([k2=((1+n)/(1-n))^2 * 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 G4[i]:coeff(area,cos((2*i+1)*sigma)),
|
|
if expand(area-sum(G4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0
|
|
then error("left over terms in G4"),
|
|
'done)$
|
|
|
|
codeG4(maxpow):=block([tab1:" ",nn:maxpow,c],
|
|
c:0,
|
|
for m:0 thru nn-1 do block(
|
|
[q:jtaylor(G4[m],n,eps,nn-1), linel:1200],
|
|
for j:m thru nn-1 do (
|
|
print(concat(tab1,"G4x(",c,"+1) = ",
|
|
string(horner(coeff(q,eps,j))),";")),
|
|
c:c+1)
|
|
),
|
|
'done)$
|
|
dispG4(ord):=(ord:ord-1,for i:0 thru ord do
|
|
block([tt:jtaylor(G4[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("G4[",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 ""))))))$
|
|
|
|
maxpow:6$
|
|
(computeG4(maxpow),
|
|
print(""),
|
|
codeG4(maxpow),
|
|
print(""),
|
|
dispG4(maxpow))$
|