The Museum of HP Calculators


Transformation of coordinates and Precession for the HP-41

This program is Copyright © 2002-2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview:
 

 1) 3 Subroutines

    a) Rectangular-Spherical conversion
    b) Spherical-Rectangular conversion
    c) A very useful subroutine: "EE"

 2) Equatorial & Ecliptic Coordinates  ( new )

    a) Equatorial >>> Ecliptic
    b) Ecliptic >>> Equatorial

 3) Equatorial & Azimuthal Coordinates  ( new )

    a) Equatorial >>> Azimuthal
    b) Azimuthal >>> Equatorial

 4) Galactic Coordinates  ( new )

    a) Equatorial >>> Galactic
    b) Galactic >>> Equatorial

 5) Precession   ( updated )   ( new precession formulae are used - cf references [4] & [5] )

    a) Equatorial coordinates
    b) Ecliptic coordinates, program#1
    c) Ecliptic coordinates, program#2
 

1)  3 Subroutines
 

    a) Rectangular-Spherical conversion:
 

                                      x = r cos b cos l
-We have the formulae:   y = r cos b sin l
                                      z = r sin b

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

-However, the results can be obtained more easily by the R-P and P-R functions:
 T-register is saved and no data register is used!
 

01  LBL "R-S"
02  R-P
03  X<>Y
04  RDN
05  R-P
06  R^
07  X<>Y
08  END

( 16 bytes / SIZE 000 )
 
 
      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
 

   b) Spherical-Rectangular conversion:
 

01  LBL "S-R"
02  X<>Y
03  RDN
04  P-R
05  R^
06  X<>Y
07  P-R
08  END

( 16 bytes / SIZE 000 )
 
 
      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

- "R-S" and "S-R" work in all angular modes
 

   c) A very useful subroutine: "EE"
 

-Many transformations 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

- But once again, P-R and R-P lead to a shorter and faster algorithm.
 

01  LBL "EE"
02  1
03  XEQ "S-R"
04  RDN
05  R-P
06  X<> Z
07  ST- Y
08  X<> Z
09  P-R
10  R^
11  XEQ "R-S"
12  RDN
13  END

( 31 bytes / SIZE 000 )
 

-If "EE" is useful for you but "R-S" and "S-R" are not, here is another version of this program that doesn't need any subroutine:
 

01  LBL "EE"
02  1
03  X<>Y
04  RDN
05  P-R
06  R^
07  X<>Y
08  P-R
09  RDN
10  R-P
11  X<> Z
12  ST- Y
13  X<> Z
14  P-R
15  R^
16  R-P
17  X<>Y
18  RDN
19  R-P
20  X<> T
21  END

( 32 bytes / SIZE 000 )
 
 
      STACK      INPUTS      OUTPUTS
           Z            e            e
           Y         decl            b
           X          RA            l

  where  e = obliquity of the ecliptic

Example:     if right-ascension = RA = 116.328942 , declination = decl = 28.026183 and  e = 23.4392911

    23.4392911  ENTER^
    28.026183    ENTER^
  116.328942    XEQ "EE"  >>>> l  =  113.215630   RDN   b  =  6.684170
 

-Like "R-S" and "S-R" , "EE"  works in all angular modes.
-But here are more easy-to-use programs:
 

2) Equatorial & Ecliptic Coordinates
 

    a) Equatorial >>> Ecliptic
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h        R01 & R02: temp
Flag:   F01   Set flag F01 if you are using apparent coordinates
                  Clear flag F01 ---------------   mean     -----------
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "OBL"  ( cf "Nutation - Obliquity of the Ecliptic - Sidereal time for the HP-41" )
                        "NUT"  ( idem )  if SF 01  ;  "EE"
 

01  LBL "EQ-ECL"
02  STO 01
03  RDN
04  STO 02
05  X<>Y
06  XEQ "J0"
07  365250
08  /
09  STO 00
10  XEQ "OBL"
11  LASTX
12  RCL 02
13  HR
14  RCL 01
15  HR
16  15
17  *
18  XEQ "EE"
19  X<>Y
20  HMS
21  X<>Y
22  360
23  MOD
24  HMS
25  END

( 54 bytes SIZE 003 )
 
 
        STACK         INPUTS       OUTPUTS
             Z   YYYY.MNDDdd          e or em
             Y       decl ( . ' " )         b ( . ' " )
             X      RA (hh.mnss)         l ( . ' " )

  where   l  = ecliptic longitude , b  = ecliptic latitude
 

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

          CF 01
    2134.0404  ENTER^
        25.1249  ENTER^
        12.3456  XEQ "EQ-ECL"  >>>>  l  =  17713'44"69    RDN    b  =  2627'18"71
 

Example2:     On  2134/04/04 at 0h  if  RA = 12h34m56s  and  decl =  2512'49"  ( apparent coordinates )

          SF 01
    2134.0404  ENTER^
        25.1249  ENTER^
        12.3456  XEQ "EQ-ECL"  >>>>  l  =  17713'41"39    RDN    b  =  2627'18"39
 
 

    b) Ecliptic >>> Equatorial
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h        R01 & R02: temp
Flag:   F01   Set flag F01 if you are using apparent coordinates
                  Clear flag F01 ---------------   mean     -----------
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "OBL"  ( cf "Nutation - Obliquity of the Ecliptic - Sidereal time for the HP-41" )
                        "NUT"  ( idem )  if SF 01 ;  "EE"
 

01  LBL "ECL-EQ"
02  STO 01
03  RDN
04  STO 02
05  X<>Y
06  XEQ "J0"
07  365250
08  /
09  STO 00
10  XEQ "OBL"
11  LASTX
12  CHS
13  RCL 02
14  HR
15  RCL 01
16  HR
17  XEQ "EE"
18  X<>Y
19  HMS
20  X<>Y
21  15
22  /
23  24
24  MOD
25  HMS
26  END

( 54 bytes / SIZE 003 )
 
 
        STACK         INPUTS       OUTPUTS
             Z   YYYY.MNDDdd          e or em
             Y       decl ( . ' " )         b ( . ' " )
             X      RA (hh.mnss)         l ( . ' " )

  where   l  = ecliptic longitude , b  = ecliptic latitude
 

Example1:     On  2134/04/04 at 0h  if  l  =  17713'44"69  and b  =  2627'18"71  ( mean coordinates )

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

Example2:     On  2134/04/04 at 0h  if  l  =  17713'41"39  and b  =  2627'18"39  ( apparent coordinates )

          SF 01
    2134.0404  ENTER^
        26.271839  ENTER^
      177.134139  XEQ "EQ-ECL"  >>>>  RA = 12h34m56s00    RDN    decl =  2512'49"00
 

Note:    If you know the mean equatorial coordinates and you want to find the apparent equatorial coordinates,

  a)  Execute "EQ-ECL" with CF 01 to calculate the mean ecliptic coordinates
  b)  Add the nutation in longitude to the mean ecliptic longitude
  c)  Execute "ECL-EQ" with SF 01
 

3) Equatorial & Azimuthal Coordinates
 

    a) Equatorial >>> Azimuthal
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h        R01 & R05: temp   R03 = local sidereal time, R04 = hour angle
Flag:   F01   Set flag F01 if you are using apparent coordinates
                  Clear flag F01 ---------------   mean     -----------
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "OBL"  ( cf "Nutation - Obliquity of the Ecliptic - Sidereal time for the HP-41" )
                        "ST"      ( idem )
                        "NUT"  ( idem )  if SF 01 ;  "EE"
 

01  LBL "EQ-AZ"
02  STO 01
03  RDN
04  STO 02
05  RDN
06  XEQ "ST"
07  -77.0356              The East longitude of the US Naval Observatory.  Change this line according to your location or store your longitude in a data register.
08  HR
09  15
10  /
11  HMS
12  HMS+
13  STO 03
14  RCL 01
15  HMS-
16  STO 04
17  HR
18  15
19  *
20  38.55172              The latitude of the US Naval Observatory.  Change this line according to your location or store your latitude in a data register.
21  HR
22  RCL 02
23  HR
24  90
25  ST- Z
26  R^
27  -
28  XEQ "EE"
29  CHS
30  90
31  +
32  X<>Y
33  HMS
34  X<>Y
35  360                              Lines 35 to 42 are not really essential:
36  MOD                           they are only useful to give an azimuth
37  180                              between -180 and +180
38  X<Y?
39  X=0?
40  CLX
41  ST+ X
42  -
43  HMS
44  END

