元の 2 次 ODE は次のとおりです。
x'' - 2 * omega * y' - omega ** 2 * x = - mue * (x + pi2 * r12) / np.sqrt((x + pi2 * r12) ** 2 + y ** 2) ** 3 - mum * (x - pi1 * r12) / np.sqrt((x - pi1 * r12) ** 2 + y ** 2)
y'' + 2 * omega * x' - omega **2 * y = - mue * y / np.sqrt((x + pi2 * r12) ** 2 + y ** 2) ** 3 - mum * y / np.sqrt((x - pi1 * r12) ** 2 + y ** 2)
z'' = 0
ODE を解くために使用したコードは次のとおりですが、最初に 2 つの最初の順序に分割しました。
61 行目のモジュールが呼び出せないというエラーが表示されます。
61行目はu = odeint(deriv, u0, dt)
#!/usr/bin/env python
import numpy as np
import scipy.integrate as odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
me = 5.974 * 10 ** (24) # mass of the earth
mm = 7.348 * 10 ** (22) # mass of the moon
G = 6.67259 * 10 ** (-20) # gravitational parameter
re = 6378.0 # radius of the earth in km
rm = 1737.0 # radius of the moon in km
r12 = 384400.0 # distance between the CoM of the earth and moon
M = me + mm
pi1 = me / M
pi2 = mm / M
mue = 398600.0 # gravitational parameter of earth km^3/sec^2
mum = G * mm # grav param of the moon
mu = mue + mum
omega = np.sqrt(mu / r12 ** 3)
nu = 0.0 # flight path angle
x = 327156.0 # x location where the moon's SOI effects the spacecraft
y = 33050.0 # y location
vbo = 10.85 # velocity at burnout
gamma = -141.868 * np.pi / 180 # angle in radians of true anomaly
vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
# velocity of the bo in the x direction
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
# velocity of the bo in the y direction
xrel = (re + 300.0) * np.cos(gamma)
# spacecraft x location relative to the earth
yrel = (re + 300.0) * np.sin(gamma)
# r0 = [xrel, yrel, 0]
# v0 = [vx, vy, 0]
u0 = [xrel, yrel, 0, vx, vy, 0]
def deriv(u, dt):
n1 = -((mue * (u[0] + pi2 * r12) / np.sqrt((u[0] + pi2 * r12) ** 2
+ u[1] ** 2) ** 3)
- (mum * (u[0] - pi1 * r12) / np.sqrt((u[0] - pi1 * r12) ** 2
+ u[1] ** 2) ** 3))
n2 = -((mue * u[1] / np.sqrt((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3)
- (mum * u[1] / np.sqrt((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3))
return [u[3], # dotu[0] = u[3]
u[4], # dotu[1] = u[4]
u[5], # dotu[2] = u[5]
2 * omega * u[5] + omega ** 2 * u[0] + n1, # dotu[3] = that
omega ** 2 * u[1] - 2 * omega * u[4] + n2, # dotu[4] = that
0] # dotu[5] = 0
dt = np.arange(0.0, 250000.0, .1)
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = u.T
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z)
plt.show()