/*
Compute the series expansions for the ellipsoidal geodesic problem.

Copyright (c) Charles Karney (2009-2015) <karney@alum.mit.edu> and
licensed under the MIT/X11 License.  For more information, see
https://geographiclib.sourceforge.io/

References:

   Charles F. F. Karney,
   Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
   https://doi.org/10.1007/s00190-012-0578-z
   Addenda: https://geographiclib.sourceforge.io/geod-addenda.html

There are 4 sections in this file

(1) Functions to compute the expansions
(2) Functions to print C++ code
(3) Functions to display the results
(4) Calls to the above.

Edit the section at the end, to modify what is done.  As distributed
this code computes the 8th order series.  This takes about 10 secs.  If
you want to compute accurate geodesics using geodesic.mac, then you need
alse to uncomment the last line of this file so that the series get
saved as geodNN.lsp.

To run the code, start Maxima and enter

  load("geod.mac")$
*/

/* EXPANSIONS FOR INTEGRALS */

taylordepth:5$
ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$
jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
    ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$

/*

Express

    I1 = integrate( sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma )

as a series

    A1 * ( sigma + sum(C1[l] * sin(2*l*sigma), l, 1, maxpow) )

valid for k2 small.  It is convenient to write k2 = 4 * eps / (1 - eps)^2
and to expand (1 - eps) * I1 retaining terms up to order eps^maxpow
in A1 and C1[l].  This leads to a series where half the terms drop out.

*/

computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
  sintegrand:sqrt(1+k2*sin(sigma)^2),
/* Multiplicative factor 1/(1-eps) */
  sintegrandexp:ataylor(
      (1-eps)*subst([k2=4*eps/(1-eps)^2],sintegrand),
      eps,maxpow),
  s:trigreduce(integrate(sintegrandexp,sigma)),
  s:s-subst(sigma=0,s),
  A1:expand(subst(sigma=2*%pi,s)/(2*%pi)),
  tau1:ataylor(s/A1,eps,maxpow),
  for i:1 thru maxpow do C1[i]:coeff(tau1,sin(2*i*sigma)),
  if expand(tau1-sigma-sum(C1[i]*sin(2*i*sigma),i,1,maxpow)) # 0
  then error("left over terms in B1"),
  A1:A1/(1-eps),
  'done)$

/*

Write

    tau1 = sigma + sum(C1[l] * sin(2*l*sigma), l, 1, maxpow)

and revert this to obtain

    sigma = tau1 + sum(C1p[l] * sin(2*tau1), l, 1, maxpow)

retaining terms up to order eps^maxpow in tp[l].

Write

  tau = sigma + B1(sigma)
  sigma = tau + B1p(tau)
  B1(sigma) = sum(C1[l] * sin(2*l*sigma), l, 1, inf)
  B1p(tau) = sum(C1p[l] * sin(2*tau), l, 1, inf)

Then the Lagrange Inversion Theorem

   J. L. Lagrange, Nouvelle methode pour resoudre les equations
   litterales par le moyen des series, Mem. de l'Acad. Roy. des Sciences
   de Berlin 24, 251-326 (1768, publ. 1770), Sec. 16,
   https://books.google.com/books?id=YywPAAAAIAAJ&pg=PA25

gives

  B1p(tau) = sum( (-1)^n/n! * diff( B1(tau)^n, tau, n-1 ),
                         n, 1, inf)

Call this after computeI1(maxpow)$

*/

