Modelling lunar motion in Python
I'm trying to solve the equations of motion of the Moon in the Moon-Sun-Earth system (equations 1139-1141 here:http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node130.html)
My initial conditions were $ \boldsymbol{r} = $, where a is the Lunar semi-major axis. The initial velocities were $\boldsymbol{\dot{r}} = \omega a <0, cos\phi, sin\phi>$, where $\omega$ is the Lunar angular velocity and $\phi$ is the inclination of the Lunar orbit from the ecliptic plane.
This is my code:
```python
import matplotlib.pyplot as plt import numpy as np import scipy.integrate as spi
Data
Moon orbit semi-major axis
a = 384.748e6 # m
Moon angular velocity around the Earth (27 days)
omega = 2.6934e-6 # s^-1
omega = 0.2327 # days^-1
Earth orbit semi-major axis
a_dash = 149.6e9 # m
Earth angular velocity around the Sun (365 days)
omega_dash = 1.9923e-7 # s^-1
omega_dash = 0.0172 # days^-1
Initial conditions
at t = 0,
[x, xdot, y, ydot, z, zdot] - same as the motion vector
0.<PHONE_NUMBER> rad = 5.15 deg
vy = omegaanp.cos(0.<PHONE_NUMBER>) vz = omegaanp.sin(0.<PHONE_NUMBER>)
init = [a, 0, 0, vy, 0, vz]
Equations of motion
Returns the derivatives of x,y,z (u,v,w) and u,v,w in a vector [xdot, udot, ydot, vdot, zdot, wdot]
r - states
rdot = np.zeros(6)
def motion(t, r): x = r[0] u = r[1] # xdot or velocity in x direction
y = r[2]
v = r[3] # ydot
z = r[4]
w = r[5] # zdot
# the derivatives with respect to time
# [xdot, udot, ydot, vdot, zdot, wdot] = [u, udot, v, vdot, w, wdot]
# xdot
rdot[0] = u
# udot =
rdot[1] = 2*omega*v + (omega**2+omega_dash**2/2)*x - omega**2 * x *a**3/(np.sqrt(x**2+y**2+z**2)**3) + 3/2*omega_dash**2*(x*np.cos(2*(omega-omega_dash)*t) - y*np.sin(2*(omega-omega_dash)*t) )
# ydot
rdot[2] = v
# vdot
rdot[3] = -2*omega*u + (omega**2+omega_dash**2/2)*y - omega**2 * y *a**3/(np.sqrt(x**2+y**2+z**2)**3) - 3/2*omega_dash**2*( x*np.sin(2*(omega-omega_dash)*t) + y*np.cos(2*(omega-omega_dash)*t) )
# zdot
rdot[4] = w
# wdot
rdot[5] = -omega_dash**2*z - omega**2 * z *a**3/(np.sqrt(x**2+y**2+z**2)**3)
return rdot
Solve system of ODEs
t_start = 0. t_end = 5000 t_step = 1 ode = spi.ode(motion)
BDF method suited to stiff systems of ODEs
ode.set_integrator('vode',nsteps=500,method='bdf', rtol=0.1) ode.set_initial_value(init,t_start)
ts = [] soln = []
while ode.successful() and ode. 380