Let's add an approximations to take into account some of the General Relativity (GR) effects — at least for bodies orbiting close to the massive Sun — and start to look at $J_2$ the lowest order multipole term for a body's gravitational potential beyond the monopole term $-GM/r$.
Newton:
The acceleration of a body in the gravitation field of another body of standard gravitational parameter $GM$ can be written:
$$\mathbf{a_{Newton}} = -GM \frac{\mathbf{r}}{|r|^3},$$
where $r$ is the vector from the body $M$ to the body who's acceleration is being calculated. Remember that the acceleration of each body depends only on the mass of the other body, even though the force depends on both masses, because the first mass cancels out by $a=F/m$.
General Relativity (approximate):
Although I'm not familliar with GR, I'm going to recommend an equation that seems to work well and seems to be supported by several links. It is an approximate relativistic correction to Newtonian gravity that is used in orbital mechanics simulations. You will see various forms in the following links, most with additional terms not shown here:
- https://physics.stackexchange.com/q/313146/83380
- Eq. 1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf
- Eq. 27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
- Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
- also see David Hammen's discussion/advice in this answer.
- Eq. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S (click "print") see this answer
The following approximation should be added to the Newtonian term:
$$\mathbf{a_{GR}} = GM \frac{1}{c^2 |r|^3}\left(4 GM \frac{\mathbf{r}}{|r|} - (\mathbf{v} \cdot \mathbf{v}) \mathbf{r} + 4 (\mathbf{r} \cdot \mathbf{v}) \mathbf{v} \right),$$
Oblateness ($J_2$ only, approximate):
I'm just using the math from Wikipedia's article on the Geopotential Model with a very important-to-remember approximation; I am assuming that the oblateness is in the plane of the ecliptic — that the rotational axis of the oblate body is in the $\mathbf{\hat{z}}$ direction, perpendicular to the ecliptic. Don't forget that this is an approximation! The full vector equations are messier than this, I'll try to come back and update this once I'm sure I've got it correct. In the mean time, here is an approximation:
$$\mathbf{r} = x \mathbf{\hat{x}} + y \mathbf{\hat{y}} + z \mathbf{\hat{z}} $$
$$a_x = J_2 \frac{x}{|r|^7} (6z^2 - 1.5(x^2+y^2)) $$
$$a_y = J_2 \frac{y}{|r|^7} (6z^2 - 1.5(x^2+y^2)) $$
$$a_z = J_2 \frac{z}{|r|^7} (3z^2 - 4.5(x^2+y^2)) $$
The following should be added to the Newtonian term:
$$\mathbf{a_{J2}} = a_x \mathbf{\hat{x}} + a_y \mathbf{\hat{y}} + a_z \mathbf{\hat{z}} $$
Tidal Forces:
It gets even more complicated when looking at terms that involve non-spherical mass distributions in both bodies at the same time, whether they are static or not. At this point it's probably necessary to hit the books.
Here's a test run. I've compared to downloaded data from JPL's Horizons. For the outer planets I use the Horizons data for each planet's barycenter, which makes sure it's the velocity for the center of mass of the planet plus all of its moons. I haven't added the correction to the planet's masses, but that's a much smaller effect since it only affects the movement of other, distant bodies.
For the Earth data, make sure to download the Earth's geocenter and the Moon separately (not the Earth-Moon barycenter). For the outer planets remember to download the barycenters.
I've integrated for a year, and everything is within about one kilometer of the Horizons data except for Earth's Moon. That's not a surprise considering all the intimate tidal effects between these two. Adding Earth's $J_2$ to the potential felt by the Moon only helps partially, it's really not the right way to do it, especially considering that the Earth's axis (and therefore oblateness) is so far away from the ecliptic.
So this is overall an illustration that the more physics you put in, the closer you can get to a really serious JPL simulation! This is the absolute distance between the simulated positions here and the Horizons output from 2017-01-01 00:00 to 2018-01-01 00:00. Following that is the Python script I used to generate it.
def deriv_Newton_Only(X, t):
x, v = X.reshape(2, -1)
xs, vs = x.reshape(-1, 3), v.reshape(-1, 3)
things = zip(bodies, xs, vs)
accs, vels = [], []
for a, xa, va in things:
acc_a = np.zeros(3)
for b, xb, vb in things:
if b != a:
r = xa - xb
acc_a += -b.GM * r * ((r**2).sum())**-1.5
accs.append(acc_a)
vels.append(va)
return np.hstack((np.hstack(vels), np.hstack(accs)))
def deriv_sunGRJ2EarthJ2(X, t):
x, v = X.reshape(2, -1)
xs, vs = x.reshape(-1, 3), v.reshape(-1, 3)
things = zip(bodies, xs, vs)
accs, vels = [], []
for a, xa, va in things:
acc_a = np.zeros(3)
for b, xb, vb in things:
if b != a:
r = xa - xb
acc_a += -b.GM * r * ((r**2).sum())**-1.5
if a.do_SunGR and not a.name == 'Sun':
a.flag0 = True
r = xa - xs[0]
v = va - vs[0]
rsq = (r**2).sum()
rm3 = rsq**-1.5
rm1 = rsq**-0.5
# https://physics.stackexchange.com/q/313146/83380
# Eq. 1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf
# Eq. 27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
# Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
# Eq. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S
term_0 = Sun.GM / (clight**2) * rm3
term_1 = 4.*Sun.GM * r * rm1
term_2 = -np.dot(v, v) * r
term_3 = 4.*np.dot(r, v) * v
accGR = term_0*(term_1 + term_2 + term_3)
acc_a += accGR
if a.do_SunJ2 and not a.name == 'Sun':
a.flag1 = True
r = xa - xs[0] # position relative to Sun
x, y, z = r
xsq, ysq, zsq = r**2
rsq = (r**2).sum()
rm7 = rsq**-3.5
# https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))
accJ2 = J2s * np.hstack((accJ2x, accJ2y, accJ2z))
acc_a += accJ2
if a.do_EarthJ2 and not a.name == 'Earth':
a.flag2 = True
r = xa - xs[3] # position relative to Earth
x, y, z = r
xsq, ysq, zsq = r**2
rsq = (r**2).sum()
rm7 = rsq**-3.5
# https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))
accJ2 = J2e * np.hstack((accJ2x, accJ2y, accJ2z))
acc_a += accJ2
accs.append(acc_a)
vels.append(va)
return np.hstack((np.hstack(vels), np.hstack(accs)))
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
names = ['Sun', 'Mercury', 'Venus',
'Earth', 'Moon', 'Mars',
'Ceres', 'Pallas', 'Vesta',
'Jupiter', 'Saturn', 'Uranus',
'Neptune']
GMsDE430 = [1.32712440040944E+20, 2.203178E+13, 3.24858592E+14,
3.98600435436E+14, 4.902800066E+12, 4.2828375214E+13,
6.28093938E+10, 1.3923011E+10, 1.7288009E+10,
1.267127648E+17, 3.79405852E+16, 5.7945486E+15,
6.83652719958E+15 ] # https://ipnpr.jpl.nasa.gov/progress_report/42-178/178C.pdf
# for masses also see ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck
# https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/
# https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de436s.bsp.lbl
# https://astronomy.stackexchange.com/questions/13488/where-can-i-find-visualize-planets-stars-moons-etc-positions
# https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.cmt
# ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck
GMs = GMsDE430
clight = 2.9979E+08 # m/s
halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]
# J2 values
J2_sun = 2.110608853272684E-07 # unitless
R_sun = 6.96E+08 # meters
J2s = J2_sun * (GMs[0] * R_sun**2) # is there a minus sign?
J2_earth = 1.08262545E-03 # unitless
R_earth = 6378136.3 # meters
J2e = J2_earth * (GMs[3] * R_earth**2) # is there a minus sign?
JDs, positions, velocities, linez = [], [], [], []
use_outer_barycenters = True
for name in names:
fname = name + ' horizons_results.txt'
if use_outer_barycenters:
if name in ['Jupiter', 'Saturn', 'Uranus', 'Neptune']:
fname = name + ' barycenter horizons_results.txt'
with open(fname, 'r') as infile:
lines = infile.read().splitlines()
iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]
# print name, iSOE, iEOE, lines[iSOE], lines[iEOE]
lines = lines[iSOE+1:iEOE]
lines = [line.split(',') for line in lines]
JD = np.array([float(line[0]) for line in lines])
pos = np.array([[float(item) for item in line[2:5]] for line in lines])
vel = np.array([[float(item) for item in line[5:8]] for line in lines])
JDs.append(JD)
positions.append(pos * 1000.) # km to m
velocities.append(vel * 1000.) # km/s to m/s
linez.append(lines)
JD = JDs[0] # assume they are identical
class Body(object):
def __init__(self, name):
self.name = name
bodies = []
for name, GM, pos, vel in zip(names, GMs, positions, velocities):
body = Body(name)
bodies.append(body)
body.GM = GM
body.daily_positions = pos
body.daily_velocities = vel
body.initial_position = pos[0]
body.initial_velocity = vel[0]
x0s = np.hstack([b.initial_position for b in bodies])
v0s = np.hstack([b.initial_velocity for b in bodies])
X0 = np.hstack((x0s, v0s))
ndays = 365
nperday = 144
time = np.arange(0, ndays*24*3600+1, 24*3600./nperday)
days = time[::nperday]/(24*3600.)
for body in bodies:
body.do_SunGR = False
body.do_SunJ2 = False
body.do_EarthJ2 = False
body.flag0 = False
body.flag1 = False
body.flag2 = False
Sun, Mercury, Venus, Earth, Moon, Mars = bodies[:6]
Ceres, Pallas, Vesta = bodies[6:9]
Jupiter, Saturn, Uranus, Neptune = bodies[9:]
Mercury.do_SunGR = True
Venus.do_SunGR = True
Earth.do_SunGR = True
Moon.do_SunGR = True
Mars.do_SunGR = True
Ceres.do_SunGR = True
Pallas.do_SunGR = True
Vesta.do_SunGR = True
Mercury.do_SunJ2 = True
Moon.do_EarthJ2 = True
rtol = 1E-12 # this is important!!!
answer, info = ODEint(deriv_sunGRJ2EarthJ2, X0, time,
rtol = rtol, full_output=True)
print answer.shape
nbodies = len(bodies)
xs, vs = answer.T.reshape(2, nbodies, 3, -1)
for body, x, v in zip(bodies, xs, vs):
body.x = x
body.v = v
body.x_daily = body.x[:, ::nperday]
body.v_daily = body.v[:, ::nperday]
body.difference = np.sqrt(((body.x_daily - body.daily_positions.T)**2).sum(axis=0))
if True:
for body in bodies[:6]:
print body.name, body.flag0, body.flag1, body.flag2
if True:
plt.figure()
for i, body in enumerate(bodies[:12]): # Sorry Neptune!!!
plt.subplot(4, 3, i+1)
plt.plot(days, 0.001*body.difference)
plt.title(body.name, fontsize=14)
plt.xlim(0, 365)
plt.suptitle("calc vs JPL Horizons (km vs days)", fontsize=16)
plt.show()