1/*
2Compute the series expansion for the great ellipse area.
3
4Copyright (c) Charles Karney (2014) <charles@karney.com> and licensed
5under the MIT/X11 License.  For more information, see
6https://geographiclib.sourceforge.io/
7
8Area of great ellipse quad
9
10area from equator to phi
11dA = dlambda*b^2*sin(phi)/2*
12       (1/(1-e^2*sin(phi)^2) + atanh(e*sin(phi))/(e*sin(phi)))
13Total area = 2*pi*(a^2 + b^2*atanh(e)/e)
14
15convert to beta using
16  sin(phi)^2 = sin(beta)^2/(1-e^2*cos(beta)^2)
17
18dA = dlambda*sin(beta)/2*
19  ( a^2*sqrt(1-e^2*cos(beta)^2)
20  + b^2*atanh(e*sin(beta)/sqrt(1-e^2*cos(beta)^2))
21          /(e*sin(beta))) )
22
23subst for great ellipse
24  sin(beta) = cos(gamma0)*sin(sigma)
25  dlambda = dsigma * sin(gamma0) / (1 - cos(gamma0)^2*sin(sigma)^2)
26  (1-e^2*cos(beta)^2) = (1-e^2)*(1+k^2*sin(sigma)^2)
27  k^2 = ep^2*cos(gamma0)^2
28  e*sin(beta) = sqrt(1-e^2)*k*sin(sigma)
29
30dA = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*sqrt(1+ep^2)*
31    ( (1+k^2*sin(sigma)^2)
32    + atanh(k*sin(sigma)/sqrt(1+k^2*sin(sigma)^2)) / (k*sin(sigma))) )
33    / (ep^2 - k^2*sin(sigma)^2)
34
35Spherical case radius c, c^2 = a^2/2 + b^2/2*atanh(e)/e
36dA0 = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*
37    ( a^2 + b^2*atanh(e)/e)
38    / (1 - cos(gamma0)^2*sin(sigma)^2)
39   = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*sqrt(1+ep^2)
40    ( sqrt(1+ep^2) + atanh(ep/sqrt(1+ep^2))/ep )
41    / (ep^2 - k^2*sin(sigma)^2)
42
43dA-dA0 = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*sqrt(1+ep^2)*
44    -( ( sqrt(1+ep^2)
45        + atanh(ep/sqrt(1+ep^2))/ep ) -
46       ( sqrt(1+k^2*sin(sigma)^2)
47        + atanh(k*sin(sigma)/sqrt(1+k^2*sin(sigma)^2)) / (k*sin(sigma)) ) ) /
48     / (ep^2 - k^2*sin(sigma)^2)
49
50atanh(y/sqrt(1+y^2)) = asinh(y)
51
52dA-dA0 = dsigma*sin(gamma0)*cos(gamma0)*sin(sigma)/2*e^2*a^2*
53    - ( ( sqrt(1+ep^2) + asinh(ep)/ep ) -
54        ( sqrt(1+k^2*sin(sigma)^2)
55         + asinh(k*sin(sigma)) / (k*sin(sigma)) ) ) /
56      / (ep^2 - k^2*sin(sigma)^2)
57
58r(x) = sqrt(1+x) + asinh(sqrt(x))/sqrt(x)
59dA-dA0 = e^2*a^2*dsigma*sin(gamma0)*cos(gamma0)*
60    - ( r(ep^2) - r(k^2*sin(sigma)^2)) /
61      (   ep^2  -   k^2*sin(sigma)^2 ) *
62      sin(sigma)/2*sqrt(1+ep^2)*
63
64subst
65  ep^2=4*n/(1-n)^2 -- second eccentricity in terms of third flattening
66  ellipse semi axes = [a, a * sqrt(1-e^2*cos(gamma0)^2)]
67  third flattening for ellipsoid
68  eps = (1 - sqrt(1-e^2*cos(gamma0)^2)) / (1 + sqrt(1-e^2*cos(gamma0)^2))
69  e^2*cos(gamma0)^2 = 4*eps/(1+eps)^2 -- 1st ecc in terms of 3rd flattening
70  k2=((1+n)/(1-n))^2 * 4*eps/(1+eps)^2
71
72Taylor expand in n and eps, integrate, trigreduce
73*/
74
75taylordepth:5$
76jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
77    ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
78
79/* Great ellipse via r */
80computeG4(maxpow):=block([int,r,intexp,area, x,ep2,k2],
81  maxpow:maxpow-1,
82  r : sqrt(1+x) + asinh(sqrt(x))/sqrt(x),
83  int:-(rf(ep2) - rf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2)
84  * sin(sigma)/2*sqrt(1+ep2),
85  int:subst([rf(ep2)=subst([x=ep2],r),
86    rf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],r)],
87    int),
88  int:subst([abs(sin(sigma))=sin(sigma)],int),
89  int:subst([k2=((1+n)/(1-n))^2 * 4*eps/(1+eps)^2,ep2=4*n/(1-n)^2],int),
90  intexp:jtaylor(int,n,eps,maxpow),
91  area:trigreduce(integrate(intexp,sigma)),
92  area:expand(area-subst(sigma=%pi/2,area)),
93  for i:0 thru maxpow do G4[i]:coeff(area,cos((2*i+1)*sigma)),
94  if expand(area-sum(G4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0
95  then error("left over terms in G4"),
96  'done)$
97
98codeG4(maxpow):=block([tab1:"  ",nn:maxpow,c],
99  c:0,
100  for m:0 thru nn-1 do block(
101    [q:jtaylor(G4[m],n,eps,nn-1), linel:1200],
102    for j:m thru nn-1 do (
103      print(concat(tab1,"G4x(",c,"+1) = ",
104          string(horner(coeff(q,eps,j))),";")),
105      c:c+1)
106    ),
107  'done)$
108dispG4(ord):=(ord:ord-1,for i:0 thru ord do
109block([tt:jtaylor(G4[i],n,eps,ord),
110  ttt,t4,linel:1200],
111  for j:i thru ord do (
112    ttt:coeff(tt,eps,j),
113    if ttt # 0 then block([a:taylor(ttt+n^(ord+1),n,0,ord+1),paren,s],
114      paren : is(length(a) > 2),
115      s:if j = i then concat("G4[",i,"] = ") else "        ",
116      if subst([n=1],part(a,1)) > 0 then s:concat(s,"+ ")
117      else (s:concat(s,"- "), a:-a),
118      if paren then s:concat(s,"("),
119      for k:1 thru length(a)-1 do block([term:part(a,k),nn],
120        nn:subst([n=1],term),
121        term:term/nn,
122        if nn > 0 and k > 1 then s:concat(s," + ")
123        else if nn < 0 then (s:concat(s," - "), nn:-nn),
124        if lopow(term,n) = 0 then s:concat(s,string(nn))
125        else (
126          if nn # 1 then s:concat(s,string(nn)," * "),
127          s:concat(s,string(term))
128          )),
129      if paren then s:concat(s,")"),
130      if j>0 then s:concat(s," * ", string(eps^j)),
131      print(concat(s,if j = ord then ";" else ""))))))$
132
133maxpow:6$
134(computeG4(maxpow),
135  print(""),
136  codeG4(maxpow),
137  print(""),
138  dispG4(maxpow))$
139