( 84 bytes / SIZE 005 )
 
 
        STACK         INPUTS       OUTPUTS
             T    YYYY.MNDD             /
             Z    HH.MNSS (UT)             /
             Y       decl ( . ' " )    altitude ( . ' " )
             X      RA (hh.mnss)    azimuth ( . ' " )

-Azimuths are measured clockwise ( westwards ) from South
-If you reckon the azimuths clockwise from North, as the navigators prefer, replace the + ( line 31 ) with a  -
 

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

         CF 01
    2005.1212  ENTER^
        20.5129  ENTER^
        60.2137  ENTER^
          7.4116  XEQ "EQ-AZ"  >>>>   Azimuth = -16903'28"33   RDN   Altitude = 1055'57"90

Example2:    On  2005/12/12  at  20h51m29s (UT)  we have  R.A. = 7h41m16s  &  decl = 6021'37"   ( apparent coordinates )

         SF 01
    2005.1212  ENTER^
        20.5129  ENTER^
        60.2137  ENTER^
          7.4116  XEQ "EQ-AZ"  >>>>   Azimuth = -16903'29"87   RDN   Altitude = 1055'57"43
 

-If you want to take the refraction into account, cf for instance "Atmospheric Refraction for the HP-41"
 

    b) Azimuthal >>> Equatorial
 

Data Registers:    R00 = number of millenia since 2000/01/01  0h        R01 & R05: temp   R03 = local sidereal time
Flag:   F01   Set flag F01 if you are using apparent coordinates
                  Clear flag F01 ---------------   mean     -----------
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "OBL"  ( cf "Nutation - Obliquity of the Ecliptic - Sidereal time for the HP-41" )
                        "ST"      ( idem )
                        "NUT"  ( idem )  if SF 01 ;  "EE"
 

01  LBL "EQ-AZ"
02  STO 01
03  RDN
04  STO 02
05  RDN
06  XEQ "ST"
07  -77.0356              The East longitude of the US Naval Observatory.  Change this line according to your location or store your longitude in a data register.
08  HR
09  15
10  /
11  HMS
12  HMS+
13  STO 03
14  38.55172              The latitude of the US Naval Observatory.  Change this line according to your location or store your latitude in a data register.
15  CHS
16  HR
17  RCL 02
18  HR
19  90
20  ST+ Z                   Add   CHS  after line 20 if you measure the azimuths from North.
21  RCL 01
22  HR
23  -
24  XEQ "EE"
25  90
26  -
27  X<>Y
28  HMS
29  X<>Y
30  15
31  /
32  RCL 03
33  HR
34  +
35  24
36  MOD
37  HMS
38  END

( 74 bytes / SIZE 004 )
 
 
        STACK         INPUTS       OUTPUTS
             T    YYYY.MNDD             /
             Z    HH.MNSS (UT)             /
             Y     altitude ( . ' " )       decl ( . ' " )
             X     azimuth ( . ' " )      RA (hh.mnss)

Example1:    On  2005/12/12  at  20h51m29s (UT)  we have  Azimuth = 19056'31"67   &   Altitude = 1055'57"90   ( mean coordinates ),
                      compute the equatorial coordinates at the US Naval Observatory.

         CF 01
    2005.1212      ENTER^
        20.5129      ENTER^
        10.555790  ENTER^
      190.563167  XEQ "AZ-EQ"  >>>>   R.A. = 7h41m16s00   RDN   decl = 6021'37"00

Example2:    On  2005/12/12  at  20h51m29s (UT)  we have  Azimuth = 19056'30"13   &  Altitude = 1055'57"43   ( apparent coordinates )

         SF 01
    2005.1212      ENTER^
        20.5129      ENTER^
        10.555743  ENTER^
      190.563013  XEQ "EQ-AZ"  >>>>   R.A. = 7h41m16s00   RDN   decl = 6021'37"00
 

4) Galactic Coordinates
 

    a) Equatorial >>> Galactic
 

Data Registers:  /
Flags:  /
Subroutine:   "EE"
 

01  LBL "EQ-GAL"
02  DEG
03  HR
04  15
05  *
06  77.140702
07  +
08  X<>Y
09  HR
10  62.871732
11  X<> Z
12  XEQ "EE"
13  X<>Y
14  HMS
15  X<>Y
16  32.93192
17  +
18  360
19  MOD
20  HMS
21  END

( 62 bytes / SIZE 000 )
 
 
        STACK         INPUTS       OUTPUTS
             Y       decl ( . ' " )         b ( . ' " )
             X      RA (hh.mnss)         l ( . ' " )

  where   l  = galactic longitude , b  = galactic latitude

Example:    The coordinates of Procyon are  Right Ascension = 7h39m18s1  ,  Declination = 513'30"  for the equinox of  J2000.0

   5.1330    ENTER^
   7.39181  XEQ "EQ-GAL"  >>>>  l  = 21342'08"19   RDN   b  =  1301'10"16

-In this program, right ascensions and declinations must be referred to J2000.0  ( Julian day = 2451545 )
-If  R.A. & decl are referred to B1950.0 ( Julian day = 2433282.4235 ),  replace lines 06-10-16  with  77.75  ;  62.6  ;  33  respectively.
 

    b) Galactic >>> Equatorial
 

Data Registers:  /
Flags:  /
Subroutine:   "EE"
 

01  LBL "GAL-EQ"
02  DEG
03  HR
04  32.93192
05  -
06  X<>Y
07  HR
08  62.871732
09  CHS
10  X<> Z
11  XEQ "EE"
12  X<>Y
13  HMS
14  X<>Y
15  77.140702
16  -
17  15
18  /
19  24
20  MOD
21  HMS
22  END

( 62 bytes / SIZE 000 )
 
 
        STACK         INPUTS       OUTPUTS
             Y         b ( . ' " )       decl ( . ' " )
             X         ( . ' " )      RA ( hh.mnss )

  where   l  = galactic longitude , b  = galactic latitude

Example:    l  = 21342'08"19   &   b  =  1301'10"16   for the equinox of  J2000.0

      13.011016    ENTER^
    213.420819    XEQ "GAL-EQ"  >>>>  RA  = = 7h39m18s10   RDN   Declination = 513'30"00

-Right ascensions and declinations are referred to J2000.0  ( Julian day = 2451545 )
-If  R.A. & decl are referred to B1950.0 ( Julian day = 2433282.4235 ),  replace lines 04-08-15  with  33 ;  62.6  ;  77.75  respectively.
 
 

5) Precession
 

    a) Equatorial coordinates
 

Data Registers:    R00 thru R03: temp
Flag:   F00
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "EE"
 

  01  LBL "PREQ"
  02  DEG
  03  HR
  04  STO 01
  05  RDN
  06  HR
  07  STO 02
  08  15
  09  ST* 01
  10  R^
  11  STO 00
  12  R^
  13  SF 00
  14  LBL 00
  15  X<> 00
  16  XEQ "J0"
  17  .5
  18  -
  19  365250
  20  /
  21  17
  22  RCL Y
  23  9
  24  *
  25  +
  26  *
  27  5005
  28  -
  29  *
  30  8301
  31  -
  32  *
  33  6405787
  34  -
  35  *
  36  CHS
  37  736
  38  +
  39  STO 03
  40  CLX
  41  8
  42  *
  43  79
  44  +
  45  *
  46  5075
  47  -
  48  *
  49  30354
  50  -
  51  *
  52  6405770
  53  -
  54  *
  55  736
  56  +
  57   E6
  58  ST/ 03
  59  /
  60  FC? 00
  61  X<> 03
  62  90
  63  -
  64  ST+ 01
  65  CLX
  66  5
  67  +
  68  *
  69  4
  70  *
  71  11617
  72  +
  73  *
  74  11930
  75  +
  76  *
  77   E6
  78  /
  79  5.5672
  80  -
  81  *
  82  FC? 00
  83  CHS
  84  RCL 02
  85  RCL 01
  86  XEQ "EE"
  87  90
  88  +
  89  RCL 03
  90  -
  91  STO 01
  92  X<>Y
  93  STO 02
  94  FS?C 00
  95  GTO 00                   ( theoretically a three-byte GTO )
  96  HMS
  97  X<>Y
  98  15
  99  /
100  24
101  MOD
102  HMS
103  END

( 187 bytes / SIZE 004 )
 
 
        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 = -1621'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 =  -1652'22"05
 

    b) Ecliptic coordinates, program#1
 

-One could use similar expressions for the ecliptic coordinates ( see the next program below ).
-But the following program converts ecliptic coordinates into equatorial coordinates,
  then "PREQ" is called and finally, "EQ-ECL" produces the required results.
-This method has at least one advantage: it saves bytes!
-However, it's also slower than direct formulae.
 
 

Data Registers:    R00 thru R04: temp
Flags:   F00 & F01
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "ECL-EQ"  "PREQ"  "EQ-ECL"  "EE"
 

01  LBL "PREC"
02  CF 01
03  R^
04  STO 03
05  X<> T
06  STO 04
07  RDN
08  XEQ "ECL-EQ"
09  RCL 03
10  RCL 04
11  R^
12  R^
13  XEQ "PREQ"
14  X<>Y
15  RCL 04
16  X<> Z
17  XEQ "EQ-ECL"
18  END

( 49 bytes / SIZE 005 )
 
 
        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 = 9830'58"32  ;  b1 =  -3939'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=  10557'44"15   RDN b1=  -3935'19"15
 

    c) Ecliptic coordinates, program#2
 

