-ASTRO2010   Manual


 Overview
 

-The first 3 functions ( KG , OB , KPL ) are M-Code routines.
-The other ones are focal programs.

-The origin of time is 2000/01/01 0h TT ( 12 hours before J2000.0 )
-Flag F04 must be cleared for Gregorian dates
-Flag F04 must be set for Julian dates.
-The longitudes are measured positively Eastward from the meridian of Greenwich
 
 
 
XROM  Function  Desciption
 06,00
 06,01
 06,02
 06,03
 06,04
 06,05
 06,06
 06,07
 06,08
 06,09
 06,10
 06,11
 06,12
 06,13
 06,14
 06,15
 06,16
 06,17
 06,18
 06,19
 06,20
 06,21
 06,22
 06,23
 06,24
 06,25
 06,26
 06,27
 06,28
 06,29
 06,30
 06,31
 06,32
 06,33
 06,34
 06,35
 06,36
 06,37
 06,38
 06,39
 06,40
 06,41
 06,42
 06,43
 06,44
 06,45
 06,46
 06,47
 06;48
 06,49
 06,50
 06,51
 06,52
 06,53
 06,54
 06,55
 06,56
 06,57
 06,58
 06,59
 06,60
 06,61
 06,62
 06,63
 -ASTRO2010
  KG
  OB
  KPL
 "J2000"
 "DT"
 "EE"
 "R-S"
 "S-R"
 "MST"
 -EPHEMERIS
 "EPH"
 "SUN"
 "MO"
 "ME"
 "VE"
 "MA"
 "JU"
 "SA"
 "UR"
 "NE"
 "PL"
 "N"
 "K"
 "O"
 "L"
 -SATELLITES
 "SATJSU"
 -COMETS
 "PHEM"
 "GAUSS"
 -TRANSCOOR
 "EQ-ECL"
 "ECL-EQ"
 "EQ-AZ"
 "AZ-EQ"
 "PREC"
 "PREQ"
 "IOO"
-REFRACTION
 "H-H0"
 "H0-H"
 "REFR"
-MISCELLAN
 "BOSS"
 "BSTAR"
 "COSMO"
 "DELTAT"
 "DT-KY"
 "EASTER"
 "ECLIPSE"
 "FIX"
 "GNB+1"
 "GNBP"
 "IFMD"
 "KY-DT"
 "LST"
 "MNODE"
 "NUT"
 "OBL"
 "PHASE"
 "RTS"
 "STYLM"
 "ZODIAC"
 Section Header
 Gaussian Gravitational Constant
 Obliquity of the Ecliptic on 2000/01/01 0h
 Kepler's Equation
 Date >>> Number of days since 2000/01/01 0h
 Number of days since 2000/01/01 0h >>> Date
 Equatorial <> Ecliptic ( subroutine )
 Rectangular >>> Spherical
 Spherical >>> Rectangular
 Mean Sidereal Time at Greenwich
 Section Header
 Initialization
 The Sun
 The Moon
 Mercury
 Venus
 Mars
 Jupiter
 Saturn
 Uranus
 Neptune
 Pluto
 "L" "N" "O" "K"
 are 4 subroutines that are
 used to compute the position of celestial bodies.
 "N" may be useful for itself, however.
 Section Header
 Satellites of Jupiter(4), Saturn(7) & Uranus(5)
 Section Header
 Parabolic, Hyperbolic, Elliptic Motions
 Preliminary Orbit Determination by Gauss' Method
 Section Header
 Equatorial > Ecliptic
 Ecliptic > Equatorial
 Equatorial > Azimuthal
 Azimuthal > Equatorial
 Precession ( ecliptic coordinates )
 Precession ( equatorial coordinates )
 Orbital elements from one equinox to another one
 Section Header
 True Altitude > Apparent Altitude ( mean conditions
 Apparent Altitude > True Altitude   of temp & press )
 A more general refraction program
 Section Header
 Barycenter of the Solar System
 Binary Stars
 Redshift > Cosmological Distances
 Calculation of Delta T = TT-UT for a given year
 Julian/Gregorian Calendar > Hindu Calendar
 Date of Easter
 Solar & Lunar Eclipses
 Positional fix from the heights of 2 stars
 Gravitational (n+1)-body Problem
 Gravitational n-body Problem
 Illuminated Fraction of the Moon's Disk
 Hindu Calendar > Julian/Gregorian Calendar
 Local Solar Time
 Passage of the Moon through the Nodes
 Nutation ( simplified formulae )
 Obliquity of the Ecliptic
 Phases of the Moon
 Rising, Transit, Setting
 Sidereal & Tropical years and Lunar Months
 Date & Time for a given longitude of the Sun

 
 

Gaussian Gravitational Constant
 

 KG  simply returns  k = 0.01720209895  in X-register.
 

Obliquity of the Ecliptic on 2000/01/01 0h
 

 Likewise,  OB  returns  23.4392796  in X-register,
 which is the mean obliquity of the ecliptic ( expressed in degrees ) on 2000/01/01 0h  TT
 

 Kepler's Equation
 

 This M-Code routine uses Newton's method to solve Kepler's equation:

      M = E - e sin E     if  0 <= e < 1       ( elliptic case )
      M = e sinh E - E   if      e > 1        ( hyperbolic case )

               E = eccentric anomaly,
 where    M = mean anomaly
               e = eccenticity
 
 
      STACK        INPUTS      OUTPUTS
           Y             e            e
           X            M            E
           L             /            M

Examples:

  e  = 1.2     ENTER^
 M = 0.4     XEQ "KPL"  returns   E = 0.987853767

  e  = 0.9     ENTER^
 M =  7°      XEQ "KPL"  returns   E = 40.46693453°

Notes:

-When e < 1 , the angles are expressed in degrees, whatever the angular mode is.
-When e < 1 , first take M modulo 360 to get a better accuracy.
-In some pathological cases - like M = 7° , e = 0.999 - more than 40 iterations are needed before the solution ( 52°27026153 ) is found.
 

 Date >>> Number of days since 2000/01/01 0h
 

-"J2000"  takes a date expressed in YMD format ( Gregorian Calendar if CF 04 , Julian Calendar if SF 04 )
  and calculate the number of days since 2000/01/01 at 0h. ( i-e the Julian Day number minus 2451544.5 )
-Z & T-registers are lost.
-Y register is saved in Y & Z registers, and the day of the week is returned in T-register ( Sunday = 0 , Monday = 1 , Tuesday = 2 , .... , Saturday = 6 ).
 

Data Registers: /
Flag:  F04             SF04 for the Julian Calendar
                             CF04 for the Gregorian Calendar
Subroutines: /
 
 
 
         STACK        INPUTS       OUTPUTS
             T             /           dow
             Z             /             Y
             Y            Y             Y
             X  YYYY.MNDDdd             N

Examples:

  CF 04  2134.0404 ( Gregorian )  XEQ "J2000"   >>>     N =    49036   RDN   RDN   gives  0 = Sunday
  SF 04  1234.0428     ( Julian )     XEQ "J2000"   >>>     N =  -279651  RDN   RDN   gives  5 = Friday

Notes:

-The days may have decimals
-You can choose the Julian Calendar after 1582/10/15
  and you can choose the Gregorian Calendar before 1582/10/04
-This program works for all dates: positive, 0 or negative.
-The "B.C." years are counted astronomically:
     year  0 = 1 B.C.
     year -1 = 2 B.C. ...etc...
 

 Number of days since 2000/01/01 0h >>> Date
 

 "DT"  performs the reciprocal conversion.
 

Data Registers: /
Flag:  F04             SF04 for the Julian Calendar
                             CF04 for the Gregorian Calendar
Subroutines: /
 
 
 
        STACK         INPUTS       OUTPUTS
            T              /             Y
            Z              /             Y
            Y             Y             Y
            X             N  YYYY.MNDDdd

Examples:

  CF 04  N =    49036   XEQ "DT"  >>>   2134.0404
  SF 04  N = -279651   XEQ "DT"  >>>   1234.0428

Note:

-"J2000" & "DT" are separated by a RTN, so you can easily use them to convert a date in the Julian Calendar into a date in the Gregorian Calendar.
-For instance:   SF 04   1234.0428  XEQ "J2000"  >>>  -279651   CF 04  R/S  >>>  1234.0505
-Thus  1234 April 4th Julian = 1234 May 5th Gregorian.
 

 Equatorial <> Ecliptic ( subroutine )
 

-Many transformations of coordinates use the same type of formulae which appear in the equatorial-ecliptic conversion, namely:

                       sin  b = cos e  sin d - sin e  cos d  sin a
             cos cos l = cos d  cos a
             cos b  sin  l = sin e  sin d  +  cos e  cos d  sin a

-So, "EE" is extensively used in the "TRANSCOOR" section.
 

Data Registers: /
Flags:  /
Subroutines: /
 
 
      STACK      INPUTS      OUTPUTS
           Z            e            e
           Y            d            b
           X            a            l

Example:     if  a = 116°328942 ,  d = 28°026183 and  e = 23°4392911 ( in  DEG mode )

    23.4392911  ENTER^
    28.026183    ENTER^
  116.328942    XEQ "EE"  >>>> l  =  113°215630   RDN   b  =  6°684170

Note:

-Like "R-S" and "S-R" , "EE"  works in all angular modes.
 

 Rectangular >>> Spherical
 

                                     x = r cos b cos l
-We use the formulae:   y = r cos b sin l        with  the built-in  R-P  function.
                                    z = r sin b

   where    x , y , z = rectangular coordinates,    r  = ( x2 + y2 + z2 )1/2  ,   l  = longitude ,  b = latitude
 

Data Registers: /
Flags:  /
Subroutines: /
 
 
      STACK       INPUTS     OUTPUTS
           T            T            T
           Z            z            b
           Y            y            l
           X            x            r
           L            /      (x2+y2)1/2

 

Example:     x = 3 ; y = 4 ; z = -7   find the spherical coordinates ( in DEG mode )

    -7  ENTER^
     4  ENTER^
     3  XEQ "R-S"    r  =   8.602325267
             RDN           =  53.13010235°
             RDN         b  = -54.46232221°
 

 "R-S"  works in all angular modes.
 

 Spherical >>> Rectangular
 

-The same formulas are used with P-R
 

Data Registers: /
Flags:  /
Subroutines: /
 
 
      STACK       INPUTS     OUTPUTS
           T            T             T
           Z            b             z
           Y            l             y
           X            r             x
           L            /      (x2+y2)1/2

Example:     r = 10 ; l = 124° ; b = 37°   find    x ; y ; z   ( in DEG mode )

     37  ENTER^
   124  ENTER^
    10  XEQ "S-R"    x = -4.465913097
                    RDN   y =  6.620988446
                    RDN   z =  6.018150232

-"S-R" works in all angular modes
 

 Mean Sidereal Time at Greenwich
 

 "MST" computes the mean sidereal time at Greenwich

    T = the number of millenia since 2000/01/01 0h is stored into R00
 

Data Register:       R00 = Number of millenia since 2000/01/01 0h
Flag: /
Subroutine:  "J2000"
 
 
       STACK        INPUTS      OUTPUTS
           Y   YYYY.MNDD             /
           X     hh.mnss (UT)  HH.MNSS (MST)

Example:   Find the mean sidereal time at Greenwich on 2010/07/16 at 7h41m UT

   2010.0716  ENTER^
         7.41      XEQ "MST"  >>>>   MST =  3h17m10s

  and  R00 = T = 0.01053886417
 

 Initialization
 

  "EPH"  stores your longitude in R16  measured positively eastwards from the meridian of Greenwich.
                        your latitude   in R17  positive = North , negative = South
the value of DeltaT = TT-UT  in R18  expressed in seconds, ( in 2010 , TT-UT is about 66 seconds )

  "EPH" also calculates  MST in R19  and stores the number of millenia since 2000/01/01 0h TT in R00

  Suppose we are at the US Naval Observatory at Washington ( D.C. )

      Longitude = 77°03'56" W = -77°03'56"
       Latitude  =  38°55'17" N = +38°55'17"

  It is 7h41m UT , on July 16th 2010
  And we take TT-UT = 66 seconds
 

  XEQ "EPH"  the HP41 displays "LONG^LAT"

    -77.0356   ENTER^
      38.5517  R/S              next  PROMPT   "TT-UT=SEC?"

          66       R/S              next  PROMPT   "DATE^TIME?"

     2010.0716  ENTER^
           7.41      R/S           the HP-41 returns  the date in Y & the time in X

>>>  R16-R17-R18-R19 are initialized, R00 = 0.01053886626  ( DeltaT has been added to R00 after executing "MST" ) and the HP-41 is set in DEG mode.
 

 We suppose these values remain unchanged in all the examples below
 

-We are now at the beginning of "SUN"

