1/*
2Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
3
4Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program.
5
6Use, modification and distribution is subject to the Boost Software License,
7Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
8http://www.boost.org/LICENSE_1_0.txt)
9
10This file is converted from GeographicLib, https://geographiclib.sourceforge.io
11GeographicLib is originally written by Charles Karney.
12
13Author: Charles Karney (2008-2017)
14
15Last updated version of GeographicLib: 1.49
16
17Original copyright notice:
18
19Copyright (c) Charles Karney (2009-2015) <charles@karney.com> and
20licensed under the MIT/X11 License. For more information, see
21https://geographiclib.sourceforge.io
22
23Compute the series expansions for the ellipsoidal geodesic problem.
24
25References:
26
27   Charles F. F. Karney,
28   Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
29   https://doi.org/10.1007/s00190-012-0578-z
30   Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
31
32The code below contains minor modifications to conform with
33Boost Geometry style guidelines.
34
35To run the code, start Maxima and enter
36
37  load("geod.mac")$
38*/
39
40taylordepth:5$
41ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$
42jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
43ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
44
45computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
46  sintegrand:sqrt(1+k2*sin(sigma)^2),
47  sintegrandexp:ataylor(
48      (1-eps)*subst([k2=4*eps/(1-eps)^2],sintegrand),
49      eps,maxpow),
50  s:trigreduce(integrate(sintegrandexp,sigma)),
51  s:s-subst(sigma=0,s),
52  A1:expand(subst(sigma=2*%pi,s)/(2*%pi)),
53  tau1:ataylor(s/A1,eps,maxpow),
54  for i:1 thru maxpow do C1[i]:coeff(tau1,sin(2*i*sigma)),
55  if expand(tau1-sigma-sum(C1[i]*sin(2*i*sigma),i,1,maxpow)) # 0
56  then error("left over terms in B1"),
57  A1:A1/(1-eps),
58  'done)$
59
60codeA1(maxpow):=block([tab2:"    ",tab3:"        "],
61print("// The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
62static inline CT evaluate_series_A1(CT eps) {
63    CT eps2 = math::sqr(eps);
64    CT t;
65    switch (SeriesOrder/2) {"),
66  for n:0 thru entier(maxpow/2) do block([
67    q:horner(ataylor(subst([eps=sqrt(eps2)],A1*(1-eps)-1),eps2,n)),
68    linel:1200],
69    print(concat(tab2,"case ",string(n),":")),
70    print(concat(tab3,"t = ",string(q),";")),
71    print(concat(tab3,"break;"))),
72  print("    }
73    return (t + eps) / (1 - eps);
74}"),
75'done)$
76
77computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
78  sintegrand:1/sqrt(1+k2*sin(sigma)^2),
79  sintegrandexp:ataylor(
80      (1+eps)*subst([k2=4*eps/(1-eps)^2],sintegrand),
81      eps,maxpow),
82  s:trigreduce(integrate(sintegrandexp,sigma)),
83  s:s-subst(sigma=0,s),
84  A2:expand(subst(sigma=2*%pi,s)/(2*%pi)),
85  tau1:ataylor(s/A2,eps,maxpow),
86  for i:1 thru maxpow do C2[i]:coeff(tau1,sin(2*i*sigma)),
87  if expand(tau1-sigma-sum(C2[i]*sin(2*i*sigma),i,1,maxpow)) # 0
88  then error("left over terms in B2"),
89  A2:A2/(1+eps),
90  'done)$
91
92codeA2(maxpow):=block([tab2:"    ",tab3:"        "],
93print("// The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
94CT evaluate_series_A2(CT const& eps)
95{
96    CT const eps2 = math::sqr(eps);
97    CT t;
98    switch (SeriesOrder/2) {"),
99  for n:0 thru entier(maxpow/2) do block([
100    q:horner(ataylor(subst([eps=sqrt(eps2)],A2*(1+eps)-1),eps2,n)),
101    linel:1200],
102    print(concat(tab2,"case ",string(n),":")),
103    print(concat(tab3,"t = ",string(q),";")),
104    print(concat(tab3,"break;"))),
105  print("    }
106    return (t - eps) / (1 + eps);
107}"),
108'done)$
109
110computeI3(maxpow):=block([int,intexp,dlam,eta,del,eps,nu,f,z,n],
111  maxpow:maxpow-1,
112  int:subst([k2=4*eps/(1-eps)^2],
113    (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma)^2))),
114  int:subst([f=2*n/(1+n)],int),
115  intexp:jtaylor(int,n,eps,maxpow),
116  dlam:trigreduce(integrate(intexp,sigma)),
117  dlam:dlam-subst(sigma=0,dlam),
118  A3:expand(subst(sigma=2*%pi,dlam)/(2*%pi)),
119  eta:jtaylor(dlam/A3,n,eps,maxpow),
120  A3:jtaylor(A3,n,eps,maxpow),
121  for i:1 thru maxpow do C3[i]:coeff(eta,sin(2*i*sigma)),
122  if expand(eta-sigma-sum(C3[i]*sin(2*i*sigma),i,1,maxpow)) # 0
123  then error("left over terms in B3"),
124  'done)$
125
126codeA3(maxpow):=block([tab2:"    ",tab3:"        "],
127print("// The scale factor A3 = mean value of (d/dsigma)I3
128static inline void evaluate_series_A3(CT const& n, CT c[])
129{
130    switch (SeriesOrder) {"),
131  for nn:0 thru maxpow do block(
132    [q:if nn=0 then 0 else
133    jtaylor(subst([n=n],A3),n,eps,nn-1),
134    linel:1200],
135    print(concat(tab2,"case ",string(nn),":")),
136    for i : 0 thru nn-1 do
137    print(concat(tab3,"c[",i,"] = ",
138        string(horner(coeff(q,eps,i))),";")),
139    print(concat(tab3,"break;"))),
140  print("    }
141}"),
142'done)$
143
144codeC1(maxpow):=block([tab2:"    ",tab3:"        "],
145  print("// The coefficients C1[l] in the Fourier expansion of B1
146static inline evaluate_coeffs_C1(CT eps, CT c[]) {
147    CT eps2 = math::sqr(eps);
148    CT d = eps;
149    switch (SeriesOrder) {"),
150  for n:0 thru maxpow do (
151    print(concat(tab2,"case ",string(n),":")),
152    for m:1 thru n do block([q:d*horner(
153        subst([eps=sqrt(eps2)],ataylor(C1[m],eps,n)/eps^m)),
154      linel:1200],
155      if m>1 then print(concat(tab3,"d *= eps;")),
156      print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
157    print(concat(tab3,"break;"))),
158  print("    }
159}"),
160'done)$
161
162revertI1(maxpow):=block([tau,eps,tauacc:1,sigacc:0],
163  for n:1 thru maxpow do (
164    tauacc:trigreduce(ataylor(
165          -sum(C1[j]*sin(2*j*tau),j,1,maxpow-n+1)*tauacc/n,
166          eps,maxpow)),
167    sigacc:sigacc+expand(diff(tauacc,tau,n-1))),
168  for i:1 thru maxpow do C1p[i]:coeff(sigacc,sin(2*i*tau)),
169  if expand(sigacc-sum(C1p[i]*sin(2*i*tau),i,1,maxpow)) # 0
170  then error("left over terms in B1p"),
171  'done)$
172
173codeC1p(maxpow):=block([tab2:"    ",tab3:"        "],
174  print("// The coefficients C1p[l] in the Fourier expansion of B1p
175static inline evaluate_coeffs_C1p(CT eps, CT c[])
176{
177    CT const eps2 = math::sqr(eps);
178    CT d = eps;
179    switch (SeriesOrder) {"),
180  for n:0 thru maxpow do (
181    print(concat(tab2,"case ",string(n),":")),
182    for m:1 thru n do block([q:d*horner(
183        subst([eps=sqrt(eps2)],ataylor(C1p[m],eps,n)/eps^m)),
184      linel:1200],
185      if m>1 then print(concat(tab3,"d *= eps;")),
186      print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
187    print(concat(tab3,"break;"))),
188  print("    }
189}"),
190'done)$
191
192codeC2(maxpow):=block([tab2:"    ",tab3:"        "],
193print("// The coefficients C2[l] in the Fourier expansion of B2
194static inline void evaluate_coeffs_C2(CT const& eps, CT c[])
195{
196    CT const eps2 = math::sqr(eps);
197    CT d = eps;
198    switch (SeriesOrder) {"),
199  for n:0 thru maxpow do (
200    print(concat(tab2,"case ",string(n),":")),
201    for m:1 thru n do block([q:d*horner(
202        subst([eps=sqrt(eps2)],ataylor(C2[m],eps,n)/eps^m)),
203      linel:1200],
204      if m>1 then print(concat(tab3,"d *= eps;")),
205      print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
206    print(concat(tab3,"break;"))),
207print("    }
208}"),
209'done)$
210
211codeC3(maxpow):=block([tab2:"    ",tab3:"        "],
212print("// The coefficients C3[l] in the Fourier expansion of B3
213static inline void evaluate_coeffs_C3(CT const& n, CT c[])
214{
215    const CT n2 = math::sqr(n);
216    switch (SeriesOrder) {"),
217  for nn:0 thru maxpow do block([c],
218    print(concat(tab2,"case ",string(nn),":")),
219    c:0,
220    for m:1 thru nn-1 do block(
221      [q:if nn = 0 then 0 else
222      jtaylor(subst([n=n],C3[m]),_n,eps,nn-1),
223      linel:1200],
224      for j:m thru nn-1 do (
225        print(concat(tab3,"c[",c,"] = ",
226            string(horner(coeff(q,eps,j))),";")),
227        c:c+1)
228    ),
229    print(concat(tab3,"break;"))),
230  print("    }
231}"),
232'done)$
233
234maxpow:8$
235
236computeI1(maxpow)$
237computeI2(maxpow)$
238computeI3(maxpow)$
239
240revertI1(maxpow)$
241codeA1(maxpow)$
242codeA2(maxpow)$
243codeA3(maxpow)$
244
245codeC1(maxpow)$
246codeC2(maxpow)$
247codeC3(maxpow)$
248
249codeC1p(maxpow)$
250