/*
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))$