Precession

Precession for the HP-41

Overview

1°) Precession Parameters p_A & e_A
2°) Precession Parameters Q_A & P_A
3°) PolynomialApproximations from 8000 B.C. to A.D. 12000
4°) Transformation of Coordinates

      a) Ecliptic Coordinates
      b) Equatorial<>Ecliptic Conversion
      c) Equatorial Coordinates

-Several programs are listed in "Transformation of Coordinates for the HP-41" to take the precession into account.
-They use formulae given in reference [2]
-But these formulas are accurate only a few centuries from J2000.

-The following routines employ the formulae given in reference [1] which are valid 200 thousands of years from J2000
-The accuracy is the same as IAU2006 near J2000 but the errors can reach a few arcminutes 200000 years before or after J2000.
-So, we'll get a much better precision though the execution time will be increased too...

1°) Precession Parameters p_A & e_A

-"PAEA" takes the time expressed in 10000 years from J2000 ( 2000/01/01 12h TT ) in register R00
and returns the mean obliquity of the ecliptic e_A in Y-register and the general precession in longitude p_A in X-register, both in decimal degrees.

Data Registers: • R00 = T = time from J2000 ( Register R00 is to be initialized before executing "PAEA" )

unit = 10000 years = 3652500 days ( -20 < T < +20 )

R01 = p_A R02 = e_A , R03-R04: temp
Flags: /
Subroutines: /

01 LBL "PAEA"
02 DEG
03 RCL 00
04 360
05 *
06 STO 04
07 2.03
08 /
09 STO 03
10 86.021161
11 +
12 COS
13 47604764
14 *
15 STO 01
16 RCL 03
17 4.429265
18 -
19 COS
20 10841074
21 *
22 STO 02
23 RCL 04
24 2.77
25 /
26 STO 03
27 10.754166
28 +
29 COS
30 54765036
31 *
32 ST+ 01
33 RCL 03
34 75.298375
35 -
36 COS

37 19377033
38 *
39 ST+ 02
40 RCL 04
41 3.06
42 /
43 STO 03
44 19.460619
45 +
46 COS
47 63817823
48 *
49 ST- 01
50 RCL 03
51 69.576035
52 -
53 COS
54 22485751
55 *
56 ST- 02
57 RCL 04
58 40.43
59 /
60 STO 03
61 32.946452
62 +
63 COS
64 123082089
65 *
66 ST+ 01
67 RCL 03
68 77.972005
69 -
70 COS
71 84131084
72 *

73 ST- 02
74 RCL 04
75 2.8892
76 /
77 STO 03
78 8.939928
79 -
80 SIN
81 280606764
82 *
83 ST+ 01
84 RCL 03
85 9.088533
86 -
87 COS
88 99471091
89 *
90 ST- 02
91 RCL 04
92 4.1715
93 /
94 STO 03
95 19.648463
96 -
97 COS
98 346046807
99 *
100 ST+ 01
101 RCL 03
102 64.673899
103 +
104 COS
105 58513547
106 *
107 ST- 02
108 RCL 04

109 4.029
110 /
111 STO 03
112 45.952
113 +
114 COS
115 342696586
116 *
117 ST- 01
118 RCL 03
119 3.4661186
120 +
121 SIN
122 247556156
123 *
124 ST- 02
125 RCL 04
126 5.3722
127 /
128 STO 03
129 40.8256734
130 +
131 COS
132 533629325
133 *
134 ST+ 01
135 RCL 03
136 49.72366
137 -
138 COS
139 163051596
140 *
141 ST+ 02
142 RCL 04
143 3.9615
144 /

145 STO 03
146 8.0052016
147 +
148 COS
149 897273066
150 *
151 ST- 01
152 RCL 03
153 73.9666746
154 -
155 COS
156 249224635
157 *
158 ST- 02
159 RCL 04
160 4.099
161 /
162 STO 03
163 22.3842941
164 -
165 COS
166 2075345042
167 *
168 ST- 01
169 RCL 03
170 66.1437076
171 +
172 COS
173 517770288
174 *
175 ST+ 02
176 RCL 00
177 19742583
178 CHS
179 RCL 00
180 75278

