Geodesy for the HP48S/SX/G/GX


-This zip file contains:

   GEODESY48  a binary file ( #28455d / 13315 bytes )  which works with an HP48S/SX/G/GX
   GEODESY50  a binary file ( #32617d / 13315 bytes )  which works with an HP49/50G
   GEODESY      an txt file which works with all HP48/49/50

-All are directories.
 

Notes:

-Longitudes are positive East.
-The azimuths are measured clockwise from North.

-For the spheroidal model of the Earth, I've used  a = semi-major axis = 6378136.61 m anf  f = flattening = 1/298.256421
-If you prefer the WGS84 constants, place 1/298.257223563 in 'f' and 6378137 in 'a'

-For the triaxial ellipsoid, the equatorial major axis is in the direction   14°92911W
  and the semiaxes are 6378171.27379 m , 6378101.94621 m , 6356751.86801 m
-Of course, the last decimals are somewhat arbitrary...
-If you prefer, store other values in the variables A B C.
 
 
 Functions  Description
 AXI3
   GRV3
     L0
     A
     B
     C
     GABC
     GRV
   JACOBI
    XYZ->UV
    UV->XYZ
    EQ
    BETA
    U1
    V1
    U2
    V2
    N
    DIST
      N
      DGT
    APPROX
      GGRD
      GDTD
      GECD
      DDD
      COEF
      N
    A
    B
    C
    L0
    GGR->XYZ
    GDT->XYZ
    GL16
    S->R
 DGT
 DGT2
 DG2
 DGT3
 FWD2
 FWD3
 GRV
 LOX0
 LOX3
 GM
 w
  f
 a		  Directory = Triaxial Ellipsoid
   Directory = Gravity
     Longitude 0 = -14°92911
     A = equatorial semimajor axis = 6378171.27379 m
     B = equatorial semiminor axis = 6378101.94621 m
     C = semiminor axis = 6356751.86801 m
     GABC = Gravity at the points (A,0,0) (0,B,0) (0,0,C)
     GRV = Gravity on a triaxial ellipsoid
  Directory = Jacobi's Method
   Rectangular Coord. -> Jacobi's Coordinates
   Jacobi's Coordinates -> Rectangular Coordinates
   Jacobi's Equation
   Parameter "beta"
   Coordinate U of the 1st point
   Coordinate V of the 1st point
   Coordinate U of the 2nd point
   Coordinate V of the 2nd point
   Number of subintervals for Gaussian Integration Formula
   Directory = Rigorous Geodesic Distance
     Number of subintervals for Newton-Cotes 7 pts Formula
     Rigorous Geodesic Distance on a Triaxial Ellipsoid
   Directory = Approximate Method
     Geographic Coord. -> Geodesic Distance
     Geodetic Coord. -> Geodesic Distance
     Geocentric Coord. -> Geodesic Distance
     Subroutine
     Coefficient = 1000 
     Number of subintervals for Gaussian Integration Formula
   A = equatorial semimajor axis = 6378171.27379 m
   B = equatorial semiminor axis = 6378101.94621 m
   C = semiminor axis = 6356751.86801 m
   Longitude0 = -14°92911
   Geographic Coordinates -> Rectangular Coordinates
   Geographic Coordinates -> Rectangular Coordinates
   Gauss-Legendre 16-point Formula
   Spherical -> Rectangular Conversion
 Geodesic Distance, Andoyer's 1st order formula
 Geodesic Distance, Andoyer's 2nd order formula
 Geodesic Dist, Andoyer's 2nd order formula ( 1st order Az )
 Geodesic Distance, Vincenty's formula
 Forward Problem, Andoyer's 2nd order formula
 Forward Problem, Vincenty's formula
 Gravity on an Ellipsoid of Revolution - rigorous formula
 Loxodromy on a Sphere R = 6371 km
 Loxodromy on an Ellipsoid of revolution, 3rd order formula
 Earth Gravitational Constant = 3.9860044188 E14 SI
 Earth Angular Velocity = 0.00007292115 rd/s
 Earth Flattening = 1/298.256421
 Earth Equatorial Radius = 6378136.61 m

 

Geodesic Distance - Andoyer's Formula ( 1st order )
 

Formulae:

-Let                    L = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
-We calculate     S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

-Then,                D = a.{ 2µ + f.[ (3R-µ)/(1-S)  sin2 F  cos2 G  - (3R+µ)/S  sin2 G  cos2 F ] }
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")             /
       Level 3     Lat1 (° ' ")             /
       Level 2    Long2 (° ' ")             /
       Level 1     Lat2 (° ' ")     Dist ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N

   -77.0356    ENTER
   38.55172   ENTER
   2.20138     ENTER
   48.50112   DGT  >>>>   D = 6181.6153 km

-In this example, the error is about  6 meters.
-If the points are not ( nearly ) antipodal, the relative error is less than  10 -5 and the maximum error is of the order of 60 meters in absolute value:
 with ( 0° , -40° ) ( 120° , 40° )   error = 66 meters
-However, errors can reach several kilometers for nearly antipodal points.
 

Geodesic Distance - Andoyer's Formulae  ( 2nd order for the Distance & the Azimuths )
 

-Second order formulae - given in reference [4] - are used to calculate the distances & the azimuths
-They employed the parametric latitudes
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")             /
       Level 3     Lat1 (° ' ")     Az2 (° ' ")
       Level 2    Long2 (° ' ")     Az1 (° ' ")
       Level 1     Lat2 (° ' ")     Dist ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N

   -77.0356    ENTER
   38.55172   ENTER
   2.20138     ENTER
   48.50112   DGT2  >>>>  D  = 6181.62143 km
                                           Az1 =  51°47'36"81 = forward azimuth
                                           Az2 = -68°09'58"96 = back azimuth

Warning:

-The formulae that compute the azimuths do not work well if the difference in longitudes is very close to +/-90° !!!
 

Geodesic Distance - Andoyer's Formula ( 2nd order but 1st order for the Azimuths )
 

-DG2 uses a formula which is also second-order in the flattening f:
 

-With   L = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
              S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

 we define  A = ( sin2 G  cos2 F )/S + ( sin2 F  cos2 G )/( 1-S )
      and      B = ( sin2 G  cos2 F )/S - ( sin2 F  cos2 G )/( 1-S )

-The geodesic distance is then   D = 2a.µ ( 1 + eps )

  where   eps =    (f/2) ( -A - 3B R/µ )
                       + (f2/4) { { - ( µ/R + 6R/µ ).B + [ -15/4 + ( 2S - 1 ) ( µ/R + (15/4) R/µ ) ] A + 4 - ( 2S - 1 ) µ/R } A
                                      - [ (15/2) ( 2S - 1 ).B  R/µ - ( µ/R + 6R/µ ) ] B }

 >>> The forward azimuth Az1 in the direction from P1 towards P2 is obtained by the following formulae:    Az1 = Az1sph + dAz

      where  dAz = (f/2).( cos2 F - sin2 G ).[ ( 1+µ/R ).( sin G cos F )/S + ( 1-µ/R ).( sin F cos G )/(1-S) ].( sin 2L )

       and    Tan Az1sph = cos b2 sin 2L / ( cos b1 sin b2 - sin b1 cos b2 cos 2L )

  Az1sph is the forward azimuth in a spherical model of the Earth,
  dAz is a correction that takes the Earth's flattening f into account.

-The back azimuth ( from P2 to P1 ) is obtained by replacing L by -L and exchanging b1 & b2 in the formulas above.
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")             /
       Level 3     Lat1 (° ' ")     Az2 (° ' ")
       Level 2    Long2 (° ' ")     Az1 (° ' ")
       Level 1     Lat2 (° ' ")     Dist ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N

   -77.0356    ENTER
   38.55172   ENTER
   2.20138     ENTER
   48.50112   DG2  >>>>   D  = 6181.62145 km
                                         Az1 =  51°47'36"76 = forward azimuth
                                         Az2 = -68°09'59"26 = back azimuth

Notes:

-"DG2" gives accurate results if the 2 points are not nearly antipodal.
-Az is obtained with an accuracy of the order of a few arcseconds ( the error is smaller than 1" in this example ).
 

Geodesic Distance - Vincenty's Formula
 

-Let   L =  L2 - L1  = L0   ;   U = Arctan((1-f).tan b1)  ;   V = Arctan((1-f).tan b2)

-We compute recursively:      Ln+1 = L + (1-C).f.(sin µ).{ c + C.sin c [ cos 2d + C.cos c ( -1 + 2.cos2 2d ) ] }

   where      sin c = [ ( cos V sin Ln )2 + ( cos U sin V - sin U cos V cos Ln )2 ]1/2
                 cos c = sin U sin V + cos U cos V cos Ln
                       µ = Arcsin [ ( cos U cos V sin Ln )/sin c ]  ;  µ = azimuth of the geodesic at the equator
                      C = (f/16).(cos2 µ ).[ 4 + f.(4-3.cos2 µ ) ]
              cos 2d = cos c - 2.sin U sin V / cos2 µ

   until   Ln+1 = Ln = lambda

-Then,     D = a.A.(1-f).(c-D)

    with     A = 1 + (u/16384).{ 4096 + u.[-768 + u.(320-175.u)] }      ;       B = (u/1024).{ 256 + u.[ -128 + u.(74-47.u)] }
               D = (B.sin c).{ cos 2d + (B/4). [ ( cos c ).( -1 + 2.cos2 2d ) - (B/6).(cos 2d).(-3+4.sin2 c ).(-3+4.cos2 2d )] }
                u =  f(2-f)/(1-f)2 . cos2 µ

-The azimuths  Az1 &  Az2 are finally calculated by the formulas:

       Tan Az1 = [ cos V sin (lambda) ] / [ cos U sin V - sin U cos V cos (lambda) ]
       Tan Az2 = [ -cos U sin (lambda) ] / [ sin U cos V - cos U sin V cos (lambda) ]
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")             /
       Level 3     Lat1 (° ' ")     Az2 (° ' ")
       Level 2    Long2 (° ' ")     Az1 (° ' ")
       Level 1     Lat2 (° ' ")     Dist ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N

   -77.0356    ENTER
   38.55172   ENTER
   2.20138     ENTER
   48.50112   DGT3  >>>>   D  = 6181.62143367 km
                                           Az1 =  51°47'36"81 = forward azimuth
                                           Az2 = -68°09'58"97 = back azimuth

Notes:

-DGT3 will still work if the points are nearly antipodal: an intermediate point is used
 with the built-in function ROOT to get the correct geodesic distance ( tagged Dist2: ).

-In this case, however, long execution time is to be expected,
 the azimuths might be computed less accurately
 and the last executed function is KILL.
 

Direct Problem - Andoyer's 2nd order Formulae
 

-Given the coordinates ( L1 , b1 ) of the first point  P1 , the forward azimuth Az1 and the length D of the geodesic,
 "FWD" computes the longitude, the latitude ( L2 , b2 ) of the second point  P2 and the backward azimuth Az2

-The formulas are given in the reference [4]
 
 
      STACK        INPUTS      OUTPUTS
      Level 4       L1 ( ° ' " )             /
      Level 3       b1 ( ° ' " )      Az2 ( ° ' " )
      Level 2     Az1 ( ° ' " )       L2 ( ° ' " )
      Level 1        D ( km )       b2 ( ° ' " )

 where   Az1  is the forward azimuth in the direction from the first point P1 towards the second point P2

Example:     L1 = 10°30'     b1 = 49°41'     Az1 = 12°24'     D = 16000 km

   10.30   ENTER
   49.41   ENTER
   12.24   ENTER
  16000   FWD2  >>>>     b2 =   -14°06'40"75
                                         L2 = -177°03'07"98
                                        Az2 =    -8°15'03"68

Notes:

-Andoyer's formulas have one advantage: they are non-iterative!
-In most cases, their accuracy is satisfactory, all the more that the Earth's irregularities cannot be taken into account.
-If, however, a better precision is required, this can be achieved via Vincenty's formulas below.
 

Direct Problem - Vincenty's Formulae
 

-Given the coordinates ( L1 , b1 ) of the first point  P1 , the forward azimuth Az1 and the length D of the geodesic,
 "FWD" computes the longitude and the latitude ( L2 , b2 ) of the second point  P2

-The formulas are given in the reference [3]
 
 
      STACK        INPUTS      OUTPUTS
      Level 4       L1 ( ° ' " )             /
      Level 3       b1 ( ° ' " )             /
      Level 2     Az1 ( ° ' " )       L2 ( ° ' " )
      Level 1        D ( km )       b2 ( ° ' " )

 where   Az1  is the forward azimuth in the direction from the first point P1 towards the second point P2

Example:     L1 = 10°30'    b1 = 49°41'     Az1 = 12°24'     D = 16000 km

   10.30   ENTER
   49.41   ENTER
   12.24   ENTER
  16000   FWD3  >>>>     b2 =   -14°06'40"75
                                         L2 = -177°03'07"98
 

Gravity on an Ellipsoid of Revolution
 

-Let  b = geographic ( geodetic ) latitude and  h = altitude
-The gravity g above an ellipsoid of revolution may be computed by the rigorous formulae ( cf reference [6] ):

             g  =  ( p2+ q2 )1/2

 where   p =  { - [ GM + ( omega )2 a2 E ( q' / 2q0 ) ( sin2 F - 1/3 ) ] / ( u2 + E2 ) + ( omega )2 u cos2 F } / w
  and     q =   [ a2 ( q / q0 ) ( u2 + E2 ) -1/2 - ( u2 + E2 )1/2 ] ( omega )2 ( sin F cos F ) / w

            u  =  { [ ( x2 + y2 + z2 - E2 ) / 2 ] [ 1 + { 1 + 4 E2 z2 / ( x2 + y2 + z2 - E2 )2 }1/2 ] }1/2
  with    E =   a ( 2 f - f 2 )1/2  ,   sin F = z / u   ,   cos F = ( x2 + y2 )1/2 / ( u2 + E2 )1/2
            w =  [ ( u2 + E2 sin2 F )1/2 / ( u2 + E2 )1/2 ]1/2

            q  = (1/2) [ ( 1 + 3 u2 / E2 ) Atan ( E / u )  - 3 E / u ]                                                            ( x2 + y2 )1/2  = ( N + h ) cos b
  and    q0 = (1/2) [ ( 1 + 3 b' 2 / E2 ) Atan ( E / b' ) - 3 E / b' ]         ( b' = polar radius )                          z            = [ ( b'/a )2 N + h ] sin b
            q' =   3 ( 1 + u2 / E2 ) [ 1 - ( u / E ) Atan ( E / u ) ] - 1                                                                 N          = a ( 1 - e2 sin2 b ) -1/2
 

-If, however, we use ATAN to calculate  q , q0 and q' , the results will not be accurate because of cancellation of leading digits.
-So, the following program employs ascending series to compute:

       Atan t + 3 [ ( Atan t ) / t - 1 ] / t      =  4 t3  Sumk=0,1,2,.....  ( k + 1 ) ( - t2 )k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]
    3 ( 1 + 1 / t2 ) [ 1 - ( Atan t ) / t ] - 1  =  6 t2  Sumk=0,1,2,.....  ( - t2 )k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]

-Since the argument t never exceeds 0.0821 for the Earth, only a few terms must be calculated.
 
 
      STACK        INPUTS      OUTPUTS
       Level 2         b ( ° ' " )             /
       Level 1          h ( m )       g ( m/s2 )

Example1:     b = 38°55'17"2     h = 23456 m

  38.55172   ENTER
    23456      GRV  >>>>   g = 9.728751601 m/s2
 

Example2:     b = 38°55'17"2     h = 12345678 m

  38.55172   ENTER
 12345678   GRV  >>>>   g = 1.078713338 m/s2
 

Loxodromy on a Sphere
 

-A loxodromic curve is a line on the surface of the Earth which cuts all meridians at the same angle µ  ( the azimuth )
-The 2 following programs compute D = the length of these curves and the azimuth µ
 

Formulae:        L2 - L1 =  ( Tan µ ). Ln [ ( Tan ( 45° + b2/2 ) ) / (Tan  ( 45° + b1/2 ) ) ]              Li = longitudes ;  bi = latitudes
                                 D =  R.( b2 - b1 )/Cos µ                                                                              R = mean radius of the Earth = 6371 km
 
 
      STACK        INPUTS      OUTPUTS
       Level 4       L1 ( ° ' " )             /
       Level 3       b1 ( ° ' " )             /
       Level 2       L2 ( ° ' " )        µ ( ° ' " )
       Level 1       b2 ( ° ' " )        D ( km )

Example:

 U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W =  -77°03'56"0   ;   b1 = 38°55'17"2 N =  +38°55'17"2
         The Observatoire de Paris:                           L2 =   2°20'13"8 E  =   +2°20'13"8    ;   b2 = 48°50'11"2 N =  +48°50'11"2

   -77.03560  ENTER
    38.55172  ENTER
      2.20138  ENTER
    48.50112  LOX0  >>>>   D = 6436.5499 km
                                              µ = 80°08'14"
 

Loxodromy on an Ellipsoid of revolution, 3rd order formula
 

-Higher accuracy is obtained if we take the Earth's flattening into account.
 

Formulae:

    L2 - L1 =  ( Tan µ ).{ Ln [ ( Tan ( 45° + b2/2 ) ) / (Tan  ( 45° + b1/2 ) ) ] - (e/2).Ln [ ( 1 + e.sin b2 ).( 1 - e.sin b1 )/( 1 - e.sin b2 )/( 1 + e.sin b1 ) ] }
            D =  [ a.( 1 - e2 )/ cos µ ] §b1b2  ( 1 - e2 sin2 t ) -3/2 dt     ( § = Integral )

-One could use any integration formula but 'LOX3' evaluates this integral by a series expansion up to the terms in  e6
  ( the first neglected term is proportional to  e8 )

            D = [ a.( 1 - e2 )/ cos µ ] [ ( 1 + 3e2/4 + 45e4/64 + 175e6/256 ).( b2 - b1 ) - ( 3e2/8 + 15e4/32 + 525e6/1024 ).( sin 2b2 - sin 2b1 )
                   + ( 15e4/256 + 105e6/1024 ).( sin 4b2 - sin 4b1 ) - ( 35e6/3072 ).( sin 6b2 - sin 6b1 ) ] + .......................

-However, we cannot use this formula if   cos µ = 0 , in this case ( i-e  if  b1 = b2 ) we have:

           D = a.( cos b ).( L2 - L1 ) / ( 1 - e2 sin2 b ) -1/2    the loxodromic line is the parallel of latitude b = b1 = b2
 
 
      STACK        INPUTS      OUTPUTS
       Level 4       L1 ( ° ' " )             /
       Level 3       b1 ( ° ' " )             /
       Level 2       L2 ( ° ' " )        µ ( ° ' " )
       Level 1       b2 ( ° ' " )        D ( km )

Example:

 U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W =  -77°03'56"0   ;   b1 = 38°55'17"2 N =  +38°55'17"2
         The Observatoire de Paris:                           L2 =   2°20'13"8 E  =   +2°20'13"8    ;   b2 = 48°50'11"2 N =  +48°50'11"2

   -77.03560  ENTER
    38.55172  ENTER
      2.20138  ENTER
    48.50112  LOX3  >>>>   D = 6453.389608 km
                                              µ = 80°10'15"31
 

Gravity at the 3 points (A,0,0) (0,B,0) (0,0,C)
 

-Somigliana's formula works on an ellipsoid of revolution.
-In reference [11], M. Caputo gives a generalization of this 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.

-But to apply this formula, we first have to calculate the gravity at the points of intersection of the axes and the ellipsoid itself

-The formulas are given in reference [11]  ( warning: M. Caputo denotes  a1 = b , a2 = a , a3 = c ). With small differences in the notations:

  Let  E2 = (b2 - c2) / c2  and   n =  (a2 - b2) / b2

    Ai,j = A'i,j + n A"i,j

   A'1,1 = A'1,2 = A'2,1 = A'2,2 = [ 0.75 / (b2 - c2)5/2 ] [ Arc tan E - (E/3) (5.E2 + 3) / (1+E2)2 ]
   A"1,1 = [ 5 / 16 / (b2 - c2)7/2 ] [  -Arc tan E + ( E/15/(1+E2) ) (20 - (5-13.E4)/(1+E2)2 ) ]         A"2,1 = A"1,2 = 3 A"1,1 ;  A"2,2 = 5 A"1,1

   A'1,3 = A'2,3 = [ 3 / (b2 - c2)5/2 ] [  -Arc tan E + ( E/3/(1+E2) ) (2E2+3) ]
   A"1,3 = [ 15 / 8 / (b2 - c2)7/2 ] [ Arc tan E - (E/30) ( 25 + (5-9.E4) / (1+E2)2 ) ]                        A"2,3 = 3 A"1,1

-Then, we must calculate  k1 & k2  with

    D.K1 / w2 = - A'1,1 - n [ A'1,1 + 6 A"1,1 + (c2/b2) A'2,3 ]                       w2 = angular velocity of the Earth = 7.292115 E-5 rd/s
    D.K2 / w2 = - A'1,1 - n [ A'1,1 -  (c2/b2) A'2,3 ]

  where  D = 4 A'1,1 [ 2 A'1,1 -  (c2/b2) A'2,3 ] - 2 n (c2/b2) A'2,3 ( A'1,1 + 6 A"1,1 ) + 4 n A'1,1 [ 2 A'1,1 + 12 A"1,1  -  (c2/b2) A"2,3 ]

-Finally,  with  GM = Earth gravitational constant = 3.9860044188 E14 m3/s2

    •  ga/a  =  ( GM + 4 K2 / a2 ) / ( abc ) - 2 ( A1,2 K1 + 3 A2,2 K2 ) - w2
    •  gb/b =  ( GM + 4 K1 / b2 ) / ( abc ) - 2 ( 3 A1,1 K1 + A2,1 K2 ) - w2
    •  gc/c =       GM / ( abc ) - 2 ( A1,3 K1 + A2,3 K2 )
 

-Unfortunately, the expressions involving Arc Tan E  would produce low accuracy because of cancellation of leading digits.
-In the following program, I've used truncated series to get a better precision, namely:

   Arc tan E - (E/3) (5.E2 + 3) / (1+E2)2 =  (8/15) E5 - (8/7) E7 + (16/9) E9 - (80/33) E11 + (40/13) E13 - (56/15) E15 + ...
  -Arc tan E + ( E/15/(1+E2) ) (20 - (5-13.E4)/(1+E2)2 ) = -(16/35) E7 + (64/45) E9 - (32/11)  E11 + (64/13) E13 - (112/15) E15 + (896/85) E17 - ....

  -Arc tan E + ( E/3/(1+E2) ) (2E2+3) = (2/15) E5 - (4/21) E7 + (2/9) E9 - (8/33) E11 + (10/39) E13 - (4/15) E15 + ...
   Arc tan E - (E/30) ( 25 + (5-9.E4) / (1+E2)2 ) = -(8/105) E7 + (8/45) E9 - (16/55)  E11 + (16/39) E13 - (8/15) E15 + (56/85) E17 - ....

-For the Earth,  E2 = (b2 - c2) / c2 ~  0.0067  so these terms are sufficient. More terms should be used for larger E-values.
 
 
      STACK        INPUTS      OUTPUTS
      Level 3             /             ga
      Level 2              /             gb
      Level 1             /             gc

       where the gravities g are expressed  in  m/s2

Example:

             GABC   >>>>   ga = 9.780379417  m/s2
                                       gb = 9.780274111  m/s2
                                       gc = 9.832185874  m/s2

Note:

-"GABC" does not check that  a >= b >= c
 

Gravity on a triaxial ellipsoid
 

-At the surface of the ellipsoid ( altitude h = 0 ), the gravity at a point whose longitude = L and latitude = B is given by

   g(0)  = ( a ga Cos2 L' Cos2 B + b gb Sin2 L' Cos2 B + c gc Sin2 B ) / d

  with      d2 = a2 Cos2 L' Cos2 B + b2 Sin2 L' Cos2 B + c2 Sin2 B    and     L' = L + 14°92911

  ( The longitude of the major equatorial axis is  14°92911 W )
 

-If  h # 0, an approximate formula is used:

   g(h) = g(0) [ 1 - 2 (h/a') ( 1 + f + m - 2 f Sin2 B ) + 3 sign(h) (h2/a'2) ]       with    a' = (a+b)/2   and  m = a.b.c w2 / GM
 

    w  = 7.292115 E-5 rd/s  = angular velocity of the Earth
  GM = 3.986004419 E14 m3/s2 = Earth gravitational constant
 
 
      STACK        INPUTS      OUTPUTS
       Level 3         L ( ° ' " )            /
       Level 2         B ( ° ' " )            /
       Level 1          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         GRV  >>>>   g(67) = 9.800516274  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          GRV   >>>>   g(1706) = 9.790660011  m/s2
 

Notes:

-The longitudes are positive East.
-Latitudes & longitudes are both supposed geodetic, but using a geocentric longitude will not change the result significantly.

-The values of  ga , gb , gc   are obtained by GABC that is called by GRV

-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 spheroidal model of the Earth.

-The 3 parameters ga , gb , gc may also be chosen differently to get a better fit, but they should verify the Pizzetti theorem:

     ga/a + gb/b + gc/c = 3.GM / (abc) - 2 w2
 

Rectangular Coord. -> Jacobi's Coordinates( u , v )
 

-The rectangular coordinates of a point on the ellipsoid      x2 / a2 + y2 / b2 + z2 / c2 = 1    ( a > b > c )   are related to Jacobi's coordinates  ( u , v )  by:

    x =  [ a2 / ( a2 - c2 ) ]1/2 Cos u  ( a2 - b2 Sin2 v - c2 Cos2 v )1/2
    y =    b  Sin u  Cos v
    z =  [ c2 / ( a2 - c2 ) ]1/2 Sin v  ( a2 Sin2 u + b2 Cos2 u - c2 )1/2

-In the conversion  ( x , y , z )  >>>  ( u , v ) , we have to solve the equation:

   ( b2 - b.c ) Sin4 v + [ b.y2 + b.c - b2 - a.y2 - b.(a-c)/c z2 ] Sin2 v + b.(a-c)/c  z2 = 0
 
 
      STACK        INPUTS      OUTPUTS
       Level 3             X            /
       Level 2             Y       U ( deg )
       Level 1             Z       V ( deg ) 

    U & V are given in decimal degrees

Example:     X = 4398.916449   Y = 3822.64999964   Z = 2583.13552679  ( found in the next paragraph )

      4398.916449       ENTER
      3822.64999964   ENTER
      2583.13552679   XYZ->UV  >>>>   U = 41°               ( with a very small decimal part )
                                                                 V = 24°

Note:

 XYZ->UV  does not check that   x2 / a2 + y2 / b2 + z2 / c2 = 1
 

Jacobi's Coordinates( u , v ) -> Rectangular Coordinates
 

-The same formulas are used for this conversion
 
 
      STACK        INPUTS      OUTPUTS
       Level 3             /            X
       Level 2        U ( deg )            Y
       Level 1        V ( deg )             Z 

    U & V are given in decimal degrees

Example:     U = 41°    V = 24°

         41  ENTER
         24  UV->XYZ   >>>>    X = 4398.916449
                                                Y = 3822.64999964
                                                Z =  2583.13552679
 

Jacobi's Equation
 

-First, we have to determine the parameter  Beta  such that    ( § = Integral )

    §u1u2  ( a2 Sin2 u + b2 Cos2 u )1/2  ( a2 Sin2 u + b2 Cos2 u - c2 ) -1/2 ( (a2-b2) Sin2 u + Beta ) -1/2  du

  =  +/-  §v1v2  ( b2 Sin2 v + c2 Cos2 v )1/2  ( a2 - b2 Sin2 v - c2 Cos2 v ) -1/2 ( (b2-c2) Cos2 v - Beta ) -1/2  dv

-Then, the geodesic distance s between the 2 points may be computed by 2 formulae:

   s =  §u1u2  [ P + R ( dv/du )2 ]1/2  du          (1)

   s =  §v1v2  [ P ( du/dv )2 +R ]1/2  dv           (2)

  where the geodesic parameters are given by

  P = ( x / u )2 + ( y / u )2 + ( z / u )2
  R = ( x / v )2 + ( y / v )2 + ( z / v )2

  ( the other geodesic parameter Q = ( x / u )( x / v ) + ( y / u )( y / v ) + ( z / u )( z / v ) = 0 here )

-It yields:

   P = ( a2 / (a2-c2) ) Sin2 u ( a2 - b2 Sin2 v - c2 Cos2 v ) + b2 Cos2 u Cos2 v + ( c2 (a2-b2)2 / (a2-c2) ) Sin2 u Cos2 u Sin2 v ( a2 Sin2 u + b2 Cos2 u - c2 ) -1

   R = ( c2 / (a2-c2) ) Cos2 v ( a2 Sin2 u + b2 Cos2 u - c2 ) + b2 Sin2 u Sin2 v + ( a2 (b2-c2)2 / (a2-c2) ) Sin2 v Cos2 v Cos2 u  ( a2 - b2 Sin2 v - c2 Cos2 v ) -1
 

-It may also happen that the geodesic looks like this:

                                                          *
                                 *                                               *
                     *                                                                     *
              *                                                                         Paris
         *
  Washington
 

-In other words, the derivative  dv / du  changes its sign

-So, in this case, we have to compute   §v1v0  f(v)  dv - §v0v2  f(v)  dv

  with  f(v) = ( b2 Sin2 v + c2 Cos2 v )1/2  ( a2 - b2 Sin2 v - c2 Cos2 v ) -1/2 ( (b2-c2) Cos2 v - Beta ) -1/2

  where the denominator  ( (b2-c2) Cos2 v - Beta )  vanishes between v1 & v2  i-e  ( (b2-c2) Cos2 v0 - Beta ) = 0

-In order to remove this singularity without producing another singularity for v = 0, I've made the following change of variable:

   Cos w =   [ ( Cos2 v - K ) / ( 1 - K ) ] 1/2              with   K = beta / ( b2 - c2 )
    Sin w =  Sin v /  ( 1 - K ) 1/2

-The integral becomes:    ( b2-c2 ) -1/2 §w1w2  g(w)  dw

   with  g(w) = [ b2 - ( b2-c2 ).{ K + ( 1-K) Cos2 w } ]1/2  [ ( a2-b2 ) - ( b2-c2 ).{ K + ( 1-K) Cos2 w } ] -1/2  [ K + ( 1-K) Cos2 w } ] -1/2
 

>>> Flags F01 & F04 determine the form of the geodesic between the 2 points:
 

                       *                                                     *                                          *                            *
                 *                                                                *                             *             *                       *
            *                                                                         *                    *                      *                       *                                *
         *                                                                               *               *                                                       *                    *
       *                                                                                    *           *                                                                   *

      1   CF      4  CF                                                1 SF   4  CF             1  CF   4  SF                                1  SF    4  SF

-Of course, it's not always obvious to know the good case in advance.
-So, you'll probably have to make several tests before finding the right configuration.
 
 
      STACK        INPUTS      OUTPUTS
       Level 1              /         BETA 

    BETA is given in decimal degrees

Example1:     Calculate  BETA  for the ellipsoid  x2 / 41 + y2 / 37 + z2 / 35 = 1   and the 2 points:

  u1 = +10°          u2 = +75°         ( all in decimal degrees )
  v1 = -15°          v2 = +61°

-Store  sqrt(41)   sqrt(37)    sqrt(35)   in the variables  'A'  'B'  'C'  respectively  ( in the parent directory )
-Store 10 , -15 , 75 , 61  in the variables 'U1'  'V1'  'U2'  'V2'  respectively

-If you choose N = 1  ,  store  1  in the variable 'N'   ( number of subintervals for Gauss-Legendre 16-point formula )

>>>    1  CF  since ( u1 -  u2 ) & ( v1 -  v2 )  have the same sign  and if we try   4  CF

-With the initial guess for beta, BETA = 0.1  , store this value in the variable 'BETA'
-Go to the SOLVR menu and solve for BETA, it yields:

         BETA = 0.19865402015

-Trying  N = 2  returns the same value

-Go to the next paragraph to get the geodesic distance.
 

Example2:    Calculate BETA for the geodesic 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

-Assuming the longitudes are geocentric and the latitudes are geodetic , GGR->XYZ  &  XYZ->UV give the following Jacobi parameters:

    u1 = -62°1615552526        u2 = 33°4252270445           store these numbers in the variables 'U1' 'V1' 'U2' 'V2'
    v1 = +38°8438199514       v2 = -33°8883727534

>>>    Set flag F01    1  SF   since ( u1 -  u2 ).( v1 -  v2 ) < 0

-You will find, with  N = 1   'N'   STO    and   4  CF   that  BETA = 136050.686615

-With N = 2 , we get the same value

-Go to the next paragraph to get the geodesic distance.
 

Example3:      Calculate BETA for tthe geodesic between the U.S. Naval Observatory at Washington (D.C.)

     L1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N  and the "Observatoire de Paris"  L2 = 2°20'13"8 E ,  b2 = 48°50'11"2 N

-The usual routines GGR->XYZ  &  XYZ->UV  return the Jacobi parameters

       u1 = -62°1615552526        u2 = 17°300852295           store these numbers in the variables 'U1' 'V1' 'U2' 'V2'
       v1 = +38°8438199514       v2 = 48°8377638099

-The geodesic looks like this:
 

                                 *
                     *                  *                                which suggests   1  CF   4  SF
               *                             *
          *                                     *
       *
 

-With  N = 1 ,  1  CF   4 SF ,  it yields   BETA = 101353.551826

-  N = 2  returns  BETA = 101353.551827

 ( We can keep  N = 1 )
 

-Go to the next paragraph to get the geodesic distance.
 

Geodesic Distance on a Triaxial Ellipsoid
 

-Now,  BETA  is known, so we can try to evaluate the integral that gives  s

 DGT  uses the classical Runge-Kutta method of order 4 to solve - at each step - the differential equation,
 and the geodesic distance itself is obtained by the Newton-Côtes 7-point formula
 
 
      STACK        INPUTS      OUTPUTS
       Level 3              /        Distance
       Level 2              /            w
       Level 1              /            w'

  Where  w & w'  are 2 angles which would be equal if the precision were perfect

Example1:     Calculate the geodesic distance on the ellipsoid  x2 / 41 + y2 / 37 + z2 / 35 = 1   between the 2 points:

  u1 = +10°          u2 = +75°
  v1 = -15°          v2 = +61°

-Assuming BETA is known as explained in the paragraph above ( do not change flags 1 & 4 )
  we choose different N-values for the Newton-Côtes formula  ( number of sub-intervals )

-With  N = 8 ,  DGT  gives  8.59482326321         w = 67°1590005932       w' = 67°1589920238
-With N = 16   DGT  -----  8.59482263395         w = 67°1590005932       w' = 67°1590000556
-With N = 32   DGT  -----  8.59482258246         w = 67°1590005932       w' = 67°1590005568                in 4mn56s  ( HP48GX )

>>> So, we can trust the last result:    D = 8.594822582     rounded to 9decimals
 

Example2:    Calculate the geodesic 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

-Assuming BETA is known as explained in the paragraph above ( do not change flags 1 & 4 )
  we choose different N-values for the Newton-Côtes formula  ( number of sub-intervals )

-With  N = 8 ,  DGT  gives  12709.5646011 km         w = -52°0771446804       w' = -52°077144092
-With N = 16   DGT  -----  12709.5645850 km         w = -52°0771446804       w' = -52°0771446429                 in 2mn42s  ( HP48GX )
-With N = 32   DGT  -----  12709.5645839 km         w = -52°0771446804       w' = -52°0771446774                 in 5mn24s  ( HP48GX )

>>> So, the last result  D = 12709.5645839 km    is probably ( almost ) exact.
 

Example3:   Calculate the 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 the "Observatoire de Paris"  L2 = 2°20'13"8 E ,  b2 = 48°50'11"2 N

-Assuming BETA is known as explained in the paragraph above ( do not change flags 1 & 4 )
 we choose different N-values for the Newton-Côtes formula  ( number of sub-intervals )

-With  N = 8 ,  DGT  gives  6181.62547507 km         w = 108°085278506       w' = 108°085278435
-With N = 16   DGT  -----  6181.62547563 km         w = 108°085278506       w' = 108°085278502                 in 2mn28s  ( HP48GX )

>>> So, the last result  D = 6181.62547563 km    is probably almost exact !
 

Notes:

-Obviously, these programs are neither fast nor easy to use !
-Moreover, they will work in most cases, but not in all cases
-Fortunately, for the Earth, an approximate method may be used to get satisfactory results, except for nearly antipodal points.
 

Geodesic Distance on a Triaxial Ellipsoid - Approximate Method
 

-This program does not use Jacobi's equation:

-The center O of the ellipsoid and 2 points M & N define a plane ( at least if they are not aligned ) that intersects the ellipsoid.
  and we can compute the distance D0 of the arc of ellipse MN thus defined.
-This will not be - in general - the geodesic distance, but it will remain slightly larger if the flattenings are small.

-Better approximations are also computed:

-Let  W = the angle ( OM , ON )

-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(µ) that minimizes the integral below under the form:

          h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) + h3 Sin ( 720° µ / W )           ( in fact 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.

  where   T = ( A2 + B2 + C2 )   and

    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 are:

  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 } ]
 

-Let  D (h1,h2,h3)  the distance as a function of  h1 , h2 , h3

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

  D0 =  A = f(0,0,0)   B = f(H,0,0)   C = f(-H,0,0)    D= f(0,H,0)   E = f(0,-H,0)   F= f(0,0,H)   G = f(0,0,-H)

-Then,

    D1 = A - (B-C)2 / [ 8.(B+C-2.A) ]
    D2 = A - (B-C)2 / [ 8.(B+C-2.A) ] - (D-E)2 / [ 8.(D+E-2.A) ]
    D3 = 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) ]

-Though these formulas are quite elementary and though we simply add the different corrections for the extremum, the results are surprisingly good !

-In fact DDD is called as a subroutine by the following programs according to the types of longitudes & latitudes:
-See the paragraphs below for numerical examples.
 

Geographic Coordinates -> Geodesic Distance
 

-Here, we assume that the longitudes are geocentric and the latitudes are geodetic
-For the Gauss-Lendre 16-point formula, N = 1 is sufficient with the Earth
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")      D0 ( km )
       Level 3     Lat1 (° ' ")      D1 ( km )
       Level 2    Long2 (° ' ")      D2 ( km )
       Level 1     Lat2 (° ' ")      D3 ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N  ( assuming L is geocentric and b is geodetic )

   -77.0356    ENTER               D0 = 6181.62809236 km
   38.55172   ENTER               D1 = 6181.62550491 km
   2.20138     ENTER               D2 = 6181.62550389 km
   48.50112   GGRD    >>>>   D3 = 6181.62547937 km                                  ---Execution time = 2mn28s---     ( HP48GX )
 

-Compared to the "exact" value given by the Jacobi's method, i-e D = 6181.62547563 km,  D3-error is inferior to 4 mm !

Notes:

-The accuracy of D0 is of the same order as Andoyer's 1st order formula, i-e about 60 meters provided the points are not nearly antipodal.
-Gauss-Legendre 16-point formula could be replaced by a simpler method to get faster results.
-The result should always verify  D3 <= D2 <= D1 <= D0
-If the eccentricities are large, the precision is not so good. Moreover, N must be choosen > 1
-  COEF = 1000  seems very good for the Earth, but other choices may be better for other ellipsoids
 

Geodetic Coordinates -> Geodesic Distance
 

-In fact, I don't really know if the coordinates we find in books or on the net are geocentric or geodetic for the longitudes
-Anyway, if the longitudes and the latitudes are geodetic, use  GDTD
-For the Gauss-Lendre 16-point formula, N = 1 is sufficient with the Earth
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")      D0 ( km )
       Level 3     Lat1 (° ' ")      D1 ( km )
       Level 2    Long2 (° ' ")      D2 ( km )
       Level 1     Lat2 (° ' ")      D3 ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N  ( assuming L & b are geodetic )

   -77.0356    ENTER               D0 = 6181.56896411 km
   38.55172   ENTER               D1 = 6181.56637672 km
   2.20138     ENTER               D2 = 6181.56637570 km
   48.50112   GDTD    >>>>   D3 = 6181.56635118 km                                  ---Execution time = 2mn28s---     ( HP48GX )

-Thus, we can take geodesic distance = 6181.56635(118) km

Note:

-Logically, this result is significantly different from the value given by  GGRD
 

Geocentric Coordinates -> Geodesic Distance
 

-If the longitudes & latitudes are geocentric, use  GECD
-For the Gauss-Lendre 16-point formula, N = 1 is sufficient with the Earth
 
 
      STACK        INPUTS      OUTPUTS
       Level 4    Long1 (° ' ")      D0 ( km )
       Level 3     Lat1 (° ' ")      D1 ( km )
       Level 2    Long2 (° ' ")      D2 ( km )
       Level 1     Lat2 (° ' ")      D3 ( km )

Example:         Calculate the 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 the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N  ( assuming L & b are geocentric... which is surely not true in this case )

   -77.0356    ENTER               D0 = 6160.79419773 km
   38.55172   ENTER               D1 = 6160.79164579 km
   2.20138     ENTER               D2 = 6160.79164480 km
   48.50112   GECD    >>>>   D3 = 6160.79162057 km                                  ---Execution time = 2mn28s---     ( HP48GX )

-Thus, we can take geodesic distance = 6160.79162(057) km
 

Note:

-Still logically, this result is even more different from the value given by  GGRD
 

Geographic Coordinates -> Rectangular Coordinates
 

-If the longitude l is geocentric i-e equals the angle between the meridian plane passing through the observer and the Prime Meridian plane,
 and if the latitude b is the angle between the local vertical in the local meridian plane and the plane of the Equator,
 the rectangular coordinates x , y , z of a point P on the triaxial ellipsoid may be found as follows:

  r2 = b2 + ( a2-b2 ) / [ 1 + (a2/b2) Tan2 L ]                with   r = OQ   where  Q = point of intersection of the equator and the local meridian

  R2 = c2 + ( r2-c2 ) / [ 1 + (c2/r2) Tan2 b ]                  and  L = l - L0

    x = R Cos µ  Cos L
    y = R Cos µ  Sin L                        where    R = OP  and  Tan µ = (c2/r2) Tan b
    z = R Sin µ
 
 
      STACK        INPUTS      OUTPUTS
       Level 3             /            x
       Level 2    Long ( ° ' "  )            y
       Level 1     Lat ( ° ' "  )            z 

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

    2.20138   ENTER
   48.50112  GGR->XYZ     x = 4016.62586787 km
                                            y = 1248.44665689 km
                                            z = 4778.59664418 km

Notes:

-Strictly speaking, the perpendicular to the local meridian plane is not perpendicular to the triaxial ellipsoid,
  but the difference is negligible: both versions return almost the same z-values
-This would not be the case if the flattening were large.
 

Geodetic Coordinates -> Rectangular 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:

    x = w Cos b Cos l
    y = w ( 1 - e2 ) Cos b Sin l                       In fact   l - L0  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
 
 
      STACK        INPUTS      OUTPUTS
       Level 3             /            x
       Level 2    Long ( ° ' "  )            y
       Level 1     Lat ( ° ' "  )            z 

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

    2.20138   ENTER
   48.50112  GDT->XYZ     x = 4016.63356045 km
                                            y = 1248.42190799 km
                                            z = 4778.59664415 km
 
 

Gauss-Legendre 16-point Formula
 

 GL16 employs Gauss-Legendre 16-point formula to evaluate the integral to find Jacobi's beta parameter.
 It may also be used independently.
 
 
      STACK        INPUTS      OUTPUTS
       Level 5             /        << f >>
       Level 4        << f >>             a
       Level 3             a             b
       Level 2             b             N
       Level 1             N             I

    where  I  is an approximation of   §ab  f(x) dx
              N = the number of subintervals to apply Gauss formula
     and  << f >>  is a program that takes x in level 1 and returns f(x) in level 1

Example:      Evaluate  §03  exp ( - x2 ) dx

-With N = 1

       <<  SQ  NEG   EXP  >>   ENTER
                        0                       ENTER
                        3                       ENTER
                        1                       GL16             >>>>      0.886207348259

-With   N = 2      DROP   DROP   2      GL16             >>>>      0.886207348261

Notes:

-The built-in integrator returns  0.886207348259
-Of course, you can replace GL16 by your own program provided it has the same specifications
 

Spherical -> Rectangular Conversion
 

-This program employs the usual formulae ( for the astronomers ):

       x = r cos b cos L
       y = r cos b sin L
       z = r sin b
 
 
      STACK        INPUTS      OUTPUTS
       Level 3             r            x
       Level 2            L            y
       Level 1            b            z 

Example:     r = 7   L = 124°   b = 41°

                        7   ENTER
                      124 ENTER
                       41  S->R      >>>>   x = -2.95419769010
                                                        y =  4.37977818860
                                                        z =  4.59241320294

Notes:

- S->R  works in all angular mode
-But in DEG mode, the longitude L & the latitude b must be expressed in decimal degrees
 
 

References:

[1]  Jean Meeus - "Astronomical Algorithms" - Willmann-Bell -  ISBN 0-943396-35-2
[2]  Jacqueline Lelong-Ferrand - "Geometrie Differentielle" - Masson - ( in French )
[3]  Vincenty's formula: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
[4]  Paul D. Thomas - "Spheroidal Geodesics, Reference Systems & Local Geometry" - US Naval Oceanographic Office.
[5]  Paul D. Thomas - "Mathematical Models for Navigation Systems" - TR-182 -  US Naval Oceanographic Office.
[6]  http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
[7]  Jacobi - "De la ligne géodésique sur un ellipsoïde, et des différents usages d’une transformation analytique remarquable"
[8]  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
[9]  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

[10]  If you want to visualize the geodesics on a triaxial ellipsoid, go to this excellent page

[11]  Michèle Caputo - "Some Space Gravity Formulas and the Dimensions and the Mass of the Earth"