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