! Solar System (DE200) + Ceres: Barycentric
Coordinates
! General Relativity and Earth-Moon Structure
Corrections
! Fourth Order Symmetric Method (Runge-Kutta start).
! Ian Gatland, Georgia Institute of Technology
! Initiated 2002/04/01, last modified 2004/07/28.
!
program SolSys33
option nolet
title$ = "SolSys33"
fname$ = "SolSys33r"
open #1: name fname$, access outin, create newold,
org text
erase #1
print #1
print #1: " " & title$
print
print " " & title$
public np, m(0), g, c, q, sm, tm
dim r(0), v(0), u0(0), u1(0), u2(0), a0(0), a1(0),
a2(0)
np = 11
n = np*3
mat redim m(1:np), r(1:n), v(1:n)
mat redim u0(1:n), u1(1:n), u2(1:n), a0(1:n),
a1(1:n), a2(1:n)
call Init(nr, nt, dt, t, r(), v())
call CoM(r())
call CoM(v())
call Report(#1, 1, dt, t, r(), v())
call Prep(dt, t, n, r(), v(), u0(), u1(), u2(),
a0(), a1(), a2())
for k = 1 to nr
print
print " nt = "; nt
print " t = "; t + nt*dt
call Sym4(nt, dt, t, n, r(), v(), u0(), u1(),
u2(), a0(), a1(), a2())
call Report(#1, 0, dt, t, r(), v())
next k
close #1
print
print " Done"
end
!
sub Init(nr, nt, dt, t, r(), v())
declare public np, m(), g, c, q, sm, tm
m(1) = 1/6023600 ! mass of Mercury
m(2) = 1/408523.5 ! mass of Venus
m(3) = 1/328900.55 ! mass of Earth+Moon
m(4) = 1/3098710 ! mass of Mars
m(5) = 1/1700000000 ! mass of Ceres
m(6) = 1/1047.350 ! mass of Jupiter
m(7) = 1/3498.0 ! mass of Saturn
m(8) = 1/22960 ! mass of Uranus
m(9) = 1/19314 ! mass of Neptune
m(10) = 1/130000000 ! mass of Pluto
m(np) = 1 ! mass of Sun
g = 0.01720209895^2 ! gravitational constant
c = 173.144632721 ! speed of light
q = 5.46e-8 ! Earth+Moon quadrupole
factor
sm = 100 ! minimum separation
distance
tm = 0 ! minimum separation time
nr = 7 ! number of reports
nt = 13440 ! number of time steps
dt = 0.0625 ! time step
t = 2440400.5 ! initial time
mat read r ! initial positions, x1,
y1, z1, x2, y2, etc.
data 3.5726021254696e-01, -9.1549055285616e-02,
-8.5981004134536e-02
data 6.0824943776644e-01, -3.4913244404770e-01,
-1.9554432558022e-01
data 1.1601491704454e-01, -9.2660555805310e-01,
-4.0180626511782e-01
data -1.1468856546204e-01, -1.3283665333858e+00,
-6.0615518746928e-01
data 1.4380476396912e+00, -2.2047841358670e+00,
-1.3264690976797e+00
data -5.3842086414064e+00, -8.3124999735360e-01,
-2.2509802926003e-01
data 7.8898894267323e+00, 4.5957099267226e+00,
1.5584291663445e+00
data -1.8269891137986e+01, -1.1627330499135e+00,
-2.5037650434585e-01
data -1.6059504334173e+01, -2.3942941306015e+01,
-9.4004277295751e+00
data -3.0487996972540e+01, -8.7321653623023e-01,
8.9113520872503e+00
data 0.0000000000000e+00, 0.0000000000000e+00,
0.0000000000000e+00
mat read v ! initial velocities, vx1,
vy1, vz1, vx2, etc.
data 3.3678452045578e-03, 2.4889342837586e-02,
1.2944071597159e-02
data 1.0952419935474e-02, 1.5612506911548e-02,
6.3288764369226e-03
data 1.6811620039589e-02, 1.7431312618369e-03,
7.5597507976519e-04
data 1.4482004836578e-02, 2.3728517456873e-04,
-2.8374875686161e-04
data 8.4663246357757e-03, 4.6824117231476e-03,
4.6509329975024e-04
data 1.0923674506708e-03, -6.5232939031698e-03,
-2.8230121107231e-03
data -3.2172051412201e-03, 4.3306320894907e-03,
1.9264168192697e-03
data 2.2154446129588e-04, -3.7676549166365e-03,
-1.6532438908973e-03
data 2.6431259526341e-03, -1.5034868645846e-03,
-6.8126855659202e-04
data 3.2254176879840e-04, -3.1487599655419e-03,
-1.0801855125339e-03
data 0.0000000000000e+00, 0.0000000000000e+00,
0.0000000000000e+00
end sub
!
sub CoM(r())
declare public np, m(), g, c, q, sm, tm
mt = 0
x = 0
y = 0
z = 0
for j = 1 to np
mt = mt + m(j)
i = 3*(j-1)
x = x + m(j)*r(i+1)
y = y + m(j)*r(i+2)
z = z + m(j)*r(i+3)
next j
for j = 1 to np
i = 3*(j-1)
r(i+1) = r(i+1) - x/mt
r(i+2) = r(i+2) - y/mt
r(i+3) = r(i+3) - z/mt
next j
end sub
!
sub Acc(t, r(), v(), a())
declare public np, m(), g, c, q, sm, tm
n = np*3
for i = 1 to n
a(i) = 0
next i
for j1 = 1 to np-1
i1 = 3*(j1-1)
for j2 = j1+1 to np
i2 = 3*(j2-1)
x = r(i1+1) - r(i2+1)
y = r(i1+2) - r(i2+2)
z = r(i1+3) - r(i2+3)
s = sqr(x*x + y*y + z*z)
if s < sm then
sm = s
tm = t
end if
f = -g/(s*s*s)
if j2 <> np then
a(i1+1) = a(i1+1) + m(j2)*f*x
a(i1+2) = a(i1+2) + m(j2)*f*y
a(i1+3) = a(i1+3) + m(j2)*f*z
a(i2+1) = a(i2+1) - m(j1)*f*x
a(i2+2) = a(i2+2) - m(j1)*f*y
a(i2+3) = a(i2+3) - m(j1)*f*z
else
if j1 = 3 then f = f*(1 + q/(s*s))
vx = v(i1+1) - v(i2+1)
vy = v(i1+2) - v(i2+2)
vz = v(i1+3) - v(i2+3)
rfac = 1 + ((vx*vx + vy*vy + vz*vz) -
4*g*m(np)/s)/(c*c)
vfac = -4*(x*vx + y*vy + z*vz)/(c*c)
a(i1+1) = a(i1+1) + m(j2)*f*(x*rfac +
vx*vfac)
a(i1+2) = a(i1+2) + m(j2)*f*(y*rfac +
vy*vfac)
a(i1+3) = a(i1+3) + m(j2)*f*(z*rfac + vz*vfac)
a(i2+1) = a(i2+1) - m(j1)*f*(x*rfac + vx*vfac)
a(i2+2) = a(i2+2) - m(j1)*f*(y*rfac + vy*vfac)
a(i2+3) = a(i2+3) - m(j1)*f*(z*rfac + vz*vfac)
end if
next j2
next j1
end sub
!
sub Sym4(nt, dt, t, n, r(), v(), u0(), u1(), u2(),
a0(), a1(), a2())
if nt <= 0 then exit sub
nt3 = int(nt/3)
for j = 1 to nt3
t = t + dt
for i = 1 to n
u2(i) = u2(i) + (5*a0(i) + 2*a1(i) +
5*a2(i))*dt/4
r(i) = r(i) + u2(i)*dt
v(i) = u2(i) + (27*a0(i) - 22*a1(i) +
7*a2(i))*dt/24
next i
call Acc(t, r(), v(), a2())
t = t + dt
for i = 1 to n
u1(i) = u1(i) + (5*a2(i) + 2*a0(i) +
5*a1(i))*dt/4
r(i) = r(i) + u1(i)*dt
v(i) = u1(i) + (27*a2(i) - 22*a0(i) +
7*a1(i))*dt/24
next i
call Acc(t, r(), v(), a1())
t = t + dt
for i = 1 to n
u0(i) = u0(i) + (5*a1(i) + 2*a2(i) + 5*a0(i))*dt/4
r(i) = r(i) + u0(i)*dt
v(i) = u0(i) + (27*a1(i) - 22*a2(i) + 7*a0(i))*dt/24
next i
call Acc(t, r(), v(), a0())
next j
for j = 3*nt3+1 to nt
t = t + dt
for i = 1 to n
u2(i) = u2(i) + (5*a0(i) + 2*a1(i) +
5*a2(i))*dt/4
r(i) = r(i) + u2(i)*dt
v(i) = u2(i) + (27*a0(i) - 22*a1(i) +
7*a2(i))*dt/24
next i
call Acc(t, r(), v(), a2())
for i = 1 to n
u = u2(i)
u2(i) = u1(i)
u1(i) = u0(i)
u0(i) = u
a = a2(i)
a2(i) = a1(i)
a1(i) = a0(i)
a0(i) = a
next i
next j
for i = 1 to n
v(i) = u0(i) + (7*a0(i) + 6*a1(i) - a2(i))*dt/24
next i
end sub
!
sub Prep(dt, t, n, r(), v(), u0(), u1(), u2(), a0(),
a1(), a2())
dim rc(0), vc(0)
mat redim rc(1:n), vc(1:n)
tc = t
for i = 1 to n
rc(i) = r(i)
vc(i) = v(i)
next i
call Acc(tc, rc(), vc(), a0())
call RKstep(-dt, tc, n, rc(), vc(), u0())
call Acc(tc, rc(), vc(), a1())
call RKstep(-dt, tc, n, rc(), vc(), u1())
call Acc(tc, rc(), vc(), a2())
call RKstep(-dt, tc, n, rc(), vc(), u2())
end sub
!
sub RKstep(dt, t, n, r(), v(), u())
dim r1(0), r2(0), r3(0), r4(0)
dim v1(0), v2(0), v3(0), v4(0)
dim a1(0), a2(0), a3(0), a4(0)
mat redim r1(1:n), r2(1:n), r3(1:n), r4(1:n)
mat redim v1(1:n), v2(1:n), v3(1:n), v4(1:n)
mat redim a1(1:n), a2(1:n), a3(1:n), a4(1:n)
t1 = t
for i = 1 to n
r1(i) = r(i)
v1(i) = v(i)
next i
call Acc(t1, r1(), v1(), a1())
t2 = t + dt/2
for i = 1 to n
r2(i) = r(i) + v1(i)*dt/2
v2(i) = v(i) + a1(i)*dt/2
next i
call Acc(t2, r2(), v2(), a2())
t3 = t + dt/2
for i = 1 to n
r3(i) = r(i) + v2(i)*dt/2
v3(i) = v(i) + a2(i)*dt/2
next i
call Acc(t3, r3(), v3(), a3())
t4 = t + dt
for i = 1 to n
r4(i) = r(i) + v3(i)*dt
v4(i) = v(i) + a3(i)*dt
next i
call Acc(t4, r4(), v4(), a4())
t = t + dt
for i = 1 to n
u(i) = (v1(i) + 2*v2(i) + 2*v3(i) + v4(i))/6
r(i) = r(i) + u(i)*dt
v(i) = v(i) + (a1(i) + 2*a2(i) + 2*a3(i) +
a4(i))*dt/6
next i
end sub
!
sub Report(#1, k, dt, t, r(), v())
declare public np, m(), g, c, q, sm, tm
n = np*3
o$ = " -%.%%%%%%%%%%%%%^^^^"
if k = 1 then
print #1
print #1, using o$: "g" , "c" , "q"
print #1, using o$: g , c , q
print #1, using o$: "dt"
print #1, using o$: dt
print #1, using o$: "m"
for j = 1 to np
print #1, using o$: m(j);
if mod(j,3) = 0 or j = np then print #1
next j
print
print using o$: "g", "c", "q"
print using o$: g , c , q
print using o$: "dt"
print using o$: dt
print using o$: "m"
for j = 1 to np
print using o$: m(j);
if mod(j,3) = 0 or j = np then print
next j
end if
print #1
print #1, using o$: "t", "sm", "tm"
print #1, using o$: t , sm , tm
print #1, using o$: "rx", "ry", "rz"
for j = 1 to np
i = 3*(j-1)
print #1, using o$: r(i+1), r(i+2), r(i+3)
! print #1, using o$: r(i+1)-r(n-2), r(i+2)-r(n-1),
r(i+3)-r(n)
next j
print #1, using o$: "vx", "vy", "vz"
for j = 1 to np
i = 3*(j-1)
print #1, using o$: v(i+1), v(i+2), v(i+3)
! print #1, using o$: v(i+1)-v(n-2), v(i+2)-v(n-1),
v(i+3)-v(n)
next j
print
print using o$: "t", "sm", "tm"
print using o$: t , sm , tm
print using o$: "rx", "ry", "rz"
for j = 1 to np
i = 3*(j-1)
print using o$: r(i+1), r(i+2), r(i+3)
! print using o$: r(i+1)-r(n-2), r(i+2)-r(n-1),
r(i+3)-r(n)
next j
print using o$: "vx", "vy", "vz"
for j = 1 to np
i = 3*(j-1)
print using o$: v(i+1), v(i+2), v(i+3)
! print using o$: v(i+1)-v(n-2), v(i+2)-v(n-1),
v(i+3)-v(n)
next j
end sub