revertI1(maxpow):=block([tau,eps,tauacc:1,sigacc:0],
  for n:1 thru maxpow do (
    tauacc:trigreduce(ataylor(
          -sum(C1[j]*sin(2*j*tau),j,1,maxpow-n+1)*tauacc/n,
          eps,maxpow)),
    sigacc:sigacc+expand(diff(tauacc,tau,n-1))),
  for i:1 thru maxpow do C1p[i]:coeff(sigacc,sin(2*i*tau)),
  if expand(sigacc-sum(C1p[i]*sin(2*i*tau),i,1,maxpow)) # 0
  then error("left over terms in B1p"),
  'done)$

/*

Express

    I2 = integrate( 1/sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma )

as a series

    A2 * ( sigma + sum(C2[l] * sin(2*l*sigma), l, 1, maxpow) )

valid for k2 small.  It is convenient to write k2 = 4 * eps / (1 - eps)^2
and to expand 1/(1 - eps) * I2 retaining terms up to order eps^maxpow
in A2 and C2[l].  This leads to a series where half the terms drop out.

*/

computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
  sintegrand:1/sqrt(1+k2*sin(sigma)^2),
/* Multiplicative factor 1/(1+eps) */
  sintegrandexp:ataylor(
      (1+eps)*subst([k2=4*eps/(1-eps)^2],sintegrand),
      eps,maxpow),
  s:trigreduce(integrate(sintegrandexp,sigma)),
  s:s-subst(sigma=0,s),
  A2:expand(subst(sigma=2*%pi,s)/(2*%pi)),
  tau1:ataylor(s/A2,eps,maxpow),
  for i:1 thru maxpow do C2[i]:coeff(tau1,sin(2*i*sigma)),
  if expand(tau1-sigma-sum(C2[i]*sin(2*i*sigma),i,1,maxpow)) # 0
  then error("left over terms in B2"),
  A2:A2/(1+eps),
  'done)$

/*

Express

   I3 = integrate( (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma1)^2)), sigma1, 0, sigma )

as a series

   A3 * ( sigma + sum(C3[l] * sin(2*l*sigma), l, 1, maxpow-1) )

valid for f and k2 small.  It is convenient to write k2 = 4 * eps / (1 -
eps)^2 and f = 2*n/(1+n) and expand in eps and n.  This procedure leads
to a series where the coefficients of eps^j are terminating series in n.

*/

computeI3(maxpow):=block([int,intexp,dlam,eta,del,eps,nu,f,z,n],
  maxpow:maxpow-1,
  int:subst([k2=4*eps/(1-eps)^2],
    (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma)^2))),
  int:subst([f=2*n/(1+n)],int),
  intexp:jtaylor(int,n,eps,maxpow),
  dlam:trigreduce(integrate(intexp,sigma)),
  dlam:dlam-subst(sigma=0,dlam),
  A3:expand(subst(sigma=2*%pi,dlam)/(2*%pi)),
  eta:jtaylor(dlam/A3,n,eps,maxpow),
  A3:jtaylor(A3,n,eps,maxpow),
  for i:1 thru maxpow do C3[i]:coeff(eta,sin(2*i*sigma)),
  if expand(eta-sigma-sum(C3[i]*sin(2*i*sigma),i,1,maxpow)) # 0
  then error("left over terms in B3"),
  'done)$

/*

Express

    I4 = -integrate( (t(ep2) - t(k2*sin(sigma1)^2)) / (ep2 - k2*sin(sigma1)^2)
      * sin(sigma1)/2, sigma1, pi/2, sigma )
    where
      t(x) = sqrt(1+1/x)*asinh(sqrt(x)) + x

as a series

    sum(C4[l] * cos((2*l+1)*sigma), l, 0, maxpow-1) )

valid for ep2 and k2 small.  It is convenient to write k2 = 4 * eps /
(1 - eps)^2 and ep2 = 4 * n / (1 - n)^2 and to expand in eps and n.
This procedure leads to a series which converges even for very
eccentric ellipsoids.

*/

computeI4(maxpow):=block([int,t,intexp,area, x,ep2,k2],
  maxpow:maxpow-1,
  t : sqrt(1+1/x) * asinh(sqrt(x)) + x,
  int:-(tf(ep2) - tf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2)
  * sin(sigma)/2,
  int:subst([tf(ep2)=subst([x=ep2],t),
    tf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],t)],
    int),
  int:subst([abs(sin(sigma))=sin(sigma)],int),
  int:subst([k2=4*eps/(1-eps)^2,ep2=4*n/(1-n)^2],int),
  int:subst([abs(eps-1)=1-eps,abs(n-1)=1-n],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 C4[i]:coeff(area,cos((2*i+1)*sigma)),
  if expand(area-sum(C4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0
  then error("left over terms in I4"),
  'done)$

/* Call all of the above */
computeall():=(
  computeI1(maxpow), revertI1(maxpow),
  computeI2(maxpow), computeI3(maxpow), computeI4(maxpow))$

/* FORMAT FOR C++ */

/* If nA1, nC1, nC1p, nA2, nA3, nC3 are compile-time constants
indicating the required order, the compiler will include only the needed
code. */

codeA1(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
  Math::real Geodesic::A1m1f(real eps) {
    real
      eps2 = Math::sq(eps),
      t;
    switch (nA1_/2) {"),
  for n:0 thru entier(maxpow/2) do block([
    q:horner(ataylor(subst([eps=sqrt(eps2)],A1*(1-eps)-1),eps2,n)),
    linel:1200],
    print(concat(tab2,"case ",string(n),":")),
    print(concat(tab3,"t = ",string(q),";")),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nA1_ >= ",string(0),
      " && nA1_ <= ",string(maxpow),", \"Bad value of nA1_\");")),
  print(concat(tab3,"t = 0;")),
  print("    }
    return (t + eps) / (1 - eps);
  }"),
'done)$

codeC1(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The coefficients C1[l] in the Fourier expansion of B1
  void Geodesic::C1f(real eps, real c[]) {
    real
      eps2 = Math::sq(eps),
      d = eps;
    switch (nC1_) {"),
  for n:0 thru maxpow do (
    print(concat(tab2,"case ",string(n),":")),
    for m:1 thru n do block([q:d*horner(
        subst([eps=sqrt(eps2)],ataylor(C1[m],eps,n)/eps^m)),
      linel:1200],
      if m>1 then print(concat(tab3,"d *= eps;")),
      print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nC1_ >= ",string(0),
      " && nC1_ <= ",string(maxpow),", \"Bad value of nC1_\");")),
  print("    }
  }"),
'done)$

codeC1p(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The coefficients C1p[l] in the Fourier expansion of B1p
  void Geodesic::C1pf(real eps, real c[]) {
    real
      eps2 = Math::sq(eps),
      d = eps;
    switch (nC1p_) {"),
  for n:0 thru maxpow do (
    print(concat(tab2,"case ",string(n),":")),
    for m:1 thru n do block([q:d*horner(
        subst([eps=sqrt(eps2)],ataylor(C1p[m],eps,n)/eps^m)),
      linel:1200],
      if m>1 then print(concat(tab3,"d *= eps;")),
      print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nC1p_ >= ",string(0),
      " && nC1p_ <= ",string(maxpow),", \"Bad value of nC1p_\");")),
  print("    }
  }"),
'done)$

codeA2(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
  Math::real Geodesic::A2m1f(real eps) {
    real
      eps2 = Math::sq(eps),
      t;
    switch (nA2_/2) {"),
  for n:0 thru entier(maxpow/2) do block([
    q:horner(ataylor(subst([eps=sqrt(eps2)],A2*(1+eps)-1),eps2,n)),
    linel:1200],
    print(concat(tab2,"case ",string(n),":")),
    print(concat(tab3,"t = ",string(q),";")),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nA2_ >= ",string(0),
      " && nA2_ <= ",string(maxpow),", \"Bad value of nA2_\");")),
  print(concat(tab3,"t = 0;")),
  print("    }
    return (t - eps) / (1 + eps);
  }"),
'done)$

codeC2(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The coefficients C2[l] in the Fourier expansion of B2
  void Geodesic::C2f(real eps, real c[]) {
    real
      eps2 = Math::sq(eps),
      d = eps;
    switch (nC2_) {"),
  for n:0 thru maxpow do (
    print(concat(tab2,"case ",string(n),":")),
    for m:1 thru n do block([q:d*horner(
        subst([eps=sqrt(eps2)],ataylor(C2[m],eps,n)/eps^m)),
      linel:1200],
      if m>1 then print(concat(tab3,"d *= eps;")),
      print(concat(tab3,"c[",string(m),"] = ",string(q),";"))),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nC2_ >= ",string(0),
      " && nC2_ <= ",string(maxpow),", \"Bad value of nC2_\");")),
  print("    }
  }"),
'done)$

codeA3(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The scale factor A3 = mean value of (d/dsigma)I3
  void Geodesic::A3coeff() {
    switch (nA3_) {"),
  for nn:0 thru maxpow do block(
    [q:if nn=0 then 0 else
    jtaylor(subst([n=_n],A3),_n,eps,nn-1),
    linel:1200],
    print(concat(tab2,"case ",string(nn),":")),
    for i : 0 thru nn-1 do
    print(concat(tab3,"_A3x[",i,"] = ",
        string(horner(coeff(q,eps,i))),";")),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nA3_ >= ",string(0),
      " && nA3_ <= ",string(maxpow),", \"Bad value of nA3_\");")),
  print("    }
  }"),
'done)$

codeC3(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The coefficients C3[l] in the Fourier expansion of B3
  void Geodesic::C3coeff() {
    switch (nC3_) {"),
  for nn:0 thru maxpow do block([c],
    print(concat(tab2,"case ",string(nn),":")),
    c:0,
    for m:1 thru nn-1 do block(
      [q:if nn = 0 then 0 else
      jtaylor(subst([n=_n],C3[m]),_n,eps,nn-1),
      linel:1200],
      for j:m thru nn-1 do (
        print(concat(tab3,"_C3x[",c,"] = ",
            string(horner(coeff(q,eps,j))),";")),
        c:c+1)
    ),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nC3_ >= ",string(0),
      " && nC3_ <= ",string(maxpow),", \"Bad value of nC3_\");")),
  print("    }
  }"),
'done)$

codeC4(maxpow):=block([tab2:"    ",tab3:"      "],
  print("  // The coefficients C4[l] in the Fourier expansion of I4
  void Geodesic::C4coeff() {
    switch (nC4_) {"),
  for nn:0 thru maxpow do block([c],
    print(concat(tab2,"case ",string(nn),":")),
    c:0,
    for m:0 thru nn-1 do block(
      [q:jtaylor(subst([n=_n],C4[m]),_n,eps,nn-1),
      linel:1200],
      for j:m thru nn-1 do (
        print(concat(tab3,"_C4x[",c,"] = ",
            string(horner(coeff(q,eps,j))),";")),
        c:c+1)
    ),
    print(concat(tab3,"break;"))),
  print(concat(tab2,"default:")),
  print(concat(tab3,"static_assert(nC4_ >= ",string(0),
      " && nC4_ <= ",string(maxpow),", \"Bad value of nC4_\");")),
  print("    }
  }"),
'done)$

printcode():=(
  print(""),
  print(concat("  // Generated by Maxima on ",timedate())),
  print(""),
  codeA1(maxpow),
  print(""),
  codeC1(maxpow),
  print(""),
  codeC1p(maxpow),
  print(""),
  codeA2(maxpow),
  print(""),
  codeC2(maxpow),
  print(""),
  codeA3(maxpow),
  print(""),
  codeC3(maxpow),
  print(""),
  codeC4(maxpow))$

/* FORMAT FOR DISPLAY */

dispA1(ord):=block(
  [tt:ataylor(A1*(1-eps),eps,ord),ttt,linel:1200],
  for j:2 step 2 thru ord do (ttt:coeff(tt,eps,j),
    print(concat(if j = 2 then "A1 = (1 " else "        ",
        if ttt>0 then "+ " else "- ",
        string(abs(ttt)), " * ", string(eps^j),
        if j=ord or j = ord-1 then ") / (1 - eps);" else ""))))$

dispC1(ord):=for i:1 thru ord do
block([tt:ataylor(C1[i],eps,ord),ttt,linel:1200],
  print(),
  for j:i step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat(
        if j = i then concat("C1[",string(i),"] = ") else "        ",
        if ttt>0 then "+ " else "- ",
        string(abs(ttt)), " * ", string(eps^j),
        if j=ord or j=ord-1 then ";" else ""))))$

dispC1p(ord):=for i:1 thru ord do
block([tt:ataylor(C1p[i],eps,ord),ttt,linel:1200],
  print(),
  for j:i step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat(
        if j = i then concat("C1p[",string(i),"] = ") else "         ",
        if ttt>0 then "+ " else "- ",
        string(abs(ttt)), " * ", string(eps^j),
        if j=ord or j=ord-1 then ";" else ""))))$

dispA2(ord):=block(
  [tt:ataylor(A2*(1+eps),eps,ord),ttt,linel:1200],
  for j:2 step 2 thru ord do (ttt:coeff(tt,eps,j),
    print(concat(if j = 2 then "A2 = (1 " else "        ",
        if ttt>0 then "+ " else "- ",
        string(abs(ttt)), " * ", string(eps^j),
        if j=ord or j = ord-1 then ") / (1 + eps);" else ""))))$

dispC2(ord):=for i:1 thru ord do
block([tt:ataylor(C2[i],eps,ord),ttt,linel:1200],
  print(),
  for j:i step 2 thru ord do (ttt:coeff(tt,eps,j), print(concat(
        if j = i then concat("C2[",string(i),"] = ") else "        ",
        if ttt>0 then "+ " else "- ",
        string(abs(ttt)), " * ", string(eps^j),
        if j=ord or j=ord-1 then ";" else ""))))$

dispA3(ord):=(ord:ord-1,block(
  [tt:jtaylor(A3,n,eps,ord),ttt,t4,linel:1200,s],
  for j:1 thru ord do (ttt:expand(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=1 then "A3 = 1" 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,")"),
      s:concat(s," * ", string(eps^j)),
      print(concat(s,if j = ord then ";" else ""))))))$

dispC3(ord):=(ord:ord-1,for i:1 thru ord do
block([tt:jtaylor(C3[i],eps,n,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("C3[",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,")"),
      s:concat(s," * ", string(eps^j)),
      print(concat(s,if j = ord then ";" else ""))))))$

dispC4(ord):=(ord:ord-1,for i:0 thru ord do
block([tt:jtaylor(C4[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("C4[",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 ""))))))$

dispseries():=(
  print(""),
  print(concat("// Generated by Maxima on ",timedate())),
  print(""),
  dispA1(maxpow),
  print(""),
  dispC1(maxpow),
  print(""),
  dispC1p(maxpow),
  print(""),
  dispA2(maxpow),
  print(""),
  dispC2(maxpow),
  print(""),
  dispA3(maxpow),
  print(""),
  dispC3(maxpow),
  print(""),
  dispC4(maxpow),
  print(""))$

/* CALL THE FUNCTIONS */

/* Timings for computeall(n)
   n   time(s)
   8     8
  10    19
  12    43
  20   571
  30  6671 (111m)

For computeI4(n)
   n   time(s)
   8      1.5
  16       29
  24      215
  30      702 = 11.7 m
  32     1040 = 17.3 m
  40     3851
  48    11521
  64    69960 = 19.4 h

  */
maxpow:8$
computeall()$
/* printcode()$ */
dispseries()$

load("polyprint.mac")$
printgeod():= block([macro:if simplenum then "GEOGRAPHICLIB_GEODESIC_ORDER"
  else "GEOGRAPHICLIB_GEODESICEXACT_ORDER"],
  value1('(A1*(1-eps)-1),'eps,2,0),
  array1('C1,'eps,2,0),
  array1('C1p,'eps,2,0),
  value1('(A2*(1+eps)-1),'eps,2,0),
  array1('C2,eps,2,0),
  value2('A3,'n,'eps,1),
  array2('C3,'n,'eps,1),
  array2('C4,'n,'eps,1))$
printgeod()$
/* Save the values needed for geodesic.mac  This is commented out
here to avoid accidentally overwriting files in a user's directory. */
/* (file:concat("geod",maxpow,".lsp"), save(file, values, arrays))$ */