beginner’s question with Maxim a

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),
   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]
https://lists.sourceforge.net/lists/listinfo/maxima-discuss