Precession

# Precession for the HP-41

Overview

1°)  Precession Parameters pA & eA
2°)  Precession Parameters QA & PA
3°)  Polynomial Approximations 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 pA & eA

-"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 eA in Y-register and the general precession in longitude pA 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 = pA  R02 = eA ,  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 / eA X / pA

with  pA = general precession in longitude  and eA = obliquity of the ecliptic ( in decimal degrees )

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

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

Note:

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

2°)  Precession Parameters QA & PA

-"QAPA" takes the time expressed in 10000 years from J2000 ( 2000/01/01 12h TT ) in register R00
and returns PA = Sin pA Sin ¶A  in Y-register and QA = Sin pA 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 = QA  R02 = PA ,  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 / PA X / QA

with  PA = Sin pA Sin ¶A  &  QA = Sin pA Cos ¶A

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

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

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

Notes:

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

3°)  Polynomial Approximations 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 )

pA = 139°688785 T + 3°070652 T2 + 0°011186 T3 - 0°646288 T4  - 0°023670 T5
+ 0°050722 T6 + 0°013594 T7 - 0°001690 T8 - 0°002000 T9 - 0°000003 T10

eA = 23°4392794 - 1°301021 T - 0°000459 T2 + 0°553207 T3 - 0°021415 T4  - 0°068157 T5
- 0°000784 T6 + 0°003674 T7 + 0°001093 T8 - 0°000081 T9 - 0°000141 T10

¶A = 174°874109 - 24°111006 T + 0°339983 T2 - 0°025435 T3 - 0°017068 T4  - 0°000671 T5
+ 0°000418 T6 + 0°000029 T7 - 0°000007 T8 + 0 T9 - 0°000002 T10

pA = 1°305527 T - 0°093051 T2 - 0°034510 T3 + 0°002859 T4  + 0°000058 T5
- 0°000016 T6 + 0°000003 T7 - 0°000002 T8

-Simply store T into R00 before executing these routines and you'll get  pA or  eA or  ¶A or  pA

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 / pA or eA or ¶A or pA

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

Examples:

•  T = 0.987  STO 00

XEQ "PA"    >>>>  pA = 140°295412
XEQ "EA"   >>>>   eA =  22°605889
XEQ "PIA"   >>>>   ¶A = 151°366870
XEQ "PIa"   >>>>    pA =  1°167480

•  T = -0.987  STO 00

XEQ "PA"    >>>>  pA = -135°448671
XEQ "EA"   >>>>   eA =  24°231402
XEQ "PIA"   >>>>   ¶A = 199°012112
XEQ "PIa"   >>>>    pA =  -1°343381

Notes:

-If you compare the results given by "QAPA" followed by  R-P when  T < 0,  pA 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 b1 ( ° . ' " ) b2 ( ° . ' " ) X l1  ( ° . ' " ) l2  ( ° . ' " )

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

Example1:

-On  -8000/01/01 0h TT ( gregorian calendar, -8000 = 8001 B.C. )   l1 = 167°23'45"  and b1 = -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,   l2 = 191°53'48"69  &  b2 = -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,   l2 = 87°07'03"54  &  b2 = -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 b1 ( ° . ' " ) b2 ( ° . ' " ) X l1  ( ° . ' " ) l2  ( ° . ' " )

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

Example:

-On  -8000/01/01 0h TT ( gregorian calendar, -8000 = 8001 B.C. )   l1 = 167°23'45"  and b1 = -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,   l2 = 87°07'03"53  &  b2 = -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 eA YYYY.MNDD -eA 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 Decl1 ( ° . ' " ) Decl2 ( ° . ' " ) X RA1  ( ° . ' " ) RA2  ( ° . ' " )

where   RA = right-ascension , Decl = declination

Example:

-On  -8000/01/01 0h TT ( gregorian calendar, -8000 = 8001 B.C. )   RA1 = 12h34m56s  and  Decl1 = 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,  RA2 = 5h31m22s10  &  Decl2 = 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