I am trying to numerically solve the orbit of a space probe around the moon in a non-spherical gravitational field .
$ \Large \begin{align} \frac{d^{2} x}{dt^2} = - \frac{G M x}{(x^2 + y^2)^{3/2}} - \frac{q\ G M x^\prime}{(x^{^\prime2} + y^{^\prime2})^{3/2}} \end{align} $
$ \Large \begin{align} \frac{d^{2} y}{dt^2} = - \frac{G M y}{(x^2 + y^2)^{3/2}} - \frac{q\ G M y^\prime}{(x^{^\prime2} + y^{^\prime2})^{3/2}} \end{align} $
$ \Large \begin{align} x^\prime = x + 0.8\ R \cos \left( \frac{2 \pi t}{T_{Moon}} \right) \end{align} $
$ \Large \begin{align} y^\prime = y + 0.8\ R \sin \left( \frac{2 \pi t}{T_{Moon}} \right) \end{align} $
where q = 0.00025 and T_moon = 27.3
I have tried to implement this with the velocity Verlet algorithm and have obtained a height verus time graph for the orbit. The space probe is supposed to hit the lunar surface at one point and here is the graph of the height versus time until the probe collides with the surface :
Here are the height oscillations versus time when I set the initial height so that the probe does not collide with the lunar surface after many orbits, compared with the previous result:
Is this physically correct with the equations of motion given above? Would this really happen in reality? How can I interpret this result? Here is the code I used in python to generate this. I used the velocity verlet algorithm to generate the numerical solution to the orbit and I am not sure if the implementation is correct, any comments on my implementation are very welcome.
from IPython.display import HTML
import matplotlib.pyplot as plt
from matplotlib import animation, rc
import numpy as np
from scipy import integrate, optimize
import time
plt.rcParams['figure.figsize'] = (10, 6)
plt.rcParams['font.size'] = 16
# Gravitational constant
gg = 6.67408e-11 # m^3 s^-1 kg^-2
# Lunar mass
mass = 7.342e22 # kg
# Lunar radius
radius = 1738000 # m
# 1 day in seconds
day = 3600*24 # seconds
### Initial positions and velocities at t=0
rs = [1842280, 0] # m
vs = [0, 1634] # m/s
N=10**5
t_moon = 27.3*24*3600
ts = np.linspace(0,t_moon, N-1)
q = 0.00025
rs = np.array([1842280, 0]) # m
vs = np.array([0, 1634]) # m/s
distance = np.array([(np.linalg.norm(rs) - radius),0])
for i in ts:
x_prime = rs[0] + 0.8*radius*np.cos((2*np.pi*i)/t_moon)
y_prime = rs[1] + 0.8*radius*np.sin((2*np.pi*i)/t_moon)
r_primes = np.array([x_prime,y_prime])
acs = np.array([-((gg*mass*rs[0])/(np.linalg.norm(rs)**3))-((q*(gg*mass*r_primes[0]))/(np.linalg.norm(r_primes)**3)),-((gg*mass*rs[1])/(np.linalg.norm(rs)**3))-((q*(gg*mass*r_primes[1]))/(np.linalg.norm(r_primes)**3))])
rs = rs + vs*ts[1] + acs*((ts[1]**2)/2)
next_acs = np.array([-((gg*mass*rs[0])/(np.linalg.norm(rs)**3))-((q*(gg*mass*r_primes[0]))/(np.linalg.norm(r_primes)**3)),-((gg*mass*rs[1])/(np.linalg.norm(rs)**3))-((q*(gg*mass*r_primes[1]))/(np.linalg.norm(r_primes)**3))])
vs = vs + ((next_acs + acs)*ts[1])/2
pos = np.c_[pos,rs]
distance = np.c_[distance,[(np.linalg.norm(rs)-radius),i+ts[1]]]
if np.linalg.norm(rs) <= radius:
print("BANG")
print("time is: ", ts[len(pos[0,:])])
break
plt.xlabel("time(hrs)")
plt.ylabel("height(km)")
# plt.plot(pos[0,:], pos[1,:])
plt.plot(distance[1,:]/3600,distance[0,:]/(10**3))
