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) + k₁/6 + 2k₂/6 )
y'(t+h) = y'(t) + k₁/6 + 2k₂/3 + k₃/6 where k₁ = h.f (y) ; k₂ = h.f(y+h.y'/2+h.k₁/8) ; k₃ = h.f(y+h.y'+h.k₂/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 ):

              d²x_i /dt² = SUM_j#i   Gm_j ( x_j - x_i ) / [ ( x_i - x_j )²+ ( y_i - y_j )² + ( z_i - z_j )² ]^3/2
              d²y_i /dt² = SUM_j#i   Gm_j ( y_j - y_i ) / [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^3/2          i , j = 1 , ..... , n
              d²z_i /dt² = SUM_j#i   Gm_j ( z_j - z_i ) / [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^3/2

G = gravitational constant = 0.01720209895² = 1 / 3379.380681 m_i = 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 = t₀ = initial instant                                    R04 to R14: temp
                                      • R02 = h = step size ( in days )
                                      • R03 = N = nb of steps                       • R16 = m₁ = mass of the 1st body
                                                                                                       ................................................

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

                                      • R16+n = x₁              • R16+4n = x'₁
                                      • R17+n = y₁              • R17+4n = y'₁                       R16+7n to R15+10n = k₂ , k₃
                                      • R18+n = z₁              • R18+4n = z₁                       R16+10n to R15+13n = k₁ , k₁ - k₂
                                         .................                   ..................

                                      • R13+4n = x_n            • R13+7n = x'_n
                                      • R14+4n = y_n            • R14+7n = y'_n
                                      • R15+4n = z_n            • 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	/	t₀ + N.h

Example:

-Let's imagine a system of three stars      M₁( 2 ; 0 ;0 )       M₂ ( 0 ; 4 ; 0 )      M₃ ( 0 ; 0 ; 1 )                     unit = 1 AU
              with initial velocities                 V₁ ( 0 ; 0.03 ; 0 ) V₂ ( 0 ; 0 ; 0.01 ) V₃ ( -0.02 ; 0 ; 0 )               unit = 1 AU/day
                     and masses                        m₁ = 2 ; m₂ = 1 ; m₃ = 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:

   k₁ = h.f(x₀,y₀)
   k₂ = h.f(x₀+ h/4 , y₀+ (h/4) y'₀+ (h/32) k₁)
   k₃ = h.f(x₀+ h/2 , y₀+ (h/2) y'₀+ h (-k₁/24 + k₂/6 )
   k₄ = h.f(x₀+ 3h/4 , y₀+ (3h/4) y'₀+ h (3k₁/32 + k₂/8 + k₃/16 )
   k₅ = h.f(x₀+ h , y₀+ h y'₀+ h (3k₂/7 - k₃/14 + k₄/7 )

-Then, y = y₀ + h y'₀+ h² ( 7 k₁/90 + 4 f₂/15 + f₃/15 + 4 f₄/45 ) + O(h⁷)
and y' = y'₀ + h ( 7 f₁/90 + 16 f₂/45 + 2 f₃/15 + 16 f₄/45 + 7 f₅/90 ) + O(h⁷)

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 = t₀ = initial instant                            R04 to R18: temp     R19: unused
                                      • R02 = h = step size ( in days )
                                      • R03 = N = nb of steps                       • R21 = m₁ = mass of the 1st body
                                                                                                       ................................................

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

                                      • R21+n = x₁              • R21+4n = x'₁                       R21+7n to R20+10n = k₄
                                      • R22+n = y₁              • R22+4n = y'₁                      R21+10n to R20+13n = k₃
                                      • R23+n = z₁              • R23+4n = z₁                       R21+13n to R20+16n = k₂ , 7.k₁/90-17k₂/105
                                         .................                   ..................                         R21+16n to R20+19n = k₁ , k₅

                                      • R18+4n = x_n            • R18+7n = x'_n
                                      • R19+4n = y_n            • R19+7n = y'_n
                                      • R20+4n = z_n            • 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	/	t₀ + N.h

Example: The same used by "GM4":

-A system of three stars              M₁( 2 ; 0 ;0 )       M₂ ( 0 ; 4 ; 0 )      M₃ ( 0 ; 0 ; 1 )                     unit = 1 AU
with initial velocities                 V₁ ( 0 ; 0.03 ; 0 ) V₂ ( 0 ; 0 ; 0.01 ) V₃ ( -0.02 ; 0 ; 0 )               unit = 1 AU/day
        and masses                        m₁ = 2 ; m₂ = 1 ; m₃ = 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:

k₁ = h.f (y) y(t+h) = y(t) + h [ y' + ( k₁ + k₂ )/4 ]
k₂ = h.f (y+2.h.y'/3+2.h.k₁/9) y'(t+h) = y'(t) + k₁/4 + 3.k₂/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 = t₀ = initial instant                            R04 to R16: temp
                                      • R02 = h = step size ( in days )
                                      • R03 = N = nb of steps                       • R18 = m₁ = mass of the 1st body
                                                                                                       ................................................

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

                                      • R18+n = x₁              • R18+4n = x'₁                       R18+7n to R17+10n = k₁
                                      • R19+n = y₁              • R19+4n = y'₁
                                      • R20+n = z₁              • R20+4n = z₁
                                         .................                   ..................

                                      • R15+4n = x_n            • R15+7n = x'_n
                                      • R16+4n = y_n            • R16+7n = y'_n
                                      • R17+4n = z_n            • 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	/	t₀ + 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

t₀ = 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) + k₁/6 + 2k₂/6 )
y'(t+h) = y'(t) + k₁/6 + 2k₂/3 + k₃/6 where k₁ = h.f (y) ; k₂ = h.f(y+h.y'/2+h.k₁/8) ; k₃ = h.f(y+h.y'+h.k₂/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 ):

               d²x_i /dt² = - Gm₀ x_i/r_i³ - SUM_j Gm_j x_j/r_j³ + SUM_j#i   Gm_j ( x_j - x_i ) / r_ij³
               d²y_i /dt² = - Gm₀ y_i/r_i³ - SUM_j Gm_j y_j/r_j³ + SUM_j#i   Gm_j ( y_j - y_i ) / r_ij³                  i , j = 1 , ..... , n
               d²z_i /dt² = - Gm₀ z_i/r_i³ - SUM_j Gm_j z_j/r_j³ + SUM_j#i   Gm_j ( z_j - z_i ) / r_ij³

where r_i = ( x_i² + y_i² + z_i² )^1/2 and r_ij = [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^1/2

G = gravitational constant = 0.01720209895² = 1 / 3379.380681

m₀ = mass of the Sun
m_i = 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 = t₀ = initial instant                                                   R04 to R17: temp
                                      • R02 = h = step size ( in days )            • R19 = m₀ = mass of the Sun ( 1 for the Solar system )
                                      • R03 = N = nb of steps                       • R20 = m₁ = mass of the 1st planet
                                                                                                       ................................................

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

                                      • R20+n = x₁              • R20+4n = x'₁
                                      • R21+n = y₁              • R21+4n = y'₁                       R20+7n to R19+10n = k₂ , k₃
                                      • R22+n = z₁              • R22+4n = z₁                       R20+10n to R19+13n = k₁ , k₁ - k₂
                                         .................                   ..................

                                      • R17+4n = x_n            • R17+7n = x'_n
                                      • R18+4n = y_n            • R18+7n = y'_n
                                      • R19+4n = z_n            • 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	/	t₀ + 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 = t₀ = initial instant                                                   R04 to R19: temp
                                      • R02 = h = step size ( in days )            • R21 = m₀ = mass of the Sun ( 1 for the Solar system )
                                      • R03 = N = nb of steps                       • R22 = m₁ = mass of the 1st planet
                                                                                                       ................................................

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

                                      • R22+n = x₁              • R22+4n = x'₁
                                      • R23+n = y₁              • R23+4n = y'₁                       R22+7n to R21+10n = k₄
                                      • R24+n = z₁              • R24+4n = z₁                       R22+10n to R21+13n = k₃
                                         .................                   ..................                         R22+13n to R21+16n = k₂ , 7.k₁/90-17k₂/105

                                      • R19+4n = x_n            • R19+7n = x'_n                      R22+17n to R21+19n = k₁ , k₅
                                      • R20+4n = y_n            • R20+7n = y'_n
                                      • R21+4n = z_n            • 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	/	t₀ + 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]

dr_i'' = ( GM/r_i³ ) [ 4 ( GM/r_i ) r_i - v_i²r_i + 4 ( r_i.v_i ) v_i ] / c² c = speed of light = 173.446327 AU/day

where r_i = ( x_i , y_i , z_i ) = position-vector of the planet i r_i = distance Sun-planet i M = mass of the Sun = m₀ = 1 for our Solar system
v_i = ( x'_i , y'_i , z'_i ) = velocity-vector of the planet i v_i = norm of v_i

-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:

   k₁ = h.f (y,y')
   k₂ = h.f( y+h.y'/2+h.k₁/8 , y'+k₁/2 )
   k₃ = h.f( y+h.y'/2+h.k₁/8 , y'+k₂/2)
   k₄ = h.f( y+h.y'+h.k₃/2 , y'+k₃)

y(t+h) = y(t) + h [ y'(t) + ( k₁ + k₂ + k₃ )/6 ]
y'(t+h) = y'(t) + ( k₁ + 2.k₂ + 2.k₃ + k₄ ) / 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 = t₀ = initial instant                  R04 to R27: temp            R28: unused
                                      • R02 = h = step size ( in days )            • R30 = m₀ = mass of the Sun ( 1 for the Solar system )
                                      • R03 = N = nb of steps                       • R31 = m₁ = mass of the 1st planet
                                                                                                       ................................................

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

                                      • R31+n = x₁              • R31+4n = x'₁                      R31+7n to R30+10n = k₅ and linear combinations of other k_i
                                      • R32+n = y₁              • R32+4n = y'₁                      R31+10n to R30+13n = k₂ , k₂ - 2.k₃
                                      • R33+n = z₁              • R33+4n = z₁                       R31+13n to R30+16n = k₁
                                         .................                   ..................

                                      • R28+4n = x_n            • R28+7n = x'_n
                                      • R29+4n = y_n            • R29+7n = y'_n
                                      • R30+4n = z_n            • 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	/	t₀ + 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:

   k₁ = h f( y , y' )
   k₂ = h f( y+h.( (2/5).y'+(2/25).k₁ ), y'+(2/5) k₁ )
   k₃ = h f( y+h.( (1/4).y'+(11/256).k₁-(3/256).k₂ ), y'+(11/64) k₁+(5/64).k₂ )
   k₄ = h f( y+h.( (1/2).y'+(1/8).k₃ ), y'+(1/2) k₃ )
   k₅ = h f( y+h.( (3/4).y'+(9/256).k₁-(21/256).k₂+(3/16).k₃+(9/64).k₄ ), y'+(3/64) k₁-(15/64).k₂+(3/8).k₃+(9/16).k₄ )
   k₆ = h f( y+h.( y'+(3/7).k₂+(9/14).k₃-(6/7).k₄+(2/7).k₅ ), y'+(5/7) k₂+(6/7).k₃-(12/7).k₄+(8/7).k₅ )

y(x+h) - y(x) = h.[ y' + ( 7 k₁ + 24 k₃ + 6 k₄ + 8 k₅ ) / 90 ]
y'(x+h) - y'(x) = ( 7 k₁ + 32 k₃ + 12 k₄ + 32 k₅ + 7 k₆ ) / 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 = t₀ = initial instant                                                     R04 to R28: temp
                                      • R02 = h = step size ( in days )            • R30 = m₀ = mass of the Sun ( 1 for the Solar system )
                                      • R03 = N = nb of steps                       • R31 = m₁ = mass of the 1st planet
                                                                                                       ................................................

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

                                      • R31+n = x₁              • R31+4n = x'₁                      R31+7n to R30+10n = k₅ and linear combination of other k_i
                                      • R32+n = y₁              • R32+4n = y'₁                      R31+10n to R30+13n = k₄ and linear combination of other k_i
                                      • R33+n = z₁              • R33+4n = z₁                       R31+13n to R30+16n = k₃
                                         .................                   ..................                        R31+16n to R30+19n = k₂
                                                                                                                       R31+19n to R30+22n = k₁ , k₆
                                      • R28+4n = x_n            • R28+7n = x'_n
                                      • R29+4n = y_n            • R29+7n = y'_n
                                      • R30+4n = z_n            • 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	/	t₀ + 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]:

r_i'' = Sum_j#i ( Gm_j/r_ij³ ) { 1 - ( 4/c² ) Sum_k#i Gm_k/r_ik - ( 1/c² ) Sum_k#j Gm_k/r_jk + ( v_i² / c² ) + 2 ( v_j² / c² ) + ( 4/c² ) v_i.v_j
- ( 3/c² ) [ ( r_ij.v_j ) / r_ij ]² + ( 1/c² ) r_ij . r_j" }
+ ( 1/c² ) Sum_j#i ( Gm_j/r_ij³ ) [ r_ij. ( 4.v_i - 3.v_j ) ] v_ij + ( 7/(2.c²) ) Sum_j#i ( Gm_j/r_ij ) r_j"

-Since the acceleration r_i'' appears in the right-hand side of the formula - but multiplied by ( 1/c² ) - 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 = t₀ = initial instant
                                      • R03 = N = nb of steps                       • R41 = m₁ = mass of the 1st celestial body
                                                                                                       ................................................

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

                                      • R41+n = x₁              • R41+4n = x'₁                      R41+7n to R40+10n = k₃ , k₄
                                      • R42+n = y₁              • R42+4n = y'₁                      R41+10n to R40+13n = k₂ , 2.k₃ - k₂
                                      • R43+n = z₁              • R43+4n = z₁                       R41+13n to R40+16n = k₁
                                         .................                   ..................
                                                                                                                       R41+16.n to R40+19.n = r_i'' ( Newtonian )
                                      • R38+4n = x_n            • R38+7n = x'_n                     R41+19.n to R40+20.n = Sum_j#i Gm_j/r_ij
                                      • R39+4n = y_n            • R39+7n = y'_n
                                      • R40+4n = z_n            • 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	/	t₀ + 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 = t₀ = initial instant
                                      • R03 = N = nb of steps                       • R41 = m₁ = mass of the 1st celestial body
                                                                                                       ................................................

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

                                      • R41+n = x₁              • R41+4n = x'₁
                                      • R42+n = y₁              • R42+4n = y'₁                      R41+7n to R40+10n = k₅ and linear combination of other k_i
                                      • R43+n = z₁              • R43+4n = z₁                       R41+10n to R40+13n = k₄ and linear combination of other k_i
                                         .................                   ..................                        R41+13n to R40+16n = k₃
                                                                                                                       R41+16n to R40+19n = k₂
                                      • R38+4n = x_n            • R38+7n = x'_n                     R41+19n to R40+22n = k₁ , k₆
                                      • R39+4n = y_n            • R39+7n = y'_n                     R41+22.n to R40+25.n = r_i'' ( Newtonian )
                                      • R40+4n = z_n            • R40+7n = z'_n                     R41+25.n to R40+26.n = Sum_j#i Gm_j/r_ij
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	/	t₀ + 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