ELLIPSOID Module


 Overview
 

-A few programs that deal with ellipsoids, hyperellipsoids and elliptic or ellipsoidal functions.

-For the Earth, the geodesic distance is calculated with Vincenty's formula by "DGT2"
-A triaxial ellipsoid model of the Earth is also used with "DGT3" & "LOX3" for the geodesics & loxodromes,
 and with "GRV3" for the Earth gravity.

-The longitudes are positive East.
 
 
 
XROM  Function  Desciption
 21,00
 21,01
 21,02
 21,03
 21,04
 21,05
 21,06
 21,07
 21,08
 21,09
 21,10
 21,11
 21,12
 21,13
 21,14
 21,15
 21,16
 21,17
 21,18
 21,19
 21,20
 21,21
 21,22
 21,23
 21,24		 -ELLIPSOID
 abcL
 GX7
 GW7
 "SNX"
 "1/G"
 X#Y??
 "LB-XYZ"
 "XYZ-LB"
 "UV-XYZ"
 "XYZ-UV"
 "DPE"
 "ELKGM"
 "DGT2"
 "DGT3"
 "LOX3"
 "GRV3"
 "SAE"
 "SAHE"
 "QHS"
 "EWF"
 "EWF2"
 "RF"
 "RG"
 "RJ"		  Section Header
 Semi-axes of the Ellipsoidal Earth
 Abscissas of Gauss-Legendre 7-point formula
 Weights of Gauss-Legendre 7-point formula
 Jacobian Elliptic Function sn ( x | k^2 ) k < 1
 Reciprocal of Gamma Function
 Skips the next line if X almost equal to Y
 Conversion Geodetic Coord. -> Cartesian Coord.
 Conversion Cartesian Coord. -> Geodetic Coord.
 Conversion Ellipsoidal Coord. -> Cartesian Coord.
 Conversion Cartesian Coord. -> Ellipsoidal Coord.
 Distance from a Point to an Ellipsoid
 Gaussian & Mean Curvature on an Ellipsoid
 Geodesic Distance ( Vincenty Formula )
 Geodesic Distance on a Triaxial Ellipsoid
 Loxodrome on a Triaxial Ellipsoid
 Gravity on a Triaxial Ellipsoid
 Surface Area of an Ellipsoid
 Surface Area of a Hyperellipsoid
 Quadratic Hypersurface
 Ellipsoidal Wave Differential Equation ( RK6 )
 Ellipsoidal Wave Function ( Series )
 Carlson Elliptic Integral of the 1st kind
 Carlson Symmetric Elliptic Integral of the 2nd kind
 Carlson Elliptic Integral of the 3rd kind

 

Semi-axes of the Ellipsoidal Earth
 

 abcL is an M-Code routine that places in registers X Y Z T the semi-axes a b c of the Earth and the longitude of the equatorial major axis
 
 
      STACK        INPUTS      OUTPUTS
           T             /      -14.92911
           Z             /    6356.751868
           Y             /    6378.101946
           X             /    6378.171274

 I've chosen values given in reference [1] for the mean equatorial radius R = 6378.13661 km  and for the Earth flattening  f = 1 / 298.256421
 The difference a - b is about 69 km, so the decimals in a & b are somewhat arbitrary...

 a b c are expressed in kilometers
 The longitude L = -14.92911  is expressed in decimal degrees.     ( longitudes > 0 East )

 The previous content of the stack is overwritten.
 

Abscissas & Weights of Gauss-Legendre 7-point formula
 

-Two programs need to evaluate integrals and GX7 & GW7 returns the coefficients:

   1  XEQ "GX7"  ->  0.4058451514             the 1st abscissa is zero, so I didn't coded it
   2  XEQ "GX7"  ->  0.7415311856
   3  XEQ "GX7"  ->  0.9491079123

   0  XEQ "GW7"  ->  0.4179591837
   1  XEQ "GW7"  ->  0.3818300505
   2  XEQ "GW7"  ->  0.2797053915
   3  XEQ "GW7"  ->  0.1294849662
 

Jacobian Elliptic Function Sn ( x | k^2 ) k < 1
 

-The ellipsoidal wave function W(x) involves Sn(x)
-So I had to write a specific program to calculate this function.
 
 
      STACK        INPUTS      OUTPUTS
           Y           k^2             /
           X             x    Sn ( x | k^2 )

 
Example:

    0.8   ENTER^
    0.7   XEQ "SNX"  >>>>  0.612365484                ---Execution time = 7s---

Notes:

-Though it will not work if k = 1, it will with k2 = .9999999999

    .9999999999  ENTER^
             0.7            R/S        >>>>     0.604367777   which is actually  Tanh 0.7

-With k = 0 and x = 0.7  again,  it yields  0.644217687 = Cos 0.7

-Registers R19 thru R27 are used
 

Reciprocal of Gamma Function
 

-The program "SAHE" which calculates the surface area of a hyperellipsoid calls "1/G" as a subroutine
-"1/G" employs a continued fraction to evaluate 1 / Gamma(x)
 
 
 
      STACK        INPUTS      OUTPUTS
           X             x    1 / Gamma(x)

 
Example:

   PI  XEQ "1/G"  >>>>   0.437055717                                ---Execution time = 3s---
   -3      R/S         >>>>    0                                                   ---Execution time = 4s---
-4.16    R/S         >>>>   -4.696326589                                ---Execution time = 5s---
 

X#Y??  Skips the next line in a program if X almost equal to Y
 

-An M-Code routine to avoid infinite loops when the last results indefinitely alternates between 2 closed values.
-No check for alpha data
 

Conversion Geodetic Coordinates <> Cartesian Coordinates
 

-The formulae between the geodetic longitude l &  latitude b and the cartesian coordinates x , y , z
  of a point on a triaxial ellipsoid are given in references [2] & [3]

    x = w Cos b Cos l
    y = w ( 1 - e2 ) Cos b Sin l                       In fact   l + 14°92911  instead of l
    z = w ( 1 - e' 2 ) Sin b

  where    e2 = 1 - b2/a2    ,    e' 2 = 1 - c2/a2    ,     w = a / ( 1 - e' 2 Sin2b - e2 Cos2 b Sin2l ) 1/2
 

