# hp41programs

Gravitation2 The Gravitational n-body problem(2) for the HP-41

Overview

1°)  Newtonian Gravitation - Inertial Frame of Reference

a) 4th-order Runge-Kutta Formula
b) 6th-order Runge-Kutta Formula
c) Up to n=30 ? 3rd-order RKN !

2°)  Newtonian Gravitation - Heliocentric Coordinates

a) 4th-order Runge-Kutta Formula
b) 6th-order Runge-Kutta Formula

3°)  General Relativity - Solar Systems

a) 4th-order Runge-Kutta Formula
b) 5th-order Runge-Kutta Formula

4°)  General Relativity - Barycentric Coordinates

a) 4th-order Runge-Kutta Formula
b) 5th-order Runge-Kutta Formula

-The programs listed in paragraphs 1a & 2a are simply improvements of 2 similar programs listed in"The Gravitational n-body Problems(1) for the HP41"
-They are shorter, faster and the calculations are performed so that the routines can solve the problems up to n = 23 or 24.
but only 18 planets if we take the relativistic perturbations into account (§3).

-Runge-Kutta-Nystrom formulas of order 6 are more accurate but require more data registers. So "GM6" can solve 13-body problems.

-On the other hand, we can also use a third-order Runge-Kutta-Nystrom formula:
-With a few tricks, the 30-body gravitational problem may be integrated by a HP41.

1°)  Newtonian Gravitation - Inertial Frame of Reference

a) 4th-order Runge-Kutta Formula

-"GM4" employs a 4th-order Runge-Kutta formula to calculate the positions of n planets in a Solar system
-Since the differential equations are of the type y" = f(x,y) , we can use a 3-stage method:

y(t+h) = y(t) + h ( y'(t) + k1/6 + 2k2/6 )
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)

-The equations are ( with x , y , z = rectangular coordinates in  AU of the planets and  x' , y' , z' = the speeds of the planets in AU/day ):

d2xi /dt2 = SUMj#i   Gmj ( xj - xi ) / [ ( xi - xj )2+ ( yi - yj )2 + ( zi - zj )2 ]3/2
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

G = gravitational constant = 0.017202098952  = 1 / 3379.380681               mi = mass of the body i

Data Registers:              R00 = G.h           ( Registers R01-R02-R03 & R15 thru R15+7n are to be initialized before executing "GM4" )

•  R01 = n = nb of celestial bodies         •  R15 = t0 = initial instant                                    R04 to R14: temp
•  R02 = h = step size  ( in days )
•  R03 = N = nb of steps                       •  R16 = m1 = mass of the 1st body
................................................

•  R15+n = mn = mass of the n-th body

•  R16+n = x1              •  R16+4n = x'1
•  R17+n = y1              •  R17+4n = y'1                       R16+7n to R15+10n = k2 , k3
•  R18+n = z1              •  R18+4n = z1                       R16+10n to R15+13n = k1 ,  k1 - k2
.................                   ..................

•  R13+4n = xn            •  R13+7n = x'n
•  R14+4n = yn            •  R14+7n = y'n
•  R15+4n = zn            •  R15+7n = z'n
Flags: /
Subroutines: /

-Lines 29-189 are three-byte GTOs

 01 LBL "GM4"  02 RCL 03  03 STO 14  04 RCL 02  05 3379.380681   06 /  07 STO 00  08 15.015  09 RCL 01  10 +  11 STO 05  12 LASTX  13 .1  14 %  15 X<>Y  16 3  17 *  18 STO T  19 +  20 +  21 STO 04  22 INT  23 +  24 STO 11  25 +  26 STO 12  27 + 28 STO 13  29 GTO 01  30 LBL 10  31 RCL 04  32 STO 06  33 LBL 03  34 RCL 04  35 STO 08  36 RCL 05  37 STO 07  38 RCL 10  39 ENTER  40 CLX  41 STO IND Y  42 DSE Y  43 STO IND Y  44 DSE Y  45 STO IND Y  46 LBL 04  47 RCL IND 08  48 RCL IND 06  49 -  50 ENTER  51 X^2  52 DSE 06  53 DSE 08  54 RCL IND 08 55 RCL IND 06  56 -  57 STO 09  58 X^2  59 +  60 DSE 06  61 DSE 08  62 RCL IND 08  63 RCL IND 06  64 -  65 STO T  66 X^2  67 +  68 ENTER  69 SQRT  70 *  71 X=0?  72 SIGN  73 RCL 00  74 X<>Y  75 /  76 RCL IND 07  77 *  78 ST* 09  79 ST* T  80 *  81 ST+ IND 10 82 DSE 10  83 RCL 09  84 ST+ IND 10  85 DSE 10  86 R^  87 ST+ IND 10  88 CLX  89 SIGN  90 ST- 07  91 ST+ X  92 ST+ 06  93 ST+ 10  94 DSE 08  95 GTO 04  96 ST- 06  97 PI  98 INT  99 ST- 10 100 DSE 06 101 GTO 03 102 RCL 13 103 STO 06 104 RCL 11 105 STO 07 106 RCL 04 107 STO 08 108 RTN 109 LBL 01 110 RCL 13 111 STO 10 112 XEQ 10 113 LBL 05 114 RCL IND 06 115 4 116 / 117 RCL IND 07 118 + 119 2 120 / 121 RCL 02 122 * 123 ST+ IND 08 124 CLX 125 SIGN 126 ST- 06 127 ST- 07 128 DSE 08 129 GTO 05 130 XEQ 10 131 RCL 12 132 STO 09 133 LBL 07 134 RCL IND 09 135 RCL IND 06 136 4 137 / 138 - 139 RCL IND 07 140 + 141 RCL 02 142 * 143 2 144 / 145 ST+ IND 08 146 RCL IND 06 147 RCL IND 09 148 ST- IND 06 149 4 150 * 151 + 152 6 153 / 154 ST+ IND 07 155 CLX 156 SIGN 157 ST- 06 158 ST- 07 159 ST- 09 160 DSE 08 161 GTO 07 162 RCL 12 163 STO 10 164 XEQ 10 165 RCL 12 166 STO 09 167 LBL 08 168 RCL IND 06 169 RCL 02 170 * 171 6 172 / 173 ST+ IND 08 174 RCL IND 09 175 LASTX 176 / 177 ST+ IND 07 178 CLX 179 SIGN 180 ST- 06 181 ST- 07 182 ST- 09 183 DSE 08 184 GTO 08 185 RCL 02 186 ST+ 15 187 VIEW 15 188 DSE 14 189 GTO 01 190 RCL 15 191 END

( 292 bytes / SIZE 16+13.n )

 STACK INPUT OUTPUT X / t0 + N.h

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/day
and masses                        m1 = 2  ;  m2 = 1  ;  m3 = 3                                                           unit = 1 Solar mass

-Calculate the positions and velocities 10 days later.

3     bodies                    >>>                                        3     STO 01
step size   h = 10 days   >>>                                       10    STO 02
number of steps:    1      >>>                                        1     STO 03

initial time     0  STO 15

then, we store the 3 masses and the 18 position and velocity coordinates as shown below  and  XEQ "GM4"

>>>  the new positions and velocities have replaced the previous ones in registers R19 thru R36  ( next to last column )       Execution time = 1mn26s

 Register input h = 10 ; N=1 h = 5 ; N=2 m R16=m1 R17=m2 R18=m3 2    1    3 undisturbed undisturbed p  o  s. R19=x1 R20=y1 R21=z1 -------- R22=x2 R23=y2 R24=z2 -------- R25=x3 R26=y3 R27=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  . R28=x'1 R29=y'1 R30=z'1 -------- R31=x'2 R32=y'2 R33=z'2 -------- R34=x'3 R35=y'3 R36=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

-With  h = 5 days & N = 2 steps, we get the results listed above in the last column.
-Simply press R/S to continue with the same h & N

Note:

-"GM4" can solve a gravitational 23-body problem.

b) 6th-order Runge-Kutta Formula

-The following program employs Albrecht method:

k1 = h.f(x0,y0)
k2 = h.f(x0+ h/4 , y0+ (h/4) y'0+ (h/32) k1)
k3 = h.f(x0+ h/2 , y0+ (h/2) y'0+ h (-k1/24 + k2/6 )
k4 = h.f(x0+ 3h/4 , y0+ (3h/4) y'0+ h (3k1/32 + k2/8 + k3/16 )
k5 = h.f(x0+ h , y0+ h y'0+ h (3k2/7 - k3/14 + k4/7 )

-Then,    y = y0 + h y'0+ h2 ( 7 k1/90 + 4 f2/15 + f3/15 + 4 f4/45 ) + O(h7)
and      y' = y'0 + h ( 7 f1/90 + 16 f2/45 + 2 f3/15 + 16 f4/45 + 7 f5/90 ) + O(h7)

Data Registers:              R00 = G.h           ( Registers R01-R02-R03 & R20 thru R20+7n are to be initialized before executing "GM6" )

•  R01 = n = nb of celestial bodies         •  R20 = t0 = initial instant                            R04 to R18: temp     R19: unused
•  R02 = h = step size  ( in days )
•  R03 = N = nb of steps                       •  R21 = m1 = mass of the 1st body
................................................

•  R20+n = mn = mass of the n-th body

•  R21+n = x1              •  R21+4n = x'1                       R21+7n to R20+10n = k4
•  R22+n = y1              •  R22+4n = y'1                      R21+10n to R20+13n = k3
•  R23+n = z1              •  R23+4n = z1                       R21+13n to R20+16n = k2 , 7.k1/90-17k2/105
.................                   ..................                         R21+16n to R20+19n = k1 , k5

•  R18+4n = xn            •  R18+7n = x'n
•  R19+4n = yn            •  R19+7n = y'n
•  R20+4n = zn            •  R20+7n = z'n
Flags: /
Subroutines: /

-Lines 33-309 are three-byte GTOs

 01 LBL "GM6"  02 RCL 03  03 STO 18  04 RCL 02  05 3379.380681   06 /  07 STO 00  08 20.02  09 RCL 01  10 +  11 STO 05  12 LASTX  13 .1  14 %  15 X<>Y  16 3  17 *  18 STO T  19 +  20 +  21 STO 04  22 INT  23 +  24 STO 11  25 +  26 STO 15  27 +  28 STO 14  29 +  30 STO 13  31 +  32 STO 12  33 GTO 01  34 LBL 10  35 RCL 04  36 STO 06  37 LBL 03  38 RCL 04  39 STO 08  40 RCL 05  41 STO 07  42 RCL 10  43 ENTER  44 CLX 45 STO IND Y  46 DSE Y  47 STO IND Y  48 DSE Y  49 STO IND Y  50 LBL 04  51 RCL IND 08  52 RCL IND 06  53 -  54 ENTER  55 X^2  56 DSE 06  57 DSE 08  58 RCL IND 08  59 RCL IND 06  60 -  61 STO 09  62 X^2  63 +  64 DSE 06  65 DSE 08  66 RCL IND 08  67 RCL IND 06  68 -  69 STO T  70 X^2  71 +  72 ENTER  73 SQRT  74 *  75 X=0?  76 SIGN  77 RCL 00  78 X<>Y  79 /  80 RCL IND 07  81 *  82 ST* 09  83 ST* T  84 *  85 ST+ IND 10  86 DSE 10  87 RCL 09  88 ST+ IND 10 89 DSE 10  90 R^  91 ST+ IND 10  92 CLX  93 SIGN  94 ST- 07  95 ST+ X  96 ST+ 06  97 ST+ 10  98 DSE 08  99 GTO 04 100 ST- 06 101 PI 102 INT 103 ST- 10 104 DSE 06 105 GTO 03 106 RCL 12 107 STO 06 108 RCL 11 109 STO 07 110 RCL 04 111 STO 08 112 RTN 113 LBL 01 114 RCL 12 115 STO 10 116 XEQ 10 117 LBL 05 118 RCL IND 06 119 8 120 / 121 RCL IND 07 122 + 123 RCL 02 124 * 125 4 126 / 127 ST+ IND 08 128 CLX 129 SIGN 130 ST- 06 131 ST- 07 132 DSE 08 133 GTO 05 134 XEQ 10 135 RCL 13 136 STO 09 137 LBL 06 138 RCL IND 09 139 6 140 / 141 RCL IND 06 142 7 143 * 144 96 145 / 146 - 147 RCL IND 07 148 4 149 / 150 + 151 RCL 02 152 * 153 ST+ IND 08 154 CLX 155 SIGN 156 ST- 06 157 ST- 07 158 ST- 09 159 DSE 08 160 GTO 06 161 XEQ 10 162 RCL 13 163 STO 09 164 RCL 14 165 STO 16 166 LBL 07 167 RCL IND 06 168 13 169 * 170 24 171 / 172 RCL IND 09 173 6 174 / 175 - 176 RCL IND 16 177 4 178 / 179 + 180 RCL IND 07 181 + 182 RCL 02 183 * 184 4 185 / 186 ST+ IND 08 187 CLX 188 SIGN 189 ST- 06 190 ST- 07 191 ST- 09 192 ST- 16 193 DSE 08 194 GTO 07 195 XEQ 10 196 RCL 13 197 STO 09 198 RCL 14 199 STO 16 200 RCL 15 201 STO 17 202 LBL 08 203 RCL IND 09 204 17 205 * 206 56 207 / 208 RCL IND 06 209 3 210 * 211 32 212 / 213 - 214 RCL IND 16 215 15 216 * 217 112 218 / 219 - 220 RCL IND 17 221 7 222 / 223 + 224 RCL IND 07 225 4 226 / 227 + 228 RCL 02 229 * 230 ST+ IND 08 231 RCL IND 06 232 7 233 * 234 90 235 / 236 ST+ IND 07 237 RCL IND 09 238 17 239 * 240 105 241 / 242 - 243 X<> IND 09 244 16 245 * 246 45 247 / 248 ST+ IND 07 249 CLX 250 SIGN 251 ST- 06 252 ST- 07 253 ST- 09 254 ST- 16 255 ST- 17 256 DSE 08 257 GTO 08 258 RCL 12 259 STO 10 260 XEQ 10 261 RCL 13 262 STO 09 263 RCL 14 264 STO 10 265 RCL 15 266 STO 16 267 LBL 12 268 RCL IND 10 269 87 270 * 271 RCL IND 16 272 34 273 * 274 - 275 630 276 / 277 RCL IND 09 278 + 279 RCL 02 280 * 281 ST+ IND 08 282 RCL IND 10 283 12 284 * 285 RCL IND 16 286 32 287 * 288 + 289 RCL IND 06 290 7 291 * 292 + 293 90 294 / 295 ST+ IND 07 296 CLX 297 SIGN 298 ST- 06 299 ST- 07 300 ST- 09 301 ST- 10 302 ST- 16 303 DSE 08 304 GTO 12 305 RCL 02 306 ST+ 20 307 VIEW 20 308 DSE 18 309 GTO 01 310 RCL 20 311 END

( 474 bytes / SIZE 21+19.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:   The same used by "GM4":

-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/day
and masses                        m1 = 2  ;  m2 = 1  ;  m3 = 3                                                           unit = 1 Solar mass

-Calculate the positions and velocities 10 days later.

3     bodies                   >>>                                        3     STO 01
step size   h = 10 days   >>>                                       10    STO 02
number of steps:    1      >>>                                        1     STO 03

initial time     0  STO 19

then, we store the 3 masses and the 18 position and velocity coordinates as shown below  and  XEQ "GM6"

>>>  the new positions and velocities have replaced the previous ones in registers R24 thru R41  ( next to last column )       Execution time = 1mn52s

 Register input h = 10 ; N=1 h = 5 ; N=2 m R21=m1 R22=m2 R23=m3 2    1    3 undisturbed undisturbed p  o  s. R24=x1 R25=y1 R26=z1 -------- R27=x2 R28=y2 R29=z2 -------- R30=x3 R31=y3 R32=z3 2    0    0 -----    0    4    0 -----    0    0    1 1.992077587 0.300333550 0.003673675  -------------- 0.000661669 3.996080573 0.100603412  -------------- -0.194938948 0.001084109 0.997349746 1.992077588     0.300333550     0.003673676      --------------     0.000661669     3.996080575     0.100603412     --------------    -0.194938948     0.001084109     0.997349746 v  e  l  . R33=x'1 R34=y'1 R35=z'1 -------- R36=x'2 R37=y'2 R38=z'2 -------- R39=x'3 R40=y'3 R41=z'3 0  0.03    0 -----    0    0  0.01 ----- -0.02    0    0 -0.001550083 0.030038158 0.000706684  -------------- 0.000132598 -0.000790385 0.010117549  -------------- -0.019010811  0.000238023 -0.000510306 -0.001550083     0.030038158     0.000706684     --------------     0.000132598    -0.000790385     0.010117549     --------------    -0.019010811     0.000238023    -0.000510306

-With  h = 5 days & N = 2 steps, we get the results listed above in the last column.
-Simply press R/S to continue with the same h & N

Note:

-"GM6" may solve a gravitational 13-body problem.

c) Up to n=30 ? 3rd-order RKN !

-"GM3" uses a Runge-Kutta-Nystrom formula of order 3 with only 2 stages:

k1 = h.f (y)                                             y(t+h) = y(t) + h [ y' + ( k1 + k2 )/4 ]
k2 = h.f (y+2.h.y'/3+2.h.k1/9)                 y'(t+h) = y'(t) + k1/4 + 3.k2/4

Data Registers:              R00 = G.h           ( Registers R01-R02-R03 & R17 thru R17+7n are to be initialized before executing "GM3" )

•  R01 = n = nb of celestial bodies         •  R17 = t0 = initial instant                            R04 to R16: temp
•  R02 = h = step size  ( in days )
•  R03 = N = nb of steps                       •  R18 = m1 = mass of the 1st body
................................................

•  R17+n = mn = mass of the n-th body

•  R18+n = x1              •  R18+4n = x'1                       R18+7n to R17+10n = k1
•  R19+n = y1              •  R19+4n = y'1
•  R20+n = z1              •  R20+4n = z1
.................                   ..................

•  R15+4n = xn            •  R15+7n = x'n
•  R16+4n = yn            •  R16+7n = y'n
•  R17+4n = zn            •  R17+7n = z'n
Flag:  F10
Subroutines: /

-Lines 27-189 are three-byte GTOs

 01 LBL "GM3"  02 RCL 03  03 STO 16  04 RCL 02  05 3379.380681   06 /  07 STO 00  08 17.017  09 RCL 01  10 +  11 STO 05  12 LASTX  13 .1  14 %  15 X<>Y  16 3  17 *  18 STO T  19 +  20 +  21 STO 04  22 INT  23 +  24 STO 11  25 +  26 STO 12  27 GTO 01 28 LBL 10  29 RCL 04  30 STO 06  31 RCL 12  32 STO 10  33 RCL 11  34 STO 13  35 LBL 03  36 RCL 04  37 STO 08  38 RCL 05  39 STO 07  40 FC? 10  41 GTO 04  42 RCL 10  43 ENTER  44 CLX  45 STO IND Y  46 DSE Y  47 STO IND Y  48 DSE Y  49 STO IND Y  50 LBL 04  51 RCL IND 08  52 RCL IND 06  53 -  54 ENTER 55 X^2  56 DSE 06  57 DSE 08  58 RCL IND 08  59 RCL IND 06  60 -  61 STO 09  62 X^2  63 +  64 DSE 06  65 DSE 08  66 RCL IND 08  67 RCL IND 06  68 -  69 STO T  70 X^2  71 +  72 ENTER  73 SQRT  74 *  75 X=0?  76 SIGN  77 RCL 00  78 X<>Y  79 /  80 RCL IND 07  81 * 82 ST* 09  83 ST* Z  84 *  85 RCL X  86 RCL 14  87 *  88 ST+ IND 13  89 CLX  90 RCL 15  91 *  92 ST+ IND 10  93 DSE 10  94 DSE 13  95 CLX  96 RCL 09  97 RCL 14  98 *  99 ST+ IND 13 100 X<> 09 101 RCL 15 102 * 103 ST+ IND 10 104 DSE 10 105 DSE 13 106 CLX 107 RCL 14 108 * 109 ST+ IND 13 110 CLX 111 RCL 15 112 * 113 ST+ IND 10 114 CLX 115 SIGN 116 ST- 07 117 ST+ X 118 ST+ 06 119 ST+ 10 120 ST+ 13 121 DSE 08 122 GTO 04 123 ST- 06 124 PI 125 INT 126 ST- 10 127 ST- 13 128 DSE 06 129 GTO 03 130 RCL 12 131 STO 06 132 RCL 11 133 STO 07 134 RCL 04 135 STO 08 136 RTN 137 LBL 01 138 4 139 1/X 140 STO 14 141 SIGN 142 STO 15 143 SF 10 144 XEQ 10 145 LBL 05 146 RCL IND 06 147 RCL IND 07 148 12 149 * 150 + 151 RCL 02 152 * 153 18 154 / 155 ST+ IND 08 156 CLX 157 SIGN 158 ST- 06 159 ST- 07 160 DSE 08 161 GTO 05 162 .75 163 STO 14 164 CLX 165 STO 15 166 CF 10 167 XEQ 10 168 LBL 06 169 RCL IND 07 170 RCL IND 06 171 6 172 / 173 - 174 RCL 02 175 * 176 3 177 / 178 ST+ IND 08 179 CLX 180 SIGN 181 ST- 06 182 ST- 07 183 DSE 08 184 GTO 06 185 RCL 02 186 ST+ 17 187 VIEW 17 188 DSE 16 189 GTO 01 190 RCL 17 191 END

( 295 bytes / SIZE 18+10.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:    Once again the same one:

 Register input h = 1 ; N=10 h = 0.5 ; N=20 m R18=m1 R19=m2 R20=m3 2    1    3 undisturbed undisturbed p  o  s. R21=x1 R22=y1 R23=z1 ------- R24=x2 R25=y2 R26=z2 -------- R27=x3 R28=y3 R29=z3 2    0    0 -----    0    4    0 -----    0    0    1 1.992077583 0.300333564 0.003673681  -------------- 0.000661670 3.996080572 0.100603411  -------------- -.194938947 0.001084100 0.997349742 1.992077587     0.300333552     0.003673676      --------------     0.000661670     3.996080575     0.100603412     --------------    -0.194938947     0.001084108     0.997349745 v  e  l  . R30=x'1 R31=y'1 R32=z'1 -------- R33=x'2 R34=y'2 R35=z'2 -------- R36=x'3 R37=y'3 R38=z'3 0  0.03    0 -----    0    0  0.01 ----- -0.02    0    0 -0.001550083 0.030038158 0.000706684  -------------- 0.000132598 -0.000790385 0.010117549  -------------- -0.019010811  0.000238023 -0.000510306 -0.001550083     0.030038158     0.000706684     --------------     0.000132598    -0.000790385     0.010117549     --------------    -0.019010811     0.000238023    -0.000510306

t0 =  0  STO 17

n =   3  STO 01
h =   1  STO 02
N = 10 STO 03   XEQ "GM3"  >>>  X = 10

R21 thru R38 contain the next to last column above

-With h = 0.5  STO 02  & N = 20  STO 03 , the HP41 returns the last column.

2°)  Newtonian Gravitation - Heliocentric Coordinates

a) 4th-order Runge-Kutta Formula

-Like "GM4" , "PLN4" employs a 4th-order Runge-Kutta formula to calculate the positions of n planets in a Solar system
-Since the differential equations are still of the type y" = f(x,y) , we can use the same 3-stage method:

y(t+h) = y(t) + h ( y'(t) + k1/6 + 2k2/6 )
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)

-These equations are ( with x , y , z = heliocentric coordinates in  AU of the planets and  x' , y' , z' = the speeds of the planets in AU/day ):

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

G = gravitational constant = 0.017202098952  = 1 / 3379.380681

m0 = mass of the Sun
mi = mass of the planet i

Data Registers:              R00 = G.h           ( Registers R01-R02-R03 & R18 thru R19+7n are to be initialized before executing "PLN4" )

•  R01 = n = nb of planets                      •  R18 = t0 = initial instant                                                   R04 to R17: temp
•  R02 = h = step size ( in days )            •  R19 = m0 = mass of the Sun ( 1 for the Solar system )
•  R03 = N = nb of steps                       •  R20 = m1 = mass of the 1st planet
................................................

•  R19+n = mn = mass of the n-th planet

•  R20+n = x1              •  R20+4n = x'1
•  R21+n = y1              •  R21+4n = y'1                       R20+7n to R19+10n = k2 , k3
•  R22+n = z1              •  R22+4n = z1                       R20+10n to R19+13n = k1 , k1 - k2
.................                   ..................

•  R17+4n = xn            •  R17+7n = x'n
•  R18+4n = yn            •  R18+7n = y'n
•  R19+4n = zn            •  R19+7n = z'n
Flags: /
Subroutines: /

-Lines 29-187-275 are three-byte GTOs

 01 LBL "PLN4"  02 RCL 03  03 STO 17  04 RCL 02  05 3379.380681   06 /  07 STO 00  08 19.019  09 RCL 01  10 +  11 STO 05  12 LASTX  13 .1  14 %  15 X<>Y  16 3  17 *  18 STO T  19 +  20 +  21 STO 04  22 INT  23 +  24 STO 11  25 +  26 STO 12  27 +  28 STO 13  29 GTO 01  30 LBL 10  31 RCL 04  32 STO 06  33 STO 09  34 CLX  35 STO 14  36 STO 15  37 STO 16  38 RCL 05  39 STO 07  40 LBL 13 41 RCL IND 09  42 ENTER  43 X^2  44 DSE 09  45 RCL IND 09  46 STO T  47 X^2  48 +  49 DSE 09  50 RCL IND 09  51 STO 08  52 X^2  53 +  54 ENTER  55 SQRT  56 *  57 RCL IND 07  58 X<>Y  59 /  60 ST* 08  61 ST* Z  62 *  63 ST- 16  64 X<>Y  65 ST- 15  66 RCL 08  67 ST- 14  68 CLX  69 SIGN  70 ST- 07  71 DSE 09  72 GTO 13  73 RCL 01  74 STO 07  75 RCL 10  76 STO 08  77 LBL 12  78 RCL 16  79 RCL 00  80 * 81 STO IND 08  82 DSE 08  83 RCL 15  84 LASTX  85 *  86 STO IND 08  87 DSE 08  88 RCL 14  89 LASTX  90 *  91 STO IND 08  92 CLX  93 SIGN  94 ST- 08  95 DSE 07  96 GTO 12  97 LBL 03  98 RCL 04  99 STO 08 100 RCL 05 101 STO 07 102 LBL 04 103 RCL IND 08 104 RCL IND 06 105 - 106 ENTER 107 X^2 108 DSE 06 109 DSE 08 110 RCL IND 08 111 RCL IND 06 112 - 113 STO 09 114 X^2 115 + 116 DSE 06 117 DSE 08 118 RCL IND 08 119 RCL IND 06 120 - 121 STO T 122 X^2 123 + 124 ENTER 125 SQRT 126 * 127 X=0? 128 SIGN 129 RCL 00 130 X<>Y 131 / 132 RCL IND 07 133 * 134 ST* 09 135 ST* T 136 * 137 ST+ IND 10 138 DSE 10 139 RCL 09 140 ST+ IND 10 141 DSE 10 142 R^ 143 ST+ IND 10 144 CLX 145 SIGN 146 ST- 07 147 ST+ X 148 ST+ 06 149 ST+ 10 150 DSE 08 151 GTO 04 152 RCL IND 06 153 ENTER 154 X^2 155 DSE 06 156 RCL IND 06 157 STO T 158 X^2 159 + 160 DSE 06 161 RCL IND 06 162 STO 09 163 X^2 164 + 165 ENTER 166 SQRT 167 * 168 RCL 00 169 X<>Y 170 / 171 RCL 19 172 * 173 ST* 09 174 ST* Z 175 * 176 ST- IND 10 177 DSE 10 178 X<>Y 179 ST- IND 10 180 DSE 10 181 RCL 09 182 ST- IND 10 183 CLX 184 SIGN 185 ST- 10 186 DSE 06 187 GTO 03 188 RCL 13 189 STO 06 190 RCL 11 191 STO 07 192 RCL 04 193 STO 08 194 RTN 195 LBL 01 196 RCL 13 197 STO 10 198 XEQ 10 199 LBL 05 200 RCL IND 06 201 4 202 / 203 RCL IND 07 204 + 205 2 206 / 207 RCL 02 208 * 209 ST+ IND 08 210 CLX 211 SIGN 212 ST- 06 213 ST- 07 214 DSE 08 215 GTO 05 216 XEQ 10 217 RCL 12 218 STO 09 219 LBL 07 220 RCL IND 09 221 RCL IND 06 222 4 223 / 224 - 225 RCL IND 07 226 + 227 RCL 02 228 * 229 2 230 / 231 ST+ IND 08 232 RCL IND 06 233 RCL IND 09 234 ST- IND 06 235 4 236 * 237 + 238 6 239 / 240 ST+ IND 07 241 CLX 242 SIGN 243 ST- 06 244 ST- 07 245 ST- 09 246 DSE 08 247 GTO 07 248 RCL 12 249 STO 10 250 XEQ 10 251 RCL 12 252 STO 09 253 LBL 08 254 RCL IND 06 255 RCL 02 256 * 257 6 258 / 259 ST+ IND 08 260 RCL IND 09 261 LASTX 262 / 263 ST+ IND 07 264 CLX 265 SIGN 266 ST- 06 267 ST- 07 268 ST- 09 269 DSE 08 270 GTO 08 271 RCL 02 272 ST+ 18 273 VIEW 18 274 DSE 17 275 GTO 01 276 RCL 18 277 END

( 415 bytes / SIZE 20+13.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:   Given the masses & initial positions & velocities in the array below, calculate the positions & velocities 10 days later

n =  2   STO 01     and if we choose stepsize = h  = 10 days & N = 1 step:
h = 10  STO 02
N =  1   STO 03

-Initialize registers R19 thru R33:   3  STO 19   2 STO 20  .........  0.01  STO 33
-If we take initial instant = 0          0  STO 18

 Registers input h = 10 , N=1 h = 5 , N = 2 m R19=m0 R20=m1 R21=m2 3    2    1 unchanged unchanged p  o  s. R22=x1 R23=y1 R24=z1 -------- R25=x2 R26=y2 R27=z2 2    0   -1 ----    0    4   -1 2.187016538  0.299249966 -0.993675929  --------------  0.195600614  3.994996700 -0.896746283 2.187016532   0.299249475  -0.993676059   --------------   0.195600615   3.994996481  -0.896746330 v  e  l. R28=x'1 R29=y'1 R30=z'1 -------- R31=x'2 R32=y'2 R33=z'2 0.02 0.03   0 ---- 0.02   0 0.01 0.017460717  0.029800137  0.001216996  --------------  0.019143404 -0.001028406  0.010627856 0.017460727    0.029800135    0.001216990    -------------    0.019143408   -0.001028407    0.010627854

XEQ "PLN4"  >>>>  X = 10 = R18                                           ---Execution time = 63s---

-And we have in registers R22 to R33 the next to last column of the array above.
-With h = 5 days & N = 2 , the HP41 returns the last comlum above

Notes:

-To continue with the same h and N , simply press  R/S
-With n = 9 & N = 1 , execution time = 10m59s

-The maximum n-value is 23 ( so 24-body problem ) -> SIZE 319

b) 6th-order Runge-Kutta Formula

-"PLN6" employs the same formula as "GM6" ( Albrecht method )

Data Registers:              R00 = G.h           ( Registers R01-R02-R03 & R20 thru R21+7n are to be initialized before executing "PLN6" )

•  R01 = n = nb of planets                      •  R20 = t0 = initial instant                                                   R04 to R19: temp
•  R02 = h = step size ( in days )            •  R21 = m0 = mass of the Sun ( 1 for the Solar system )
•  R03 = N = nb of steps                       •  R22 = m1 = mass of the 1st planet
................................................

•  R21+n = mn = mass of the n-th planet

•  R22+n = x1              •  R22+4n = x'1
•  R23+n = y1              •  R23+4n = y'1                       R22+7n to R21+10n = k4
•  R24+n = z1              •  R24+4n = z1                       R22+10n to R21+13n = k3
.................                   ..................                         R22+13n to R21+16n = k2 , 7.k1/90-17k2/105

•  R19+4n = xn            •  R19+7n = x'n                      R22+17n to R21+19n = k1 , k5
•  R20+4n = yn            •  R20+7n = y'n
•  R21+4n = zn            •  R21+7n = z'n
Flags: /
Subroutines: /

-Lines 33-191-395 are three-byte GTOs

 01 LBL "PLN6"  02 RCL 03  03 STO 19  04 RCL 02  05 3379.380681   06 /  07 STO 00  08 21.021  09 RCL 01  10 +  11 STO 05  12 LASTX  13 .1  14 %  15 X<>Y  16 3  17 *  18 STO T  19 +  20 +  21 STO 04  22 INT  23 +  24 STO 11  25 +  26 STO 15  27 +  28 STO 14  29 +  30 STO 13  31 +  32 STO 12  33 GTO 01  34 LBL 10  35 RCL 04  36 STO 06  37 STO 09  38 CLX  39 STO 16  40 STO 17  41 STO 18  42 RCL 05  43 STO 07  44 LBL 13  45 RCL IND 09  46 ENTER  47 X^2  48 DSE 09  49 RCL IND 09  50 STO T  51 X^2  52 +  53 DSE 09  54 RCL IND 09  55 STO 08  56 X^2  57 + 58 ENTER  59 SQRT  60 *  61 RCL IND 07   62 X<>Y  63 /  64 ST* 08  65 ST* Z  66 *  67 ST- 18  68 X<>Y  69 ST- 17  70 RCL 08  71 ST- 16  72 CLX  73 SIGN  74 ST- 07  75 DSE 09  76 GTO 13  77 RCL 01  78 STO 07  79 RCL 10  80 STO 08  81 LBL 14  82 RCL 18  83 RCL 00  84 *  85 STO IND 08  86 DSE 08  87 RCL 17  88 LASTX  89 *  90 STO IND 08  91 DSE 08  92 RCL 16  93 LASTX  94 *  95 STO IND 08  96 CLX  97 SIGN  98 ST- 08  99 DSE 07 100 GTO 14 101 LBL 03 102 RCL 04 103 STO 08 104 RCL 05 105 STO 07 106 LBL 04 107 RCL IND 08 108 RCL IND 06 109 - 110 ENTER 111 X^2 112 DSE 06 113 DSE 08 114 RCL IND 08 115 RCL IND 06 116 - 117 STO 09 118 X^2 119 + 120 DSE 06 121 DSE 08 122 RCL IND 08 123 RCL IND 06 124 - 125 STO T 126 X^2 127 + 128 ENTER 129 SQRT 130 * 131 X=0? 132 SIGN 133 RCL 00 134 X<>Y 135 / 136 RCL IND 07 137 * 138 ST* 09 139 ST* T 140 * 141 ST+ IND 10 142 DSE 10 143 RCL 09 144 ST+ IND 10 145 DSE 10 146 R^ 147 ST+ IND 10 148 CLX 149 SIGN 150 ST- 07 151 ST+ X 152 ST+ 06 153 ST+ 10 154 DSE 08 155 GTO 04 156 RCL IND 06 157 ENTER 158 X^2 159 DSE 06 160 RCL IND 06 161 STO T 162 X^2 163 + 164 DSE 06 165 RCL IND 06 166 STO 09 167 X^2 168 + 169 ENTER 170 SQRT 171 * 172 RCL 00 173 X<>Y 174 / 175 RCL 21 176 * 177 ST* 09 178 ST* Z 179 * 180 ST- IND 10 181 DSE 10 182 X<>Y 183 ST- IND 10 184 DSE 10 185 RCL 09 186 ST- IND 10 187 CLX 188 SIGN 189 ST- 10 190 DSE 06 191 GTO 03 192 RCL 12 193 STO 06 194 RCL 11 195 STO 07 196 RCL 04 197 STO 08 198 RTN 199 LBL 01 200 RCL 12 201 STO 10 202 XEQ 10 203 LBL 05 204 RCL IND 06 205 8 206 / 207 RCL IND 07 208 + 209 RCL 02 210 * 211 4 212 / 213 ST+ IND 08 214 CLX 215 SIGN 216 ST- 06 217 ST- 07 218 DSE 08 219 GTO 05 220 XEQ 10 221 RCL 13 222 STO 09 223 LBL 06 224 RCL IND 09 225 6 226 / 227 RCL IND 06 228 7 229 * 230 96 231 / 232 - 233 RCL IND 07 234 4 235 / 236 + 237 RCL 02 238 * 239 ST+ IND 08 240 CLX 241 SIGN 242 ST- 06 243 ST- 07 244 ST- 09 245 DSE 08 246 GTO 06 247 XEQ 10 248 RCL 13 249 STO 09 250 RCL 14 251 STO 16 252 LBL 07 253 RCL IND 06 254 13 255 * 256 24 257 / 258 RCL IND 09 259 6 260 / 261 - 262 RCL IND 16 263 4 264 / 265 + 266 RCL IND 07 267 + 268 RCL 02 269 * 270 4 271 / 272 ST+ IND 08 273 CLX 274 SIGN 275 ST- 06 276 ST- 07 277 ST- 09 278 ST- 16 279 DSE 08 280 GTO 07 281 XEQ 10 282 RCL 13 283 STO 09 284 RCL 14 285 STO 16 286 RCL 15 287 STO 17 288 LBL 08 289 RCL IND 09 290 17 291 * 292 56 293 / 294 RCL IND 06 295 3 296 * 297 32 298 / 299 - 300 RCL IND 16 301 15 302 * 303 112 304 / 305 - 306 RCL IND 17 307 7 308 / 309 + 310 RCL IND 07 311 4 312 / 313 + 314 RCL 02 315 * 316 ST+ IND 08 317 RCL IND 06 318 7 319 * 320 90 321 / 322 ST+ IND 07 323 RCL IND 09 324 17 325 * 326 105 327 / 328 - 329 X<> IND 09 330 16 331 * 332 45 333 / 334 ST+ IND 07 335 CLX 336 SIGN 337 ST- 06 338 ST- 07 339 ST- 09 340 ST- 16 341 ST- 17 342 DSE 08 343 GTO 08 344 RCL 12 345 STO 10 346 XEQ 10 347 RCL 13 348 STO 09 349 RCL 14 350 STO 10 351 RCL 15 352 STO 16 353 LBL 12 354 RCL IND 10 355 87 356 * 357 RCL IND 16 358 34 359 * 360 - 361 630 362 / 363 RCL IND 09 364 + 365 RCL 02 366 * 367 ST+ IND 08 368 RCL IND 10 369 12 370 * 371 RCL IND 16 372 32 373 * 374 + 375 RCL IND 06 376 7 377 * 378 + 379 90 380 / 381 ST+ IND 07 382 CLX 383 SIGN 384 ST- 06 385 ST- 07 386 ST- 09 387 ST- 10 388 ST- 16 389 DSE 08 390 GTO 12 391 RCL 02 392 ST+ 20 393 VIEW 20 394 DSE 19 395 GTO 01 396 RCL 20 397 END

( 600 bytes / SIZE 22+19.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:   Given the masses & initial positions & velocities in the array below, calculate the positions & velocities 10 days later

n =  2   STO 01     and if we choose stepsize = h  = 10 days & N = 1 step:
h = 10  STO 02
N =  1   STO 03

-Initialize registers R22 thru R36:   3  STO 21  2 STO 22  .........  0.01  STO 35
-If we take initial instant = 0          0  STO 20

 Registers input h = 10 , N=1 h = 5 , N = 2 m R21=m0 R22=m1 R23=m2 3    2    1 unchanged unchanged p  o  s. R24=x1 R25=y1 R26=z1 -------- R27=x2 R28=y2 R29=z2 2    0   -1 ----    0    4   -1 2.187016535  0.299249441 -0.993676070  --------------  0.195600617  3.994996465 -0.896746334 2.187016534   0.299249441  -0.993676070   --------------   0.195600617   3.994996464  -0.896746334 v  e  l. R30=x'1 R31=y'1 R32=z'1 -------- R33=x'2 R34=y'2 R35=z'2 0.02 0.03   0 ---- 0.02   0 0.01 0.017460728  0.029800135  0.001216990  --------------  0.019143409 -0.001028408  0.010627854 0.017460728    0.029800135    0.001216990    -------------    0.019143409   -0.001028408    0.010627854

XEQ "PLN6"  >>>>  X = 10 = R20                                           ---Execution time = 122s---

-And we have in registers R24 to R35 the next to last column of the array above.
-With h = 5 days & N = 2 , the HP41 returns the last comlum above.
-Since the results are - almost - equal, they are -almost - exact ( rounded to 9 decimals ).

Notes:

-To continue with the same h and N , simply press  R/S
-With n = 9 & N = 1 , execution time = 19m49s

-"PLN6" can solve a Solar system with 15 planets.

3°)  General Relativity - Solar Systems

a) 4th-order Runge-Kutta Formula

-The equations of the Einsteinian gravitation are much more complicated, but for a Solar system,
the masses of the planets are usually very small, compared to the mass of the Sun.
-So, we can use a simpler way to calculate the relativistic perturbations, given in reference [2]

dri'' = ( GM/ri3 ) [ 4 ( GM/ri ) ri - vi2ri + 4 ( ri.vi ) vi ] / c2              c = speed of light = 173.446327  AU/day

where  ri =  ( xi , yi , zi )  = position-vector of the planet i     ri = distance Sun-planet i      M = mass of the Sun = m0 = 1 for our Solar system
vi = ( x'i , y'i , z'i ) = velocity-vector of the planet i     vi = norm of vi

-Since the system of differential equations explicitly depends on the speeds of the planets, we cannot use the same Runge-Kutta formula.
-"PLN4G" employs a similar - but 4-stage - Runge-Kutta formula of order 4:

k1 = h.f (y,y')
k2 = h.f( y+h.y'/2+h.k1/8 , y'+k1/2 )
k3 = h.f( y+h.y'/2+h.k1/8 , y'+k2/2)
k4 = h.f( y+h.y'+h.k3/2 , y'+k3)

y(t+h) = y(t) + h [ y'(t) + ( k1 + k2 + k3 )/6 ]
y'(t+h) = y'(t) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6

Data Registers:              R00 = G           ( Registers R01-R02-R03 & R29 thru R30+7n are to be initialized before executing "PL4G" )

•  R01 = n = nb of planets                      •  R29 = t0 = initial instant                  R04 to R27: temp            R28: unused
•  R02 = h = step size ( in days )            •  R30 = m0 = mass of the Sun ( 1 for the Solar system )
•  R03 = N = nb of steps                       •  R31 = m1 = mass of the 1st planet
................................................

•  R30+n = mn = mass of the n-th planet

•  R31+n = x1              •  R31+4n = x'1                      R31+7n  to  R30+10n = k5 and linear combinations of other ki
•  R32+n = y1              •  R32+4n = y'1                      R31+10n to R30+13n = k2 ,  k2 - 2.k3
•  R33+n = z1              •  R33+4n = z1                       R31+13n to R30+16n = k1
.................                   ..................

•  R28+4n = xn            •  R28+7n = x'n
•  R29+4n = yn            •  R29+7n = y'n
•  R30+4n = zn            •  R30+7n = z'n
Flags: /
Subroutines: /

-Lines 48-268-392 are three-byte GTOs

 01 LBL "PL4G"  02 RCL 03  03 STO 16  04 3379.380681   05 1/X  06 STO 00  07 RCL 02  08 *  09 STO 27  10 RCL 30  11 *  12 STO 15  13 29979.06382  14 STO 23  15 30.03  16 RCL 01  17 +  18 STO 05  19 LASTX  20 .1  21 %  22 X<>Y  23 3  24 *  25 STO T  26 +  27 +  28 STO 04  29 INT  30 .1  31 %  32 +  33 +  34 STO 11  35 INT  36 +  37 STO 12  38 +  39 STO 13  40 +  41 STO 14  42 2  43 STO 25  44 4  45 STO 26  46 +  47 STO 24  48 GTO 01  49 LBL 10  50 RCL 04  51 STO 06  52 STO 09  53 RCL 11  54 STO 20  55 CLX  56 STO 17 57 STO 18  58 STO 19  59 RCL 05  60 STO 07  61 LBL 13  62 RCL IND 09  63 ENTER  64 X^2  65 DSE 09  66 RCL IND 09  67 STO T  68 X^2  69 +  70 DSE 09  71 RCL IND 09  72 STO 08  73 X^2  74 +  75 ENTER  76 SQRT  77 *  78 RCL IND 07  79 X<>Y  80 /  81 ST* 08  82 ST* Z  83 *  84 ST- 19  85 X<>Y  86 ST- 18  87 RCL 08  88 ST- 17  89 CLX  90 SIGN  91 ST- 07  92 DSE 09  93 GTO 13  94 RCL 01  95 STO 07  96 RCL 10  97 STO 08  98 LBL 12  99 RCL 19 100 RCL 27 101 * 102 STO IND 08 103 DSE 08 104 RCL 18 105 LASTX 106 * 107 STO IND 08 108 DSE 08 109 RCL 17 110 LASTX 111 * 112 STO IND 08 113 CLX 114 SIGN 115 ST- 08 116 DSE 07 117 GTO 12 118 LBL 03 119 RCL IND 06 120 STO 09 121 X^2 122 DSE 06 123 RCL IND 06 124 STO 08 125 X^2 126 + 127 DSE 06 128 RCL IND 06 129 STO 07 130 X^2 131 + 132 ENTER 133 SQRT 134 STO 21 135 * 136 RCL 15 137 X<>Y 138 / 139 STO 22 140 RCL IND 20 141 STO 19 142 X^2 143 DSE 20 144 RCL IND 20 145 STO 18 146 X^2 147 + 148 DSE 20 149 RCL IND 20 150 STO 17 151 X^2 152 + 153 RCL 00 154 RCL 26 155 * 156 RCL 21 157 / 158 - 159 RCL 07 160 RCL 17 161 * 162 RCL 08 163 RCL 18 164 * 165 + 166 RCL 09 167 RCL 19 168 * 169 + 170 RCL 26 171 * 172 RCL 23 173 ST/ Z 174 / 175 ST* 17 176 ST* 18 177 ST* 19 178 CLX 179 SIGN 180 ST- 20 181 + 182 ST* 07 183 ST* 08 184 RCL 09 185 * 186 RCL 19 187 - 188 RCL 22 189 STO 09 190 * 191 ST- IND 10 192 DSE 10 193 RCL 08 194 RCL 18 195 - 196 RCL 09 197 * 198 ST- IND 10 199 DSE 10 200 RCL 07 201 RCL 17 202 - 203 RCL 09 204 * 205 ST- IND 10 206 RCL 25 207 ST+ 06 208 ST+ 10 209 RCL 04 210 STO 08 211 RCL 05 212 STO 07 213 LBL 04 214 RCL IND 08 215 RCL IND 06 216 - 217 ENTER 218 X^2 219 DSE 06 220 DSE 08 221 RCL IND 08 222 RCL IND 06 223 - 224 STO 09 225 X^2 226 + 227 DSE 06 228 DSE 08 229 RCL IND 08 230 RCL IND 06 231 - 232 STO T 233 X^2 234 + 235 ENTER 236 SQRT 237 * 238 X=0? 239 SIGN 240 RCL 27 241 X<>Y 242 / 243 RCL IND 07 244 * 245 ST* 09 246 ST* T 247 * 248 ST+ IND 10 249 DSE 10 250 RCL 09 251 ST+ IND 10 252 DSE 10 253 R^ 254 ST+ IND 10 255 CLX 256 SIGN 257 ST- 07 258 RCL 25 259 ST+ 06 260 ST+ 10 261 DSE 08 262 GTO 04 263 ST- 06 264 PI 265 INT 266 ST- 10 267 DSE 06 268 GTO 03 269 RCL 14 270 STO 06 271 RCL 11 272 STO 07 273 RCL 04 274 STO 08 275 RTN 276 LBL 01 277 RCL 14 278 STO 10 279 XEQ 10 280 LBL 05 281 RCL IND 07 282 RCL IND 06 283 RCL 26 284 / 285 + 286 RCL 02 287 * 288 RCL 25 289 / 290 ST+ IND 08 291 RCL IND 06 292 LASTX 293 / 294 ST+ IND 07 295 CLX 296 SIGN 297 ST- 06 298 ST- 07 299 DSE 08 300 GTO 05 301 XEQ 10 302 RCL 13 303 STO 09 304 LBL 07 305 RCL IND 09 306 RCL IND 06 307 - 308 RCL 25 309 / 310 ST+ IND 07 311 CLX 312 SIGN 313 ST- 06 314 ST- 09 315 DSE 07 316 GTO 07 317 XEQ 10 318 RCL 13 319 STO 09 320 RCL 12 321 STO 20 322 LBL 08 323 RCL IND 20 324 RCL IND 09 325 RCL 25 326 / 327 - 328 RCL IND 06 329 RCL 26 330 / 331 - 332 RCL IND 07 333 + 334 RCL 02 335 * 336 RCL 25 337 / 338 ST+ IND 08 339 RCL IND 09 340 LASTX 341 / 342 RCL IND 20 343 ST- IND 09 344 ST- IND 09 345 - 346 ST- IND 07 347 CLX 348 SIGN 349 ST- 06 350 ST- 07 351 ST- 09 352 ST- 20 353 DSE 08 354 GTO 08 355 RCL 12 356 STO 10 357 XEQ 10 358 RCL 13 359 STO 09 360 RCL 12 361 STO 10 362 LBL 09 363 RCL IND 06 364 RCL IND 09 365 + 366 RCL 02 367 * 368 RCL 24 369 / 370 ST+ IND 08 371 RCL IND 06 372 RCL IND 09 373 ST+ X 374 + 375 RCL IND 10 376 + 377 RCL 24 378 / 379 ST+ IND 07 380 CLX 381 SIGN 382 ST- 06 383 ST- 07 384 ST- 09 385 ST- 10 386 DSE 08 387 GTO 09 388 RCL 02 389 ST+ 29 390 VIEW 29 391 DSE 16 392 GTO 01 393 RCL 29 394 END

( 611 bytes / SIZE 31+16.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:   Integrate the Solar system with h = 0.1 day , starting with the initial values in the 1st column below 2019/11/13 0h TDB

9   STO 01
0.1  STO 02
320  STO 03   to obtain the positions 32 days later ( 2019/12/15 0h TDB )

-Store the initial values of positions & velocities into R40 thru R93 ( use for instance the following routine )

 01 LBL "INIPLN"  02 6023682.156  03 1/X  04 STO 31  05 408523.7187  06 1/X  07 STO 32  08 328900.5598   09 1/X  10 STO 33  11 3098703.590  12 1/X  13 STO 34  14 1047.348625  15 1/X  16 STO 35  17 3497.901768  18 1/X  19 STO 36  20 22902.98161          21 1/X  22 STO 37  23 19412.25978   24 1/X  25 STO 38  26 135836683.8   27 1/X  28 STO 39 29 169.3287041 E-3  30 STO 40  31 236.7984713 E-3  32 STO 41  33 108.9436218 E-3  34 STO 42  35 206.8552787 E-3  36 STO 43  37 -631.2899640 E-3  38 STO 44  39 -297.1381466 E-3  40 STO 45  41 635.7893837 E-3  42 STO 46  43 695.9559735 E-3  44 STO 47  45 301.6942285 E-3  46 STO 48  47 -1604.907774 E-3  48 STO 49  49 -301.5483654 E-3  50 STO 50  51 -94.99990761 E-3  52 STO 51  53 161.3509664 E-3  54 STO 52  55 -4817.408291 E-3  56 STO 53 57 -2068.801992 E-3  58 STO 54  59 3557.619174 E-3  60 STO 55  61 -8621.479227 E-3  62 STO 56  63 -3714.330343 E-3  64 STO 57  65 16336.48164 E-3  66 STO 58  67 10370.59623 E-3  68 STO 59  69 4311.051187 E-3  70 STO 60  71 29210.34115 E-3  72 STO 61  73 -5766.105481 E-3  74 STO 62  75 -3087.429294 E-3  76 STO 63  77 12830.55991 E-3  78 STO 64  79 -28665.41178 E-3  80 STO 65  81 -12811.23541 E-3  82 STO 66  83 -29.19594877 E-3  84 STO 67 85 13.50298685 E-3  86 STO 68  87 10.23959860 E-3  88 STO 69  89 19.25542713 E-3  90 STO 70  91 5.625334541 E-3  92 STO 71  93 1.312762462 E-3  94 STO 72  95 -13.46613287 E-3  96 STO 73  97 10.08021408 E-3  98 STO 74  99 4.369793699 E-3 100 STO 75 101 3.212151368 E-3 102 STO 76 103 -11.36579203 E-3 104 STO 77 105 -5.299891568 E-3 106 STO 78 107 7.459099614 E-3 108 STO 79 109 .6076288015 E-3 110 STO 80 111 .7886915195 E-4 112 STO 81 113 4.916558968 E-3 114 STO 82 115 1.891205960 E-3 116 STO 83 117 .5695514479 E-3 118 STO 84 119 -2.249753958 E-3 120 STO 85 121 2.790004087 E-3 122 STO 86 123 1.253624408 E-3 124 STO 87 125 .6706339629 E-3 126 STO 88 127 2.860964834 E-3 128 STO 89 129 1.154306180 E-3 130 STO 90 131 2.985297449 E-3 132 STO 91 133 .8432457347 E-3 134 STO 92 135 -.6338481038 E-3 136 STO 93 137 1 138 STO 30 139 END

-Then   XEQ "PL4G"    >>>>    X = 32                        ---Execution time =8m39s---        with V41 in turbo mode

-Execution time is about 2 seconds with free42binary and about 86h02mn with a real HP41.

-Here are the results:

 2019/11/13 0h        Initial Values 2019/12/15 0h             DE436 2019/12/15 0h         HP41  PLN4 2019/12/15 0h         HP41  PL4G 2019/12/15 0h       free42  PL4G M 169.3287041 E-3 E  236.7984713 E-3 R  108.9436218 E-3 -0.362228680      -0.226192417      -0.083282304 -0.362228535      -0.226192480      -0.083282351 -0.362228673      -0.226192421      -0.083282306 -0.362228680      -0.226192417      -0.083282304 V 206.8552787 E-3 E -631.2899640 E-3 N -297.1381466 E-3 0.669481274      -0.241877222      -0.151193495 0.669481269      -0.241877213      -0.151193491 0.669481278      -0.241877219      -0.151193496 0.669481274      -0.241877222      -0.151193495 E  635.7893837 E-3 M 695.9559735 E-3 B  301.6942285 E-3 0.129773202       0.895203527       0.388069088 0.129773212       0.895203528       0.388069090 0.129773204       0.895203533       0.388069090 0.129773209       0.895203532       0.388069091 M-1604.907774E-3 A -301.5483654E-3 R -94.99990761E-3 -1.447205476      -0.650423094      -0.259275988 -1.447205487      -0.650423104      -0.259275978 -1.447205486      -0.650423103      -0.259275978 -1.447205476      -0.650423094      -0.259275988 J  161.3509664 E-3 U -4817.408291E-3 P -2068.801992 E-3 0.399788668      -4.792899165      -2.064100995 0.399788667      -4.792899161      -2.064100993 0.399788667      -4.792899161      -2.064100993 0.399788668      -4.792899165      -2.064100995 S  3557.619174E-3 A -8621.479227E-3 T -3714.330343E-3 3.714404180      -8.559663453      -3.695545437 3.714404185      -8.559663447      -3.695545434 3.714404184      -8.559663445      -3.695545434 3.714404180      -8.559663453      -3.695545436 U  16336.48164 E-3 R  10370.59623 E-3 A  4311.051187 E-3 16.26417108       10.45967901       4.351085196 16.26417110       10.45967902       4.351085197 16.26417110       10.45967902       4.351085198 16.26417109       10.45967902       4.351085195 N 29210.34115 E-3 E -5766.105481 E-3 P -3087.429294 E-3 29.23163521      -5.674517011      -3.050471876 29.23163520      -5.674517009      -3.050471876 29.23163520      -5.674517008      -3.050471876 29.23163521      -5.674517011      -3.050471876 P 12830.55991 E-3 L-28665.41178 E-3 U-12811.23541 E-3 12.92603848      -28.63831132      -12.83146648 12.92603849      -28.63831131      -12.83146645 12.92603849      -28.63831131      -12.83146645 12.92603848      -28.63831132      -12.83146648

-As you can see, the results are almost perfect with free42, except for the Earth-Moon barycenter
-The small errors given by the HP41 are due to roundoff-errors because the HP41 works with 10 digits.

-The differences between PLN4 & PLNG are the most important for Mercury:
-Logical since the main effect of relativistic perturbation ( advance of the perihelion ) occurs for Mercury.

Notes:

-To continue with the same h and N , simply press  R/S
-With n = 9 & N = 1 , execution time = 16m08s

-Using the same program and the same initial values, we get the following positions 1000 days after 2019/11/13 0h TDB i-e 2022/09/08 0h TDB
( I've used free42binary )

 2022/08/09 0h             DE436 2022/08/09 0h         free42  PLN4 2022/08/09 0h        free42  PL4G M     -.358841135 E     -.232656389 R     -.087091439 -0.358838544      -0.232661808      -0.087094602 -0.358841136      -0.232656384      -0.087091436 V      .027401134  E      .656557789  N      .293690656 0.027399505       0.656557746       0.293690739 0.027401123       0.656557789       0.293690657 E      .730692795 M     -.644947005  B     -.279584371 0.730692351      -0.644947426      -0.279584592 0.730691632      -0.644948176      -0.279584917 M     1.387877753  A      .162635812  R      .037151340 1.387877685       0.162636252       0.037151544 1.387877749       0.162635820       0.037151344 J      4.956737203 U     -.042241507 P     -.138761018 4.956737191      -0.042241491      -0.138761011 4.956737215      -0.042241518      -0.138761024 S      7.713325685 A    -5.568694720 T    -2.632229528 7.713325688      -5.568694718      -2.632229527 7.713325690      -5.568694722      -2.632229529 U      13.786845002 R      12.949153815 A       5.476298203 13.78684501       12.949153818       5.476298202 13.786845011       12.949153818       5.476298202 N      29.71561523  E      -2.874349013  P      -1.916321904 29.715615232      -2.874349013      -1.916321905 29.715615232      -2.874349014      -1.916321905 P      15.761636163  L     -27.712190899  U     -13.394699107 15.761636169      -27.712190896      -13.394699104 15.761636169      -27.712190896      -13.394699104

Notes:

-For an interval of 1000 days, "PLN4" seems enough for Jupiter, Saturn, Uranus, Neptune & Pluto.
-The precision is rather good, except for the Earth-Moon barycenter
( the Moon influences the EMB motion, but it's not taken into account in these routines )
-With an HP41, there would be probably a little more roundoff-errors

-We can also use these routines with 2 planets Earth & Moon instead of EMB:

-Here is a routine to initialize R30 to R100 ( positions & speeds at 2019/11/13 0h TDB )
-The position of the Earth is registers R47-R48-R49 and its velocity in R77-R78-R79
-The position of the Moon is registers R68-R69-R70 and its velocity in R98-R99-R100.

 01 LBL "INIPL2"  02 1  03 STO 30  04 6023682.156  05 1/X  06 STO 31  07 408523.7187  08 1/X  09 STO 32  10 332946.0488          11 1/X  12 STO 33  13 3098703.590   14 1/X  15 STO 34  16 1047.348625  17 1/X  18 STO 35  19 3497.901768  20 1/X  21 STO 36  22 22902.98161   23 1/X  24 STO 37  25 19412.25978   26 1/X  27 STO 38  28 135836683.8  29 1/X  30 STO 39  31 27068703.24 32 1/X  33 STO 40  34 169.3287041 E-3  35 STO 41  36 236.7984713 E-3  37 STO 42  38 108.9436218 E-3  39 STO 43  40 206.8552787 E-3  41 STO 44  42 -631.2899640 E-3  43 STO 45  44 -297.1381466 E-3  45 STO 46  46 635.7711661 E-3  47 STO 47  48 695.9312501 E-3  49 STO 48  50 301.6856667 E-3  51 STO 49  52 -1604.907774 E-3  53 STO 50  54 -301.5483654 E-3  55 STO 51  56 -94.99990761 E-3  57 STO 52  58 161.3509664 E-3  59 STO 53  60 -4817.408291 E-3  61 STO 54  62 -2068.801992 E-3 63 STO 55  64 3557.619174 E-3  65 STO 56  66 -8621.479227 E-3  67 STO 57  68 -3714.330343 E-3  69 STO 58  70 16336.48164 E-3  71 STO 59  72 10370.59623 E-3  73 STO 60  74 4311.051187 E-3  75 STO 61  76 29210.34115 E-3  77 STO 62  78 -5766.105481 E-3  79 STO 63  80 -3087.429294 E-3  81 STO 64  82 12830.55991 E-3  83 STO 65  84 -28665.41178 E-3  85 STO 66  86 -12811.23541 E-3  87 STO 67  88 .6372704910  89 STO 68  90 .6979659403  91 STO 69  92 .3023903046  93 STO 70 94 -29.19594877 E-3  95 STO 71  96 13.50298685 E-3  97 STO 72  98 10.23959860 E-3  99 STO 73 100 19.25542713 E-3 101 STO 74 102 5.625334541 E-3 103 STO 75 104 1.312762462 E-3 105 STO 76 106 -13.46021537 E-3 107 STO 77 108 10.07688356 E-3 109 STO 78 110 4.367839045 E-3 111 STO 79 112 3.212151368 E-3 113 STO 80 114 -11.36579203 E-3 115 STO 81 116 -5.299891568 E-3 117 STO 82 118 7.459099614 E-3 119 STO 83 120 .6076288015 E-3 121 STO 84 122 .7886915195 E-4 123 STO 85 124 4.916558968 E-3 125 STO 86 126 1.891205960 E-3 127 STO 87 128 .5695514479 E-3 129 STO 88 130 -2.249753958 E-3 131 STO 89 132 2.790004087 E-3 133 STO 90 134 1.253624408 E-3 135 STO 91 136 .6706339629 E-3 137 STO 92 138 2.860964834 E-3 139 STO 93 140 1.154306180 E-3 141 STO 94 142 2.985297449 E-3 143 STO 95 144 .8432457347 E-3 145 STO 96 146 -.6338481038 E-3 147 STO 97 148 -13.94722891 E-3 149 STO 98 150 10.35098743 E-3 151 STO 99 152 100 153 4.528708175 E-3 154 STO IND Y 155 END

-Here are the results with stepsize  h = 0.1 day ( h = 0.05 day gives approximately the same accuracy )

XEQ "INIPL2"

0     STO 29
10    STO 01
0.1    STO 02
10000  STO 03   XEQ "PL4G"   >>>>   1000                                       ---Execution time = 57s---        with free42binary

( To use "PLN4" instead of "PL4G",

XEQ "INIPL2"    30.019071  XEQ "REGMOVE"

10     STO 01
0.1    STO 02
10000  STO 03   XEQ "PLN4"   )

 2022/08/09 0h             DE436 2022/08/09 0h         free42  PLN4 2022/08/09 0h        free42  PL4G M     -.358841135 E     -.232656389 R     -.087091439 -0.358838545      -0.232661806      -0.087094601 -0.358841137      -0.232656382      -0.087091435 V      .027401134  E      .656557789  N      .293690656 0.027399505       0.656557746       0.293690739 0.027401124       0.656557789       0.293690656 E      .730691272 A     -.644920826 R     -.279571034 0.730692018      -0.644920056      -0.279570700 0.730691299      -0.644920805      -0.279571025 M     1.387877753  A      .162635812  R      .037151340 1.387877685       0.162636252       0.037151544 1.387877749       0.162635820       0.037151344 J      4.956737203 U     -.042241507 P     -.138761018 4.956737191      -0.042241491      -0.138761011 4.956737215      -0.042241518      -0.138761024 S      7.713325685  A    -5.568694720 T    -2.632229528 7.713325688      -5.568694718      -2.632229527 7.713325690      -5.568694722      -2.632229529 U      13.786845002 R      12.949153815 A       5.476298203 13.786845010       12.949153818       5.476298202 13.786845011       12.949153818       5.476298202 N      29.71561523  E      -2.874349013  P      -1.916321904 29.715615232      -2.874349013      -1.916321905 29.715615232      -2.874349014      -1.916321905 P      15.761636163  L     -27.712190899  U     -13.394699107 15.761636169      -27.712190896      -13.394699104 15.761636169      -27.712190896      -13.394699104 M     0.730816620 O    -0.647075360 O    -0.280668649 0.730816733       -0.647074645       -0.280668329 0.730816017       -0.647075394       -0.280668654

-If we calculate the heliocentric coordinates of the Earth-Moon barycenter from the values above, it yields:

x =  0.730692814                                                x =  0.730692795
y = -0.644946984       whereas  DE436  gives     y = -0.644947005
z = -0.279584362                                               z = -0.279584371

-The errors were about 10^(-6)  AU  with EMB as a single planet.
-They are about 2 10^(-8)  AU  with the Earth & the Moon positions computed as 2 planets.

-If you execute PL4G from the last results, with h = - 0.1 day,
you get the initial values with an error of about 2 E-9 AU for Mercury
and even smaller errors for the other planets !

-Since "PL4G" can solve a Solar system with 18 planets - and "PLN4" with 23 planets - you can add several asteroids !
( of course much more with free42 )

b) 5th-order Runge-Kutta Formula

-The following 6-stage method is used:

k1 = h f( y , y' )
k2 = h f( y+h.( (2/5).y'+(2/25).k1 ), y'+(2/5) k1 )
k3 = h f( y+h.( (1/4).y'+(11/256).k1-(3/256).k2 ), y'+(11/64) k1+(5/64).k2 )
k4 = h f( y+h.( (1/2).y'+(1/8).k3 ), y'+(1/2) k3 )
k5 = h f( y+h.( (3/4).y'+(9/256).k1-(21/256).k2+(3/16).k3+(9/64).k4 ), y'+(3/64) k1-(15/64).k2+(3/8).k3+(9/16).k4 )
k6 = h f( y+h.( y'+(3/7).k2+(9/14).k3-(6/7).k4+(2/7).k5 ), y'+(5/7) k2+(6/7).k3-(12/7).k4+(8/7).k5 )

y(x+h) - y(x) = h.[ y' + ( 7 k1 + 24 k3 + 6 k4 + 8 k5 ) / 90 ]
y'(x+h) - y'(x) = ( 7 k1 + 32 k3 + 12 k4 + 32 k5 + 7 k6 ) / 90

Data Registers:              R00 = G            ( Registers R01-R02-R03 & R29 thru R30+7n are to be initialized before executing "PL5G" )

•  R01 = n = nb of planets                      •  R29 = t0 = initial instant                                                     R04 to R28: temp
•  R02 = h = step size ( in days )            •  R30 = m0 = mass of the Sun ( 1 for the Solar system )
•  R03 = N = nb of steps                       •  R31 = m1 = mass of the 1st planet
................................................

•  R30+n = mn = mass of the n-th planet

•  R31+n = x1              •  R31+4n = x'1                      R31+7n  to  R30+10n = k5 and linear combination of other ki
•  R32+n = y1              •  R32+4n = y'1                      R31+10n to R30+13n = k4 and linear combination of other ki
•  R33+n = z1              •  R33+4n = z1                       R31+13n to R30+16n = k3
.................                   ..................                        R31+16n to R30+19n = k2
R31+19n to R30+22n = k1 , k6
•  R28+4n = xn            •  R28+7n = x'n
•  R29+4n = yn            •  R29+7n = y'n
•  R30+4n = zn            •  R30+7n = z'n
Flags: /
Subroutines: /

-Lines 46-266-569-602 are three-byte GTOs

 01 LBL "PL5G"  02 RCL 03  03 STO 16  04 3379.380681   05 1/X  06 STO 00  07 RCL 02  08 *  09 STO 24  10 RCL 30  11 *  12 STO 15  13 29979.06382  14 STO 23  15 30.03  16 RCL 01  17 +  18 STO 05  19 LASTX  20 .1  21 %  22 X<>Y  23 3  24 *  25 STO T  26 +  27 +  28 STO 04  29 INT  30 +  31 STO 11  32 +  33 STO 27  34 +  35 STO 28  36 +  37 STO 12  38 +  39 STO 13  40 +  41 STO 14  42 2  43 STO 25  44 X^2  45 STO 26  46 GTO 01  47 LBL 10  48 RCL 04  49 STO 06  50 STO 09  51 RCL 11  52 STO 20  53 CLX  54 STO 17  55 STO 18  56 STO 19  57 RCL 05  58 STO 07  59 LBL 13  60 RCL IND 09  61 ENTER  62 X^2  63 DSE 09  64 RCL IND 09  65 STO T  66 X^2  67 +  68 DSE 09  69 RCL IND 09  70 STO 08  71 X^2  72 +  73 ENTER  74 SQRT  75 *  76 RCL IND 07  77 X<>Y  78 /  79 ST* 08  80 ST* Z  81 *  82 ST- 19  83 X<>Y  84 ST- 18  85 RCL 08  86 ST- 17 87 CLX  88 SIGN  89 ST- 07  90 DSE 09  91 GTO 13  92 RCL 01  93 STO 07  94 RCL 10  95 STO 08  96 LBL 12  97 RCL 19  98 RCL 24  99 * 100 STO IND 08 101 DSE 08 102 RCL 18 103 LASTX 104 * 105 STO IND 08 106 DSE 08 107 RCL 17 108 LASTX 109 * 110 STO IND 08 111 CLX 112 SIGN 113 ST- 08 114 DSE 07 115 GTO 12 116 LBL 03 117 RCL IND 06 118 STO 09 119 X^2 120 DSE 06 121 RCL IND 06 122 STO 08 123 X^2 124 + 125 DSE 06 126 RCL IND 06 127 STO 07 128 X^2 129 + 130 ENTER 131 SQRT 132 STO 21 133 * 134 RCL 15 135 X<>Y 136 / 137 STO 22 138 RCL IND 20 139 STO 19 140 X^2 141 DSE 20 142 RCL IND 20 143 STO 18 144 X^2 145 + 146 DSE 20 147 RCL IND 20 148 STO 17 149 X^2 150 + 151 RCL 00 152 RCL 26 153 * 154 RCL 21 155 / 156 - 157 RCL 07 158 RCL 17 159 * 160 RCL 08 161 RCL 18 162 * 163 + 164 RCL 09 165 RCL 19 166 * 167 + 168 RCL 26 169 * 170 RCL 23 171 ST/ Z 172 / 173 ST* 17 174 ST* 18 175 ST* 19 176 CLX 177 SIGN 178 ST- 20 179 + 180 ST* 07 181 ST* 08 182 RCL 09 183 * 184 RCL 19 185 - 186 RCL 22 187 STO 09 188 * 189 ST- IND 10 190 DSE 10 191 RCL 08 192 RCL 18 193 - 194 RCL 09 195 * 196 ST- IND 10 197 DSE 10 198 RCL 07 199 RCL 17 200 - 201 RCL 09 202 * 203 ST- IND 10 204 RCL 25 205 ST+ 06 206 ST+ 10 207 RCL 04 208 STO 08 209 RCL 05 210 STO 07 211 LBL 04 212 RCL IND 08 213 RCL IND 06 214 - 215 ENTER 216 X^2 217 DSE 06 218 DSE 08 219 RCL IND 08 220 RCL IND 06 221 - 222 STO 09 223 X^2 224 + 225 DSE 06 226 DSE 08 227 RCL IND 08 228 RCL IND 06 229 - 230 STO T 231 X^2 232 + 233 ENTER 234 SQRT 235 * 236 X=0? 237 SIGN 238 RCL 24 239 X<>Y 240 / 241 RCL IND 07 242 * 243 ST* 09 244 ST* T 245 * 246 ST+ IND 10 247 DSE 10 248 RCL 09 249 ST+ IND 10 250 DSE 10 251 R^ 252 ST+ IND 10 253 CLX 254 SIGN 255 ST- 07 256 RCL 25 257 ST+ 06 258 ST+ 10 259 DSE 08 260 GTO 04 261 ST- 06 262 PI 263 INT 264 ST- 10 265 DSE 06 266 GTO 03 267 RCL 14 268 STO 06 269 RCL 11 270 STO 07 271 RCL 04 272 STO 08 273 RTN 274 LBL 01 275 RCL 14 276 STO 10 277 XEQ 10 278 LBL 05 279 RCL IND 07 280 RCL IND 06 281 .4 282 ST* Z 283 * 284 ST+ IND 07 285 5 286 / 287 + 288 RCL 02 289 * 290 ST+ IND 08 291 CLX 292 SIGN 293 ST- 06 294 ST- 07 295 DSE 08 296 GTO 05 297 XEQ 10 298 RCL 13 299 STO 09 300 LBL 07 301 RCL IND 06 302 147 303 * 304 25 305 / 306 RCL IND 09 307 3 308 * 309 - 310 256 311 / 312 RCL IND 07 313 .15 314 * 315 - 316 RCL 02 317 * 318 ST+ IND 08 319 RCL IND 09 320 25 321 * 322 RCL IND 06 323 73 324 * 325 - 326 320 327 / 328 ST+ IND 07 329 CLX 330 SIGN 331 ST- 06 332 ST- 09 333 ST- 07 334 DSE 08 335 GTO 07 336 XEQ 10 337 RCL 13 338 STO 09 339 RCL 12 340 STO 17 341 LBL 08 342 RCL IND 07 343 RCL 26 344 / 345 RCL IND 06 346 11 347 * 348 RCL IND 09 349 + 350 128 351 / 352 - 353 RCL IND 17 354 8 355 / 356 + 357 RCL 02 358 * 359 ST+ IND 08 360 RCL IND 17 361 RCL 25 362 / 363 RCL IND 06 364 11 365 * 366 RCL IND 09 367 5 368 * 369 + 370 64 371 / 372 - 373 ST+ IND 07 374 CLX 375 SIGN 376 ST- 06 377 ST- 07 378 ST- 09 379 ST- 17 380 DSE 08 381 GTO 08 382 XEQ 10 383 RCL 13 384 STO 09 385 RCL 12 386 STO 17 387 RCL 28 388 STO 18 389 LBL 09 390 RCL IND 06 391 9 392 * 393 RCL IND 09 394 21 395 * 396 - 397 RCL 26 398 / 399 RCL IND 18 400 9 401 * 402 + 403 RCL 26 404 / 405 RCL IND 17 406 - 407 RCL 26 408 / 409 RCL IND 07 410 + 411 RCL 26 412 / 413 RCL 02 414 * 415 ST+ IND 08 416 RCL IND 06 417 3 418 * 419 RCL IND 09 420 15 421 * 422 - 423 RCL 26 424 / 425 RCL IND 18 426 9 427 * 428 + 429 RCL 25 430 / 431 RCL IND 17 432 - 433 8 434 / 435 ST+ IND 07 436 CLX 437 SIGN 438 ST- 06 439 ST- 09 440 ST- 17 441 ST- 18 442 ST- 07 443 DSE 08 444 GTO 09 445 XEQ 10 446 RCL 13 447 STO 09 448 RCL 12 449 STO 17 450 RCL 28 451 STO 18 452 RCL 27 453 STO 19 454 LBL 14 455 RCL IND 09 456 255 457 * 458 RCL IND 17 459 162 460 * 461 + 462 RCL IND 18 463 510 464 * 465 - 466 7 467 / 468 RCL IND 06 469 3 470 * 471 - 472 16 473 / 474 RCL IND 07 475 + 476 RCL 26 477 / 478 RCL IND 19 479 ST+ X 480 7 481 / 482 STO 20 483 + 484 RCL 02 485 * 486 ST+ IND 08 487 RCL IND 09 488 425 489 * 490 RCL IND 17 491 216 492 * 493 + 494 RCL IND 18 495 1020 496 * 497 - 498 7 499 / 500 RCL IND 06 501 3 502 * 503 - 504 64 505 / 506 RCL 20 507 4 508 * 509 + 510 ST+ IND 07 511 RCL IND 18 512 582 513 * 514 RCL IND 17 515 158 516 * 517 - 518 RCL IND 19 519 248 520 * 521 - 522 45 523 / 524 RCL IND 09 525 5 526 * 527 - 528 7 529 / 530 RCL IND 06 531 LASTX 532 * 533 90 534 / 535 STO 20 536 + 537 X<> IND 19 538 62 539 * 540 CHS 541 RCL IND 18 542 291 543 * 544 + 545 RCL IND 17 546 118.5 547 * 548 - 549 45 550 / 551 RCL IND 09 552 3 553 * 554 - 555 7 556 / 557 RCL 20 558 + 559 STO IND 18 560 CLX 561 SIGN 562 ST- 06 563 ST- 09 564 ST- 17 565 ST- 18 566 ST- 07 567 ST- 19 568 DSE 08 569 GTO 14 570 RCL 14 571 STO 10 572 XEQ 10 573 RCL 27 574 STO 10 575 RCL 28 576 STO 09 577 LBL 00 578 RCL IND 09 579 RCL 02 580 * 581 ST+ IND 08 582 RCL IND 06 583 7 584 * 585 90 586 / 587 RCL IND 10 588 + 589 ST+ IND 07 590 CLX 591 SIGN 592 ST- 06 593 ST- 09 594 ST- 07 595 ST- 10 596 DSE 08 597 GTO 00 598 RCL 02 599 ST+ 29 600 VIEW 29 601 DSE 16 602 GTO 01 603 RCL 29 604 END

( 942 bytes / SIZE 31+22.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:   Integrate the Solar system with h = 0.2 day , starting with the initial values in the 1st column below 2019/11/13 0h TDB

9   STO 01
0.2  STO 02
160  STO 03   to obtain the positions 32 days later ( 2019/12/15 0h TDB )

( cf paragraph 3°)a) for the initialization routine )

-Then   XEQ "PL5G"    >>>>    X = 32                                        ---Execution time = 7m21s---              With V41 in turbo mode

-Here are the results ( in R40 to R66 )

 2019/11/13 0h        Initial Values 2019/12/15 0h             DE436 2019/12/15 0h         HP41  PL5G M 169.3287041 E-3 E  236.7984713 E-3 R  108.9436218 E-3 -0.362228680      -0.226192417      -0.083282304 -0.362228683      -0.226192417      -0.083282305 V 206.8552787 E-3 E -631.2899640 E-3 N -297.1381466 E-3 0.669481274      -0.241877222      -0.151193495 0.669481273      -0.241877222      -0.151193495 E  635.7893837 E-3 M 695.9559735 E-3 B  301.6942285 E-3 0.129773202       0.895203527       0.388069088 0.129773211       0.895203532       0.388069091 M-1604.907774E-3 A -301.5483654E-3 R -94.99990761E-3 -1.447205476      -0.650423094      -0.259275988 -1.447205467      -0.650423098      -0.259275995 J  161.3509664 E-3 U -4817.408291E-3 P -2068.801992 E-3 0.399788668      -4.792899165      -2.064100995 0.399788668      -4.792899163      -2.064100997 S  3557.619174E-3 A -8621.479227E-3 T -3714.330343E-3 3.714404180      -8.559663453      -3.695545437 3.714404190      -8.559663444      -3.695545433 U  16336.48164 E-3 R  10370.59623 E-3 A  4311.051187 E-3 16.26417108       10.45967901       4.351085196 16.26417103       10.45967900       4.351085196 N 29210.34115 E-3 E -5766.105481 E-3 P -3087.429294 E-3 29.23163521      -5.674517011      -3.050471876 29.23163520      -5.674517012      -3.050471878 P 12830.55991 E-3 L-28665.41178 E-3 U-12811.23541 E-3 12.92603848      -28.63831132      -12.83146648 12.92603846      -28.63831130      -12.83146649

Notes:

-With h = 0.2 , the precision is almost as good as with h = 0.1 & "GM4G"
-These results were obtained with V41.
-With free42, the roundoff-errors are negligible

-The maximum n-value is 13 with an HP41 ( 14-body problem )

4°)  General Relativity - Barycentric Coordinates

a) 4th-order Runge-Kutta Formula

-"GM4G" employs the same Runge-Kutta formula of order 4 as "PLN4G"

-In "PL4G" the relativistic perturbations from the mass of the Sun are taken into account.
-Here, the relativistic corrections from the masses of the other celestial bodies are also taken into account .

-So "GM4G" employs the formula given in reference [4]:

ri'' = Sumj#i ( Gmj/rij3 ) { 1 - ( 4/c2 ) Sumk#i Gmk/rik - ( 1/c2 ) Sumk#j Gmk/rjk + ( vi2 / c2 ) + 2 ( vj2 / c2 ) + ( 4/c2 ) vi.vj
- ( 3/c2 ) [ ( rij.vj ) / rij ]2 + ( 1/c2rij . rj" }
+ ( 1/c2 ) Sumj#i ( Gmj/rij3 ) [ rij. ( 4.vi - 3.vj ) ] vij + ( 7/(2.c2) ) Sumj#i ( Gmj/rij ) rj"

-Since the acceleration ri'' appears in the right-hand side of the formula - but multiplied by ( 1/c2 ) - the Newtonian acceleration may be used there.

Data Registers:              R00 = G           ( Registers R01-R02-R03 & R30 thru R30+7n are to be initialized before executing "GM4G" )

•  R01 = n = nb of celestial bodies                                                          R04 to R36&R39: temp    R37-R38: unused
•  R02 = h = step size ( in days )            •  R40 = t0 = initial instant
•  R03 = N = nb of steps                       •  R41 = m1 = mass of the 1st celestial body
................................................

•  R40+n = mn = mass of the n-th celestial body

•  R41+n = x1              •  R41+4n = x'1                      R41+7n  to  R40+10n = k3 , k4
•  R42+n = y1              •  R42+4n = y'1                      R41+10n to R40+13n = k2 , 2.k3 - k2
•  R43+n = z1              •  R43+4n = z1                       R41+13n to R40+16n = k1
.................                   ..................
R41+16.n to R40+19.n = ri''  ( Newtonian )
•  R38+4n = xn            •  R38+7n = x'n                     R41+19.n to R40+20.n = Sumj#i Gmj/rij
•  R39+4n = yn            •  R39+7n = y'n
•  R40+4n = zn            •  R40+7n = z'n
Flags: /
Subroutines: /

-Lines 41-134-377-385-504 are three-byte GTOs

 01 LBL "GM4G"  02 RCL 03  03 STO 39  04 3379.380681   05 1/X  06 STO 00  07 29979.06382  08 STO 04  09 40.04  10 RCL 01  11 +  12 STO 17  13 LASTX  14 .1  15 %  16 X<>Y  17 3  18 *  19 STO T  20 +  21 +  22 STO 18  23 INT  24 .1  25 %  26 +  27 +  28 STO 19  29 INT  30 +  31 STO 20  32 +  33 STO 21  34 +  35 STO 22  36 +  37 STO 23  38 RCL 01  39 +  40 STO 24  41 GTO 01  42 LBL 10  43 RCL 18  44 STO 06  45 STO 26  46 RCL 23  47 STO 07  48 RCL 24  49 STO 10  50 LBL 00  51 RCL 18  52 STO 08  53 RCL 17  54 STO 09  55 RCL 07  56 ENTER  57 CLX  58 STO IND Y  59 DSE Y  60 STO IND Y  61 DSE Y  62 STO IND Y  63 STO IND 10  64 RCL IND 06  65 STO 15  66 DSE 06  67 RCL IND 06  68 STO 14  69 DSE 06  70 RCL IND 06  71 STO 13  72 LBL 02 73 RCL IND 08  74 RCL 15  75 -  76 ENTER  77 X^2  78 DSE 08  79 RCL IND 08  80 RCL 14  81 -  82 STO 11  83 X^2  84 +  85 DSE 08  86 RCL IND 08  87 RCL 13  88 -  89 STO T  90 X^2  91 +  92 ENTER  93 SQRT  94 STO 12  95 *  96 X=0?  97 SIGN  98 RCL 00  99 X<>Y 100 / 101 RCL IND 09 102 * 103 ST* 11 104 ST* T 105 * 106 ST+ IND 07 107 DSE 07 108 RCL 11 109 ST+ IND 07 110 DSE 07 111 R^ 112 ST+ IND 07 113 RCL 00 114 RCL IND 09 115 * 116 RCL 12 117 X#0? 118 1/X 119 * 120 ST+ IND 10 121 CLX 122 SIGN 123 ST- 09 124 ST+ X 125 ST+ 07 126 DSE 08 127 GTO 02 128 PI 129 INT 130 ST- 07 131 SIGN 132 ST- 10 133 DSE 06 134 GTO 00 135 RCL 19 136 STO 27 137 RCL 24 138 STO 25 139 LBL 03 140 RCL 16 141 ENTER 142 CLX 143 STO IND Y 144 DSE Y 145 STO IND Y 146 DSE Y 147 STO IND Y 148 RCL IND 26 149 STO 10 150 DSE 26 151 RCL IND 26 152 STO 09 153 DSE 26 154 RCL IND 26 155 STO 08 156 RCL IND 27 157 STO 13 158 X^2 159 DSE 27 160 RCL IND 27 161 STO 12 162 X^2 163 + 164 DSE 27 165 RCL IND 27 166 STO 11 167 X^2 168 + 169 STO 14 170 RCL 18 171 STO 28 172 RCL 19 173 STO 29 174 RCL 17 175 STO 30 176 RCL 23 177 STO 31 178 RCL 24 179 STO 32 180 LBL 04 181 RCL IND 29 182 STO 07 183 X^2 184 DSE 29 185 RCL IND 29 186 STO 06 187 X^2 188 + 189 DSE 29 190 RCL IND 29 191 STO 05 192 X^2 193 + 194 RCL 05 195 RCL 11 196 * 197 RCL 06 198 RCL 12 199 * 200 + 201 RCL 07 202 RCL 13 203 * 204 + 205 ST+ X 206 - 207 ST+ X 208 RCL 14 209 + 210 RCL IND 28 211 RCL 10 212 - 213 STO 35 214 RCL 07 215 * 216 DSE 28 217 RCL IND 28 218 RCL 09 219 - 220 STO 34 221 RCL 06 222 * 223 + 224 DSE 28 225 RCL IND 28 226 RCL 08 227 - 228 STO 33 229 RCL 05 230 * 231 + 232 X^2 233 RCL 33 234 X^2 235 RCL 34 236 X^2 237 + 238 RCL 35 239 X^2 240 + 241 X#0? 242 1/X 243 STO 36 244 * 245 1.5 246 * 247 - 248 RCL 35 249 RCL IND 31 250 * 251 RCL 34 252 DSE 31 253 RCL IND 31 254 * 255 + 256 RCL 33 257 DSE 31 258 RCL IND 31 259 * 260 + 261 2 262 ST+ 31 263 / 264 + 265 RCL IND 32 266 - 267 RCL IND 25 268 4 269 * 270 - 271 RCL 04 272 / 273 1 274 + 275 STO 15 276 RCL 05 277 RCL 11 278 - 279 STO 05 280 3 281 * 282 RCL 11 283 - 284 RCL 33 285 * 286 RCL 06 287 RCL 12 288 - 289 STO 06 290 3 291 * 292 RCL 12 293 - 294 RCL 34 295 * 296 + 297 RCL 07 298 RCL 13 299 - 300 STO 07 301 3 302 * 303 RCL 13 304 - 305 RCL 35 306 * 307 + 308 RCL 04 309 / 310 ST* 05 311 ST* 06 312 ST* 07 313 RCL 15 314 ST* 33 315 ST* 34 316 RCL 35 317 * 318 RCL 07 319 - 320 RCL 36 321 * 322 3.5 323 RCL 04 324 / 325 STO 07 326 RCL IND 31 327 * 328 + 329 RCL 00 330 RCL 02 331 * 332 RCL IND 30 333 * 334 RCL 36 335 SQRT 336 * 337 STO 15 338 * 339 ST+ IND 16 340 DSE 16 341 DSE 31 342 RCL 34 343 RCL 06 344 - 345 RCL 36 346 * 347 RCL 07 348 RCL IND 31 349 * 350 + 351 RCL 15 352 * 353 ST+ IND 16 354 DSE 16 355 DSE 31 356 RCL 33 357 RCL 05 358 - 359 RCL 36 360 * 361 RCL 07 362 RCL IND 31 363 * 364 + 365 RCL 15 366 * 367 ST+ IND 16 368 CLX 369 SIGN 370 ST- 29 371 ST- 30 372 ST- 31 373 ST- 32 374 ST+ X 375 ST+ 16 376 DSE 28 377 GTO 04 378 SIGN 379 ST- 25 380 ST- 27 381 PI 382 INT 383 ST- 16 384 DSE 26 385 GTO 03 386 RCL 22 387 STO 06 388 RCL 19 389 STO 12 390 RCL 18 391 STO 13 392 RTN 393 LBL 01 394 RCL 22 395 STO 16 396 XEQ 10 397 LBL 05 398 RCL IND 12 399 RCL IND 06 400 4 401 / 402 + 403 RCL 02 404 * 405 2 406 / 407 ST+ IND 13 408 RCL IND 06 409 LASTX 410 / 411 ST+ IND 12 412 CLX 413 SIGN 414 ST- 06 415 ST- 12 416 DSE 13 417 GTO 05 418 XEQ 10 419 RCL 21 420 STO 09 421 LBL 07 422 RCL IND 09 423 RCL IND 06 424 - 425 2 426 / 427 ST+ IND 12 428 CLX 429 SIGN 430 ST- 06 431 ST- 09 432 DSE 12 433 GTO 07 434 XEQ 10 435 RCL 21 436 STO 09 437 RCL 20 438 STO 10 439 LBL 08 440 RCL IND 12 441 RCL IND 06 442 4 443 / 444 - 445 RCL IND 10 446 ST+ X 447 RCL IND 09 448 - 449 STO IND 09 450 2 451 / 452 ST+ IND 12 453 + 454 RCL 02 455 * 456 2 457 / 458 ST+ IND 13 459 CLX 460 SIGN 461 ST- 06 462 ST- 09 463 ST- 10 464 ST- 12 465 DSE 13 466 GTO 08 467 RCL 20 468 STO 16 469 XEQ 10 470 RCL 21 471 STO 09 472 RCL 20 473 STO 10 474 LBL 09 475 RCL IND 06 476 RCL IND 09 477 - 478 RCL 02 479 * 480 6 481 / 482 ST+ IND 13 483 RCL IND 06 484 RCL IND 09 485 ST+ X 486 - 487 RCL IND 10 488 + 489 6 490 / 491 ST+ IND 12 492 CLX 493 SIGN 494 ST- 06 495 ST- 09 496 ST- 10 497 ST- 12 498 DSE 13 499 GTO 09 500 RCL 02 501 ST+ 40 502 VIEW 40 503 DSE 39 504 GTO 01 505 RCL 40 506 END

( 771 bytes / SIZE 41+20.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example1:   Integrate the Solar system with h = 0.1 day , starting with the initial values in the 1st column below 2019/11/13 0h TDB

10     STO 01
0.1    STO 02
10000  STO 03   to obtain the positions 1000 days later ( 2022/08/09 0h TDB )

( Execute for instance the following routine )

 01 LBL "INGM"  02 1  03 STO 31  04 6023682.156  05 1/X  06 STO 32  07 408523.7187  08 1/X  09 STO 33  10 328900.5598  11 1/X  12 STO 34  13 3098703.59  14 1/X  15 STO 35  16 1047.348625   17 1/X  18 STO 36  19 3497.901768   20 1/X  21 STO 37  22 22902.98161  23 1/X  24 STO 38  25 19412.25978         26 1/X  27 STO 39  28 135836683.8  29 1/X  30 STO 40 31 -0.003386657126  32 STO 41  33 0.006899049618  34 STO 42  35 0.003003849853  36 STO 43  37 0.1659420469  38 STO 44  39 0.2436975209  40 STO 45  41 0.1119474717  42 STO 46  43 0.2034686215  44 STO 47  45 -.6243909144  46 STO 48  47 -0.2941342968  48 STO 49  49 0.6324027266  50 STO 50  51 0.7028550231  52 STO 51  53 0.3046980783  54 STO 52  55 -1.608294431  56 STO 53  57 -0.2946493158  58 STO 54  59 -0.09199605775   60 STO 55 61 0.1579643093  62 STO 56  63 -4.810509241  64 STO 57  65 -2.065798143  66 STO 58  67 3.554232517  68 STO 59  69 -8.614580178  70 STO 60  71 -3.711326494  72 STO 61  73 16.33309498  74 STO 62  75 10.37749528  76 STO 63  77 4.314055037  78 STO 64  79 29.2069545  80 STO 65  81 -5.759206431  82 STO 66  83 -3.084425444  84 STO 67  85 12.82717325  86 STO 68  87 -28.65851273  88 STO 69  89 -12.80823156   90 STO 70 91 -8.454840562 E-6  92 STO 71  93 -1.431104541 E-6  94 STO 72  95 -3.683195062 E-7  96 STO 73  97 -0.02920440361  98 STO 74  99 0.01350155574 100 STO 75 101 0.01023923028 102 STO 76 103 0.01924697229 104 STO 77 105 0.005623903436 106 STO 78 107 0.001312394143 108 STO 79 109 -0.01347458771 110 STO 80 111 0.01007878298 112 STO 81 113 0.00436942538 114 STO 82 115 0.003203696527 116 STO 83 117 -0.01136722313 118 STO 84 119 -0.005300259888 120 STO 85 121 0.007450644773 122 STO 86 123 6.061976969 E-4 124 STO 87 125 7.850083244 E-5 126 STO 88 127 0.004908104128 128 STO 89 129 0.001889774855 130 STO 90 131 5.691831284 E-4 132 STO 91 133 -0.002258208799 134 STO 92 135 0.002788572982 136 STO 93 137 0.001253256089 138 STO 94 139 6.621791223 E-4 140 STO 95 141 0.002859533729 142 STO 96 143 0.00115393786 144 STO 97 145 0.002976842608 146 STO 98 147 8.418146302 E-4 148 STO 99 149 100 150 -6.342164234 E-4 151 STO IND Y 152 END

CLX  STO 40   XEQ "GM4G"   and you get the following results ( in about 213 seconds with free42binary ) in registers R51 to R80:

 2013/11/13 0h             DE436 2022/08/09 0h              DE436 2022/08/09 0h        free42  GM4G S    -0.003386657126 U     0.006899049618 N     0.003003849853 -0.009061149270      0.001213932393      0.0007438465588 -0.0090611537       0.0012139328       0.0007438472 M     0.1659420469 E      0.2436975209 R      0.1119474717 -0.3679022839      -0.2314424565      -0.0863475922 -0.3679022823      -0.2314424638      -0.0863475964 V      0.2034686215  E     -0.6243909144  N    -0.2941342968 0.01833998496       0.6577717218       0.2944345025 0.0183399722       0.6577717220       0.2944345037 E      0.6324027266 M     0.7028550231 B      0.3046980783 0.7216316454      -0.6437330727      -0.2788405241 0.7216304734      -0.6437342480      -0.2788410724 M     -1.608294431 A     -0.2946493158 R     -0.09199605775 1.3788166042       0.1638497441       0.03789518642 1.3788165924       0.1638497657       0.0378951967 J       0.1579643093 U     -4.810509241 P     -2.065798143 4.9476760534      -0.04102757478      -0.1380171712 4.9476760604      -0.0410275855      -0.1380171765 S      3.554232517 A    -8.614580178 T    -3.711326494 7.7042645359      -5.5674807876      -2.6314856814 7.7042645376      -5.5674807902      -2.6314856827 U      16.33309498 R      10.37749528 A      4.314055037 13.777783853       12.950367747       5.4770420499 13.7777838535       12.9503677508       5.4770420500 N      29.20695450 E      -5.759206431 P      -3.084425444 29.706554081      -2.8731350804      -1.9155780574 29.706554086      -2.8731350808      -1.9155780579 P      12.82717325 L     -28.65851273 U     -12.80823156 15.752575014      -27.710976967      -13.393955261 15.752575012      -27.710976963      -13.393955256

-The accuracy is rather good - except for the Earth-Moon barycenter !
-"GM4G" may also be used with Earth & Moon separately.

Notes:

-Execution time is about 70 minutes per iteration with a real HP41, so I've not tried to do 10000 iterations !
-The roundoff-errors would be also much larger...

Example2:      3 stars , let's calculate the positions 100 days after the input data below:

 Register input h =0.02 N=5000 m R41=m1 R42=m2 R43=m3 3    2    1 unchanged p  o  s. R44=x1 R45=y1 R46=z1 -------- R47=x2 R48=y2 R49=z2 -------- R50=x3 R51=y3 R52=z3 0    0    0 -----    1    0    0 -----   -2    0    0 -0.113356390  1.196226886 -0.398742295  --------------  0.839927385  1.026872301 -0.342290767  -------------- -1.339785645 -5.642425267  1.880808422 v  e  l  . R53=x'1 R54=y'1 R55=z'1 -------- R56=x'2 R57=y'2 R58=z'2 -------- R59=x'3 R60=y'3 R61=z'3 0    0    0 -----    0  0.03 -0.01 -----    0 -0.06  0.02 -0.005571436 -0.001076303  0.000358768  --------------  0.004102293  0.028530254 -0.009510085  --------------  0.008509721 -0.053831599  0.017943866

-These values are obtained with free42.
-The results in the last column above are the same ( rounded to 9 decimals ) with h = 0.01 , h = 0.02 , h = 0.04
-With h = 0.1 day, the differences are - at most - about E-9 AU

Notes:

-The maximum n-value is 13 with HP41/V41 but if you replace line 09 by 38.038 & lines 501 to 505
by  ST+38  VIEW 38   DSE 37   GTO 01  RCL 38,  the maximum n-value becomes 14 ( with SIZE 319 )

b) 5th-order Runge-Kutta Formula

-"GM5G" employs the same formula as "PL5G"

Data Registers:              R00 = G            ( Registers R01-R02-R03 & R30 thru R30+7n are to be initialized before executing "GM5G" )

•  R01 = n = nb of celestial bodies                                                          R04 to R39: temp
•  R02 = h = step size ( in days )            •  R40 = t0 = initial instant
•  R03 = N = nb of steps                       •  R41 = m1 = mass of the 1st celestial body
................................................

•  R40+n = mn = mass of the n-th celestial body

•  R41+n = x1              •  R41+4n = x'1
•  R42+n = y1              •  R42+4n = y'1                      R41+7n  to R40+10n = k5 and linear combination of other ki
•  R43+n = z1              •  R43+4n = z1                       R41+10n to R40+13n = k4 and linear combination of other ki
.................                   ..................                        R41+13n to R40+16n = k3
R41+16n to R40+19n = k2
•  R38+4n = xn            •  R38+7n = x'n                     R41+19n to R40+22n = k1 , k6
•  R39+4n = yn            •  R39+7n = y'n                     R41+22.n to R40+25.n = ri''  ( Newtonian )
•  R40+4n = zn            •  R40+7n = z'n                     R41+25.n to R40+26.n = Sumj#i Gmj/rij
Flags: /
Subroutines: /

-Lines 41-134-377-385-688-721 are three-byte GTOs

 01 LBL "GM5G"  02 RCL 03  03 STO 37  04 3379.380681   05 1/X  06 STO 00  07 29979.06382  08 STO 04  09 40.04  10 RCL 01  11 +  12 STO 17  13 LASTX  14 .1  15 %  16 X<>Y  17 3  18 *  19 STO T  20 +  21 +  22 STO 18  23 INT  24 +  25 STO 19  26 +  27 STO 38  28 +  29 STO 39  30 +  31 STO 20  32 +  33 STO 21  34 +  35 STO 22  36 +  37 STO 23  38 RCL 01  39 +  40 STO 24  41 GTO 01  42 LBL 10  43 RCL 18  44 STO 06  45 STO 26  46 RCL 23  47 STO 07  48 RCL 24  49 STO 10  50 LBL 00  51 RCL 18  52 STO 08  53 RCL 17  54 STO 09  55 RCL 07  56 ENTER  57 CLX  58 STO IND Y  59 DSE Y  60 STO IND Y  61 DSE Y  62 STO IND Y  63 STO IND 10  64 RCL IND 06  65 STO 15  66 DSE 06  67 RCL IND 06  68 STO 14  69 DSE 06  70 RCL IND 06  71 STO 13  72 LBL 02  73 RCL IND 08  74 RCL 15  75 -  76 ENTER  77 X^2  78 DSE 08  79 RCL IND 08  80 RCL 14  81 -  82 STO 11  83 X^2  84 +  85 DSE 08  86 RCL IND 08  87 RCL 13  88 -  89 STO T  90 X^2  91 +  92 ENTER  93 SQRT  94 STO 12  95 *  96 X=0?  97 SIGN  98 RCL 00  99 X<>Y 100 / 101 RCL IND 09 102 * 103 ST* 11 104 ST* T 105 * 106 ST+ IND 07 107 DSE 07 108 RCL 11 109 ST+ IND 07 110 DSE 07 111 R^ 112 ST+ IND 07 113 RCL 00 114 RCL IND 09 115 * 116 RCL 12 117 X#0? 118 1/X 119 * 120 ST+ IND 10 121 CLX 122 SIGN 123 ST- 09 124 ST+ X 125 ST+ 07 126 DSE 08 127 GTO 02 128 PI 129 INT 130 ST- 07 131 SIGN 132 ST- 10 133 DSE 06 134 GTO 00 135 RCL 19 136 STO 27 137 RCL 24 138 STO 25 139 LBL 03 140 RCL 16 141 ENTER 142 CLX 143 STO IND Y 144 DSE Y 145 STO IND Y 146 DSE Y 147 STO IND Y 148 RCL IND 26 149 STO 10 150 DSE 26 151 RCL IND 26 152 STO 09 153 DSE 26 154 RCL IND 26 155 STO 08 156 RCL IND 27 157 STO 13 158 X^2 159 DSE 27 160 RCL IND 27 161 STO 12 162 X^2 163 + 164 DSE 27 165 RCL IND 27 166 STO 11 167 X^2 168 + 169 STO 14 170 RCL 18 171 STO 28 172 RCL 19 173 STO 29 174 RCL 17 175 STO 30 176 RCL 23 177 STO 31 178 RCL 24 179 STO 32 180 LBL 04 181 RCL IND 29 182 STO 07 183 X^2 184 DSE 29 185 RCL IND 29 186 STO 06 187 X^2 188 + 189 DSE 29 190 RCL IND 29 191 STO 05 192 X^2 193 + 194 RCL 05 195 RCL 11 196 * 197 RCL 06 198 RCL 12 199 * 200 + 201 RCL 07 202 RCL 13 203 * 204 + 205 ST+ X 206 - 207 ST+ X 208 RCL 14 209 + 210 RCL IND 28 211 RCL 10 212 - 213 STO 35 214 RCL 07 215 * 216 DSE 28 217 RCL IND 28 218 RCL 09 219 - 220 STO 34 221 RCL 06 222 * 223 + 224 DSE 28 225 RCL IND 28 226 RCL 08 227 - 228 STO 33 229 RCL 05 230 * 231 + 232 X^2 233 RCL 33 234 X^2 235 RCL 34 236 X^2 237 + 238 RCL 35 239 X^2 240 + 241 X#0? 242 1/X 243 STO 36 244 * 245 1.5 246 * 247 - 248 RCL 35 249 RCL IND 31 250 * 251 RCL 34 252 DSE 31 253 RCL IND 31 254 * 255 + 256 RCL 33 257 DSE 31 258 RCL IND 31 259 * 260 + 261 2 262 ST+ 31 263 / 264 + 265 RCL IND 32 266 - 267 RCL IND 25 268 4 269 * 270 - 271 RCL 04 272 / 273 1 274 + 275 STO 15 276 RCL 05 277 RCL 11 278 - 279 STO 05 280 3 281 * 282 RCL 11 283 - 284 RCL 33 285 * 286 RCL 06 287 RCL 12 288 - 289 STO 06 290 3 291 * 292 RCL 12 293 - 294 RCL 34 295 * 296 + 297 RCL 07 298 RCL 13 299 - 300 STO 07 301 3 302 * 303 RCL 13 304 - 305 RCL 35 306 * 307 + 308 RCL 04 309 / 310 ST* 05 311 ST* 06 312 ST* 07 313 RCL 15 314 ST* 33 315 ST* 34 316 RCL 35 317 * 318 RCL 07 319 - 320 RCL 36 321 * 322 3.5 323 RCL 04 324 / 325 STO 07 326 RCL IND 31 327 * 328 + 329 RCL 00 330 RCL 02 331 * 332 RCL IND 30 333 * 334 RCL 36 335 SQRT 336 * 337 STO 15 338 * 339 ST+ IND 16 340 DSE 16 341 DSE 31 342 RCL 34 343 RCL 06 344 - 345 RCL 36 346 * 347 RCL 07 348 RCL IND 31 349 * 350 + 351 RCL 15 352 * 353 ST+ IND 16 354 DSE 16 355 DSE 31 356 RCL 33 357 RCL 05 358 - 359 RCL 36 360 * 361 RCL 07 362 RCL IND 31 363 * 364 + 365 RCL 15 366 * 367 ST+ IND 16 368 CLX 369 SIGN 370 ST- 29 371 ST- 30 372 ST- 31 373 ST- 32 374 ST+ X 375 ST+ 16 376 DSE 28 377 GTO 04 378 SIGN 379 ST- 25 380 ST- 27 381 PI 382 INT 383 ST- 16 384 DSE 26 385 GTO 03 386 RCL 22 387 STO 06 388 RCL 19 389 STO 12 390 RCL 18 391 STO 13 392 RTN 393 LBL 01 394 RCL 22 395 STO 16 396 XEQ 10 397 LBL 05 398 RCL IND 12 399 RCL IND 06 400 .4 401 ST* Z 402 * 403 ST+ IND 12 404 5 405 / 406 + 407 RCL 02 408 * 409 ST+ IND 13 410 CLX 411 SIGN 412 ST- 06 413 ST- 12 414 DSE 13 415 GTO 05 416 XEQ 10 417 RCL 21 418 STO 09 419 LBL 07 420 RCL IND 06 421 147 422 * 423 25 424 / 425 RCL IND 09 426 3 427 * 428 - 429 256 430 / 431 RCL IND 12 432 .15 433 * 434 - 435 RCL 02 436 * 437 ST+ IND 13 438 RCL IND 09 439 25 440 * 441 RCL IND 06 442 73 443 * 444 - 445 320 446 / 447 ST+ IND 12 448 CLX 449 SIGN 450 ST- 06 451 ST- 09 452 ST- 12 453 DSE 13 454 GTO 07 455 XEQ 10 456 RCL 21 457 STO 09 458 RCL 20 459 STO 10 460 LBL 08 461 RCL IND 12 462 4 463 / 464 RCL IND 06 465 11 466 * 467 RCL IND 09 468 + 469 128 470 / 471 - 472 RCL IND 10 473 8 474 / 475 + 476 RCL 02 477 * 478 ST+ IND 13 479 RCL IND 10 480 2 481 / 482 RCL IND 06 483 11 484 * 485 RCL IND 09 486 5 487 * 488 + 489 64 490 / 491 - 492 ST+ IND 12 493 CLX 494 SIGN 495 ST- 06 496 ST- 09 497 ST- 10 498 ST- 12 499 DSE 13 500 GTO 08 501 XEQ 10 502 RCL 21 503 STO 09 504 RCL 20 505 STO 10 506 RCL 39 507 STO 11 508 LBL 09 509 RCL IND 06 510 9 511 * 512 RCL IND 09 513 21 514 * 515 - 516 4 517 / 518 RCL IND 11 519 9 520 * 521 + 522 4 523 / 524 RCL IND 10 525 - 526 4 527 / 528 RCL IND 12 529 + 530 4 531 / 532 RCL 02 533 * 534 ST+ IND 13 535 RCL IND 06 536 3 537 * 538 RCL IND 09 539 15 540 * 541 - 542 4 543 / 544 RCL IND 11 545 9 546 * 547 + 548 2 549 / 550 RCL IND 10 551 - 552 8 553 / 554 ST+ IND 12 555 CLX 556 SIGN 557 ST- 06 558 ST- 09 559 ST- 10 560 ST- 11 561 ST- 12 562 DSE 13 563 GTO 09 564 XEQ 10 565 RCL 21 566 STO 09 567 RCL 20 568 STO 10 569 RCL 39 570 STO 11 571 RCL 38 572 STO 14 573 LBL 14 574 RCL IND 09 575 255 576 * 577 RCL IND 10 578 162 579 * 580 + 581 RCL IND 11 582 510 583 * 584 - 585 7 586 / 587 RCL IND 06 588 3 589 * 590 - 591 16 592 / 593 RCL IND 12 594 + 595 4 596 / 597 RCL IND 14 598 ST+ X 599 7 600 / 601 STO 15 602 + 603 RCL 02 604 * 605 ST+ IND 13 606 RCL IND 09 607 425 608 * 609 RCL IND 10 610 216 611 * 612 + 613 RCL IND 11 614 1020 615 * 616 - 617 7 618 / 619 RCL IND 06 620 3 621 * 622 - 623 64 624 / 625 RCL 15 626 4 627 * 628 + 629 ST+ IND 12 630 RCL IND 11 631 582 632 * 633 RCL IND 10 634 158 635 * 636 - 637 RCL IND 14 638 248 639 * 640 - 641 45 642 / 643 RCL IND 09 644 5 645 * 646 - 647 7 648 / 649 RCL IND 06 650 LASTX 651 * 652 90 653 / 654 STO 15 655 + 656 X<> IND 14 657 62 658 * 659 CHS 660 RCL IND 11 661 291 662 * 663 + 664 RCL IND 10 665 118.5 666 * 667 - 668 45 669 / 670 RCL IND 09 671 3 672 * 673 - 674 7 675 / 676 RCL 15 677 + 678 STO IND 11 679 CLX 680 SIGN 681 ST- 06 682 ST- 09 683 ST- 10 684 ST- 11 685 ST- 12 686 ST- 14 687 DSE 13 688 GTO 14 689 RCL 22 690 STO 16 691 XEQ 10 692 RCL 38 693 STO 14 694 RCL 39 695 STO 11 696 LBL 13 697 RCL IND 11 698 RCL 02 699 * 700 ST+ IND 13 701 RCL IND 06 702 7 703 * 704 90 705 / 706 RCL IND 14 707 + 708 ST+ IND 12 709 CLX 710 SIGN 711 ST- 06 712 ST- 11 713 ST- 12 714 ST- 14 715 DSE 13 716 GTO 13 717 RCL 02 718 ST+ 40 719 VIEW 40 720 DSE 37 721 GTO 01 722 RCL 40 723 END

( 1104 bytes / SIZE 41+26.n )

 STACK INPUT OUTPUT X / t0 + N.h

Example:      The same 3 stars , let's calculate the positions 100 days after the input data below:

 Register input h =0.2 N=500 m R41=m1 R42=m2 R43=m3 3    2    1 unchanged p  o  s. R44=x1 R45=y1 R46=z1 -------- R47=x2 R48=y2 R49=z2 -------- R50=x3 R51=y3 R52=z3 0    0    0 -----    1    0    0 -----   -2    0    0 -0.113356390  1.196226886 -0.398742295  --------------  0.839927385  1.026872301 -0.342290767  -------------- -1.339785645 -5.642425267  1.880808422 v  e  l  . R53=x'1 R54=y'1 R55=z'1 -------- R56=x'2 R57=y'2 R58=z'2 -------- R59=x'3 R60=y'3 R61=z'3 0    0    0 -----    0  0.03 -0.01 -----    0 -0.06  0.02 -0.005571436 -0.001076303  0.000358768  --------------  0.004102293  0.028530254 -0.009510085  --------------  0.008509721 -0.053831599  0.017943866

Notes:

-So, we get the same accuracy with a larger stepsize.
-These values are obtained with free42.
-However, with an HP41, the roundoff-errors are of course less small.

-The maximum n-value = 10 with HP41/V41

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  J-F Lestrade - Perturbations Relativistes des Orbites Planétaires dans la métrique de Schwarzschid... - Astronomy & Astrophysics 100, 143-155 ( 1981 )
[3]  J-F Lestrade & P Bretagnon - Perturbations Relativistes pour l'ensemble des planètes - Astronomy & Astrophysics 105, 42-52 ( 1982 )
[4]  William M. Folkner, James G. Williams, Dale H. Boggs, Ryan S. Park, and Petr Kuchynka - "The Planetary and Lunar Ephemerides DE430 and DE431"
[5]  ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii
[6]  P.W. Sharp & J.M. Fine - A contrast of direct and transformed Nystrom pairs
[7]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5