-All the programs "SUN" "MO" "ME" "VE" "MA" "JU" "SA" "UR" "NE" "PL"  return
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

 R03 = geocentric longitude in decimal degrees = l
 R04 = geocentric latitude    ------------------  = b
 R05 = distance to the Earth in Astronomical Units = r ( except for the Moon: R05 = the Moon's parallax in sexagesimal degrees )
 R06 = Right Ascension in  hh.mnss = R.A
 R07 = Declination in          ° ' ''       = Decl
 R08 = heliocentric longitude in decimal degrees = L  ( except for the Sun and the Moon )
 R09 = heliocentric latitude   ------------------- = B   ---------------------------------
 R10 = radius vector in Astronomical Units        = R    ---------------------------------
 R11 = Azimuth in ° ' ''  reckoned clockwise from the South  = Az
 R12 = True Altitude in ° ' '' = h
 R13 = Apparent Altitude in ° ' ''  ( after correction for refraction ) = h0
 R14 = Elongation from the Sun = Elg

 R01 = X & R02 = Y are the rectangular coordinates of the Sun:
"SUN" may be executed first to get the exact geocentric or topocentric coordinates of the planets.

 The effects of aberration and nutation are neglected: we obtain mean coordinates.
 The correction for parallax is neglected too, except for the Moon's altitude ( it may reach 1 degree )

 These ephemerides have an accuracy of about 0.01° in the heliocentric longitudes over the interval 1000-3000.
 ( except for Pluto: its coordinates are obtained with this precision over 1880-2110 only )
 The Sun's coordinates are more accurate. However, errors in the geocentric coordinates are increased when the distance between the Earth and the planet
  is small and in this case, for Venus and Mars, the error in the obtained geocentric longitude can reach 1 arcminute, perhaps a little more.
 All these coordinates are referred to the mean ecliptic and equinox of the date.

Notes:

-If R00 is correctly initialized but R16-R17-R18-R19 are not, the topocentric coordinates will be wrong,
  but the right ascension, declination, heliocentric & geocentric longitudes & latitudes and the distances from the Sun or the Earth will be exact !
 

 The Sun
 

-If  CF 01, all the coordinates described above are computed
-If  SF 01, the stack registers X & Y contain the Sun's radius vector & its longitude, without disturbing R01, R02, ....
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

   CF 01  XEQ "SUN"  yields   Az = 216°27'26"  = R11
                                    RDN     h  =  -20°37'58" = R12
                                    RDN     h0 = -20°37'58" = R13        ( when h0 < 0 , h0 is meaningless )
                                    RDN   Elg = 0  ... logical !

  R03 = l = 113°692               R06 = R.A. = 7h42m14s
  R04 = b = 0                         R07 = Decl = 21°21'37"
  R05 = r = 1.016434 AU

and the rectangular coordinates of the Sun:  R01 = X = -0.408421  &  R02 = Y = 0.930760

  SF 01  XEQ "SUN"  returns  X =  1.016434  &  Y = 113°692
 

 The Moon
 

-If you've just executed "SUN" with CF 01, simply press R/S instead of XEQ "MO"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "MO"  >>>>   Az = 145°13'58"  = R11
                             RDN     h  =  -49°32'56" = R12
                             RDN     h0 = -49°32'56" = R13        ( when h0 < 0 , h0 is meaningless )
                             RDN   Elg =   62°383

  R03 = l = 175°957                             R06 = R.A. =11h37m06s
  R04 = b = -5°073                               R07 = Decl = -3°02'57"
  R05 = parallax = 0°59'48"
 

 Mercury
 
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "ME"  >>>>   Az = 198°51'53"  = R11
                            RDN     h  =  -30°14'17"  = R12
                            RDN     h0 = -30°14'17"  = R13        ( when h0 < 0 , h0 is meaningless )
                            RDN   Elg =   18°528

  R03 = l = 132°157                               R06 = R.A. = 9h00m19s                         R08 = L = 186°609
  R04 = b = 1°556                                  R07 = Decl = 18°38'33"                          R09 = B =    4°687
  R05 = r = 1.19459 AU                                                                                        R10 = R = 0.39700 AU
 

 Venus
 

-If you've just XEQ "ME" , simply press  R/S  instead of   XEQ "VE"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "VE"  >>>>   Az = 171°34'59"  = R11
                            RDN     h  =  -40°29'28" = R12
                            RDN     h0 = -40°29'28" = R13        ( when h0 < 0 , h0 is meaningless )
                            RDN   Elg =   42°846

  R03 = l = 156°525                               R06 = R.A. = 10h34m53s                         R08 = L = 229°255
  R04 = b = 1°184                                  R07 = Decl = 10°12'59"                            R09 = B =   1°570
  R05 = r = 0.96043 AU                                                                                          R10 = R = 0.72394 AU
 

 Mars
 

-If you've just XEQ "VE" , simply press  R/S  instead of   XEQ "MA"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "MA"  >>>>   Az = 151°20'12"  = R11
                            RDN     h  =  -43°12'11"  = R12
                            RDN     h0 = -43°12'11"  = R13        ( when h0 < 0 , h0 is meaningless )
                            RDN   Elg =   58°114

  R03 = l = 171°803                               R06 = R.A. = 11h30m57s                         R08 = L = 204°193
  R04 = b = 0°676                                  R07 = Decl =  3°52'18"                             R09 = B =   0°796
  R05 = r = 1.8976 AU                                                                                            R10 = R = 1.6113 AU
 

 Jupiter
 

-If you've just XEQ "MA" , simply press  R/S  instead of   XEQ "JU"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "JU"   >>>>   Az = 315°49'15"  = R11
                            RDN     h  =   41°41'30" = R12
                            RDN     h0 =  41°42'35" = R13       ...  at last a planet above the horizon !
                            RDN   Elg =   110°370

  R03 = l =   3°316                                  R06 = R.A. =  0h14m20s                          R08 = L =  -7°751
  R04 = b = -1°361                                  R07 = Decl =  0°04'11"                             R09 = B = -1°238
  R05 = r =  4.519 AU                                                                                               R10 = R =  4.965 AU
 

 Saturn
 

-If you've just XEQ "JU" , simply press  R/S  instead of   XEQ "SA"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "SA"   >>>>   Az = 140°54'06"  = R11
                            RDN     h  =  -41°17'24"  = R12
                            RDN     h0 =  -41°17'24" = R13        ( when h0 < 0 , h0 is meaningless )
                            RDN   Elg =   65°925

  R03 = l = 179°597                                  R06 = R.A. =  12h02m08s                          R08 = L =  185°186
  R04 = b =  2°268                                    R07 = Decl =   2°14'28"                             R09 = B =    2°356
  R05 = r =  9.904 AU                                                                                                 R10 = R =  9.535 AU
 

 Uranus
 

-If you've just XEQ "SA" , simply press  R/S  instead of   XEQ "UR"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "UR"   >>>>   Az = 319°25'07"  = R11
                            RDN     h  =   42°40'41"   = R12
                            RDN     h0 =   42°41'43"  = R13
                            RDN   Elg =  113°153

  R03 = l =   0°537                                  R06 = R.A. =   0h03m11s                           R08 = L =  -2°129
  R04 = b = -0°765                                  R07 = Decl =  -0°29'18"                             R09 = B =  -0°749
  R05 = r =  19.672 AU                                                                                              R10 = R =  20.093 AU
 

 Neptune
 

-If you've just XEQ "UR" , simply press  R/S  instead of   XEQ "NE"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "NE"   >>>>   Az =  2°07'56"   = R11
                            RDN     h  =  38°30'15"  = R12
                            RDN     h0 =  38°31'27" = R13
                            RDN   Elg =  145°508

  R03 = l =  -31°819                                  R06 = R.A. =  22h02m04s                           R08 = L =  -32°918
  R04 = b = -0°473                                    R07 = Decl = -12°32'59"                             R09 = B =  -0°460
  R05 = r =  29.176 AU                                                                                                 R10 = R =  30.019 AU
 

 Pluto
 

-If you've just XEQ "NE" , simply press  R/S  instead of   XEQ "PL"
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

         XEQ "PL"   >>>>   Az =  55°29'06"   = R11
                            RDN     h  =  10°53'39"  = R12
                            RDN     h0 =  10°58'25" = R13
                            RDN   Elg =  159°289

  R03 = l =  -86°413                                   R06 = R.A. =  18h15m03s                           R08 = L =  274°218
  R04 = b =   5°081                                    R07 = Decl = -18°18'35"                             R09 = B =   4°929
  R05 = r =  30.912 AU                                                                                                  R10 = R =  31.865 AU
 

 A Subroutine: "N"
 

-This subroutine takes  e  in Y-register and  M  in X-register  ( M = mean anomaly , e = eccentricity < 1 here )
-And it returns  v  in Y-register  and   r / q  in X-register  where

        v = true anomaly
        r = radius vector
       q = perihelion distance

-This is used, for instance, to study the motion of a comet, an asteroid.... on an elliptic orbit.
 
 
      STACK        INPUTS      OUTPUTS
           Y          e < 1            v
           X            M          r/q

Example:

  e  = 0.9     ENTER^
 M =  7°      XEQ "N"  >>>>   r / q = 3.152974273
                                   X<>Y     v   = 116°2027755

Note:

-The HP-41 must be in DEG mode.
 

 Satellites of Jupiter, Saturn and Uranus
 

-This program allows you to get the configuration of
  4 satellites of Jupiter ( Io - Europa - Ganymede - Callisto )
  7 satellites of Saturn ( Mimas - Enceladus - Tethys - Dione - Rhea -Titan - Hyperion )
  5 satellites of Uranus ( Ariel - Umbriel -Titania - Oberon - Miranda )
-The x-axis coincides with the equator of the planet.
-The center of the Planet is the origin and x , y are measured in equatorial radii.
 

                                     y ( North )
                                      |
                                      |
                                      |
( East ) ------------Jup/Sat/Ura------------------ x   ( West )
                                      |
                                      |
                                ( South )
 
 

Data Registers:  R00 thru R20 are used for temporay data storage and when the program stops:

                   Mimas  -  Enceladus  -   Tethys  -    Dione    -     Rhea    -    Titan   -  Hyperion

                             R01 = x1        R03 = x2         R05 = x3        R07 = x4         R09 = x5       R11 = x6         R13 = x7
                             R02 = y1        R04 = y2         R06 = y3        R08 = y4         R10 = y5       R12 = y6         R14 = y7

      or              Io      -    Europa   - Ganymede - Callisto

      or            Ariel   -    Umbriel  -   Titania    -   Oberon  -  Miranda

  >>>>     R19 =  - sin DE  where DE is the planetocentric declination of the Earth.

Flags:     F01  F02  F03  F04 F05  F06  F07
              -Flag nn  is set when the distance Earth-Satellite n  is shorter than the distance Earth-Planet   ( 0 < nn < 8 )
 

Subroutines:   "MST"  "J2000"  "SUN"  "JU"  "SA"  "UR"
 
 
         STACK         INPUTS       OUTPUTS
             Z    YYYY.MNDD              /
             Y    HH.MNSS(TT)              y1
             X       1 or 2 or 3              x1

        X = 1  for the Satellites of Jupiter
        X = 2  ------------------  Saturn
        X = 3  ------------------  Uranus

Example:      On  2009/07/16 at 12h41m  TT
 

     Jovian Satellites

     2009.0716  ENTER^
         12.41      ENTER^
             1         XEQ "SATJSU"     >>>>    x1 = -1.075     X<>Y     y1 = -0.053        and in registers  R01 thru R08:

                Io         -    Europa    -    Ganymede   -    Callisto

                  x1 = -1.075       x2 =  9.261         x3 = 13.120         x4 = -12.053         DE = 0°52
                  y1 = -0.053       y2 =  0.019         y3 = -0.059          y4 =   0.223

-Flags F01 and F03 are set:  Io & Ganymede are closer to the Earth than Jupiter.
 

     Saturnian Satellites

     2009.0716  ENTER^
         12.41      ENTER^
             2            R/S             >>>>    x1 = -2.426     X<>Y     y1 = -0.058        and in registers  R01 thru R14:

                   Mimas     -   Enceladus  -     Tethys    -     Dione     -    Rhea      -     Titan     -    Hyperion

                  x1 = -2.426       x2 = -0.945         x3 =  2.034      x4 = -5.336     x5 = -8.749       x6 = -11.495     x7 = -10.289 DE = -2°63
                  y1 = -0.058       y2 = -0.176         y3 = -0.112      y4 =  0.150     y5 = -0.039       y6 =  0.626        y7 =   0.764

-Flags   F04  F05  F06  F07  are set, whence  Dione, Rhea, Titan and Hyperion are closer to the Earth than Saturn.
 

     Uranian Satellites

     2009.0716  ENTER^
         12.41      ENTER^
             3            R/S            >>>>    x1 = 6.052     X<>Y     y1 = -0.667        and in registers  R01 thru R10:

            Ariel    -    Umbriel    -    Titania    -     Oberon   -   Miranda

                  x1 =  6.052       x2 = -7.880       x3 = 5.296        x4 = 22.755     x5 = 4.557            DE = 8°76
                  y1 = -0.667      y2 =  1.039       y3 = -2.483        y4 =  0.265     y5 = 0.449

-Flags  F01 & F03 are set, so  Ariel & Titania are closer to the Earth than Uranus.
 

Notes:

-The accuracy is of order of a few hundredths of the planet's radius.
-Hyperion's coordinates are less accurate than those of the other Saturnian Satellites.
-For Hyperion, the series converge slowly and several terms should be added to get more accurate results.
-The complete series may be downloaded from the server   ftp://ftp.imcce.fr
-For the Uranian Satellites, the error in the y-values may reach 0.1
 

 Parabolic, Hyperbolic, Elliptic Motions
 

 Parabolic Motion
 

-Let  s = Tan v/2  where  v = true anomaly,
 s  is the unique solution of the cubic equation:   s3 + 3s - W = 0
 where  W = (3k/sqrt(2)) (t-T) q -3/2   and  k = 0.01720209895 = Gaussian gravitational constant

-Then, the radius vector is obtained by  r = q(1+s2)
 

 Hyperbolic Motion
 

-For hyperbolic orbits, Kepler's equation is:    k(t-T) ((e-1)/q)3/2  = e Sinh E - E
-Newton's method is used in the M-Code routine KPL.
-Then, the true anomaly v is obtained by   Tan v/2 = ((e+1)/(e-1))1/2  Tanh E/2
  and the radius vector by   r = q(1+e)/(1+e cos v)
 

 Elliptic Motion
 

-We solve Kepler's equation   M = E - e sin E  ( in radians )     where  M = k(t-T) ((1-e)/q)3/2  = mean anomaly and  E = excentric anomaly
-This equation is solved by Newton's method and the M-code routine KPL.
-Then, the true anomaly v is obtained by   Tan v/2 = ((e+1)/(1-e))1/2  Tan E/2
  and the radius vector is   r = q(1+e)/(1+e cos v)
 

Data Registers:       R00 = Number of millenia since 2000/01/01 0h

  R01 & R02 = rectangular ecliptic coordinates of the Sun
  R03 = geocentric longitude in decimal degrees = l
  R04 = geocentric latitude    ------------------  = b
  R05 = distance to the Earth in Astronomical Units = r
  R06 = Right Ascension in  hh.mnss = R.A
  R07 = Declination in          ° ' ''       = Decl
  R08 = heliocentric longitude in decimal degrees = L
  R09 = heliocentric latitude   ------------------- = B
  R10 = radius vector in Astronomical Units        = R
  R11 = Azimuth in ° ' '' = Az     reckoned clockwise from the South
  R12 = True Altitude in ° ' '' = h
  R13 = Apparent Altitude in ° ' '' = h0 ( after correction for refraction )
  R14 = Elongation from the Sun = Elg
  R16 = Longitude ( positive East ) in ° ' "
  R17 = Latitude in ° ' "
  R18 = Delta T in seconds
  R19 = Mean sidereal time at Greenwich                        R15 & R20 are unused

   R21 = T = time of passage in perihelion - expressed in days since 2000/01/01 0h TT  ( not 12h )
   R22 = q = perihelion distance in AU
   R23 = e = eccentricity
   R24 = i = inclination in degree                                                       referred to
   R25 = omega = argument of the perihelion in degree                      the ecliptic
   R26 = OMEGA = longitude of the ascending node in degree         of the date

Flag: /
Subroutines:  "MST"  "SUN"  "N"  ......
 
 
 STACK  INPUTS                 OUTPUTS
      T       /   Elg = elongation from the Sun ( deg )
      Z       /      h0 = apparent altitude  ( ° ' '' )
      Y       /         h = true altitude  ( ° ' '' )
      X       /          Az = Azimuth  ( ° ' '' )

Example1:        Calculate the position of comet Kohler for 1977 September 29  0h UT using the following elements

    T = 1977 November 10.5659  TT         i2000      =   48°7131
    q =  0.990662                                 omega2000   = 163°4788
    e =  1                                            OMEGA2000 = 182°1660

                                                                                         i       =    48°7160
-First, we reduce the elements to 1977/11/29  it gives   omega   =  163°4798           use for instance  "IOO"
                                                                                  OMEGA = 181°8582

-We have  T = -8086.4341 days  from  2000/01/01 0h

   XEQ "PHEM"   >>>>      the HP-41 displays  "T^Q^E?"

     -8086.4341   ENTER^
       0.990662    ENTER^
             1              R/S        the HP-41 displays  "I^O^OO?"

       48.7160    ENTER^
      163.4798   ENTER^
      181.8582      R/S        the HP-41 displays  "LONG^LAT?"  and,  if we were at the USNO:

       -77.0356    ENTER^
         38.5517      R/S           next  PROMPT   "TT-UT=SEC?"

            49            R/S           next  PROMPT   "DATE^TIME?"

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

     1977.0929   ENTER^
             0             R/S            Az =  78°07'22"   = R11

                          RDN             h  =  46°36'14"  = R12
                          RDN             h0 = 46°37'08"  = R13
                          RDN            Elg =  62°505     = R14

  R03 = l =  -121°944                                  R06 = R.A. =  16h19m10s                           R08 = L =  302°797
  R04 = b =   40°960                                    R07 = Decl =  20°16'22"                             R09 = B =   44°329
  R05 = r =     1.306 AU                                                                                                   R10 = R =   1.225 AU

>>>  To calculate its position for another date, place the date & time ( UT ) in Y & X-registers and  R/S

  For instance   1977.1212  ENTER^
                              20            R/S       yields   Az = 334°24'25"  ... etc ...
 

Example2:         Calculate the position of comet C/2007 T1 for 2008 January 01  6h  UT using the following elements
 

    T = 2007 December 12.49731  TT                i2000      = 117°649041
    q =  0.969480                                          omega2000   = 233°671201
    e =  1.000785                                        OMEGA2000 = 111°418623

                                                                            i       =  117°64857
-Reduce the elements to 2008/01/01.25          omega   =  233°67226             use for instance  "IOO"
                                                                     OMEGA = 111°53086

-We have  T = 2902.49731 days  from  2000/01/01 0h

   XEQ "PHEM"   >>>>      the HP-41 displays  "T^Q^E?"

     2902.49731  ENTER^
       0.969480    ENTER^
       1.000785       R/S        the HP-41 displays  "I^O^OO?"

      117.64857   ENTER^
      233.67226   ENTER^
      111.53086      R/S        the HP-41 displays  "LONG^LAT?"  and,  if we were at the USNO:

       -77.0356    ENTER^
         38.5517      R/S           next  PROMPT   "TT-UT=SEC?"

            65            R/S           next  PROMPT   "DATE^TIME?"

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

     2008.0101   ENTER^
             6             R/S            Az = 320°13'55"   = R11

                          RDN             h  = -59°25'36"  = R12
                          RDN             h0 = -59°25'36" = R13
                          RDN            Elg =  39°158     = R14

  R03 = l =  -99°228                                  R06 = R.A. =  17h02m53s                           R08 = L =  219°595
  R04 = b = -34°699                                  R07 = Decl = -57°40'52"                             R09 = B =  -61°144
  R05 = r =     1.582 AU                                                                                                R10 = R =   1.029 AU

>>>  To calculate its position for another date, place the date & time ( UT ) in Y & X-registers and  R/S

  For instance   2008.0131  ENTER^
                              23            R/S       yields   Az = 357°55'18"  ... etc ...
 

Example3:     Calculate the position of comet C/2007 K6 for 2007 December 01  0h  UT  using the following elements
 

    T = 2007 July 01.47533  TT                 i2000      =  105°063204
    q =  3.432968                                 omega2000   = 337°140230
    e =  0.984585                               OMEGA2000 = 298°075386

                                                                                 i       =  105°06377
-We reduce the elements to 2007/12/01  it gives   omega   =  337°13934    use for instance  "IOO"
                                                                          OMEGA = 298°18570

   T = 2007 July 01.47533 = 2738.47533 days from  2000/01/01 0h

   XEQ "PHEM"   >>>>      the HP-41 displays  "T^Q^E?"

     2738.47533  ENTER^
       3.432968    ENTER^
       0.984585       R/S        the HP-41 displays  "I^O^OO?"

      105.06377   ENTER^
      337.13934   ENTER^
      298.18570      R/S        the HP-41 displays  "LONG^LAT?"  and,  if we were at the USNO:

       -77.0356    ENTER^
         38.5517      R/S           next  PROMPT   "TT-UT=SEC?"

            65            R/S           next  PROMPT   "DATE^TIME?"

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

     2007.1201   ENTER^
             0             R/S            Az =  62°27'44"   = R11

                          RDN             h  =  8°14'54"    = R12
                          RDN             h0 =  8°21'05"   = R13
                          RDN           Elg =  38°523     = R14

  R03 = l =  -73°633                                  R06 = R.A. =  19h07m27s                           R08 = L =  295°894
  R04 = b =   7°068                                    R07 = Decl = -15°25'02"                             R09 = B =   8°451
  R05 = r =    4.426 AU                                                                                                  R10 = R =   3.706 AU

>>>  To calculate its position for another date, place the date & time ( UT ) in Y & X-registers and  R/S

  For instance   2007.1207  ENTER^
                           23.4523      R/S       yields   Az = 64°46'26"  ... etc ...

Notes:

-When the HP-41 displays "T^Q^E?" , simply press R/S without any digit entry if the 6 elements are already stored in R21 thru R26
-This is the case if you've just used "GAUSS" to find the orbital elements of a new comet.
-If you have to calculate several positions, i , omega , OMEGA may be regarded as constant for several days or weeks
 ( though "IOO" is relatively easy to use... )
-The program uses the ecliptic and equator of the date and that could be criticized because this is not an inertial frame.
-However, the resulting errors are usually negligible - provided the intervals of time remain of the order of a few days or weeks or even months.
 

 Preliminary Orbit Determination by Gauss' Method
 

-Gauss' method computes approximate values of  T , q , e , i , omega , OMEGA  from 3 observations of a comet.
-This program also uses the Herrick-Gibbs formulae which are 4th-order accurate in time.
-The method fails if the inclination i = 0 or is very small.

-At an instant  t1  the right ascension of the comet =  RA1  its declination = d1 and the rectangular ecliptic coordinates of the Sun are  X1 & Y1
-At an instant  t2  the right ascension of the comet =  RA2  its declination = d2 and the rectangular ecliptic coordinates of the Sun are  X2 & Y2
-At an instant  t3  the right ascension of the comet =  RA3  its declination = d3 and the rectangular ecliptic coordinates of the Sun are  X3 & Y3

-We must have   t1 <  t2 <  t3

 >>> The results will be more accurate if  t3 - t2 is equal ( or nearly equal ) to   t2 - t1  ( this value is usually of the order of a few days )

---------------------------------------------------------FORMULAE-------------------------------------------------------------------------

-Let

           li  =  cos RAi cos di
           mi = sin obl sin di + cos obl cos di sin RAi    where   obl is the obliquity of the ecliptic,   i = 1 , 2 , 3
           ni = cos obl sin di - sin obl cos di sin RAi

-Let Di = distance Earth-Comet at the instant ti , ui ( li , mi , ni ) , -Ri = position vector of the Earth , ri = position vector of the Comet

    ri = Di ui  - Ri        (  || ui || = 1  )

-If we neglect the planetary perturbations, the 3 position-vectors ri  lie in the same plane, so there are coefficients c1 , c3 such that

   r2 = c1 r1 + c3 r3   and likewise for the velocity at the second instant:        V2 = -d1 r1 + d2 r2 + d3 r3

-One can prove that, approximately:     ci = Ai + ( 1 + B2 / r23 )  Bi / r23   ( i = 1 , 3 )    with   B2 = [ k2(t3 - t1)2/12 ] ( 1 + A1 A3 )

      and   di = Gi + Hi / ri3  ( i = 1 , 2 , 3 )

     where    A1 = (t3 - t2)/(t3 - t1)  ,   B1 = k2A1 [ (t3 - t1)2 - (t3 - t2)2 ]/6      with   k =  0.01720209895 = Gaussian gravitational constant
                  A3 = (t2 - t1)/(t3 - t1)  ,   B3 = k2A3 [ (t3 - t1)2 - (t2 - t1)2 ]/6

      and      G1 = (t3 - t2)/[(t3 - t1)(t2 - t1)]       G3 = (t2 - t1)/[(t3 - t1)(t3 - t2)]        G2 = G1 - G3
                 H1 = k2(t3 - t2)/12                           H3 = k2(t2 - t1)/12                            H2 = H1 - H3

-Let   A( A1 X1 - X2 + A3 X3 , A1 Y1 - Y2 + A3 Y3 ,  A1 Z1 - Z2 + A3 Z3 )
 and   B( B1 X1 + B3 X3 , B1 Y1 + B3 Y3 ,  B1 Z1 + B3 Z3 )                               ( all the Zi's ~ 0 if we use the rectangular ecliptic coordinates of the Sun )

         P1 = A.( u2 x u3 )/D          Q1 = B.( u2 x u3 )/D
         P2 = A.( u1 x u3 )/D          Q2 = B.( u1 x u3 )/D          where      D = u1.( u2 x u3 )           ( . = dot product , x = cross product )
         P3 = A.( u1 x u2 )/D          Q3 = B.( u1 x u2 )/D

-The distance Earth-Comet at the central instant is found by solving the equation

        D2 = P2 + ( 1 + B2 / r23 ) Q2 / r23   with    r2 =  D22 - 2 D2 ( l2 X2 + m2 Y2 ) + X22 + Y22

-Then,  D1 = [ P1 + ( 1 + B2 / r23 ) Q1 / r23 ] / [ A1 + ( 1 + B2 / r23 ) B1 / r23 ]
  and    D3 = [ P3 + ( 1 + B2 / r23 ) Q3 / r23 ] / [ A3 + ( 1 + B2 / r23 ) B3 / r23 ]

  and the rectangular ecliptic coordinates of the comet at the 3 instants are:

           xi = li Di - Xi
           yi = mi Di - Yi        i = 1 , 2 , 3         so we know the 3 vectors ri's   and the velocity V2 = -d1 r1 + d2r2 + d3 r3   may be computed
           zi = ni Di

-We calculate the cross-product  h = r2 x V2 and then:

-The parameter  p = h2/k2
-The inclination  i = Acos hz/h = Atan2 (hx2+hy2)1/2/hz                          you can also use C = r1 x r3  instead of  h  to get i and OMEGA

-The longitude of the ascending node  OMEGA = Atan2 hx/(-hy)

   ( vi + omega ) is obtained by

    ri Sin ( vi + omega ) = zi / Sin i
    ri Cos ( vi + omega ) = xi Cos OMEGA + yi Sin OMEGA        whence   (v3 - v1)/2

  ( v + omega ) may also be computed by the formula:

    v + omega = sign(z) Acos [ ( x Cos OMEGA + y Sin OMEGA ) / r ]   with perhaps less accurate results in some cases.

-The relations:   2e Sin (v3 + v1)/2 = ( p/r1 - p/r3 ) / Sin (v3 - v1)/2                   give  e and (v3 + v1)/2  whence v1 and omega
                        2e Cos (v3 + v1)/2 = ( p/r1 + p/r3 - 2 ) / Cos (v3 - v1)/2

-And the time of passage in perihelion is computed via Kepler's equation.
 

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

-Registers R01 thru R15 are to be initialized.
-However, if you set F00, the HP-41 displays several PROMPTs and calculate the rectangular coordinates of the Sun  X & Y
 

Data Registers:              R00 & R16 to R42: temp         ( Registers R01 thru R15 are to be initialized )

                                        R01 = t1                 R06 = t2                 R11 = t3                  expressed in days since 2000/01/01  0h TT
                                        R02 = RA1             R07 = RA2             R12 = RA3             in hh.mnss
                                        R03 = d1                R08 = d2                R13 = d3                 in °. ' ''
                                        R04 = X1               R09 = X2               R14 = X3                 in Astronomical Units
                                        R05 = Y1               R10 = Y2               R15 = Y3                 ---------------------

    >>>> When the program stops,         R21 = T   ( in days from 2000/01/01  0h )
                                                             R22 = q    ( in Astronomical Units )
                                                             R23 = e
                                                             R24 =  i     ( in degrees )
                                                             R25 = omega  ( in degrees )
                                                             R26 = OMEGA  ( in degrees )
                                                             R27 = n = mean motion  ( in degrees/day ) - virtual for a hyperbolic orbit.
Flag:  F00
Subroutines:   "SUN"  "MST"
 
 
      STACK         INPUT      OUTPUT
           T            /             i
           Z            /             e
           Y            /             q
           X            D             T

   where  D  is an estimation of the distance Comet-Earth ( in AU ) at the 2nd instant ( take for instance D = 1 )

Example:    You observe a comet on 2007/07/01 0h , 2007/07/05 0h , 2007/07/09 0h  TT  the positions were:

   Dates        2007/07/01 0h                   2007/07/05 0h                   2007/07/09 0h
     RA           14h27m26                         14h16m34s                        14h06m38s
    Decl          -39°30'56"                         -38°44'07"                         -37°52'59"

   SF 00   1   XEQ "GAUSS"   >>>    "DATE^TIME1?"

     2007.0701   ENTER^
             0             R/S              >>>     "RA^DECL1?"

       14.2726     ENTER^
      -39.3056        R/S             >>>      "DATE^TIME2?"

     2007.0705   ENTER^
             0             R/S              >>>     "RA^DECL2?"

       14.1634     ENTER^
      -38.4407        R/S             >>>      "DATE^TIME3?"

     2007.0709   ENTER^
             0             R/S              >>>     "RA^DECL3?"

       14.0638     ENTER^
      -37.5259        R/S             >>>      T = 2817.84871 = R21  ( in days from 2000/01/01 0h TT )     Activate the turbo mode...
                                               RDN     q = 0.70439 = R22
                                               RDN     e = 0.76153 = R23
                                               RDN     i  = 9°70768 = R24

  and   R25 = omega = -1°22487
          R26 = OMEGA = 3°49104         Remark that  T = 2007/09/18.84871    ( use "DT" to get this result )
          R27 = n = 0.19415 deg/day

Notes:

-"GAUSS" solves an equation by the secant method.
-This equation may have several solutions ( 2 or sometimes 3 )
-If the solution that is returned is not satisfactory, XEQ "GAUSS" again with CF 00 and another guess in X-register.

-Sometimes, the algorithm converges to a negative - meaningless - distance D = Comet-Earth.
-In this case, it stops and displays "D=?"
-Place a different initial guess ( say 3 ) and press  R/S

-Gauss' method is very sensitive to observational errors, so take the average of successive results to increase the accuracy.

-"GAUSS" works for elliptic, hyperbolic and parabolic orbits, though the case e = 1 exactly will be very rare...
-The orbital elements are automatically stored in the proper registers to be used by "PHEM"

-If you execute "PHEM" with the orbital elements found above ( take TT-UT = 0 ),
  you'll find - almost - the same right ascensions & declinations
-The largest errors occur at the 2nd instant 2007/07/05  0h:

  "PHEM"  returns for this date:   RA = 14h16m34s029  &  Decl = -38°44'07"173   in registers  R06 & R07
                      instead of              RA = 14h16m34s       &   Decl = -38°44'07"
 

 Equatorial > Ecliptic
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h        R11 & R12: temp
Flags:  /
Subroutines:   "J2000"   "OBL"   "EE"
 
 
        STACK         INPUTS       OUTPUTS
             Z   YYYY.MNDDdd             em
             Y       decl ( ° . ' " )         b ( ° . ' " )
             X      RA (hh.mnss)         l ( ° . ' " )

  where   l  = ecliptic longitude , b  = ecliptic latitude , em = mean obliquity of the ecliptic
 

Example:     On  2134/04/04 at 0h  if  Right-Ascension = 12h34m56s  and  Declination =  25°12'49"  ( mean coordinates )

    2134.0404  ENTER^
        25.1249  ENTER^
        12.3456  XEQ "EQ-ECL"  >>>>  l  =  177°13'44"69
                                                  RDN   b  =  26°27'18"71
 

 Ecliptic > Equatorial
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h        R11 & R12: temp
Flags:  /
Subroutines:   "J2000"   "OBL"   "EE"
 
 
        STACK         INPUTS       OUTPUTS
             Z   YYYY.MNDDdd             -em
             Y        b ( ° . ' " )       Decl ( ° . ' " )
             X         l ( ° . ' " )      RA (hh.mnss)

  where   l  = ecliptic longitude , b  = ecliptic latitude , em = mean obliquity of the ecliptic

Example:     On  2134/04/04 at 0h  if  l  =  177°13'44"69  and b  =  26°27'18"71  ( mean coordinates )

    2134.0404      ENTER^
        26.271871  ENTER^
      177.134469  XEQ "ECL-EQ"  >>>>    RA = 12h34m56s00
                                                      RDN    decl =  25°12'49"00
 

 Equatorial > Azimuthal
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h      ( Registers R16 & R17 are to be initialized before executing "EQ-AZ" )
                               R11 thru R15: temp

                            R16 = Longitude ( positive East ) in ° ' "
                            R17 = Latitude in ° ' "
Flags: /
Subroutines:    "J2000"   "OBL"   "MST"   "EE"
 
 
        STACK         INPUTS       OUTPUTS
             T    YYYY.MNDD             /
             Z    HH.MNSS (UT)          h0 ( ° . ' " )
             Y       decl ( ° . ' " )          h  ( ° . ' " )
             X      RA (hh.mnss)         Az ( ° . ' " )

 Azimuths are measured clockwise ( westwards ) from South
 h  = Altitude
 h0 = Apparent Altitude ( corrected for refraction )

Example:    On  2005/12/12  at  20h51m29s (UT)  we have  R.A. = 7h41m16s  &  decl = 60°21'37"   ( mean coordinates ),
                    Compute the azimuthal coordinates at the US Naval Observatory.

       Longitude =  77°03'56" W = -77°03'56"  STO 16
         Latitude  =  38°55'17" N = +38°55'17"  STO 17

    2005.1212  ENTER^
        20.5129  ENTER^
        60.2137  ENTER^
          7.4116  XEQ "EQ-AZ"  >>>>   Az =  190°56'31"
                                                 RDN    h  =    10°55'58"
                                                 RDN   h0 =     11°00'47"
 

 Azimuthal > Equatorial
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h      ( Registers R16 & R17 are to be initialized before executing "AZ-EQ" )
                               R11 thru R15: temp

                            R16 = Longitude ( positive East ) in ° ' "
                            R17 = Latitude in ° ' "
Flags: /
Subroutines:    "J2000"   "OBL"   "MST"   "EE"
 
 
        STACK         INPUTS       OUTPUTS
             T    YYYY.MNDD             /
             Z    HH.MNSS (UT)             /
             Y     altitude ( ° . ' " )       decl ( ° . ' " )
             X     azimuth ( ° . ' " )      RA (hh.mnss)

-Azimuths are measured clockwise ( westwards ) from South

Example:    On  2005/12/12  at  20h51m29s (UT)  we have  Azimuth = 190°56'31"   &   Altitude = 10°55'58"   ( mean coordinates ),
                    Compute the equatorial coordinates at the US Naval Observatory.

       Longitude =  77°03'56" W = -77°03'56"  STO 16
         Latitude  =  38°55'17" N = +38°55'17"  STO 17

    2005.1212    ENTER^
        20.5129    ENTER^
        10.5558    ENTER^
      190.5631    XEQ "AZ-EQ"  >>>>   R.A. = 7h41m16s
                                                   RDN    decl =  60°21'37"
 

 Precession ( ecliptic coordinates )
 

Data Registers:    R00 & R11 thru R14: temp
Flag:   F00
Subroutines:  "J2000"   "EE"
 
 
        STACK         INPUTS       OUTPUTS
             T   YYYY.MNDDdd1             /
             Z   YYYY.MNDDdd2             /
             Y        b1 ( ° . ' " )         b2 ( ° . ' " )
             X        l1  ( ° . ' " )         l2  ( ° . ' " )

  where   l  = ecliptic longitudes , b  = ecliptic latitudes

Example:   The mean ecliptic coordinates of Sirius on 1600/04/04 0h are  l1 = 98°30'58"32  ;  b1 =  -39°39'17"79
                   Calculate the ecliptic coordinates on 2134/12/12 0h

   1600.0404    ENTER^
   2134.1212    ENTER^
  -39.391779   ENTER^
    98.305832   XEQ "PREC"  >>>> l2 =  105°57'44"15
                                               RDN   b2 =  -39°35'19"15
 

-The polynomial coefficients of the precession angles ( expressed in degrees ) are rounded to 6 decimals, with T expressed in millenia from J2000.
-So, the accuracy is of the order of 0.01 arcsecond over the time span [ 1000 ; 3000 ]
 

 Precession ( equatorial coordinates )
 

Data Registers:    R00 & R11 thru R14: temp
Flag:   F00
Subroutines:  "J2000"   "EE"   "PREC"
 
 
        STACK         INPUTS       OUTPUTS
             T   YYYY.MNDDdd1             /
             Z   YYYY.MNDDdd2             /
             Y     decl1 ( ° . ' " )      decl2 ( ° . ' " )
             X     RA1 (hh.mnss)      RA2 (hh.mnss)

Example:   The mean coordinates of Sirius on 1600/04/04 0h are  AR = 6h27m17s88 ;  Decl = -16°21'56"34
                   Calculate the equatorial coordinates on 2134/12/12 0h

   1600.0404    ENTER^
   2134.1212    ENTER^
  -16.215634   ENTER^
     6.271788   XEQ "PREQ"  >>>>   RA =  6h51m10s79
                                              RDN   Decl = -16°52'22"05
 

 Orbital elements from one equinox to another one
 

Data Registers:             R00 = Date2 = YYYY.MNDDdd2              ( Register R00 is to be initialized before executing "IOO" )

                                         R01 thru R07: temp

Flag:  F00 is cleared at the end
Subroutines:    "J2000"    "EE"
 
 
        STACK        INPUTS      OUTPUTS
             T  YYYY.MNDDdd1              i2
             Z        OMEGA1       OMEGA2
             Y         omega1         omega2
             X             i1              i2

Example:    We have   i1 = 12°789   omega1 = 49°345   OMEGA1 = 166°234  on  1600/01/01
                     Compute  i2 , omega2 , OMEGA2  on 2900/12/12

   2900.1212  STO 00

   1600.0101  ENTER^
      166.234   ENTER^
       49.345    ENTER^
       12.789    XEQ "IOOM"  >>>>         i2        =   12.619940°
                                              RDN      omega2   =   49.370109°
                                              RDN    OMEGA2 = 184.401887°
 
 

 True Altitude > Apparent Altitude ( mean conditions of temp & press )
 

-The light coming from a star is curved by the atmosphere,  so that its apparent altitude differs from its true altitude.
-Knowing the refraction R allows to convert the apparent altitude h0 and the true altitude h of a given star:   h = h0 - R
-"H-H0" & "H0-H" use data from the Pulkovo Refraction Tables.

   Temperature:                                +15°Celsius
   Pressure:                                    1013.25 mbar
   Light wave-length:                         0.590 µm
   Partial pressure of water vapor:          0   ( dry air )
   Latitude:                                           45°
   Observer's altitude:                            0   ( i-e at sea-level )
 

 Formula:      R ~ (1°/62.6)/Tan( h + 5.459/( h + 19.272/( h + 6.942 ) ) )       where h is expressed in degrees

-Errors are smaller than  0"8  over the whole range [ -0°32'58" , 90° ]
 

Data Registers: /
Flags: /
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           Z             /             h
           Y             /             R
           X             h             h0

-All angles expressed in ° ' " ( but Z-output in degrees and decimals )

Example:        With  h = 1°30'00"

  1.30  XEQ "H-H0"  >>>>  h0 = 1°48'38"6
                                 RDN    R = 0°18'38"6
 

 Apparent Altitude > True Altitude  ( mean conditions of temp & press )
 

 Formula:       R ~ (1°/62.83)/Tan( h0 + 4.208/( h0 + 14.978/( h0 + 5.906 ) ) )       where h0 is expressed in degrees

-Errors are smaller than  0"34  over the whole range [ 0° , 90° ]
 

Data Registers: /
Flags: /
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           Z             /             h0
           Y             /             R
           X             h0             h

-All angles expressed in ° ' " ( but Z-output in degrees and decimals )

Examples:     With  h0 = 1°30'00"  ;   h0 = 27°  ;   h0 = 0

  1.30  XEQ "H0-H"   >>>>   h =   1°09'42"6    RDN    R = 0°20'17"4
    27   XEQ "H0-H"   >>>>   h = 26°58'08"3    RDN    R = 0°01'51"7
     0    XEQ "H0-H"   >>>>   h = -0°32'58"0     RDN   R =  0°32'58"0
 

 A more general refraction program
 

 "REFR" is a generalization of  "H0-H"

 Formulae:   The standard refraction is computed by  R0 ~ (57"298)/Tan( h0 + 4.208/( h0 + 14.978/( h0 + 5.906 ) ) )

   R = R0 ( P/1013.25 ) [ 286.68 / ( t + 271.68 ) ] [ 1 - f / 6579 - ( f/426 )2 ] [ 1 + { (15-t) + (15-t)2 / 196 } exp(-h0/7) / (1+h0) / 157 ]
                                                     x [ 1 + (P-1013.25) exp(-0.4h0)/12181 ] [ 1 - ( f + f2/16 ) exp(-0.18h0) / (1+h0) / 5198 ]
 where

       t  = Temperature
       P = Pressure
       f  = Humidity

Data Registers:    R00 = h0                                                          ( Registers R01 thru R03 are to be initialized before executing "REFR" )
                             R01 = t     ( in °C )           ( -10 <= t <= +30 )
                             R02 = P    ( in mbar )      ( 700 <= P <= 1100 )
                             R03 = f     ( in mbar )          ( 0 <= f <= 20 )
Flags: /
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           Y             /             R
           X             h0             h

   ( All angles are expressed in ° ' " )

Example:         t = 0°C ;  P = 900 mbar  ;  f = 12 mbar

      0  STO 01    900  STO 02    12  STO 03

-If   h0 =  0°                     0      XEQ "REFR"  >>>>   h = -0°33'32"       X<>Y    R = 0°33'32"
-If   h0 =  10°23'45"    10.2345       R/S          >>>>   h = 10°19'02"8    X<>Y    R = 0°04'42"2
-If   h0 =  49°12'34"    49.1234       R/S          >>>>   h = 49°11'47"9    X<>Y    R = 0°00'46"1
 

 Barycenter of the Solar System
 

-This program calculates the rectangular heliocentric ecliptic coordinates X , Y of the center of mass of the Solar system ( we make the approximation Z = 0 )
-It takes a date YYYY.MNDD and returns the distance G between the Solar system barycenter B and the center of the Sun S.

-Only the 4 giant planets are taken into account: the other planets have a negligible effect.
-The positions of these planets are not computed very accurately: the terms in e2 are neglected ( e = eccentricity ),
  and only one term in the perturbations is taken into account.
 

Data Registers:   R00 = T in millennia since 2000/01/01 0h TT
                              R01 = X
                              R02 = Y        R03 to R11: temp
Flags: /
Subroutine:    "J2000"
 
 
      STACK        INPUTS      OUTPUTS
           Y             /             µ
           X  YYYY.MNDD             G

   where   G = distance Barycenter-Center of the Sun  ( unit = solar radius )
     and     µ  = ( Sx , SB )

Example:

   2008.0824   XEQ "BOSS"    >>>>    G =     1.0761
                                                X<>Y    µ =  -74°0782

  Since G > 1 , the Solar system barycenter is currently outside the Sun's globe.
  G will be smaller than one in April 2010

  We have also   X = R01 =  0.2952
                         Y = R02 = -1.0348

Other Checkings:    on  1983/03/15    G = 2.098   ( exact result - almost the maximum possible value )
                                  on  1990/04/23    G = 0.069   ( exact value = 0.066 )
                                  on  2130/03/10    G = 0.016   ( exact value = 0.021 )

Note:

  Errors seem smaller than 0.01 over the time span [ 1000 , 3000 ]
  µ may be inaccurate if G is very small
 

 Binary Stars
 

-This program calculates the polar coordinates ( r , µ ) of a binary star  ( S1 ; S2 )
 and the eccentricity e' of the apparent orbit, as seen from the Earth.
 

                   North
                       |
            S2       |
              \        |
                \  µ  |
               r  \    |                                     r  =  the angular distance = S1S2
                    \  |                                     µ  =  the position angle = ( S1N , S2N )
                      \|
                      S1
 

Data Registers:       R07-R08: temp

                               R00 = P = period of revolution ( expressed in years )
                               R01 = T = time of perihelion passage ( given as a year and decimals )
                               R02 = a = semimajor axis ( in arcseconds )
                               R03 = e = eccentricity
                               R04 = i  = inclination of the orbit
                               R05 = omega = longitude of the periastron.
                               R06 = OMEGA = position angle of the ascending node

Flag:  F00        if  CF 00  you have to initialize  R00 thru R06
                         if  SF 00   the HP-41 displays several PROMPTs
Subroutine:   "N"
 
 
      STACK        INPUTS      OUTPUTS
           Z             /            e'
           Y            /            µ
           X       YYYY.yy             r

Example:   The elements of  Gamma Virginis are:

       P = 168.68 years
       T = 2005.13                     i       = 148.
       a = 3.697"                   omega   = 256.
       e = 0.885                  OMEGA =  36.

-Calculate  r  and  µ  for the epoch  2010.25    ( ~ 2010/04/01 )

    SF 00
   2010.25  XEQ "BSTAR"  >>>>   " P^T^A^E?"

    168.68    ENTER^
   2005.13   ENTER^
     3.697     ENTER^
     0.885        R/S               >>>>    " I^O^OO?"

     148.0    ENTER^
     256.5    ENTER^
      36.9         R/S              >>>>       r =  1.544"
                                          RDN       µ =  19.66°
                                          RDN       e' = 0.844

-Likewise   2015    R/S      >>>>       r =  2.352"
                                          RDN       µ =  4.92°
                                        ( RDN       e' = 0.844 )
 

 Redshift > Cosmological Distances
 

-Knowing the cosmological redshift z of a galaxy, and the values of the current parameters  (Omega)mat , (Omega)lambda , (Omega)rad
  "COSMO" computes several "distances", the recessional velocity, the scale factor and the age of the Universe.

        D   =  light-travel time distance ( in light-years ) = light-travel time ( in years )
        D0 = comoving distance ( in light-years )
        DL = luminosity-distance ( in light-years )
        VR = recessional velocity ( speed of light = 1 )
         t0  = age of the Universe ( in years )
         k   = curvature    ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )

Formulae:

     t0 = (c/H0)     §01     y. [  (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad  ] -1/2 dy
    D  = (c/H0)  §y(em)1  y. [  (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad  ] -1/2 dy           y(em) = y at the instant of emission
    D0 = (c/H0)  §y(em)1      [  (Omega)lambda.y4 + ( 1-(Omega)tot ).y2 + (Omega)mat.y + (Omega)rad  ] -1/2 dy               z + 1 = R0/R = 1/yem

  where  (Omega)tot = (Omega)mat + (Omega)lambda + (Omega)rad

                                                                F(x) = Sinh(x)  if  k = -1
     DL = R0 ( z + 1 ) F(D0/R0)    where    F(x) =   x         if  k = 0
                                                               F(x) = Sin(x)    if  k = +1

     VR = H0.D0

-These integrals are calculated by Gauss-Legendre 3-point formula ( the change of variable y = u2  facilitates the numerical integration )
-The number of subintervals is doubled until 2 consecutive rounded values are equal.
-Thus, the accuracy is controlled by the display settings  ( avoid FIX 9 and SCI 9 , it could lead to infinite loops ! )
 
 

Data Registers:   R00 thru R17                           ( Registers R11-R12-R13-R14 are to be initialized before executing "COSMO" )

         R11 = (Omega)mat       density parameter ( matter )
         R12 = (Omega)lambda  cosmological constant
         R13 = (Omega)rad       radiation parameter
         R14 =  1/H0 = inverse of Hubble's constant expressed in years  ( for instance 13.7 109 )

 -when the program stops:   R00 = k       ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
                                              R01 = D ,  R02 = D0 , R03 = DL , R04 = VR , R05 = t0 , R06 = R0 ( current scale factor )

Flags:  F00  set flag F00  if you don't want to (re-)calculate the age of the Universe t0
             F10

Subroutines:  /
 
 
      STACK        INPUTS      OUTPUTS
           T              /            VR
           Z              /            DL
           Y              /            D0
           X              z            D
           L              /            t0

Example1:      With    (Omega)mat = 0.044 ,  (Omega)lambda  = 0.519 ,  (Omega)rad  = 0.000049  &   1/H0 =  13.7 109

   0.044  STO 11         0.000049   STO 13
   0.519  STO 12         13.7  E9     STO 14     and, for instance,  FIX 4

     z = 7  XEQ "COSMO"  >>>>   D  =  1.4046  E10  l-y     = R01
                                            RDN    D0 = 3.4012  E10  l-y     = R02
                                            RDN    DL = 4.1177  E11  l-y     = R03
                                            RDN    VR = 2.4826                   = R04         ( unit c = 1 )
                                          LASTX   t0  = 1.5507  E10  years = R05

   and   R00 = k = -1  ( hyperbolic space )
           R06 =  R0 = 2.0725  E10  light-years

-For another redshift - in the same Universe - SF 00, place z in X-register and  R/S
 

Example2:      With    (Omega)mat = 0.271 ,  (Omega)lambda  = 0.732 ,  (Omega)rad  = 0.000049  &   1/H0 =  13.7 109
                       ( these parameters are suggested by WMAP = NASA's Wilkinson Microwave Anisotropy Probe )

   0.271  STO 11         0.000049   STO 13
   0.732  STO 12         13.7  E9     STO 14     and, for instance,  FIX 7

     z = 7  XEQ "COSMO"  >>>>   D  =  1.2823287  E10  l-y     = R01
                                            RDN    D0 = 2.8607721  E10  l-y     = R02
                                            RDN    DL = 2.2835499  E11  l-y     = R03
                                            RDN    VR = 2.0881542                   = R04         ( unit c = 1 )
                                          LASTX   t0  = 1.3596706  E10  years = R05

   and   R00 = k = +1  ( spherical space )
           R06 =  R0 = 2.4810862  E11   light-years

Notes:

-Actually, WMAP suggests that the cosmological parameters lead to a flat space  ( k = 0 ), but this is still uncertain...
-In most cases, (Omega)rad  has a negligible effect, except in the very early age of the Universe.
- t0 cannot be computed if there is no Big-Bang.
-Set F00 in this case.

-The value 1/H0 = 13771721429 years would produce - almost - exactly:  Hubble's constant =  H0 = 71 km/s/Mpc
-These numbers are linked by the relation:  ( 1/H0 ) ( years ) = 977.79222143 109 / H0 ( km/s/Mpc )

-We have supposed that our Universe is isotropic & homogeneous, which is only an approximation.
-See reference [7] for a more general point of view.
 

 Calculation of DeltaT = TT-UT for a given year
 

-This routine computes the difference between Terrestrial Time ( or Ephemeris Time ) and Universal Time for a given year y.
-It uses several polynomial approximations, all expressed in seconds:

           y < 2             >>>>    DeltaT ~   32  ( y - 1820 )2  10 -4 - 20
     2 < y < 1600       >>>>    DeltaT ~   1587 - 556 t + 66.7 t2 + 0.33 t3 - 0.523 t4 - t5/291 + t6/533   with   t = y/100 - 10
 1600 < y < 1800     >>>>    DeltaT ~   8 - 16 t + 96 t2 - 75 t3 - 38 t4 + 38 t5                                      with   t = y/100 - 17
 1800 < y < 1998     >>>>    2 polynomials are used.  See reference [1]  page 80.
 1998 < y < 2150     >>>>    DeltaT ~   64 + 35 t - 112 t2 + 444 t3 - 874 t4 + 46.4 t5                          with   t = y/100 - 20
       y > 2150           >>>>    DeltaT ~   32  ( y - 1820 )2  10 -4 - 20
 
 
      STACK        INPUTS      OUTPUTS
           X          year     Delta T (sec)

Examples:

 -2000   XEQ "DELTAT"  >>>>  46675.68 s
   400    XEQ "DELTAT"  >>>>   6702.33 s
  1880   XEQ "DELTAT"  >>>>    -5.42 s
  1984   XEQ "DELTAT"  >>>>     53.82 s
  2010   XEQ "DELTAT"  >>>>     66.80 s
  2150   XEQ "DELTAT"  >>>>    328.48 s

-The accuracy is of the order of 1 second between 1800 & 2010
-For other polynomial approximations, see also:  http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
 

 Julian/Gregorian Calendar > Hindu Calendar
 

-"DT-KY" & "KY-DT" use the Old Hindu Calendars
  that began on  -3101  January  23rd ( Gregorian )
                     = -3101 February 18th ( Julian )
                     =   0.0101 ( Kali Yuga )                    ( -3101 = 3102 BC )

-It is based on a mean sidereal year = 1577917500 / 4320000 = 210389 / 576 = 365 + 149/576  ~ 365.25868 days
-However, if you set flag F01, this value is replaced by a more accurate approximation, i-e 172401/472 ~ 365.2563559 days.
 ( in 2000 , one sidereal year = 365.256363004 days )

-There are 12 months which may have 30 or 31 days.
 

 Formulae:     Given a date   YYYY.MNDD  in the Gregorian or Julian calendar, we first compute the number of days N since 2000/01/01 , then:
 

  Let  s = N + 1863079.25  ,  the date  yyyy.mndd  in the Hindu Calendar is obtained as follows:

     y  = FLOOR ( 576 s / 210389 )
     m = 1 + [ FLOOR ( 6912 s / 210389 ) ] MOD 12
     d  = 1 + FLOOR [ s mod ( 210389 / 6912 ) ]
 

Data Registers:  R00-R01-R02: temp
Flag:   F01-F04   CF 01 = traditional old Hindu Calendar     CF 04 = Gregorian Calendar
                             SF 01 = "new-old" Hindu Calendar           SF 04 =  Julian Calendar
Subroutine:       "J2000"
 
 
           STACK            INPUT         OUTPUT
               X   YYYY.MNDDGr/Ju     YYYY.MNDDKy

Examples:    With CF 04  ( Gregorian Calendar )

    CF 01        Traditional Hindu Calendar

   -3101.0123   XEQ "DT-KY"  >>>>       0.0101
     1979.0716  XEQ "DT-KY"  >>>>    5080.0331
     2010.0404  XEQ "DT-KY"  >>>>    5110.1219

    SF 01        "New-Old" Hindu Calendar

   -3101.0123   XEQ "DT-KY"  >>>>       0.0101
     1979.0716  XEQ "DT-KY"  >>>>    5080.0412
     2010.0404  XEQ "DT-KY"  >>>>    5111.0101

Notes:

-In the Hindu calendar, the day begins at 6 a.m.
-So we also assume that a date written yyyy.mndd in the Julian or Gregorian calendar represents a time after 6 a.m.
-The reverse conversion is performed by  "KY-DT"
 

 Date of Easter
 

-This program takes a year in X-register and returns the date of Easter for this year.
 

Data Registers: /
Flag:    F04               CF 04 = Gregorian calendar
                                   SF 04 = Julian calendar
Subroutines: /
 
 
      STACK       INPUTS      OUTPUTS
           X        YYYY   YYYY.MNDD

Examples:

   CF 04      2010  XEQ "EASTER"  >>>>    2010.0404
   SF 04      1243          R/S               >>>>   1243.0412
 
 

 Solar & Lunar Eclipses
 

-This program calculates the dates of Lunar & Solar eclipses with an accuracy of about 5 minutes between 1000 & 3000 and probably a little more...
 

Data Registers:   R00 = 4.k               R05 = 0 or 1         R10 = 2.F                        R15 = 7.38264721
                               R01 = m                 R06:   temp          R11 = 2.M'                      R16:    temp
                               R02 = m'                R07 = E                R12 = M'+M-180°
                               R03 = gamma         R08 = M' - 180°   R13 = M'-M-180°
                               R04 = phase           R09 = M               R14 = 2M

       For a Solar eclipse, gamma =  the least "distance" from the center of the Earth to the axis of the Moon's shadow ( unit = Earth's radius )
       For a Lunar eclipse, gamma =  the least "distance" from the center of the Moon to the axis of the Earth's shadow ( unit = Earth's radius )

Flag:   F02 is set to indicate a Hybrid ( solar ) eclipse.
Subroutines:  "J2000" & "DT"
 
 
         STACK         INPUT       OUTPUTS
             T              /              m'
             Z              /              m
             Y              /        HH.MNSS
             X       yyyy.mndd     YYYY.MNDD

    where   yyyy.mndd  is an approximate date

        YYYY.MNDD  is the date of the eclipse
          HH.MNSS      is the time of the maximum eclipse

  For a Lunar Eclipse:    m = magnitude of the umbral eclipse          if  m < 0  there is no umbral eclipse     m , m' > 1  for a total eclipse,
                                       m' = magnitude of the penumbral eclipse                                                              m , m' < 1  for a partial eclipse

  For a Solar Eclipse:     m = magnitude of the eclipse                     m > 1  for a total or hybrid eclipse
                                                                                                       m < 1  for a partial eclipse
                                       m' = 0   for a partial eclipse
                                       m' = 1   for an annular eclipse
                                       m' = m  for a total or hybrid eclipse     F02 will be set for a hybrid eclipse

Examples:

   3000.1016   XEQ "ECLIPSE"   >>>>      " SOLAR"  is displayed. Press the backarrow key

                                                                    3000.1019            Flag  F02 is set, so we have a Hybrid ( Solar ) eclipse
                                                    RDN        16.0807
                                                    RDN     m = 1.0027              R03 = gamma = -0.2167
                                                    RDN     m' = 1                      R04 = phase = 0 ( New Moon )

-The exact instant is  3000/1019 at 16h10m16s and the exact m = 1.0049
 

>>>>  To search for the next eclipse, simply press  R/S
 

-In the last example, it yields:   " LUNAR"   press the backarrow key

     3000/11/04     RDN    5h44m53s    RDN    m = -0.4165    RDN    m' = 0.5983       R03 = gamma = -1.2355       R04 = 2  ( Full Moon )

-Since  m < 0 , there is no eclipse in the umbra:
  this Lunar eclipse is a penumbral one.
 

 Positional Fix from the heights of 2 Stars
 

-"FIX" calculates your longitude & latitude from the altitudes of 2 stars ( or planets ... )
-There are usually 2 solutions.
-Several PROMPTs allow you to store the data in the proper registers R01 thru R10.

-We have to solve the non-linear system:

           sin b sin d1 + cos b cos d1 cos ( ST1 - a1 + L ) = sin h1                          (  ST = mean Sidereal Time at Greenwich  ,  hi = altitude )
           sin b sin d2 + cos b cos d2 cos ( ST2 - a2 + L ) = sin h2                         (    L = longitude , b = latitude , ai = right-ascension , di = declination )

 may be solved without any iterative method:

-Let   x = sin b    After eliminating L from these 2 equations, we find that x is one solution of the quadratic equation  A x2 + B x + C = 0

                  A = tan2 d1 + tan2 d2 + sin2 ( ST2 - a2 - ST1 + a1 ) - 2 tan d1 tan d2 cos ( ST2 - a2 - ST1 + a1 )
   with         B = -2 sin h1 sin d1 / cos2 d1 - 2 sin h2 sin d2 / cos2 d2 + 2 ( sin h1 sin d2 + sin h2 sin d1 ) cos ( ST2 - a2 - ST1 + a1 ) / ( cos d1 cos d2 )
                  C = sin2 h1 / cos2 d1 + sin2 h2 / cos2 d2 - sin2 ( ST2 - a2 - ST1 + a1 ) - 2 sin h1 sin h2 cos ( ST2 - a2 - ST1 + a1 ) / ( cos d1 cos d2 )

-After solving this quadratic equation, we find the longitude L by using the R-P conversion and:

                  2 cos [ L + ( ST2 - a2 + ST1 - a1 )/2 ] cos [ ( ST2 - a2 - ST1 + a1 )/2 ] = E + F
                  2 sin  [ L + ( ST2 - a2 + ST1 - a1 )/2 ] sin  [ ( ST2 - a2 - ST1 + a1 )/2 ] =  E - F

   with        E = ( sin h1 - x sin d1 ) / [ cos d1 ( 1 - x2 )1/2 ]
   and         F = ( sin h2 - x sin d2 ) / [ cos d2 ( 1 - x2 )1/2 ]
 

Data Registers:                R00 & R11 thru R19: temp

                   R01 = Date1                   R06 = Date2              YYYY.MNDD
                   R02 = Time1                   R07 = Time2                 HH.MNSS     ( UT )
                   R03 = R.A1                    R08 = R.A2                     hh.mnss
                   R04 = Decl1                   R09 = Decl2                     ° . '  ''
                   R05 = h1                        R10 = h2                           ° . '  ''

Flags: /
Subroutines:  "MST"   "H0-H"
 
 
      STACK        INPUTS      OUTPUTS
           T             /       b2 ( ° ' " )
           Z             /       L2 ( ° ' " )
           Y             /       b1 ( ° ' " )
           X             /       L1 ( ° ' " )

   where   L = longitude ( positive East ) , b = latitude

Example:      on  2006/04/04,

    6h00m  (UT)   alpha Lyr   -  RA = 18h37m09s , Decl = 38°47'21" ,  Altitude = 29°41'59"  ( apparent )
    6h01m  (UT)  zeta UMaj  -  RA = 13h24m10s , Decl = 54°53'34" ,  Altitude = 60°29'56"  ( apparent )

    XEQ "FIX"    >>>>      "DATE^TIME1?"

    2006.0404   ENTER^
           6              R/S        "RA^DECL^H1?"

       18.3709    ENTER^
       38.4721    ENTER^
       29.4159       R/S        "DATE^TIME2?"

    2006.0404   ENTER^
         6.01           R/S        "RA^DECL^H2?"

       13.2410    ENTER^
       54.5354    ENTER^
       60.2956       R/S            It yields:

                                       L1 =  -136°36'29''
                         RDN      b1 =     77°27'26''
                         RDN      L2 =   -74°40'43'
                         RDN      b2 =     25°50'00''

-The 2 solutions are  ( L1 , b1 )  &  ( L2 , b2 )

Notes:

-"FIX" takes the refraction into account and corrects the apparent altitudes by calling "H0-H"

-If you know that you are on the Atlantic Ocean ( near the Bahamas ! ), the correct solution is ( L2 , b2 )
-Otherwise - or to increase accuracy - use another star and take the average of the results.

-The program doesn't work if a declination = +/-90°  but in this case, the system is very easy to solve:
-If, for instance,  d1 = 90° then  b = h1 and after substituting this result into the 2nd equation   sin b sin d2 + cos b cos d2 cos ( ST2 - a2 + L ) = sin h2
  ACOS  easily gives the 2 L-values.

-It doesn't really matter to take this case into account because no star has a declination exactly equal to 90° or -90°
-However, roundoff errors may be important - at least in the longitude - if the declination is about +/-89°59'
-But it still works well for the Pole Star.
 

 Gravitational n-body Problem
 

 We have to solve a system of 3n second order differential equations of the type:   d2y/dt2 = f (y)
 The time t and the first derivative y' = dy/dt  do not appear explicitly in this trajectory problem.
 So, we can use a special 4th order Runge-Kutta method, namely:

  y(t+h) = y(t) + h ( y'(t) + k1/6 + 2k2/3 )
  y'(t+h) = y'(t) + k1/6 + 2k2/3 + k3/6                 where  k1 = h.f (y)  ;  k2 = h.f(y+h.y'/2+h.k1/8)  ;  k3 = h.f(y+h.y'+h.k2/2)

Let   M1( x1, y1,z1) ; ....... ; Mn( xn, yn,zn)  be n celestial bodies  with initial velocities   V1( x'1, y'1,z'1) ; ....... ; Vn( x'n, y'n,z'n)  at an instant t
 and    m1 , ...... , mn   the n masses of these bodies.
We want to know their future positions and velocities at an instant  t + N.h    ( h = step size ; N = number of steps )
 

                                                                        d2xi /dt2 = SUMj#i   Gmj ( xj - xi ) / [ ( xi - xj )2+ ( yi - yj )2 + ( zi - zj )2 ]3/2
   So, we have to solve the system:                    d2yi /dt2 = SUMj#i   Gmj ( yj - yi ) / [ ( xi - xj )2 + ( yi - yj )2 + ( zi - zj )2 ]3/2          i , j = 1 , ..... , n
                                                                        d2zi /dt2 = SUMj#i   Gmj ( zj - zi ) / [ ( xi - xj )2 + ( yi - yj )2 + ( zi - zj )2 ]3/2

 where G is the gravitational constant.
 "GNBP" can solve the 15-body problem.
 

Data Registers:                R00 thru R26+19n  are used
 

      R00 = n = number of bodies                                   R27+n = x1                    R27+4n = x'1
      R24 = G = constant of gravitation                           R28+n = y1                    R28+4n = y'1
      R25 = h = step size                                                R29+n = z1                     R29+4n = z'1
      R26 = N = number of steps
                                                                                    ................                       ...................
      R27 = m1
      R28 = m2                                                               R24+4n = xn                  R24+7n = x'n
      ...............                                                               R25+4n = yn                  R25+7n = y'n
      R26+n =  mn                                                          R26+4n = zn                  R26+7n = z'n

Flag:   F10
Subroutines:  /

>>> Thus, the n masses, the 3n position coordinates and the 3n velocity coordinates are stored into contiguous registers, starting at R27

N.B.

  We use  G = k2  where  k = 0.01720209895 is Gaussian gravitational constant.
  The units are Astronomical Unit , Solar mass and day.
 
 
         STACK          INPUT         OUTPUT
             X              /          bbb.eee

   where  bbb.eee = control number of the solution.

Example:   Let's imagine a system of three stars      M1 ( 2 ; 0 ;0 )       M2 ( 0 ; 4 ; 0 )      M3 ( 0 ; 0 ; 1 )                     unit = 1 AU
                                          with  initial velocities      V1 ( 0 ; 0.03 ; 0 )  V2 ( 0 ; 0 ; 0.01 )  V3 ( -0.02 ; 0 ; 0 )               unit = 1 AU per day
                                                       and masses       m1 = 2  ;  m2 = 1  ;  m3 = 3                                                           unit = 1 Solar mass

-Calculate the positions and velocities 10 days later with stepsize = h = 10 , number of steps = 1
 

    XEQ "GNBP"   >>>>         " M1=?"

                 2           R/S           " M2=?"
                 1           R/S           " M3=?"
                 3           R/S           " M4=?"                     since there are only 3 stars, press R/S without any digit entry

                              R/S           " X1^Y1^Z1"

                 2        ENTER^
                 0        ENTER^
                 0           R/S           "VX1^VY1^VZ1"

                 0        ENTER^
               0.03     ENTER^
                 0           R/S            " X2^Y2^Z2"

                 0        ENTER^
                 4        ENTER^
                 0           R/S           "VX2^VY2^VZ2"

                 0        ENTER^
                 0        ENTER^
               0.01        R/S           " X3^Y3^Z3"

                 0        ENTER^
                 0        ENTER^
                 1           R/S           "VX3^VY3^VZ3"

              -0.02     ENTER^
                 0        ENTER^
                 0           R/S           " H^Nb?"                 Stepsize ENTER Number of steps

                10       ENTER^
                 1           R/S        >>>>    30.047 = R01      ( Turbo mode quite useful ! )

>>>  The new positions and velocities have replaced the previous ones in registers R30 thru R47 and we get:
 
 
Register  input   h = 10 ; N=1      h = 5 ; N=2
 m R27=m1
R28=m2
R29=m3
   2
   1
   3
  undisturbed       undisturbed

 

 p
 o
 s.

 

R30=x1
R31=y1
R32=z1
--------
R33=x2
R34=y2
R36=z2
--------
R36=x3
R37=y3
R38=z3
   2
   0
   0
-----
   0
   4
   0
-----
   0
   0
   1
1.992077590
0.300333861
0.003673761
 --------------
0.000661665
3.996080594
0.100603408
 --------------
-0.194938948
 0.001083895
 0.997349690
    1.992077585
    0.300333570
    0.003673682
     --------------
    0.000661669
    3.996080575
    0.100603412
    --------------
   -0.194938946
    0.001084095
    0.997349741
 v
 e
 l
 .
R39=x'1
R40=y'1
R41=z'1
--------
R42=x'2
R43=y'2
R44=z'2
--------
R45=x'3
R46=y'3
R47=z'3
   0
 0.03
   0
-----
   0
   0
 0.01
-----
-0.02
   0
   0
-0.001550090
 0.030038159
 0.000706688
 --------------
 0.000132598
-0.000790384
 0.010117548
 --------------
-0.019010806
 0.000238022
-0.000510308
  -0.001550083
    0.030038158
    0.000706684
    --------------
    0.000132598
   -0.000790385
    0.010117549
    --------------
   -0.019010811
    0.000238023
   -0.000510306

 

 -To estimate the accuracy of the results, we can perform the calculation with a half step size ( h = 5 days ) and N = 2 instead of 1:
   We obtain the results in the last column: errors are smaller than 0.000001 AU
- To continue, simply press  R/S  the HP-41 displays " H^Nb?"

   Press R/S without any digit entry if you keep the same stepsize & number of steps.

Notes:

-If you get NONEXISTENT during the initialization, increase the SIZE ( at least SIZE 27+19n = 084 in the example above )
-If you apply "GNBP" or "GNB+1" to the whole Solar system,  h = 1day  gives an error of about 7 10-6 AU in the position of Mercury after 88 days.

-The relativitic effects are not taken into account.
 

 Gravitational (n+1)-body Problem
 

 When computing planetary trajectories for instance, it's often better to know directly the heliocentric coordinates by choosing the Sun as the origin.
 In this case however, the reference frame is no more inertial , the equations become more complex and therefore, the program is longer.
 On the other hand, it executes faster and requires less data registers ( the 16-body problem can be solved ). The equations are now:
 

                       d2xi /dt2 =  - Gm0 xi/ri3 - SUMj Gmj xj/rj3 +  SUMj#i   Gmj ( xj - xi ) / rij3
                       d2yi /dt2 =  - Gm0 yi/ri3 - SUMj Gmj yj/rj3 +  SUMj#i   Gmj ( yj - yi ) / rij3                  i , j = 1 , ..... , n
                       d2zi /dt2 =  - Gm0 zi/ri3 - SUMj Gmj zj/rj3 +  SUMj#i   Gmj ( zj - zi ) / rij3
 

              where  ri =  ( xi2 + yi2 + zi2 ) 1/2   and   rij = [ ( xi - xj )2 + ( yi - yj )2 + ( zi - zj )2 ] 1/2
 

Data Registers:               R00 thru R30+19n  are used

      R00 = n = Nb of bodies - 1                                    R31+n = x1                    R31+4n = x'1
      R27 = G = constant of gravitation                           R32+n = y1                    R32+4n = y'1
      R28 = h = step size                                                R33+n = z1                     R33+4n = z'1
      R29 = N = number of steps
                                                                                    ................                       ...................
      R30 = m0 = mass of the origin
      R31 = m1                                                               R28+4n = xn                  R28+7n = x'n
      ...............                                                               R29+4n = yn                  R29+7n = y'n
                                                                                    R30+4n = zn                  R30+7n = z'n
     R30+n = mn

Flag:  F10
Subroutines:  /
 

>>> Thus, the n masses, the 3n position coordinates and the 3n velocity coordinates are to be stored into contiguous registers, starting at R30
 
 
         STACK          INPUT         OUTPUT
             X              /          bbb.eee

   where  bbb.eee = control number of the solution.

Example:     let's take the same 3 stars and let's choose the third star as the origin. The coordinates become:

                                    M1 ( 2 ; 0 ; -1 )    M2 ( 0 ; 4 ; -1 )   and  V1 ( 0.02 ; 0.03 ; 0 )    V2 ( 0.02 ; 0 ; 0.01 )

-Calculate the positions and velocities 10 days later with stepsize = h = 10 , number of steps = 1

    XEQ "GNB+1"  >>>>        " M0=?"

                 3           R/S           " M1=?"
                 2           R/S           " M2=?"
                 1           R/S           " M3=?"                     since there are only 3 stars, press R/S without any digit entry

                              R/S           " X1^Y1^Z1"

                 2        ENTER^
                 0        ENTER^
                -1           R/S           "VX1^VY1^VZ1"

               0.02     ENTER^
               0.03     ENTER^
                 0           R/S            " X2^Y2^Z2"

                 0        ENTER^
                 4        ENTER^
                -1           R/S           "VX2^VY2^VZ2"

               0.02     ENTER^
                 0        ENTER^
               0.01        R/S           " H^Nb?"                 Stepsize ENTER Number of steps

                10       ENTER^
                 1           R/S        >>>>    33.044 = R01      ( Turbo mode quite useful again ! )

 >>>  the new postions and velocities have replaced the previous ones in registers R33 thru R44  (  last column )
 
 
Register  input    h = 10 ; N=1
 m R30=m0
R31=m1
R32=m2
   3
   2
   1
    undisturbed
 p
 o
 s.
R33=x1
R34=y1
R35=z1
--------
R36=x2
R37=y2
R38=z2
   2
   0
  -1
-----
   0
   4
  -1
   2.187016538
   0.299249966
  -0.993675929
  --------------
   0.195600614
   3.994996700
  -0.896746283
 v
 e
 l.
R39=x'1
R40=y'1
R41=z'1
--------
R42=x'2
R43=y'2
R44=z'2
 0.02
 0.03
   0
-----
 0.02
   0
 0.01
   0.017460717
   0.029800137
   0.001216996
  --------------
   0.019143404
  -0.001028406
   0.010627856

Notes:

-If you get NONEXISTENT during the initialization, increase the SIZE ( at least SIZE 31+19n = 069 in the example above )
-If you apply "GNBP" or "GNB+1" to the whole Solar system,  h = 1day  gives an error of about 7 10-6 AU in the position of Mercury after 88 days.

-The relativitic effects are not taken into account.

Execution time:

-Here are the results of a few tests ( for 1 step ) and a real HP-41
 
     n   GNBP  GNB+1
     3   1mn52   1mn28
     5   5mn04   4mn27
    10  19mn50  18mn51

-Since 10 minutes become 1 second with V41, we can integrate the motion of the 9 major planets in a reasonable time !
 

 Illuminated Fraction of the Moon's Disk
 

-Given a date and a time expressed in TT, this short program returns the illuminated fraction of the Moon at that instant.
-The results - though not extremely accurate - are valid over the time span [1000,3000] at least.
 

Data Registers:  R00 = T in millennia since 2000/01/01 0h
                             R01 = D , R02 = M' then M' - 2D
Flags: /
Subroutine:   "J2000"
 
 
      STACK        INPUTS      OUTPUTS
           Y  YYYY.MNDD             /
           X     HH.MNSS             f

  where  f  = illuminated fraction of the Moon's disk  ( 0 < f < 1 )

Examples:    With CF 04 - Gregorian Calendar.

   2010.0716   ENTER^
       7.41         XEQ "IFMD"  >>>>   f = 0.2688  = 27%

   3000.0404   ENTER^
          0              R/S              >>>>    f = 0.4942  = 49%

   1001.0714   ENTER^
      8.3916         R/S              >>>>    f = 1.0000  = 100%   ( almost )
 

 Hindu Calendar > Julian/Gregorian Calendar
 

-See  "DT-KY"  for the definitions.
 

Data Registers:  R00-R01-R02: temp
Flag:   F01-F04   CF 01 = traditional old Hindu Calendar     CF 04 = Gregorian Calendar
                             SF 01 = "new-old" Hindu Calendar           SF 04 =  Julian Calendar
Subroutine:       "DT"
 
 
           STACK            INPUT         OUTPUT
               X    YYYY.MNDDKy   YYYY.MNDDGr/Ju

Examples:    With CF 04  ( Gregorian Calendar )

    CF 01        Traditional Hindu Calendar

        0.0101     XEQ "KY-DT"  >>>>   -3101.0123
     5080.0331  XEQ "KY-DT"  >>>>    1979.0716
     5110.1219  XEQ "KY-DT"  >>>>    2010.0404

    SF 01        "New-Old" Hindu Calendar

       0.0101      XEQ "KY-DT"  >>>>   -3101.0123
     5080.0412  XEQ "KY-DT"  >>>>    1979.0716
     5111.0101  XEQ "KY-DT"  >>>>    2010.0404

-"DT-KY" & "KY-DT" are just separated by a RTN
 

 Local Solar Time
 

-This program calculates the local solar time defined by:  LST  = 12 h. + the hour angle of the Sun.
-12h are added because it's more convenient to regard the solar day as commencing at midnight rather than at noon.
-The equation of time ( i-e the difference between apparent and mean time ) is obtained by the following formula:

     E = y.sin 2L - 2e.sin M + 4ey.sin M cos 2L - (y2/2).sin 4L - (5e2/4).sin 2M       ( in radians )

where

     y = tan2(obl/2)    obl = obliquity of the ecliptic
     L = Sun's mean longitude
     e = eccentricity of the Earth's orbit
    M = Sun's mean anomaly
 

Data Registers:  R00 thru R05
Flag: /
Subroutine:  "J2000"
 
 
             STACK           INPUTS          OUTPUTS
                  Z       YYYY.MNDD                 /
                  Y       HH.MNSS (UT)                 /
                  X      Longitude ( ° . '  " )       LST = HH.MNSS

Example:          In the U.S. Naval Observatory at Washington ( D.C.)  ( Longitude = 77°03'56"W )  on 2010/07/16 at 16h24m41s UT

      2010.0716  ENTER^
          16.2441  ENTER^
        -77.0356  XEQ "LST"   >>>>   LST = 11h10mn21s

Notes:

-Longitudes are measured positively westward from the meridian of Greenwich
-The accuracy is of the order of 1 or 2 seconds between A.D. 1000 & 3000
 

 Passage of the Moon through the Nodes
 

-This program takes an approximate date and returns the time when the Moon passes at the ascending or the descending node of its orbit ( Moon's latitude = 0 ).
-The accuracy is about 5 minutes over the interval [1900,2100]
 

Data Registers:  R00 thru R04: temp      ( R04 is an even integer for the ascending node or an odd integer for the descending node )
Flags: /
Subroutines:  "J2000" & "DT"
 
 
      STACK        INPUTS      OUTPUTS
           Z             /      HH.MNSS
           Y             /   YYYY.MNDD
           X      yyyy.mndd         +/-1

   X-output = +1  for a passage at the ascending node
   X-output = -1  for a passage at the descending node

Examples:

   2008.0825   XEQ "MNODE"  >>>>         1
                                                   RDN    2008.0912
                                                   RDN        18.2530

-The Moon passes through the ascending node on 2008/09/12 at 18h25m TT

   2008.0915          R/S               >>>>       -1
                                                   RDN   2008.0925
                                                   RDN       16.2137

-The Moon passes through the descending node on 2008/09/25 at 16h22m TT
 

 Nutation ( simplified formulae )
 

-The complete series involve more than 1000 terms, and the following program uses the first few terms only:
-The nutation in longitude Dpsi is obtained with an accuracy of  0"5
  and the nutation in obliquity De -----------------------------  0"1
 

Data Register:         R00 = T = number of millenia since 2000/01/01 0h        ( this register is to be initialized before executing "NUT" )
Flags: /
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           Y             /      De ( deg )
           X             /      Dpsi ( deg )

 where  Dpsi  is the nutation in longitude
   and     De   is the nutation in obliquity.

Example:    With  T = 1  ( 3000/01/08  0h )

  1  STO 00   XEQ "NUT"  >>>>   Dpsi = 0.00353°
                                          X<>Y    De  = -0.00198°
 

 Obliquity of the Ecliptic
 

 "OBL" uses the series given in reference [5] to compute the mean obliquity of the ecliptic em
 

Data Register:         R00 = T = number of millenia since 2000/01/01 0h       ( this register is to be initialized before executing "OBL" )
Flags:  /
Subroutines:  /
 
 
      STACK        INPUTS      OUTPUTS
           Y             /            T
           X              /     em ( °. ' " )

  where  em  is the mean obliquity

Example:   With  T = 1  ( 3000/01/08  0h )

    1  STO 00   XEQ "OBL"  >>>>  em = 23°18'35"02

-The accuracy is of the order of  0"01 over the time-span [1000,3000]
 

 Phases of the Moon
 

-The phases of the Moon are usually defined by the difference between the longitudes of the Moon and the Sun ( 0° 90° 180° 270° )
-In the following routine, we use this definition if flag F00 is set, but if flag F00 is clear,
  we use another definition i-e when the illuminated fraction of the Moon's disk equals 0 , 1/2 , 1 , 1/2
-Thus - unlike the standard definition - the instant of a Solar eclipse is identical to the instant of the New Moon
  ( likewise for a Lunar eclipse and a Full Moon )  which is perhaps more logical...
-This program is actually combined with "ECLIPSE"
-Accuracy is about 5 minutes over the time-span 1000-3000.
 

Data Registers:   R00 thru R16: temp
                               R04 = phase

Flag:   F00       SF 00 = Standard definition of phases
                         CF 00 = Non-Standard definition of phases

Subroutines:  "J2000" & "DT"
 
 
      STACK        INPUTS      OUTPUTS
           Z             /      HH.MNSS
           Y             /  YYYY.MNDD
           X      yyyy.mndd        Phase

  where  yyyy.mndd  is an approximate date

               Phase = 0  for a New Moon
               Phase = 1  for a First Quarter
               Phase = 2  for a Full Moon
               Phase = 3  for a Last Quarter

   YYYY.MNDD  &  HH.MNSS  are the date & time ( in Terrestrial Time ) of the Phase

Example:                                             CF 00                                              SF 00

   2008.0801   XEQ "PHASE"   >>>>          0                                                         0
                                                  RDN    2008.0801                                           2008.0801
                                                  RDN        10.2426                                             10.1618

   So, there is a New Moon on 2008/08/01 at 10h24m26s TT                        10h16m18s  TT   with the standard definition

-Likewise, you find a First Quarter on  2008/08/08  at  20h02m45s  TT          20h21m20s  TT   ---------------------------
                                  Full Moon   on   2008/08/16  at  21h12m40s TT          21h19m46s  TT   ---------------------------
                                 Last Quarter on   2008/08/24  at  00h08m41s  TT         23h52m48s  TT  on  2008/08/23  ----------

Note:

-LBL "PHASE" is the 1st line of the program
-LBL "ECLIPSE" is line 132
 

 Rising, Transit, Setting
 

-"RTS" calculates the times of rising, transit and setting of a celestial body on a day D
  and the geometric altitude at the moment of transit.
-The coordinates ( right ascension and declination ) must be obtained independently on days D-1, D , D+1 ( 0h UT )
 

Data Registers:  R00 and R07 thru R15 are used for temporary data storage.       ( Registers R01 to R06 are to be initialized before executing "RTS" )

                              ( on day D-1 )        ( on day D )        ( on day D+1 )         ( at 0h )

                             R01 = RA(-1)     R03 = RA(0)     R05 = RA(+1)          RA = right ascension ( hh.mnss )         R11 = Longitude
                             R02 = decl(-1)    R04 = decl(0)    R06 = decl(+1)         decl = declination  ( ° .  '  " )               R12 = Latitude

Flags:   Set Flag   F00  for the Sun
             Set Flag   F02  for the Moon    ( Clear these 2 flags for the other celestial bodies )

Subroutine:  "MST"
 
 
            STACK            INPUTS         OUTPUTS
                 T                 /   h = altitude* ( ° .  '  " )
                 Z   Date ( YYYY.MNDD )   setting(UT) (hh.mnss)
                 Y     Longitude ( ° .  '  " )   transit(UT)  (hh.mnss)
                 X     latitude     ( ° .  '  " )   rising (UT)  (hh.mnss)

 * where  h = the geometric altitude of the center of the body at the time of transit.

Example1:   Find the times of rising, transit and setting of the Sun at Seattle ( 122°19'51"W ;  47°36'23"N )  on 2005/01/27

-First, we find the right ascension and declination of the Sun on 2005/01/26 , 2005/01/27 , 2005/01/28 at 0h UT  ( cf "Astronomical Ephemeris for the HP-41" )
 and we store these values into R01 thru R06:

       20.3404  STO 01      20.3813  STO 03     20.4221  STO 05
     -18.4417  STO 02     -18.2900  STO 04    -18.1324  STO 06

then,            SF 00  CF 02

            2005.0127  ENTER^
            -122.1951  ENTER^
               47.3623  XEQ "RTS"   >>>>    Rising =  15h41m39s (UT)
                                                    RDN    Transit = 20h22m12s (UT)
                                                    RDN    Setting = 25h03m17s (UT)    ( i-e   1h03m17s  on  2005/01/28 )
       RDN  altitude of the Sun at the time of transit =  24°07'48"

-For the Moon,   CF 00   SF 02

Example2:   Same questions for Sirius.  The right ascension and declination are  6h45m23s  and  -16°43'18"  for this date.

         6.4523   STO 01   STO 03   STO 05       the coordinates may be regarded as constant.
      -16.4318   STO 02   STO 04   STO 06

then,             CF 00   CF 02

             2005.0127   ENTER^      ( or  2005.0127   RCL 11   RCL 12   XEQ "RTS" )
              -122.1951  ENTER^
                 47.3623     R/S            >>>>    Rising =  1h42m05s (UT)
                                                     RDN    Transit =  6h28m09s (UT)
                                                     RDN    Setting = 11h14m14s (UT)
           RDN  altitude of Sirius at the time of transit =  25°40'19"

Notes:

-If the body is circumpolar, there will be no rising and setting: the star will remain above ( or below ) the horizon.
-In this case however, "RTS" yields ( approximately ) the same values in X, Y, Z registers

-For example, with the pole star ( RA = 2h37m39s ; Decl = 89°17'10" )   "RTS" gives   X = Y = Z = 2h21m06s
  but since  T = altitude of the star at the instant of transit = 48°19'13"
  we can deduce that X and Z are meaningless results.
 

 Sidereal & Tropical years and Lunar Months
 

-"STYLM" computes the mean sidereal year, tropical year and lunar ( synodic ) month for a given instant.
-The following formulas are used:

    1 sidereal year = 365.2563630 + 0.0000115 T + 0.0000003 T2
    1 tropical year = 365.2421904 - 0.0006153 T - 0.0000061 T2 + 0.0002652 T3 + 0.0000254 T4 - 0.0000338 T5
    1 lunar month  =  29.53058886 + 0.0002496 T - 0.0000364 T2 + 0.0000235 T3

  where T is measured in unit of 10000 julian years from J2000.0
 

Data Registers:  R00 = sidereal year
                             R01 = tropical year
Flags: /
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           Z             /     lunar month
           Y             /     tropical year
           X          year     sidereal year

   All the results are expressed in days

Examples:

  2000   XEQ "STYLM"   >>>>   365.2563630  days
                                        RDN    365.2421904  days
                                        RDN     29.53058886  days

  5000      R/S      >>>>   365.2563665
                            RDN   365.2420125
                            RDN    29.53059608

 12000     R/S      >>>>   365.2563748
                            RDN   365.2418258
                            RDN    29.53061253

 -3000     R/S      >>>>   365.2563573
                            RDN   365.2424660
                            RDN    29.53057518

 -8000     R/S      >>>>   365.2563518
                            RDN   365.2425936
                            RDN    29.53055791

Notes:

-The input may be fractional.
-The outputs are less accurate before -8000 and after 12000 especially for the lunar months,
  because the tidal acceleration of the Moon is not known with a great precision.
 

 Date & Time for a given longitude of the Sun
 

-"ZODIAC" calculates the instant of a year when the longitude of the Sun equals a given value.
-The nutation & aberration are taken into account.
 

Data Register:   R00 & R01: temp
Flag:  F01
Subroutines:   "SUN"   "DT"
 
 
      STACK        INPUTS      OUTPUTS
           Y         YYYY     hh.mnss (TT)
           X            L   YYYY.MNDD

 where   YYYY =  year  and  L = longitude is expressed in decimal degrees

Examples:     For the same year 2010 , it yields

   2010   ENTER^
      0      XEQ "ZODIAC"   >>>>   2010.0320
                                           X<>Y    17.3202           which is the instant of the spring equinox

   2010   ENTER^
     90        R/S          >>>>      2010.0621
                                X<>Y        11.2910           which is the instant of the summer solstice

   2010   ENTER^
    180       R/S          >>>>      2010.0923
                                X<>Y        3.1257            which is the instant of the autumn equinox

   2010   ENTER^
    270       R/S          >>>>      2010.1221
                                X<>Y        23.4055          which is the instant of the winter solstice

-Of course, L is not necessarily a multiple of 90°, for instance:

   2010   ENTER^
    234       R/S          >>>>      2010.1116
                                X<>Y        11.3131

-However, for an instant before the spring equinox, subtract 360° to L.
-Otherwise, the result would concern the following year !
 

References:

[1]  Jean Meeus - "Astronomical Algorithms" - Willmann-Bell  -  ISBN 0-943396-61-1
[2]  Robin M. Green - "Spherical Astronomy" - Cambridge University Press - ISBN  0-521-31779-7
[3]  John D. Anderson - JPL Technical Report n° 32-497
[4]  Dan Boulet - "Methods of Orbit Determination for the Micro Computer" - Willmann-Bell  -  ISBN 0-943396-34-4
[5]  United States Naval Observatory, Circular n° 179 http://aa.usno.navy.mil/publications/docs/circular_179.html
[6]  Nachum Dershowitz & Edward M. Reingold - "Calendrical Calculations" - Cambridge University Press - ISBN 978-0-521-70238-6
[7]  Kantowski, Kao, Thomas - "Distance-Redshift Relations in Inhomogeneous Friedmann-Lemaitre-Robertson-Walker Cosmology" -
      The Astrophysical Journal, 545: 549-560