def deriv(X, t):
x, v = X.reshape(2, -1)
xx = x.reshape(-1, 3)
n = xx.shape[0]
accs = []
for i in range(n):
acc = np.zeros(3)
for j in range(4):
if j != i:
xxij = xx[i] - xx[j]
acc += -GM4[j] * xxij * ((xxij**2).sum())**-1.5
accs.append(acc)
accs = np.hstack(accs)
return np.hstack((v, accs))
def rotatem(X, theta):
cth, sth = [f(theta) for f in (np.cos, np.sin)]
x, y, z = X
xr = cth*x - sth*y
yr = cth*y + sth*x
return np.vstack((xr, yr, z))
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint as ODEint
GMe = 3.9860E+14 # m/s
GMm = 4.9049E+12
GMs = 1.3271E+20
GMj = 1.2669E+17
GM4 = GMe, GMm, GMs, GMj
R4 = 6378137., 1738100., 696392000., 71492000. # m
names = 'Earth', 'Moon', 'Sun', 'Jupiter'
Re, Rm, Rs, Rj = R4
vescs = [np.sqrt(2.*GM/R) for (GM, R) in zip(GM4, R4)]
vesce = vescs[0]
for name, vesc in zip(names, vescs):
print (name, vesc)
X0e = 1000. * np.array([1.3539E+07, -1.5044E+08, -7.7480E+03,
2.9164E+01, 2.5290E+00, 8.3979E-04])
X0m = 1000. * np.array([1.3481E+07, -1.5083E+08, 1.7757E+04,
3.0127E+01, 2.3628E+00, -6.1480E-02])
X0s = 1000. * np.array([9.8486E+04, 1.0333E+06, -1.3866E+04,
-1.2308E-02, 6.41628E-03, 3.0230E-04])
X0j = 1000. * np.array([-4.9851E+08, -6.3418E+08, 1.3781E+07,
1.0118E+01, -7.4520E+00, -1.9535E-01])
X0x = np.hstack([x[:3] for x in (X0e, X0m, X0s, X0j)])
X0v = np.hstack([x[3:] for x in (X0e, X0m, X0s, X0j)])
d = np.array(((-1, 0, 0), (1, 0, 0),
(0, -1, 0), (0, 1, 0),
(0, 0, -1), (0, 0, 1)), dtype=float)
xobs = (d*Re + X0e[:3]).flatten()
factor = 1.02
vobs = (d*vesce*factor + X0e[3:]).flatten()
X0x = np.hstack((X0x, xobs))
X0v = np.hstack((X0v, vobs))
X0 = np.hstack((X0x, X0v))
rs = np.sqrt((X0x.reshape(-1, 3)**2).sum(axis=0))
for name, vesc, r in zip(names, vescs, rs):
print (name, vesc, r)
times = np.arange(0, 365*24*3600, 10000)
answer, info = ODEint(deriv, X0, times, full_output=True)
n = answer.shape[0]
xall, vall = answer.T.reshape(2, -1, 3, n)
xe, ve = [thing[0] for thing in (xall, vall)]
xps, vpe = [thing[4:] for thing in (xall, vall)]
theta = np.arctan2(xe[1], xe[0])
xer = rotatem(xe, -theta)
xpsr = np.stack([rotatem(thing, -theta) for thing in xps])
if True:
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
w = 1.5E+08
x, y, z = 1E-03 * xe
ax1.plot(x, y, z, '-b', linewidth=1)
ax1.plot(x[:1], y[:1], z[:1], 'ok')
for x, y, z in 1E-03 * xps:
ax1.plot(x, y, z, linewidth=0.5)
ax1.set_xlim(-w, w)
ax1.set_ylim(-w, w)
ax1.set_zlim(-w, w)
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
w = 1.5E+07
x, y, z = 1E-03 * (xer-xer)
ax2.plot(x, y, z, '-b', linewidth=1)
ax2.plot(x[:1], y[:1], z[:1], 'ok')
for x, y, z in 1E-03 * (xpsr-xer):
ax2.plot(x, y, z, linewidth=0.5)
ax2.set_xlim(-w, w)
ax2.set_ylim(-w, w)
ax2.set_zlim(-w, w)
plt.show()
