1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
/*
Compute the series expansion for the great ellipse area.
Copyright (c) Charles Karney (2014) <karney@alum.mit.edu> 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))$
|
