Ellipsoid-Plane Intersection for the HP-41

Overview
 

 1°)  General Case
 2°)  Great Ellipse defined by 2 points on the Earth.
 

-The intersection of an ellipsoid and a plane is an ellipse ( or a point or the empty set )
-"EPI" calculates the elements of this ellipse, assuming the plane & the ellipsoid are given by

    (P) :  m x + n y + p z = q
    (E) :  x² / a² + y² / b² + z² / c² = 1

-The second program computes the elements of the great ellipse defined by 2 points on the Earth.
 

1°)  General Case
 

 "EPI" employs the method described in reference [1]
 
 
 

Data Registers:           •  R00 = q                                   ( Registers R00 thru R06 are to be initialized before executing "EPI" )

                                      •  R01 = m     •  R04 = a               R07-R08-R09 = f₀           R16-R17-R18 = Vertex1         R25-R26-R27 = Vertex4
                                      •  R02 = n      •  R05 = b               R10-R11-R12 = f₁           R19-R20-R21 = Vertex2         R28 = t₀
                                      •  R03 = p      •  R06 = c               R13-R14-R15 = f₂           R22-R23-R24 = Vertex3         R29 & R30 = semi-axes
Flags: /
Subroutines: /
 
 
 
 

 
     ( 296 bytes / SIZE 031 )
 
 

  Where A & B are the semi-axes of the ellipse.

Example:     (P) : 2 x + 3 y + 4 z = 5    (E) :  x² / 25 + y² / 16 + z² / 9 = 1

      5  STO 00

      2  STO 01       5  STO 04
      3  STO 02       4  STO 05
      4  STO 03       3  STO 06

      XEQ "EPI"  >>>>   A = 4.547809219 = R29                                                             ---Execution time = 15s---
                         X<>Y  B = 3.374481476 = R30
 

  And we find in registers R07 to R28  ( rounded to 4 decimals )

     R07-R08-R09 = ( 0.6443 ,  0.6186 , 0.4639 ) = f₀ =  center of the ellipse
     R10-R11-R12 = ( 3.7153 , -2.4769 ,  0 )        = f₁
     R13-R14-R15 = ( 1.8862 , 1.8107 , -2.3011 ) = f₂

-So any point of the ellipse may be obtained by the formula:  f₀ +  f₁ Cos t + f₂ Sin t      where  t belongs to [0°,360°[

  And the 4 vertices are in registers  R16 to R27:

     R16-R17-R18 = ( 4.7415 , -1.2447 , -0.1872 ) = V₁
     R19-R20-R21 = ( 1.4021 , 3.0561 , -1.7432 )  = V₂
     R22-R23-R24 = ( -0.1135 , -1.8190 , 2.6710 ) = V₃
     R25-R26-R27 = ( -3.4528 , 2.4818 , 1.1150 )  = V₄

-In register R28 we have the value of  t  which corresponds to the 1st vertex, i-e    t = 16°43698070
 

Notes:

-You can use LBL 10 to compute the coordinates of a point, for a given t-value.
-For example, with t = 41°  ( t in decimal degree )

    41  XEQ 10  returns  x = +4.6857
                          RDN   y = -0.0628
                          RDN   z = -1.0457

-There will be a data error message if the intersection is empty.
 

2°)  Great Ellipse defined by 2 points on the Earth
 

-Given 2 points M & N at the surface of the Earth by their geodetic coordinates, "GEL" returns the semi-major and semi-minor axes of the great ellipse
 passing by these 2 points and the center of the Earth.

-The coordinates of M & N on the great ellipse ( the new axes = the axes of the ellipse ) are also computed
 
 

Data Registers:    R00 =  t₀     R14 = 6378.172   R15 = 6378.102   R16 = 6356.752314

                          -> R01 = x_M   R03 = x_N      R05-R06-R07 = X_M Y_M Z_M   R11-R12-R13: temp    R23-R24-R25 = U = unit vector of the 1st vertex
                          -> R02 = y_M   R04 = y_N      R08-R09-R10 = X_N Y_N Z_N     R17 to R22: temp         R26-R27-R28 = V = unit vector of the 2nd vertex

                              R29 = A = semi-major axis     R30 = B = semi-minor axis

Flags: /
Subroutines: /
 
 
 

 
     ( 515 bytes / SIZE 031 )
 
 

  Where A & B are the semi-axes of the ellipse.

Example:   Calculate the elements of the great ellipse defined by the U.S. Naval Observatory at Washington (D.C.) :
 
                   L₁ = 77°03'56"0 W ;  b₁ = 38°55'17"2 N  and the Observatoire de Paris:  L₂ = 2°20'13"8 E ;  b₂ = 48°50'11"2 N

   -77.0356    ENTER^
    38.55172  ENTER^
     2.20138   ENTER^
    48.50112  XEQ "GEL"   >>>>  a' =  6378.104423 km                            ---Execution time = 24s---
                                          X<>Y  b' = 6364.794129 km
 

-The coordinates of the U.S. Naval Observatory at Washington (D.C.) are in R05-R06-R07:   ( 2322.282431 , -4392.706835 , 3985.543384 )

 and with respect to ( U , V ) are in R01 & R02:    ( -3881.641327 , -5050.374807 )

-The coordinates of the Observatoire de Paris Observatoire de Paris are in R08-R09-R10:   ( 4016.636531 , 1248.414112 , 4778.596908 )

  and with respect to ( U , V ) are in R03-R04:    ( 1976.680334 , -6051.415542 )

-The unit vector corresponding to the major axis is R23-R24-R25 = U = ( 0.186195902 , 0.982512530 , 0.000463853 )
-The unit vector corresponding to the minor axis is R26-R27-R28 = V = ( -0.602931120 , 0.114634184 , -0.789514451 )
 

Notes:

-This program employs a triaxial ellipsoidal model of the Earth with a = 6378.172 km, b = 6378.102 km , c = 6356.752314 km
 and longitude of the major axis = -14.929°
-Change the corresponding lines ( 47-44-40-07 ) if you want to use other constants.

-The longitudes are positive East.
 
 

Reference:

[1]  https://en.wikipedia.org/wiki/Ellipsoid