Re: beginner’s question with Maxi ma

Maxima doesn't  automatically express integrate(1/(1-e^2 * sin(x)^2)^(3/2),x) in terms of elliptic_e, but  hand work gives:

(%o12)  F(a,e,phi) := a*(elliptic_e(phi, e^2) - e^2 * cos(phi) * sin(phi)/sqrt(1-e^2*sin(phi)^2))$

(%i13) plot2d(F(1,9/10,z),[z,0,%pi]);
(%o13) false
<nice graph deleted>

For e near one and phi near zero,  evaluating `elliptic_e(phi, e^2) - e^2 ...` suffers from subtractive cancellation.  I don't know of a cure for that, but maybe in the context of the meridian arclength for the earth, this isn't a problem.

--Barton




________________________________
From: Stavros Macrakis <[email protected]>
Sent: Friday, August 21, 2026 10:31 PM
To: [email protected] <[email protected]>
Cc: <[email protected]> <[email protected]>
Subject: Re: [Maxima-discuss] beginner’s question with Maxima

Caution: Non-NU Email

Maybe the problem is that you declare

declare("-", alphabetic)$

But then later write -n, which now denotes the symbol named "-n" rather than the negation of n.



On Fri, Aug 21, 2026, 20:19 <[email protected]<mailto:[email protected]>> wrote:
Hello, all,

I’m a new subscriber to this mailing list and an absolute beginner with Maxima. I’m trying to use Maxima to calculate a meridian arc length on the WGS 84 ellipsoid, and I seem to be missing something essential in my understanding of Maxima. I’ve put my attempt so far into a .mac file, which is presented below. I’d appreciate receiving your input on where my knowledge has fallen short. Here’s my code in its current state:

/*
   meridian_arc_length.mac

   Calculate meridian arc lengths to zeptometer (zm) precision,
   but display them at attometer (am) precision.
*/
fpprec: 21$
fpprintprec: fpprec - 3$
/*
   Define the ⌊ ⌋ matchfix operator as an alias for floor().
*/
matchfix("⌊", "⌋")$
⌊n⌋ := floor(n)$
/*
   Define the ², ³, and ⁴ postfix operators as aliases for ^2, ^3, and ^4
   respectively, with the same precedence as ^.
*/
exponent_precedence: 140$
postfix("²", exponent_precedence, any, any)$
"²"(n) := n^2$
postfix("³", exponent_precedence, any, any)$
"³"(n) := n^3$
postfix("⁴", exponent_precedence, any, any)$
"⁴"(n) := n^4$
/*
   Define the °, ′, ″, ‴, and ⁗ postfix operators for converting arcdegrees,
   arcminutes, arcseconds, arcthirds, and arcfourths respectively to radians.
   The precedence of these operators is just below those of the * and /
   operators, so that e.g. “1/2°” is “(1/2)°” rather than “1/(2°)”.
*/
arc_precedence: 110$
π: bfloat(%pi)$
postfix("°", arc_precedence, any, any)$
"°"(ad) := ad * π/180$
postfix("′", arc_precedence, any, any)$
"′"(am) := am * π/(180*60)$
postfix("″", arc_precedence, any, any)$
"″"(as) := as * π/(180*60²)$
postfix("‴", arc_precedence, any, any)$
"‴"(aþ) := aþ * π/(180*60³)$
postfix("⁗", arc_precedence, any, any)$
"⁗"(af) := af * π/(180*60⁴)$
/*
   Define the ∏() and ∑() macros for product() and sum() respectively,
   so that their expressions are specified last rather than first.
*/
declare("∏", alphabetic)$
∏(i, i_min, i_max, expr) ::= product(expr, i, i_min, i_max)$
declare("∑", alphabetic)$
∑(i, i_min, i_max, expr) ::= sum(expr, i, i_min, i_max)$
/*
   Load the ezunits package so that the semi-major axis
   is specifiable in meters.
*/
load("ezunits")$
/*
   Use the WGS 84 ellipsoid definitions for Earth’s semi-major axis
   and third flattening.
*/
declare("-", alphabetic)$
semi-major_axis: 6378137 ` m;
first_flattening: 1/bfloat(298.257223563)$
third_flattening: first_flattening/(2 - first_flattening);
/*
   Find the meridian arc length from the equator to 45° (North or South).
   The meridian arc length should be a bit less than semi-major_axis/2
   on the WGS 84 ellipsoid.
*/
latitude: 45°;
/*
   Adapt S(φ), which is Equation (15) in the paper
   “A General Formula for Calculating Meridian Arc Length […]”
   (https://www.gsi.go.jp/common/000062452.pdf<https://urldefense.com/v3/__https://www.gsi.go.jp/common/000062452.pdf__;!!PvXuogZ4sRB2p-tU!AIOl_j0NVPHw_aOldXuGlGn_FsaxN3yf45AK9vzPwwHic2qe470vmcHMdP0vSIMeEcHzoEm6Q-kkwO8$>),
   to be a memoizing function in Maxima,
   except with a, n, and j_maximum specified as additional parameters,
   so that the function can be used with other ellipsoids also.

   QUESTION: Have I correctly adapted the paper’s Equation 15
   for Maxima? (Equation 15 has the equivalent of inf for j_maximum,
   and I wanted to ensure that low j_maximum values were working
   before trying inf with S[].)
*/
S[a, n, φ, j_maximum] := \
 a * (1 - n)² * (1 + n) \
   * ∑(j, 0, j_maximum, ∏(k, 1, j, (-n/(2*k) - n)²)) \
   * (φ + ∑(l, 1, 2*j, (sin(2*l*φ)/l)) \
          * ∏(m, 1, l, \
              (-n/(2*j + 2*(-1)^m*⌊m/2⌋) - n)^((-1)^m)));
for j_max: 0 thru 5 do S[semi-major_axis, third_flattening, latitude, j_max];
/*
   QUESTION: For each call of S[] with j_max from 0 through 5 inclusive,
   all of j through n are not completely resolved in the returned values.

   For example, with j_max = 0,

     (%i1) S[semi-major_axis, third_flattening, 45°, 0];
                                                                     2 j
                                       │-n                          │
     (%o1) (6.36740874757294245b6 ` m) │─── - 1.67922038638370455b-3│
                                       │2 k                         │
      ((2 j sin(1.57079632679489662b0 l)
                                                        m
                 -n                                   -1  l
      ((──────────────────── - 1.67922038638370455b-3)   ) )/l
            m       m
        2 -1  floor(─) + 2 j
                    2
      + 7.8539816339744831b-1)
     (%i2) _

   What have I failed to do to ensure the resolution of j through n?
*/


_______________________________________________
Maxima-discuss mailing list
[email protected]<mailto:[email protected]>
https://lists.sourceforge.net/lists/listinfo/maxima-discuss<https://urldefense.com/v3/__https://lists.sourceforge.net/lists/listinfo/maxima-discuss__;!!PvXuogZ4sRB2p-tU!AIOl_j0NVPHw_aOldXuGlGn_FsaxN3yf45AK9vzPwwHic2qe470vmcHMdP0vSIMeEcHzoEm6WTfY46k$>

_______________________________________________
Maxima-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/maxima-discuss