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 =  25°12'49"  ( mean coordinates )

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

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

SF 01
2134.0404  ENTER^
25.1249  ENTER^
12.3456  XEQ "EQ-ECL"  >>>>  l  =  177°13'41"39    RDN    b  =  26°27'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  =  177°13'44"69  and b  =  26°27'18"71  ( mean coordinates )

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

Example2:     On  2134/04/04 at 0h  if  l  =  177°13'41"39  and b  =  26°27'18"39  ( apparent coordinates )

SF 01
2134.0404  ENTER^
26.271839  ENTER^
177.134139  XEQ "EQ-ECL"  >>>>  RA = 12h34m56s00    RDN    decl =  25°12'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 = 60°21'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 = -169°03'28"33   RDN   Altitude = 10°55'57"90

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

SF 01
2005.1212  ENTER^
20.5129  ENTER^
60.2137  ENTER^
7.4116  XEQ "EQ-AZ"  >>>>   Azimuth = -169°03'29"87   RDN   Altitude = 10°55'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 = 190°56'31"67   &   Altitude = 10°55'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 = 60°21'37"00

Example2:    On  2005/12/12  at  20h51m29s (UT)  we have  Azimuth = 190°56'30"13   &  Altitude = 10°55'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 = 60°21'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 = 5°13'30"  for the equinox of  J2000.0

5.1330    ENTER^
7.39181  XEQ "EQ-GAL"  >>>>  l  = 213°42'08"19   RDN   b  =  13°01'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 l  ( ° . ' " ) RA ( hh.mnss )

where   l  = galactic longitude , b  = galactic latitude

Example:    l  = 213°42'08"19   &   b  =  13°01'10"16   for the equinox of  J2000.0

13.011016    ENTER^
213.420819    XEQ "GAL-EQ"  >>>>  RA  = = 7h39m18s10   RDN   Declination = 5°13'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 = -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

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 = 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 b1=  -39°35'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 = 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 "PREC2"  >>>> l2=  105°57'44"15   RDN b1=  -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 ]

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.