181 *
182 +
183 *
184 E9
185 STO 04
186 ST/ 01
187 /
188 140.0847779
189 +
190 *
191 2.259449203
192 +
193 ST+ 01
194 CLX
195 30556
196 *
197 CHS
198 112194
199 -
200 *
201 10067903
202 +
203 *
204 RCL 02
205 +
206 RCL 04
207 /
208 23.34116842
209 +
210 STO 02
211 RCL 01
212 END

( 651 bytes / SIZE 005 )

STACK	INPUTS	OUTPUTS
Y	/	e_A
X	/	p_A

with p_A = general precession in longitude and e_A = obliquity of the ecliptic ( in decimal degrees )

Example: T = 0.8 ( i-e 10000/03/01 12h TT )

0.8 STO 00 XEQ "PAEA" >>>> p_A = 113°4650921 ---Execution time = 32s---
X<>Y e_A = 22°65102781

Note:

-For T = 0 "PAEA" gives p_A = -3 E-9 degree instead of exactly 0 ( roundoff-errors )
-If you want to get p_A = 0 , replace line 191 by 2.259449206 , this will not change the other results significantly.

2°) Precession Parameters Q_A & P_A

-"QAPA" takes the time expressed in 10000 years from J2000 ( 2000/01/01 12h TT ) in register R00
and returns P_A = Sin p_A Sin ¶_A in Y-register and Q_A = Sin p_A Cos ¶_A in X-register.

Data Registers: • R00 = T = time from J2000 ( Register R00 is to be initialized before executing "QAPA" )

unit = 10000 years = 3652500 days ( -20 < T < +20 )

R01 = Q_A R02 = P_A , R03-R04: temp
Flags: /
Subroutines: /

01 LBL "QAPA"
02 RCL 00
03 360
04 *
05 STO 04
06 5.47
07 /
08 STO 03
09 12.750573
10 -
11 COS
12 28803536
13 *
14 STO 01
15 RCL 03
16 69.122621
17 +
18 COS
19 35964826
20 *
21 STO 02
22 RCL 04
23 8.82
24 /
25 STO 03
26 9.938168
27 +
28 COS
29 55921543

30 *
31 ST+ 01
32 RCL 03
33 64.659178
34 -
35 COS
36 56902710
37 *
38 ST- 02
39 RCL 04
40 6.22
41 /
42 STO 03
43 24.635276
44 -
45 COS
46 96811871
47 *
48 ST- 01
49 RCL 03
50 61.676072
51 +
52 COS
53 105812918
54 *
55 ST- 02
56 RCL 04
57 11.83
58 /

59 STO 03
60 11.161703
61 -
62 COS
63 52772369
64 *
65 ST- 01
66 RCL 03
67 66.952667
68 -
69 COS
70 55779703
71 *
72 ST+ 02
73 RCL 04
74 4.922
75 /
76 STO 03
77 49.766163
78 +
79 COS
80 153380984
81 *
82 ST- 01
83 RCL 03
84 42.299904
85 -
86 COS
87 155273469

88 *
89 ST+ 02
90 RCL 04
91 16.2
92 /
93 STO 03
94 37.887736
95 +
96 COS
97 140670994
98 *
99 ST+ 01
100 RCL 03
101 34.7353008
102 -
103 COS
104 208729670
105 *
106 ST- 02
107 RCL 04
108 23.09
109 /
110 STO 03
111 12.6655261
112 +
113 COS
114 696470770
115 *
116 ST+ 01

117 RCL 03
118 .4167179
119 +
120 SIN
121 654152372
122 *
123 ST- 02
124 RCL 04
125 7.0815
126 /
127 STO 03
128 7.0655556
129 +
130 SIN
131 1546147930
132 *
133 ST- 01
134 RCL 03
135 6.9380447
136 +
137 COS
138 1535340316
139 *
140 ST- 02
141 RCL 00
142 556
143 RCL 00
144 121389
145 *

146 +
147 *
148 CHS
149 32471717
150 +
151 *
152 444690639
153 -
154 ST+ 01
155 CLX
156 28056
157 *
158 803139
159 -
160 *
161 3302778
162 -
163 *
164 1625446580
165 +
166 RCL 02
167 +
168 PI
169 18 E10
170 /
171 ST* 01
172 *
173 STO 02
174 RCL 01
175 END

( 532 bytes / SIZE 005 )

STACK	INPUTS	OUTPUTS
Y	/	P_A
X	/	Q_A

with P_A = Sin p_A Sin ¶_A & Q_A = Sin p_A Cos ¶_A

Example: T = 0.8 ( i-e 10000/03/01 12h TT )

0.8 STO 00 XEQ "QAPA" >>>> Q_A = -0.015413403 ---Execution time = 26s---
X<>Y P_A = 0.006932629

from which we deduce ( with R-P ASIN ) p_A = 0°968386023 & ¶_A = 155°7827747

Notes:

-The HP-41 must be set in DEG mode.
-For T = 0, we get Q_A = 2.4 E-11 & P_A = -1.7 E-11 instead of exactly 0.

3°) PolynomialApproximations from 8000 B.C. to A.D. 12000

-If you plan to use these formulae between these 2 dates, polynomial approximations may also be used to get faster results.
-The following ones are valid for -1 < T < +1 ( I found them thanks to my HP-48 )

p_A = 139°688785 T + 3°070652 T² + 0°011186 T³ - 0°646288 T⁴ - 0°023670 T⁵
+ 0°050722 T⁶ + 0°013594 T⁷ - 0°001690 T⁸ - 0°002000 T⁹ - 0°000003 T¹⁰

e_A = 23°4392794 - 1°301021 T - 0°000459 T² + 0°553207 T³ - 0°021415 T⁴ - 0°068157 T⁵
- 0°000784 T⁶ + 0°003674 T⁷ + 0°001093 T⁸ - 0°000081 T⁹ - 0°000141 T¹⁰

¶_A = 174°874109 - 24°111006 T + 0°339983 T² - 0°025435 T³ - 0°017068 T⁴ - 0°000671 T⁵
+ 0°000418 T⁶ + 0°000029 T⁷ - 0°000007 T⁸ + 0 T⁹ - 0°000002 T¹⁰

p_A = 1°305527 T - 0°093051 T² - 0°034510 T³ + 0°002859 T⁴ + 0°000058 T⁵
- 0°000016 T⁶ + 0°000003 T⁷ - 0°000002 T⁸

-Simply store T into R00 before executing these routines and you'll get p_A or e_A or ¶_A or p_A

Data Registers: • R00 = T = time from J2000 ( Register R00 is to be initialized before executing "PA" "EA" "PIA" or "PIa" )

unit = 10000 years = 3652500 days ( -1 < T < +1 )
Flags: /
Subroutines: /

01 LBL "PA"
02 RCL 00
03 2 E3
04 RCL 00
05 3
06 *
07 +
08 *
09 CHS
10 1690
11 -
12 *
13 13594
14 +
15 *
16 50722
17 +
18 *
19 23670
20 -
21 *
22 646288
23 -
24 *
25 11186
26 +
27 *
28 3070652

29 +
30 *
31 139688785
32 +
33 *
34 E6
35 /
36 RTN
37 LBL "EA"
38 RCL 00
39 81
40 RCL 00
41 141
42 *
43 +
44 *
45 CHS
46 1093
47 +
48 *
49 3674
50 +
51 *
52 784
53 -
54 *
55 68157
56 -

57 *
58 21415
59 -
60 *
61 553207
62 +
63 *
64 459
65 -
66 *
67 1301021
68 -
69 *
70 E6
71 /
72 23.4392794
73 +
74 RTN
75 LBL "PIA"
76 RCL 00
77 RCL 00
78 RCL 00
79 2
80 *
81 *
82 CHS
83 7
84 -

85 *
86 29
87 +
88 *
89 418
90 +
91 *
92 671
93 -
94 *
95 17068
96 -
97 *
98 25435
99 -
100 *
101 339983
102 +
103 *
104 24111006
105 -
106 *
107 174874109
108 +
109 E6
110 /
111 RTN
112 LBL "PIa"

113 RCL 00
114 3
115 RCL 00
116 2
117 *
118 -
119 *
120 16
121 -
122 *
123 58
124 +
125 *
126 2859
127 +
128 *
129 34510
130 -
131 *
132 93051
133 -
134 *
135 1305527
136 +
137 *
138 E6
139 /
140 END

( 301 bytes / SIZE 001 )

STACK	INPUTS	OUTPUTS
X	/	p_A or e_A or ¶_A or p_A

All these angles in decimal degrees with -1 < T < +1 in R00

Examples:

• T = 0.987 STO 00

   XEQ "PA"    >>>> p_A = 140°295412
   XEQ "EA"   >>>>   e_A = 22°605889
XEQ "PIA"   >>>>   ¶_A = 151°366870
XEQ "PIa"   >>>>    p_A = 1°167480

• T = -0.987 STO 00

   XEQ "PA"    >>>> p_A = -135°448671
   XEQ "EA"   >>>>   e_A = 24°231402
XEQ "PIA"   >>>>   ¶_A = 199°012112
XEQ "PIa"   >>>>    p_A = -1°343381

Notes:

-If you compare the results given by "QAPA" followed by R-P when T < 0, p_A is positive and ¶_A differs from the above value by 180°
-But this does not change the outputs of the programs hereunder.

-The differences between the outputs of these routines and "QAPA" or "PAEA"are smaller than 0"01 provided | T | doesn't exceed 1
-Don't use these expressions outside the interval [-1,+1]

4°) Transformation of Coordinates

a) Ecliptic Coordinates

-"PREC" uses the routines listed in paragraphs 1 & 2 to perform the transformations.

Data Registers: R00 = T = time from J2000 , unit = 10000 years = 3652500 julian days ( -20 < T < +20 )

R01 thru R09: temp
Flag: F00
Subroutines:

   "J1" or "J2" or "J0" ( cf "Julian & Gregorian Calendars for the HP-41" )
   "QAPA" & "PAEA" ( cf §1 and §2 )
   "EE" ( cf "Transformation of coordinates for the HP-41" )

01 LBL "PREC"
02 DEG
03 HR
04 STO 05
05 X<>Y
06 HR
07 STO 06
08 R^
09 STO 07
10 R^

11 SF 00
12 LBL 00
13 X<> 07
14 XEQ "J2"
15 .5
16 -
17 3652500
18 /
19 STO 00
20 XEQ "QAPA"

21 R-P
22 ASIN
23 STO 09
24 X<>Y
25 STO 08
26 ST- 05
27 XEQ "PAEA"
28 FS? 00
29 ST- 05
30 X<> 09

31 FS? 00
32 CHS
33 RCL 06
34 RCL 05
35 XEQ "EE"
36 RCL 08
37 +
38 STO 05
39 X<>Y
40 STO 06

41 FS?C 00
42 GTO 00
43 HMS
44 X<>Y
45 RCL 09
46 +
47 360
48 MOD
49 HMS
50 END

( 93 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
T	YYYY.MNDD1	/
Z	YYYY.MNDD2	/
Y	b₁ ( ° . ' " )	b₂ ( ° . ' " )
X	l₁ ( ° . ' " )	l₂ ( ° . ' " )

where l = mean ecliptic longitudes , b = mean ecliptic latitudes

Example1:

-On -8000/01/01 0h TT ( gregorian calendar, -8000 = 8001 B.C. ) l₁ = 167°23'45" and b₁ = -12°34'56"
-Calculate the mean ecliptic longitude & latitude of this fictitious star on 200000/01/01 0h TT

   -8000.0101   ENTER^
200000.0101 ENTER^
       -12.3456   ENTER^
       167.2345 XEQ "PREC"   >>>>     191.534869                   ---Execution time = 131s---
                                                X<>Y    -12.194353

-Thus, l₂ = 191°53'48"69 & b₂ = -12°19'43"53

Example2:

-Same inputs but the second date = 12000/01/01 0h TT

   -8000.0101   ENTER^
   12000.0101 ENTER^
       -12.3456   ENTER^
       167.2345 XEQ "PREC"   >>>>       87.070354                    ---Execution time = 131s---
                                                X<>Y    -14.250176

-Thus, l₂ = 87°07'03"54 & b₂ = -14°25'01"76

Note:

-In the second example, we can use the polynomial approximations and the following program to get faster results, almost as accurate.
-This variant employs the polynomials given in paragraph 3.

Data Registers: R00 = T = time from J2000 , unit = 10000 years = 3652500 julian days ( -1 < T < +1 )

R01 thru R05: temp
Flag: F00
Subroutines:

   "J1" or "J2" or "J0" ( cf "Julian & Gregorian Calendars for the HP-41" )
   "PIA" "PIa" & "PA" ( cf §3 )
   "EE" ( cf "Transformation of coordinates for the HP-41" )

01 LBL "PREC"
02 DEG
03 HR
04 STO 01
05 X<>Y
06 HR
07 STO 02
08 R^
09 STO 03
10 R^

11 SF 00
12 LBL 00
13 X<> 03
14 XEQ "J2"
15 .5
16 -
17 3652500
18 /
19 STO 00
20 XEQ "PIA"

21 STO 04
22 ST- 01
23 XEQ "PA"
24 STO 05
25 FS? 00
26 ST- 01
27 XEQ "PIa"
28 FS? 00
29 CHS
30 RCL 02

31 RCL 01
32 XEQ "EE"
33 RCL 04
34 +
35 STO 01
36 X<>Y
37 STO 02
38 FS?C 00
39 GTO 00
40 HMS

41 X<>Y
42 RCL 05
43 +
44 360
45 MOD
46 HMS
47 END

( 90 bytes / SIZE 006 )

STACK	INPUTS	OUTPUTS
T	YYYY.MNDD1	/
Z	YYYY.MNDD2	/
Y	b₁ ( ° . ' " )	b₂ ( ° . ' " )
X	l₁ ( ° . ' " )	l₂ ( ° . ' " )

where l = mean ecliptic longitudes , b = mean ecliptic latitudes

Example:

-On -8000/01/01 0h TT ( gregorian calendar, -8000 = 8001 B.C. ) l₁ = 167°23'45" and b₁ = -12°34'56"
-Calculate the mean ecliptic longitude & latitude of this fictitious star on 12000/01/01 0h TT

   -8000.0101   ENTER^
   12000.0101 ENTER^
       -12.3456   ENTER^
       167.2345 XEQ "PREC"   >>>>      87.070353                    ---Execution time = 27s---
                                                X<>Y    -14.250176

-Thus, l₂ = 87°07'03"53 & b₂ = -14°25'01"76

Note:

-The differences with the version above are only of the order of 0"01

b) Equatorial <> Ecliptic Conversion

Data Registers: R00 = T = time from J2000 , unit = 10000 years = 3652500 julian days ( -20 < T < +20 )

R01 thru R06: temp
Flags: /
Subroutines:

"J1" or "J2" or "J0" ( cf "Julian & Gregorian Calendars for the HP-41" )
"PAEA" ( §1 )
"EE" ( cf "Transformation of coordinates for the HP-41" )

01 LBL "EQECL"
02 DEG
03 HR
04 15
05 *
06 STO 05
07 X<>Y
08 HR
09 STO 06
10 R^
11 XEQ "J2"

12 3652500
13 /
14 STO 00
15 XEQ "PAEA"
16 CLX
17 RCL 06
18 RCL 05
19 XEQ "EE"
20 X<>Y
21 HMS
22 X<>Y

23 360
24 MOD
25 HMS
26 RTN
27 LBL "ECLEQ"
28 DEG
29 HR
30 STO 05
31 RDN
32 HR
33 STO 06

34 X<>Y
35 XEQ "J2"
36 3652500
37 /
38 STO 00
39 XEQ "PAEA"
40 X<>Y
41 CHS
42 RCL 06
43 RCL 05
44 XEQ "EE"

45 X<>Y
46 HMS
47 X<>Y
48 15
49 /
50 24
51 MOD
52 HMS
53 END

( 110 bytes / SIZE 007 )

STACK	INPUT1	OUTPUT1	INPUT2	OUTPUT2
Z	YYYY.MNDD	e_A	YYYY.MNDD	-e_A
Y	decl ( ° . ' " )	b ( ° . ' " )	b ( ° . ' " )	decl ( ° . ' " )
X	RA (hh.mnss)	l ( ° . ' " )	l ( ° . ' " )	RA (hh.mnss)

where RA = right-ascension , decl = declination , l = ecliptic longitude , b = ecliptic latitude

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

10000.0707 ENTER^
41.1657 ENTER^
12.3456 XEQ "EQECL" >>>> l = 168°40'20"09 X<>Y b = 40°45'17"52 ---Execution time = 38s---

Example2:

10000.0707        ENTER^
        40.451752    ENTER^
      168.402009    XEQ "ECLEQ" >>>>   RA = 12h34m56s00    X<>Y   Decl = 41°16'57"00                ---Execution time = 38s---

Note:

-If -1 < T < +1 , lines 39-40 & 15-16 may be replaced by XEQ "EA"

c) Equatorial Coordinates

-"PREQ" calls "EQECL" then "PREC" and finally "ECLEQ"

Data Registers: R00 = T = time from J2000 , unit = 10000 years = 3652500 julian days ( -20 < T < +20 )

R01 thru R11: temp
Flag: F00
Subroutines: "EQECL" "PREC" ( 1st version ) "ECLEQ"

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

( 45 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
T	YYYY.MNDD1	/
Z	YYYY.MNDD2	/
Y	Decl₁ ( ° . ' " )	Decl₂ ( ° . ' " )
X	RA₁ ( ° . ' " )	RA₂ ( ° . ' " )

where RA = right-ascension , Decl = declination

Example:

-On -8000/01/01 0h TT ( gregorian calendar, -8000 = 8001 B.C. ) RA₁ = 12h34m56s and Decl₁ = 41°16'57"
-Calculate the right-ascension & declination of this fictitious star on 12000/01/01 0h TT

   -8000.0101   ENTER^
   12000.0101 ENTER^
        41.1657   ENTER^
        12.3456   XEQ "PREQ"   >>>>         5.312210                    ---Execution time = 3m29s---
                                                X<>Y       61.041700

-So, RA₂ = 5h31m22s10 & Decl₂ = 61°04'17"00

Notes:

-In fact, "PAEA" is called 4 times but 2 times would be enough and execution time could be reduced...

-One can use other precession parameters to get these results more easily
-They are given in reference [1] which also contains the precession matrices ... and so on ... all valid +/- 200,000 years from J2000.0

References:

[1] J. Vondrak, N. Capitaine, P. Wallace - "New precession expressions, valid for long time intervals" - Astronomy & Astrophysics 534, A22 ( 2011 )
[2] Circular n° 179 http://aa.usno.navy.mil/publications/docs/circular_179.html

hp41programs

Precession for the HP-41