-Now, we use the specific formulae for ecliptic coordinates.
 

Data Registers:    R00 thru R04: temp
Flag:   F00
Subroutines:  "J0"      ( cf "Phases of the Moon for the HP-41" )
                        "EE"
 

01  LBL "PREC2"
02  DEG
03  HR
04  STO 01
05  X<>Y
06  HR
07  STO 02
08  R^
09  STO 00
10  R^
11  SF 00
12  LBL 00
13  X<> 00
14  XEQ "J0"
15  .5
16  -
17  365250
18  /
19  133
20  CHS
21  RCL Y
22  ENTER^
23  +
24  +
25  *
26  149
27  -
28  *
29  4389
30  +
31  *
32  2410993
33  -
34  *
35  174874109
36  +
37   E6
38  /
39  STO 03
40  ST- 01
41  RDN
42  66
43  -
44  *
45  22
46  +
47  *
48  30707
49  +
50  *
51  13968878
52  +
53  *
54  CHS
55   E6
56  /
57  STO 04
58  FS? 00
59  ST+ 01
60  CLX
61  35
62  *
63  930
64  +
65  *
66  130553
67  -
68  *
69   E6
70  /
71  FC? 00
72  CHS
73  RCL02
74  RCL 01
75  XEQ "EE"
76  RCL 03
77  +
78  STO 01
79  X<>Y
80  STO 02
81  FS?C 00
82  GTO 00         ( a three-byte GTO )
83  HMS
84  X<>Y
85  RCL 04
86  -
87  360
88  MOD
89  HMS
90  END

( 169 bytes / SIZE 005 )
 
 
        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 = 9830'58"32  ;  b1 =  -3939'17"79
                   Calculate the ecliptic coordinates on 2134/12/12 0h

   1600.0404    ENTER^
   2134.1212    ENTER^
  -39.391779   ENTER^
    98.305832   XEQ "PREC2"  >>>> l2=  10557'44"15   RDN b1=  -3935'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 ]
 

References:

   [1] Astronomical Algorithms - Jean Meeus - Willmann-Bell   ISBN 0-943396-35-2
   [2] Spherical Astronomy - Robin M.Green - Cambridge University Press    ISBN  0-521-31779-7
   [3] Introduction aux Ephemerides Astronomiques - EDP Sciences     ISBN  2-86883-298-9  ( in French )
   [4] United States Naval Observatory, Circular n 179 http://aa.usno.navy.mil/publications/docs/circular_179.html
   [5] Report of the IAU division I working group on precession and the ecliptic.
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall