[doc] Move Maxima scripts for geodesics to doc/other/maxima/geod.mac

doc/other/maxima/geod.mac

@@ -0,0 +1,232 @@
+/*
+Compute the series expansions for the ellipsoidal geodesic problem.
+
+Copyright (c) Charles Karney (2009-2015) <charles@karney.com> 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
+
+The code below contains minor modifications to conform with
+Boost Geometry style guidelines.
+
+To run the code, start Maxima and enter
+
+  load("geod.mac")$
+*/
+
+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)))))$
+
+computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
+  sintegrand:sqrt(1+k2*sin(sigma)^2),
+  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)$
+
+codeA1(maxpow):=block([tab2:"    ",tab3:"        "],
+print("// The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
+static inline CT evaluate_series_A1(CT eps) {
+    CT eps2 = math::sqr(eps);
+    CT t;
+    switch (SeriesOrder/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("    }
+    return (t + eps) / (1 - eps);
+}"),
+'done)$
+
+computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
+  sintegrand:1/sqrt(1+k2*sin(sigma)^2),
+  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)$
+
+codeA2(maxpow):=block([tab2:"    ",tab3:"        "],
+print("// The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
+CT evaluate_series_A2(CT const& eps)
+{
+    CT const eps2 = math::sqr(eps);
+    CT t;
+    switch (SeriesOrder/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("    }
+    return (t - eps) / (1 + eps);
+}"),
+'done)$
+
+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)$
+
+codeA3(maxpow):=block([tab2:"    ",tab3:"        "],
+print("// The scale factor A3 = mean value of (d/dsigma)I3
+static inline void evaluate_series_A3(CT const& n, CT c[])
+{
+    switch (SeriesOrder) {"),
+  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,"c[",i,"] = ",
+        string(horner(coeff(q,eps,i))),";")),
+    print(concat(tab3,"break;"))),
+  print("    }
+}"),
+'done)$
+
+codeC1(maxpow):=block([tab2:"    ",tab3:"        "],
+  print("// The coefficients C1[l] in the Fourier expansion of B1
+static inline evaluate_coeffs_C1(CT eps, CT c[]) {
+    CT eps2 = math::sqr(eps);
+    CT d = eps;
+    switch (SeriesOrder) {"),
+  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("    }
+}"),
+'done)$
+
+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)$
+
+codeC1p(maxpow):=block([tab2:"    ",tab3:"        "],
+  print("// The coefficients C1p[l] in the Fourier expansion of B1p
+static inline evaluate_coeffs_C1p(CT eps, CT c[])
+{
+    CT const eps2 = math::sqr(eps);
+    CT d = eps;
+    switch (SeriesOrder) {"),
+  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("    }
+}"),
+'done)$
+
+codeC2(maxpow):=block([tab2:"    ",tab3:"        "],
+print("// The coefficients C2[l] in the Fourier expansion of B2
+static inline void evaluate_coeffs_C2(CT const& eps, CT c[])
+{
+    CT const eps2 = math::sqr(eps);
+    CT d = eps;
+    switch (SeriesOrder) {"),
+  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("    }
+}"),
+'done)$
+
+codeC3(maxpow):=block([tab2:"    ",tab3:"        "],
+print("// The coefficients C3[l] in the Fourier expansion of B3
+static inline void evaluate_coeffs_C3(CT const& n, CT c[])
+{
+    const CT n2 = math::sqr(n);
+    switch (SeriesOrder) {"),
+  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,"c[",c,"] = ",
+            string(horner(coeff(q,eps,j))),";")),
+        c:c+1)
+    ),
+    print(concat(tab3,"break;"))),
+  print("    }
+}"),
+'done)$
+
+maxpow:8$
+
+computeI1(maxpow)$
+computeI2(maxpow)$
+computeI3(maxpow)$
+
+revertI1(maxpow)$
+codeA1(maxpow)$
+codeA2(maxpow)$
+codeA3(maxpow)$
+
+codeC1(maxpow)$
+codeC2(maxpow)$
+codeC3(maxpow)$
+
+codeC1p(maxpow)$

include/boost/geometry/util/series_expansion.hpp

@@ -30,51 +30,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion above is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_series_A1 below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        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)))))$
-
-        computeI1(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
-          sintegrand:sqrt(1+k2*sin(sigma)^2),
-          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)$
-
-        codeA1(maxpow):=block([tab2:"    ",tab3:"        "],
-        print("// The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
-        static inline CT evaluate_series_A1(CT eps) {
-            CT eps2 = math::sqr(eps);
-            CT t;
-            switch (SeriesOrder/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("    }
-            return (t + eps) / (1 - eps);
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI1(maxpow)$
-        codeA1(maxpow)$
-        // Maxima script end
 
      To replace each number x by CT(x) the following
      script can be used:
@@ -123,52 +80,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion above is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_series_A2 below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        computeI2(maxpow):=block([sintegrand,sintegrandexp,s,sigma,tau1,k2,eps],
-          sintegrand:1/sqrt(1+k2*sin(sigma)^2),
-          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)$
-
-        codeA2(maxpow):=block([tab2:"    ",tab3:"        "],
-        print("// The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
-        CT evaluate_series_A2(CT const& eps)
-        {
-            CT const eps2 = math::sqr(eps);
-            CT t;
-            switch (SeriesOrder/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("    }
-            return (t - eps) / (1 + eps);
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI2(maxpow)$
-        codeA2(maxpow)$
-        // Maxima script end
-
-     To replace each number x by CT(x) the following
-     script can be used:
-       sed -e 's/[0-9]\+/CT(&)/g; s/\[CT/\[/g; s/)\]/\]/g;
-               s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;'
     */
     template <typename CT, std::size_t SeriesOrder>
     CT evaluate_series_A2(CT const& eps)
@@ -212,53 +125,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion above is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_coeffs_A3 below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        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)$
-
-        codeA3(maxpow):=block([tab2:"    ",tab3:"        "],
-        print("// The scale factor A3 = mean value of (d/dsigma)I3
-        static inline void evaluate_series_A3(CT const& n, CT c[])
-        {
-            switch (SeriesOrder) {"),
-          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,"c[",i,"] = ",
-                string(horner(coeff(q,eps,i))),";")),
-            print(concat(tab3,"break;"))),
-          print("    }
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI3(maxpow)$
-        codeA3(maxpow)$
-        // Maxima script end
-
-     To replace each number x by CT(x) the following
-     script can be used:
-       sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
-               s/case\sCT(/case /g; s/):/:/g'
     */
     // TODO: this produces different results that geographiclib
     template <typename CT, std::size_t SeriesOrder>
@@ -327,38 +195,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion below is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_coeffs_C1 below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        codeC1(maxpow):=block([tab2:"    ",tab3:"        "],
-          print("// The coefficients C1[l] in the Fourier expansion of B1
-        static inline evaluate_coeffs_C1(CT eps, CT c[]) {
-            CT eps2 = math::sqr(eps);
-            CT d = eps;
-            switch (SeriesOrder) {"),
-          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("    }
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI1(maxpow)$
-        codeC1(maxpow)$
-        // Maxima script end
-
-     To replace each number x by CT(x) the following
-     script can be used:
-       sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
-               s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;
-               s/eps(CT(2))/eps2/g;'
     */
     template <typename CT, std::size_t SeriesOrder>
     static inline void evaluate_coeffs_C1(CT eps, CT c[])
@@ -456,51 +294,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion below is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_coeffs_C1p below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        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)$
-
-        codeC1p(maxpow):=block([tab2:"    ",tab3:"        "],
-          print("// The coefficients C1p[l] in the Fourier expansion of B1p
-        static inline evaluate_coeffs_C1p(CT eps, CT c[])
-        {
-            CT const eps2 = math::sqr(eps);
-            CT d = eps;
-            switch (SeriesOrder) {"),
-          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("    }
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI1(maxpow)$
-        revertI1(maxpow)$
-        codeC1p(maxpow)$
-        // Maxima script end
-
-     To replace each number x by CT(x) the following
-     script can be used:
-       sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
-               s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;
-               s/eps(CT(2))/eps2/g;'
     */
     template <typename CT, std::size_t SeriesOrder>
     static inline void evaluate_coeffs_C1p(CT eps, CT c[])
@@ -598,39 +393,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion below is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_coeffs_C2 below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        codeC2(maxpow):=block([tab2:"    ",tab3:"        "],
-        print("// The coefficients C2[l] in the Fourier expansion of B2
-        static inline void evaluate_coeffs_C2(CT const& eps, CT c[])
-        {
-            CT const eps2 = math::sqr(eps);
-            CT d = eps;
-            switch (SeriesOrder) {"),
-          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("    }
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI2(maxpow)$
-        codeC2(maxpow)$
-        // Maxima script end
-
-     To replace each number x by CT(x) the following
-     script can be used:
-       sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
-               s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;
-               s/eps(CT(2))/eps2/g;'
     */
     template <typename CT, std::size_t SeriesOrder>
     static inline void evaluate_coeffs_C2(CT const& eps, CT c[])
@@ -728,43 +492,8 @@ namespace boost { namespace geometry { namespace series_expansion {
 
      The expansion below is performed in Maxima, a Computer Algebra System.
      The C++ code (that yields the function evaluate_coeffs_C3 below) is
-     generated by the following Maxima script and is based on script:
+     generated by the following Maxima script:
-     http://geographiclib.sourceforge.net/html/geod.mac
+     geometry/doc/other/maxima/geod.mac
-
-        // Maxima script begin
-        codeC3(maxpow):=block([tab2:"    ",tab3:"        "],
-        print("// The coefficients C3[l] in the Fourier expansion of B3
-        static inline void evaluate_coeffs_C3(CT const& n, CT c[])
-        {
-            const CT n2 = math::sqr(n);
-            switch (SeriesOrder) {"),
-          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,"c[",c,"] = ",
-                    string(horner(coeff(q,eps,j))),";")),
-                c:c+1)
-            ),
-            print(concat(tab3,"break;"))),
-          print("    }
-        }"),
-        'done)$
-
-        maxpow:8$
-        computeI3(maxpow)$
-        codeC3(maxpow)$
-        // Maxima script end
-
-     To replace each number x by CT(x) the following
-     script can be used:
-       sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
-               s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;
-               s/eps(CT(2))/eps2/g;'
     */
     template <typename CT, int SeriesOrder>
     void evaluate_coeffs_C3x(CT const& n, CT c[], const CT coeff[])