>>> Registers  R00 R01 R02 R03 must be initialized before executing "LB-XYZ" & "XYZ-LB"

      R00 = Longitude of the major equatorial axis
      R01 = a
      R02 = b
      R03 = c
 
 
      STACK        INPUTS      OUTPUTS
           Z             /             z
           Y       l ( ° ' " )             y
           X       b ( ° ' " )             x

 
Example:     The Observatoire de Paris l = 2°20'13"8 E ,  b = 48°50'11"2 N

     XEQ "abcL"  STO 01  RDN  STO 02  RDN  STO 03  RDN  STO 00

    2.20138   ENTER^
   48.50112  XEQ "LB-XYZ" >>>>  x = 4016.633561 km            ---Execution time = 4s---
                                              RDN   y = 1248.421908 km
                                              RDN   z = 4778.596641 km
 
 
      STACK        INPUTS      OUTPUTS
           Z             z             /
           Y             y      b ( ° ' " )
           X             x      l ( ° ' " )

 
Example:  the reverse transformation. Assuming x y z above are in X Y Z registers

                 R/S        returns            2.201380001                          ---Execution time = 2.5s---
                               X<>Y           48.50112000

Note:

-Of course, you can use your own values for a  b  c  L
 

Conversion Ellipsoidal Coordinates <> Cartesian Coordinates
 

-The conversion between Jacobi parameters ( u , v ) and cartesian coordinates is given by the formulae ( here we assume a > b > c and L = 0 ):
 

    x =  a / ( a - c )1/2 Cos u  ( a2 - b2 Sin2 v - c2 Cos2 v )1/2
    y =  b  Sin u  Cos v
    z =  c / ( a - c )1/2 Sin v  ( a2 Sin2 u + b2 Cos2 u - c2 )1/2
 
 
 
       STACK      INPUTS1    OUTPUTS1       INPUTS2     OUTPUTS2
           Z            /             z              z              v
           Y            u             y              y              v
           X            v             x              x              u

   u & v  are expressed in decimal degrees

Example:   The ellipsoid  x2 / 64 + y2 / 49 + z2 / 36 = 1              u = +120° ,   v = -70°

   8  STO 01
   7  STO 02
   6  STO 03

-Set the HP-41 in DEG mode:  XEQ "DEG"

   120   ENTER^
   -70   XEQ "UV-XYZ"  >>>>   x = -3.072524427                      ---Execution time = 3s---
                                       RDN    y =  2.073386929
                                       RDN    z = -5.247034535

-With these values in registers  X , Y , Z  ( RDN RDN )

   XEQ "XYZ-UV"  ( or simply R/S )  gives again  u = 120  RDN  v = -70 , in 4.6 seconds, with small roudoff-errors in the last decimals.

Notes:

-These programs work in all angular modes.
-"XYZ-UV" does not check that ( x , y , z ) is on the ellipsoid.
-R00 is unused
 

Distance from a Point P to an Ellipsoid
 

 "DPE" calculates this value by using the parametric coordinates of the point P' on the ellipsoid (E) such that PP' is perpendicular to (E)

    x = a Cos U Cos V
    y = b Sin U Cos V
    z = c Sin V

 This leads to a 2x2 non-linear system which is solved by Newton's method.
 The successive U-values are displayed.
 
 
      STACK        INPUTS      OUTPUTS
           T             /            z
           Z            Z            y
           Y            Y            x
           X            X            D

     Where D = distance PP'   P(X,Y,Z)  P'(x,y,z)

Example1:   (E):   x2 / 49 + y2 / 36 + z2 / 25 = 1      P(10,11,12)

    7  STO 01
    6  STO 02
    5  STO 03

    12  ENTER^
    11  ENTER^
    10  XEQ "DPE"   >>>>    D = 13.23891299                      ---Execution time = 32s---
                                RDN     x =  3.896191634
                                RDN     y =  3.511764224
                                RDN     z =  2.948002121
 

Example2:      P(0.1,0.2,0.3)

     0.3   ENTER^
     0.2   ENTER^
     0.1      R/S    >>>>        D =  4.691015000                        ---Execution time = 53s---
                              RDN     x =  0.192100176
                              RDN     y =  0.575653468
                              RDN     z =  4.975042648

Notes:

 "DPE" works well if P is outside the ellipsoid
 Otherwise, the results may be doubtful, for instance, with P(0,0,0) it yields D = 7  which is obviously wrong  ( D = 5 )

-You can store other initial values for U & V in R11 & R12  and press  XEQ 01  again
 

Gaussian & Mean Curvature on an Ellipsoid
 

 "ELKGM" returns the Gaussian Curvature KG in X-register and the mean curvature Km on a point P(x,y,z) on an ellipsoid x2 / a2 + y2 / b2 + z2 / c2 = 1

Formulae:

  KG =  1 / [ a2 b2 c2 ( x2 / a4 + y2 / b4 + z2 / c4 )2 ]

  Km =  - (1/2) Div n    where  n = ( x / a2 , y / b2 , z / c2 ) / ( x2 / a4 + y2 / b4 + z2 / c4 ) 1/2
 

  R00 is unused but a , b , c must be stored in R01  R02  R03  respectively before executing this routine.
 
 
      STACK        INPUTS      OUTPUTS
           Z            z           KG
           Y            y           Km
           X            x           KG

 
Example:   Again with (E):   (E):   x2 / 49 + y2 / 36 + z2 / 25 = 1         P ( 3.896191634 , 3.511764224 , 2.948002121 )

    7  STO 01
    6  STO 02
    5  STO 03

    2.948002121  ENTER^
    3.511764224  ENTER^
    3.896191634  XEQ "ELKGM"   >>>>  KG =  0.025631779                        ---Execution time = 3s---
                                                        RDN  Km = -0.163109816

Note:

 The minus sign for Km could be replaced by a "+" ...
 

Geodesic Distance ( Vincenty Formula )
 

 "DGT2" uses an ellipsoid of revolution to evaluate the geodesic distance between 2 points on the Earth
  The forward and backward azimuths are also returned in Y and Z-registers.

  The method is iterative and it doesn't work for nearly antipodal points.
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )     Az2 ( ° ' " )
           Z       b1 ( ° ' " )     Az2 ( ° ' " )
           Y       L2 ( ° ' " )     Az1 ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

 where   Az1  =  forward azimuth  ( from P1 to P2 )
   and     Az2  =  back azimuth  ( from P2 to P1 )

Example:

   -77.0356    ENTER^
    38.55172   ENTER^
      2.20138   ENTER^
   48.50112  XEQ "DGT2"  >>>>    D   =  6181.621427 km                      ---Execution time = 46s---
                                           RDN    Az1 =  51°47'36"8132
                                           RDN    Az2 = -68°09'58"9654

-We also have  R06 = µ = azimuth at the equator = 37.74571892°
 

Geodesic Distance on a Triaxial Ellipsoid
 

 "DGT3" uses an approximate method to evaluate the geodesic distance between 2 points on the Earth
  The error is about a few centimeters for long geodesics, except for nearly antipodal points ( do not use this program in this case )
 

-Let a point P ( on the ellipsoid ) with ( OM , OP' ) = µ , ( OP' , OP ) = h    ,    P' = orthogonal projection of P on the plane (OMN)

    OP = x u + y v + z w       where  u = OM / || OM ||   ,  v = ON / || ON ||  ,  w = u x v          ( x = cross-product )

-We have, with R = OP

    x = R Cos h Sin ( W - µ ) / Sin W
    y = R Cos h Sin µ / Sin W
    z = R Sin h / Sin W

-Writing that P(X,Y,Z)  is on the ellipsoid  X2/a2 + Y2/b2 + Z2/c2 = 1,  we find that

  R / Sin W = 1 / ( A2 + B2 + C2 )1/2   with

      A = [ x1 Cos h Sin ( W - µ ) + x2 Cos h Sin µ + ( y1z2 - y2z1 ) Sin h ] / a
      B = [ y1 Cos h Sin ( W - µ ) + y2 Cos h Sin µ + ( z1x2 - z2x1 ) Sin h ] / b                       where  u ( x1 , y1 , z1 )  &  v ( x2 , y2 , z2 )
      C = [ z1 Cos h Sin ( W - µ ) + z2 Cos h Sin µ + ( x1y2 - x2y1 ) Sin h ] / c
 

-We search an approximation of the function  h(µ) under the form:

   h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) + h3 Sin ( 540° µ / W )                 ( the first terms of a Fourier series )

-And the corresponding integral  s = §0W [ R2 Cos2 h  + R2 ( dh/dµ )2 + ( dR/dµ )2  ] 1/2 dµ

      =  ( Sin W )    §0W [ ( Cos2 h ) / T + (1/T) ( dh/dµ )2 + ( dr/dµ )2  ] 1/2 dµ    is evaluated by a Gaussian quadrature ( Gauss-Legendre 7-point formula )

  where   T = ( A2 + B2 + C2 )   with the same notations as in paragraph above:

    r = 1/sqrt(T)                  dr/dµ = r/µ + ( r/h ) ( dh/dµ )

    r = 1 / ( A2 + B2 + C2 )1/2   with

      A = [ x1 Cos h Sin ( W - µ ) + x2 Cos h Sin µ + ( y1z2 - y2z1 ) Sin h ] / a
      B = [ y1 Cos h Sin ( W - µ ) + y2 Cos h Sin µ + ( z1x2 - z2x1 ) Sin h ] / b                       where  u ( x1 , y1 , z1 )  &  v ( x2 , y2 , z2 )
      C = [ z1 Cos h Sin ( W - µ ) + z2 Cos h Sin µ + ( x1y2 - x2y1 ) Sin h ] / c

  the partial derivatives

  r/µ = ( -1/sqrt(T3) ) [ (A/a) { -x1 Cos h Cos ( W - µ ) + x2 Cos h Cos µ }
                                   + (B/b) { -y1 Cos h Cos ( W - µ ) + y2 Cos h Cos µ }
                                   + (C/c) { -z1 Cos h Cos ( W - µ ) + z2 Cos h Cos µ } ]
  and

  r/h = ( -1/sqrt(T3) ) [ (A/a) { -x1 Sin h Sin ( W - µ ) - x2 Sin h Sin µ + ( y1z2 - y2z1 ) Cos h }
                                    + (B/b) { -y1 Sin h Sin ( W - µ ) - y2 Sin h Sin µ + ( z1x2 - z2x1 ) Cos h }
                                    + (C/c) { -z1 Sin h Sin ( W - µ ) - z2 Sin h Sin µ + ( x1y2 - x2y1 ) Cos h } ]
 

-With  H = W / 1000  , we calculate the following lengths

    A = f(0,0,0)
    B = f(H,0,0)
    C = f(-H,0,0)
    D = f(0,H,0)                     f = length s corresponding to h1 = -H , 0 , +H ;  h2 = -H , 0 , +H ;  h3 = -H , 0 , +H
    E = f(0,-H,0)
    F = f(0,0,H)
    G = f(0,0,-H)
 

-Then, the extremum is estimated by the formula:

   s ~ A - (B-C)2 / [ 8.(B+C-2.A) ] - (D-E)2 / [ 8.(D+E-2.A) ] - (F-G)2 / [ 8.(F+G-2.A) ]
 
 

Flags:  F00-F01-F02-F04

   CF 00  abcL is called to store a b c L in registers R01-R02-R03-R00
   SF 00:  you have to initialize R00-R01-R02-R03

   CF 01  CF 02  3 terms in the Fourier series:
   CF 01  SF 02   2 terms in the Fourier series:
   SF 01               1 term in the Fourier series:

   SF 04   geocentric longitudes and latitudes
   CF 04  geodetic longitudes and latitudes
 
 
      STACK        INPUTS    OUTPUTS1    OUTPUTS2    OUTPUTS3
           T       l1 ( ° ' " )       D0 ( km )       D0 ( km )       D0 ( km )
           Z       b1 ( ° ' " )       D1 ( km )       D0 ( km )       D0 ( km )
           Y       l2 ( ° ' " )       D2 ( km )       D1 ( km )       D0 ( km )
           X       b2 ( ° ' " )       D3 ( km )       D2 ( km )       D1 ( km )

 
Example1:      Geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

   l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
   l2 = 18°28'35"7 E  ,  b2 = 33°56'02"5 S

    CF 01  CF 02  CF 04

  -77.0356    ENTER^
   38.55172   ENTER^
   18.28357   ENTER^
  -33.56025  XEQ "DGT3"  >>>>   D3 =  12709.47906 km = R33                        ---Execution time = 5m46s---
                                            RDN    D2 = 12709.47906 km
                                            RDN    D1 = 12709.48047 km
                                            RDN    D0 = 12709.48074 km = R21

Example2:    With ( -40° , -40° ) ( +120° , +50° ) it yields:

  •  If  CF 01  CF 02  CF 04   we get:

                     D3 =  18095.49323 km
         RDN    D2 = 18095.49323 km
         RDN    D1 = 18095.49858 km
         RDN    D0 = 18095.55363 km

  •  If  CF 01  CF 02  SF 04   we get:

                     D3 =  18099.69052 km
         RDN    D2 = 18099.69052 km
         RDN    D1 = 18099.69587 km
         RDN    D0 = 18099.75037 km

Example3:     On the ellipsoid    x2 / 41 + y2 / 37 + z2 / 35 = 1              the geocentric coordinates of the 2 points are:

    l1 =  9°17'35"55660                 l2 = 63°20'07"3683
    b1 = -8°11'55"79344                b2 = 57°46'42"9935

    0   STO 00

   41  SQRT  STO 01
   37  SQRT  STO 02
   35  SQRT  STO 03

    SF 00   CF 01  CF 02  SF 04

    9.173555660  ENTER^
  -8.115579344  ENTER^
   63.20073683  ENTER^
   57.46429935  XEQ "DGT3"  >>>>   D3 = 8.594823580
                                                 RDN    D2 = 8.594824326
                                                 RDN    D1 = 8.594849588
                                                 RDN    D0 = 8.594880192

-The exact value: about  8.594822580 whence an error of  E-6 with D3

Notes:

-On the Earth, the maximum error given by D0 is of the same order as Andoyer's first order formula for an ellipsoid of revolution i-e about 60 meters.

-For the Earth, the difference between D2 & D3 is often negligible
-So, you can set F02 without a great loss of precision.

    SF 02  ->  Exectution time = 4m05s
    SF 01  ->  Execution time =  2m29s
 

Loxodrome on a Triaxial Ellipsoid
 

-Let  r = semi-major axis of a parallel and  rm the radius of curvature of the meridian at the point P(L,B)  L = geodetic longitude B = geodetis latitude.

     r2 = b2 + ( a2 - b2 ) / [ 1 + (b2/a2) Tan2 L ]
   rm2 = c2 / { r [ 1 - ( 1 - b2/a2 ) Sin2 B ]3/2 }

-Let  a' & b'  the semi-axes of this meridian

    a' = ( a Cos B ) / [ 1 - ( 1 - c2/a2 ) Sin2 B ]1/2
    b' = ( b Cos B ) / [ 1 - ( 1 - c2/b2 ) Sin2 B ]1/2

-Let r'm = the radius of curvature of the parallel

    r'm2 = b'2 / { a' [ 1 - ( 1 - b'2/a'2 ) Sin2 L ]3/2 }

-Then  the azimuth µ of the loxodrome verifies

        rm Tan µ dB =  r'm dL
        ( Cos µ ) ds =  rm dB

-Since r & r'm depend on L and L is an unknown function of B, we could solve a differential equation to find Tan Az
-But the eccentricity of the equator is very small, and L only appears in a corrective term,
  so we can use the expression L(B) valid for an ellipsoid of revolution

     L - L1 =  ( Tan µ ).{ Ln [ ( Tan ( 45° + B/2 ) ) / (Tan  ( 45° + B1/2 ) ) ] - (e/2).Ln [ ( 1 + e.sin B ).( 1 - e.sin B1 )/( 1 - e.sin B )/( 1 + e.sin B1 ) ] }
 

Flag:  F00

   CF 00  abcL is called to store a b c L in registers R01-R02-R03-R00
   SF 00:  you have to initialize R00-R01-R02-R03
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )        µ ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

 
Example1:

   -77.03560  ENTER^
    38.55172  ENTER^
      2.20138  ENTER^
    48.50112  XEQ "LOX3"  >>>>   D = 6453.471591 km                          ---Execution time = 2m40s---
                                            X<>Y   µ  = 80°10'15"8102

Example2:

   0   ENTER^
   0   ENTER^
  70  ENTER^
   0      R/S       >>>>  D = 7792.380608 km                          ---Execution time = 14s---
                       X<>Y   µ = 90°

Example2:

   0  ENTER^
   0  ENTER^
   0  ENTER^
  80    R/S       >>>>  D = 8885.152534 km                          ---Execution time = 15s---
                       X<>Y  µ = 0°

Notes:

-Due to roundoff errors, there is a loss of significant digits if the 2 latitudes are very close but not equal.
-If the latitudes are equal , µ = +/- 90° and if the longitudes are equal , µ = 0 or 180° so the program is much faster
-The azimuths are measured clockwise from North.
 

Gravity on a Triaxial Ellipsoid
 

-In reference [5], M. Caputo gives a generalization of Somigliana's formula that also works for triaxial ellipsoids.

      x2/a2 + y2/b2 + z2/c2 = 1

-The unique assumption is that the gravity vector must always be normal to the ellipsoid:
-In other words, the ellipsoid is an equipotential surface.
 
 
      STACK        INPUTS      OUTPUTS
           Z         L ( ° ' " )          g(0)
           Y         B ( ° ' " )          g(0)
           X          h ( m )          g(h) 

 
Example1:     L =  77°03'56" W ,  B = 38°55'17"2 N ,   h = 67m                ( U.S. Naval Observatory , Washington D.C. )

 -77.0356    ENTER^
  38.55172   ENTER^
       67         XEQ "GRV3"  >>>>   g(67) = 9.800516275  m/s2                                  ---Execution time = 5s---
                                           X<>Y   g(0)  = 9.800723034  m/s2
 

Example2:     L = 116°51'50"4 W ,    B = 33°21'22"4 N    h = 1706 m                ( Mount Palomar Observatory )

 -116.51504    ENTER^
   33.21224      ENTER^
      1706            R/S     >>>>   g(1709) = 9.790660012  m/s2                                      ---Execution time = 5s---
                                    X<>Y      g(0)    = 9.795923287  m/s2

Notes:

-The longitudes are positive East.

-The formula is rigorous if h = 0 but only approximate otherwise.
-So, this program should not be used for large h-values.

-Unfortunately, the Earth gravitational field is much more complex than that of a triaxial ellipsoid,
 so this program only gives slightly better results than the classical Somigliana's formula.
 

Surface Area of an Ellipsoid
 

-"SAE" uses a Carlson elliptic integral of the second kind ( namely, RG ) and we have:

    Area = 4.PI.RG( a2b2 , a2c2 , b2c2 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             a             /
           Y             b             /
           X             c         Area

 
Example:

    2  ENTER^
    4  ENTER^
    9  XEQ "SAE"  >>>>   Area =  283.4273840                                   ---Execution time = 42s---
 

Surface Area of a Hyperellipsoid
 

-"SAHE" employs the method described in reference [6].
-With minor modifications and a few simplifications, it yields:

  A =  [ a2 ........ aN / Gamma((N+1)/2) ]   §01  f(u) du

       with     f(u) = (1-u)(N-1)/2 u -1/2  [ 1 + SUMi=2,...,N (a12/ai2/bi) ] ( b2 .....  bN ) -1/2        where        bi  =  1 - u + u (a12/ai2)
 

-To remove the singularities for u = 0 & u = 1 , we can use the change of variable  u = Sin2 v  and apply a Gauss-Legendre formula.
-In reference [6] , the authors employ another change of variable and the Romberg method.

-"SAHE" uses a Gauss-Chebyshev formula, namely:

     §01  [ (1-x) / x ]1/2  g(x) dx  =  SUMi=1,2,...,n  wi g(xi)

-Unlike many Gaussian formulae, the abscissas & weights may be easily calculated:

    xi =  Sin2 [ (2i-1) / (2n+1) ] PI/2    &   wi = [ 2.PI / (2n+1) ] Cos2 [ (2i-1) / (2n+1) ] PI/2

-So, we can choose n large enough to get an almost perfect accuracy.
 

Data Registers:           •  R00 = N < 11                      ( Registers R00 thru RNN are to be initialized before executing "SAHE" )

                                      •  R01 = a1         •  R02 = a2    .......................      •  RNN = aN
 
 
      STACK        INPUT      OUTPUT
           X             n            An

  where   n = number of points in the Gauss-Chebyshev formula
   and    An = corresponding approximation of the area

Examples:          N = 5  ,   semi-axes =  4  5  6  7  8

   5  STO 00

   4  STO 01
   5  STO 02
   6  STO 03
   7  STO 04
   8  STO 05

-With  n = 30  ,   30   XEQ "SAHE"  >>>>  A =  31971.28857           ---Execution time = 137s---
-With  n = 40  ,   40           R/S          >>>>  A =  31971.28859           ---Execution time = 181s---
 

Quadratic Hypersurface
 

-"QHS" reduces the cartesian equation of a quadratic hypersurface (QHS) in a n-dimensional space  ( n > 1 )
-So, it can be used for ellipsoids or hyperellipsoids, but also for more general quadratic surfaces.

-Given    (QHS):  a11 x12 + a22 x22 + .............. + ann xn2 + Sum i < j  ai j xi xj  +  b1 x1 + b2 x2 + ............. + bn xn =  c

    the elements  ai j ( i < j ) are gradually zeroed by the Jacobi's iterative method.

-When all these elements are smaller than 10 -8 , the quadratic form is diagonalized.
-Then several translations and/or rotations are performed to reduce the equation further.
 

Data Registers:           •  R00 = n                      ( All these registers are to be initialized before executing "QHS"  )

           I              •  R01 = a11    • Rn+1 = a1,2   • R2n    = a2,3   ..................................  • R(n2+n)/2 = an-1,n       • R(n2+n)/2+1 = b1        • R(n2+3n)/2+1 = c
          N             •  R02 = a22    • Rn+2 = a1,3   • R2n+1 = a2,4     ...................                                                    • R(n2+n)/2+2 = b2
          P                ...............         .................      ...................                                                                               ........................
          U               ...............         .................    • R3n-2 = a2,n
          T               ................      • R2n-1 = a1,n
          S              •  Rnn = ann                                                                                                                              • R(n2+3n)/2 = bn

 -----------------------------------------------------------------------------------------------------------------------------------------------------------------

                When the program stops:

                                       R00 = n
          O
          U              R01 = a'11      Rn+1 =  b'1       R2n+1 = c'
          T               R02 = a'22      Rn+2 =  b'2
          P              ...............         .................
          U
          T
          S               Rnn = a'nn      R2n =   b'n

The reduced equation is:    (QHS):  a'11 X12 + a'22 X22 + .............. + a'nn Xn2  + b'1 X1 + b'2 X2 + ............. + b'n Xn =  c'

 where   b'i  = 0  if  a'ii # 0  and  c' = 0 or 1
 
 
      STACK        INPUTS      OUTPUTS
           Y             /        bbb.eee'
           X             /          1.eee

   Where  1.eee = control number of the reduced equation
     and   bbb.eee' = -------------------  eigenvalues

Example:    (QHS):  3 x2 + 4 y2 + 7 z2 + 6 t2 + 9 x.y + 3 x.z + 7 x.t + 4 y.z + 8 y.t + 2 z.t + 9 x + 2 y + 5 z + 4 t = 10  in a 4-dimensional space.

       SIZE 016  or greater

   4  STO 00  and we store these ( n2+3n+2 )/2 = 15 coefficients as shown below:

           3     9     4      2       9       10                      R01     R05     R08     R10       R11       R15
           4     3     8               2                     in          R02     R06     R09                  R12                     respectively
           7     7                      5                                  R03     R07                             R13
           6                             4                                  R04                                         R14
 

 XEQ "QHS"  the HP-41 displays the successive greatest  | ai j |  ( which tend to zero )  and when the program stops,

         X = 1.009   &   Y = 10.013                                               ---Execution time = 2mn19s---

         R01 =  -0.209131158     R05 = 0      R09 = 1
         R02 =   0.297006169     R06 = 0
         R03 =   1.235467092     R07 = 0
         R04 =   2.716166061     R08 = 0

  So, the reduced equation is   -0.209131158 X2 +  0.297006169 Y2 + 1.235467092 Z2 + 2.716166061 T2 = 1

-And we find the eigenvalues in R10 thru R13, namely:

         -1.035428817
          1.470506591
          6.116918408
          13.44800383

Notes:

-Register P is used, so don't interrupt this program.
-For n = 7 , the execution time is of the order of  18 minutes.
-Since "QHS is executed from a module, n could reach 23.
 

Ellipsoidal Wave Differential Equation ( RK6 )
 

 "EWF" uses a Runge-Kutta-Nystrom method of order 6 to solve the differential equation:

      W''(x) = W(x) [ -h + n(n+1) k2 Sn2 (x) - q k4 Sn4 (x) ]     where  Sn (x) is a Jacobian elliptic function.
 

Data Registers:              R00 = temp           ( Registers R01 thru R09 are to be initialized before executing "EWF" )

                                      •  R01 = h                •  R05 = x0           •  R08 = H = stepsize          R10 to R14: temp         R19 to R27: used by "SNX"
                                      •  R02 = n                •  R06 = W(x0)     •  R09 = N = nb of steps     R15 to R18: unused
                                      •  R03 = k2 < 1        •  R07 = W'(x0)
                                      •  R04 = q
Flags: /
Subroutine:  "SNX"
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             /   W'(x0 + N.H)
           Y             /   W(x0 + N.H)
           X             /      x0 + N.H

 
Example:    h = 1.2    n = 1.7   k2 = 0.8   q = sqrt(2)     W(0) = 1   W'(0) = 0

    1.2  STO 01                       0  STO 05
    1.7  STO 02                       1  STO 06
    0.8  STO 03                       0  STO 07
      2   SQRT  STO 04

-If you want to compute W(1) with  H = 0.1 and therefore,  N = 10

    0.1  STO 08
    10   STO 09

    XEQ "EWF"  >>>>      X   =  1
                          RDN   W(1) =  0.640627307                                    ---Execution time = 7m15s---
                          RDN  W'(1) = -0.401079727

-Simply press R/S to continue with the same parameters, it yields:

               W(2) = 0.542730302
              W'(2) = 0.260776593

Note:

-Though "EWF" doesn't work if k2 = 1 , it will work with k2 = .9999999999 and the difference will be negligible.
 

Ellipsoidal Wave Function ( Series )
 

-"EWF2" is less general than "EWF".
-The solution is of the form  W(x) = 1 + A1 t + A2 t2 + ............. + Am tm + ................

     where  t = Sn2 (x)

-This leads to the recurrence relation

     ( 2m+1 ) (2m+2 ) Am+1 = - q k4 Am-2 + Am-1 k2 [ n(n+1) - 2 ( m-1 ) ( 2m - 1 ) ] + Am [ - h + 4 m2 ( k2 + 1 ) ]
 

Data Registers:              R00 = x                 ( Registers R01 thru R04 are to be initialized before executing "EWF" )

                                      •  R01 = h                R05 to R13: temp
                                      •  R02 = n
                                      •  R03 = k2 < 1        R10 = t = Sn2 (x)   R12 = Sum
                                      •  R04 = q
Flags: /
Subroutine:  "SNX"
 
 
 
      STACK        INPUTS      OUTPUTS
           X             x          W(x)

 
Example:             h = 1.2    n = 1.7   k2 = 0.8   q = sqrt(2)

    1.2  STO 01
    1.7  STO 02
    0.8  STO 03
      2   SQRT  STO 04    and to calculate  W(1)

      1   XEQ "EWF2"  >>>>   W(1) = 0.640627308                                    ---Execution time = 79s---
 
Notes:

-In this example, "EWF2" is faster than "EWF", but it's not always the case...
-For instance, W(2) = 0.542730301  is returned after a very long time.

-Though "EWF" doesn't work if k2 = 1 , it will work with k2 = .9999999999 and the difference will be negligible.
 

Carlson Elliptic Integral of the 1st kind
 

-Carlson has given a new definition of a standard elliptic integral of the first kind:

       RF(x;y;z) =  (1/2) §0infinity  ( ( t + x ).( t + y ).( t + z ) ) -1/2 dt           with  x , y , z  non-negative and at most one is zero

-This program uses the following property:

      RF(x;y;z) = RF((x+L)/4;(y+L)/4;(z+L)/4)  where  L = x1/2y1/2 + x1/2z1/2 + y1/2z1/2

-This transformation is performed until  x , y , z are close enough to apply  RF(x;y;z) = µ -1/2   with  µ = (x+y+z)/3     ( we have  RF(x;x;x) = x -1/2 )
 
 
      STACK       INPUTS    OUTPUTS
          Z            z           /
          Y            y           /
          X            x     RF(x;y;z)

 
Example:

    4  ENTER^
    3  ENTER^
    2  XEQ "RF"  >>>>  RF(2;3;4) = 0.584082842                                   ---Execution time = 8s---
 

Symmetric Carlson Elliptic Integral of the 2nd kind
 

  RG(x;y;z) =  (1/4)  §0infinity  ( ( t + x ).( t + y ).( t + z ) ) -1/2 .( x/(t+x) + y/(t+y) + z/(t+z) )  t.dt

And we have:  2.RG(x;y;z) = z.RF(x;y;z) - (x-z)(y-z)/3 RD(x;y;z) + ( x.y/z )1/2
 
 
      STACK       INPUTS    OUTPUTS
          Z            z           /
          Y            y           /
          X            x     RG(x;y;z)

 
Example:

    4  ENTER^
    3  ENTER^
    2  XEQ "RG"  >>>>  RG(2;3;4) = 1.725503028                                   ---Execution time = 40s---

Note:

-If  one of the arguments is zero, do not place it in  Z-register ( there would be a DATA ERROR )
 

Carlson Elliptic Integral of the 3rd kind
 

-The elliptic integral of the third kind is defined by

   RJ(x;y;z;p) =  (3/2)  §0infinity  ( t + p ) -1  ( ( t + x ).( t + y ).( t + z ) ) -1/2  dt        with  x , y , z  non-negative and at most one is zero    p > 0

-We have  RJ(x,x,x,x) = x -3/2
-"RJ" applies  RJ(x;y;z;p) = 3 Sumn=0,1,2,....,k  RF(an,bn,bn)/4n + 1/(4k+1µ 3/2)

  where   x0 = x , y0 = y , z0 = z , p0 = p        ;      xn+1 = ( xn+ Ln )/4 , yn+1 = ( yn+ Ln )/4 , zn+1 = ( zn + Ln )/4 , pn+1 = ( pn +Ln )/4

    with    Ln = xn1/2yn1/2 + xn1/2zn1/2 + yn1/2zn1/2
               an = ( pn( xn1/2 + yn1/2 + zn1/2 ) + ( xnynzn )1/2 )2    ;   bn = pn ( pn + Ln )2

    and    µ = (xk+1+yk+1+zk+1+2pk+1)/5

-The iterations stop when  xn , yn , zn , pn  are close enough to produce a good accuracy.
 
 
      STACK       INPUTS    OUTPUTS
          T          p > 0           /
          Z            z           /
          Y            y           /
          X            x     RJ(x;y;z;p)

 
Example:

     4  ENTER^
     3  ENTER^
     2  ENTER^
     1  XEQ "RJ"  >>>>  RJ(1;2;3;4) = 0.239848100                                   ---Execution time = 34s---
 
 
 
 

References:

[1]  Geodetic Constants
[2]  J. Feltens - Vector method to compute the Cartesian (X, Y, Z) to geodetic (f, l, h) transformation on a triaxial ellipsoid - J Geod (2009) 83:129–137
[3]  Marcin Ligas  - Cartesian to geodetic coordinates conversion on a triaxial ellipsoid - J Geod (2012) 86:249–256

-With many thanks to Charles Karney for the link that he sent me to get Jacobi's method:
  http://lists.maptools.org/pipermail/proj/2012-October/006448.html

[4]  If you want to visualize the geodesics on a triaxial ellipsoid, go to this excellent page
[5]  Michèle Caputo - "Some Space Gravity Formulas and the Dimensions and the Mass of the Earth"
[6]  Dunkl & Ramirez - "Computing Hyperelliptic Integrals for Surface Measure of Ellipsoids"
[7]  B.C. Carlson, "Numerical Computation of Real or Complex Elliptic Integrals"