In [17]:
#%reset
import scipy.integrate as integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
def linearDerivativesFunction(inputstate, timespan):
x, y, z, xdot, ydot, zdot = inputstate
#BL1 = 3.329168
#BL1 = 4.06107
BL1 = 0.012155092
derivs = [xdot,
ydot,
zdot,
2.0*ydot + (2.0*BL1 + 1.0)*x,
-2.0*xdot - (BL1 - 1)*y,
-BL1*z]
return derivs
def nonlinearDerivativesFunction(inputstate, timespan):
x, y, z, xdot, ydot, zdot = inputstate
# distances
r1 = np.sqrt((mu+x)**2.0 + y**2.0 + z**2.0);
r2 = np.sqrt((1-mu-x)**2.0 + y**2.0 + z**2.0);
# masses
#m1 = 1 - mu;
#m2 = mu;
#G = 1;
derivs = [xdot,
ydot,
zdot,
x + 2.0*ydot + (1 - mu)*(-mu - x)/(r1**3.0) + mu*(1 - mu - x)/(r2**3.0),
y - 2.0*xdot - (1 - mu)*y/(r1**3.0) - mu*y/(r2**3.0),
-(1 - mu)*z/(r1**3.0) - mu*z/(r2**3.0)]
return derivs
# From Sharp, A Collection of Restricted Three-Body Test Problems
# For problems 1 to 15, mu = 0.012277471 and for problems 16 to 20, mu = 0.000953875
muArray_S = np.ones(20)
muArray_S[0:15] *= 0.012277471
muArray_S[15:20] *= 0.000953875
testcases_S = np.array([[0.994000E+00, 0.0, -0.21138987966945026683E+01, 0.54367954392601899690E+01], # x, z, ydot, T
[0.994000E+00, 0.0, -0.20317326295573368357E+01, 0.11124340337266085135E+02],
[0.994000E+00, 0.0, -0.20015851063790825224E+01, 0.17065216560157962559E+02],
[0.997000E+00, 0.0, -0.16251217072210773125E+01, 0.22929723423442969481E+02],
[0.879962E+00, 0.0, -0.66647197988564140807E+00, 0.63006757422352314657E+01],
[0.879962E+00, 0.0, -0.43965281709207999128E+00, 0.12729711861022426544E+02],
[0.879962E+00, 0.0, -0.38089067106386964470E+00, 0.19138746281183026809E+02],
[0.997000E+00, 0.0, -0.18445010489730401177E+01, 0.12353901248612092736E+02],
[0.100000E+01, 0.0, -0.16018768253456252603E+01, 0.12294387796695023304E+02],
[0.100300E+01, 0.0, -0.14465123738451062297E+01, 0.12267904265603897140E+02],
[0.120000E+01, 0.0, -0.71407169828407848921E+00, 0.18337451820715063383E+02],
[0.120000E+01, 0.0, -0.67985320356540547720E+00, 0.30753758552146029263E+02],
[0.120000E+01, 0.0, -0.67153130632829144331E+00, 0.43214375227857454128E+02],
[0.120000E+01, 0.0, -0.66998291305226832207E+00, 0.55672334134347612727E+02],
[0.120000E+01, 0.0, -0.66975741517271092087E+00, 0.68127906604713772763E+02],
[-0.102745E+01, 0.0, 0.40334488290490413053E-01, 0.18371316400018903965E+03],
[-0.976680E+00, 0.0, -0.61191623926410837000E-01, 0.17733241131524483004E+03],
[-0.766650E+00, 0.0, -0.51230158665978820282E+00, 0.17660722897242937108E+03],
[-0.109137E+01, 0.0, 0.14301959822238380020E+00, 0.82949461922342093092E+02],
[-0.110137E+01, 0.0, 0.15354250908611454510E+00, 0.60952121909407746612E+02]])
# From Howell, Three-Dimensional, Periodic, 'Halo' Orbits
muArray_H = [0.04]
testcases_H = np.array([[0.723268, 0.040000, 0.198019, 1.300177*2.0]]) # x, z, ydot, T
# From Barbee, Notional Mission 4 (Earth-Moon)
muArray_B = [0.012277471]
testcases_B = np.array([[0.862307159058101, 0.0, -0.187079489569182, 2.79101343456226]]) # x, z, ydot, T
# user inputs
testcase = 16
mode = 'Sharp' # 'Sharp' 'Howell' 'Barbee'
#timespan = np.linspace(0, 1.300177, 1000.0) # start, stop, numpoints # seconds
#timespan = np.linspace(0, 5.349501906, 1000.0) # start, stop, numpoints # seconds # L2 nonlinear
#timespan = np.linspace(0, 6.303856312, 1000.0) # L3 nonlinear
#stateinit = [0.723268, 0, 0.40000, 0, 0.198019, 0]
#stateinit = [120000, 0, 120000, 0, -0.151, 0] # km and km/s
#stateinit = [0.300000161, 0, 0.0, 0, -2.536145497, 0] # periodic orbit about L2 for nonlinear EOM's (from D. Eagle)
#stateinit = [-1.600000312, 0, 0, 0, 2.066174572, 0] # periodic orbit about L3 for nonlinear EOM's (from D. Eagle)
if mode == 'Sharp':
mu = muArray_S[testcase]
timespan = np.linspace(0, testcases_S[testcase, 3], 1000.0) # Sharp
stateinit = [testcases_S[testcase, 0], 0, testcases_S[testcase,1], 0, testcases_S[testcase, 2], 0] # Sharp
elif mode == 'Howell':
mu = muArray_H[testcase]
timespan = np.linspace(0, testcases_H[testcase, 3], 2000.0) # Howell
stateinit = [testcases_H[testcase, 0], 0, testcases_H[testcase,1], 0, testcases_H[testcase, 2], 0] # Howell
elif mode == 'Barbee':
mu = muArray_B[testcase]
timespan = np.linspace(0, testcases_B[testcase, 3], 2000.0) # Barbee
stateinit = [testcases_B[testcase, 0], 0, testcases_B[testcase,1], 0, testcases_B[testcase, 2], 0] # Barbee
statesOverTime = integrate.odeint(nonlinearDerivativesFunction, stateinit, timespan)
x, y, z, xdot, ydot, zdot = statesOverTime.T
#afig, ax = plt.subplots()
#ax.plot(timespan, x)
#plt.grid(True)
#plt.show()
figXY, axXY = plt.subplots()
axXY.plot(x, y)
axXY.set_aspect('equal')
figXY.suptitle('XY Plane')
figXZ, axXZ = plt.subplots()
axXZ.plot(x, z)
axXZ.set_aspect('equal')
figXZ.suptitle('XZ Plane')
figYZ, axYZ = plt.subplots()
axYZ.plot(y, z)
axYZ.set_aspect('equal')
figYZ.suptitle('YZ Plane')
#figXY, axXDotYDot = plt.subplots()
#axXDotYDot.plot(timespan, xdot)
#axXDotYDot.plot(timespan, ydot)
fig3D = plt.figure()
ax3D = fig3D.add_subplot(111, projection='3d')
ax3D.plot(x, y, z)
#ax3D.set_aspect('equal')
Out[17]:
In [17]: