I have assignment and part of it involves getting a satellite to orbit both the Earth and the Moon at least once. How ever no matter how I set my initial condition or parameterize my position and time evolve it my satellite won't orbit the moon and comeback. Note both the Moon and Earth lie on the x axis the system is in their co-rotating frame, I am using a fourth order Runge Kutta to simulate the orbit.
Update 1: After some suggestion from the comments below my script to explain the problem a bit more when I run my code the satellite spirals out and does not orbit , I have also added a plot of the satellite
Update 2 : For those wondering where I got my acceleration term for the satellite and why it includes Coriolis and the Centrifugal force along with the gravitational force between satellite and Earth and Satellite and moon see this reference http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node123.html and chapter 7 of this http://farside.ph.utexas.edu/teaching/336k/Newton.pdf
Updated 3: Fixed an error in the position of the satellite, and a typo in the comments of the code,updated the image that corresponded to the correction in the code
from pylab import plot,show,xlabel,ylabel,cla,clf
from numpy import linspace,exp,sin,cos,zeros,array,pi ,exp
M2 = 5.972e24 # M Earth kg
M1 = 7.34767309e22 # M moon kg
G= 6.6742e-11
R=3.844e8 #m Earth-moon distance
R_e =6371e3 # Earth radius
R_m = 1737e3 # Moon radius
M_0 = ((M2 - M1)/(M2 + M1))
alpha = pi/3
omega_0 = ((G*(M2 + M1))/(R**3))**(1/2)
T_a = (2*pi)/omega_0 # System period
def dot_p(a,b):
n= len(a)
d=0
for i in range(n):
d = d+ a[i]*b[i]
return d
def abs_v(x):
return (dot_p(x,x))**(1/2)
r_1 = array([(M2*R)/(M2 + M1),0,0],float) # Position of Moon
r_2 = array([-(M1*R)/(M2 + M1),0,0],float) #Position of Earth
omega = array([0,0,omega_0],float) #Angular velocity of Earth and Moon
def F(r,t):
r1 = r[0,:] #Earth position
vr1 = r[1,:] #Earth velocity
r2 = r[2,:] #Moon position
vr2 = r[3,:] #Moon velocity
r3 = array([R*cos(t),.5*R*sin(2*t),0],float) +r[4,:] #Satellite position Figure 8 parameterization
vr3 = r[5,:] #Satellite velocity
fr1 = -G*M2*(r1 -r2)/(abs_v((r1-r2))**3)
fr2 = -G*M1*(r2 -r1)/(abs_v((r2-r1))**3)
fr3 = (G*M1*(r3 -r_1)/(abs_v((r3-r_1))**3)) -(G*M2*(r3 -r_2)/(abs_v((r3-r_2))**3)) +2*omega_0*array([vr3[1],-vr3[0],0],float)+ (omega_0**2)*array([r3[0],r3[1],0],float) #Gravitational interaction of Satellite with Earth, Moon, Coriolis, and Centrifugal forces
return array((vr1,fr1 ,vr2,fr2,vr3,fr3),float)
def RungeKutta4o(g,ksi10,vksi10,ksi20,vksi20,ksi30,vksi30,a,b): #Runge Kutta 4th order
N=1000
M = 6
h = float((b-a)/N)
t = linspace(a,b,N)
x1 = zeros((N,3),float)
x2 = zeros((N,3),float)
x3 = zeros((N,3),float)
vx1 = zeros((N,3),float)
vx2 = zeros((N,3),float)
vx3 = zeros((N,3),float)
r = array((ksi10,vksi10,ksi20,vksi20,ksi30,vksi30),float)
k1 = zeros((M,3),float)
k2 = zeros((M,3),float)
k3 = zeros((M,3),float)
k4 = zeros((M,3),float)
for i in range(N):
x1[i,:] = r[1,:]
vx1[i,:] = r[0,:]
x2[i,:] = r[3,:]
vx2[i,:] = r[2,:]
x3[i,:] = r[5,:]
vx3[i,:] = r[4,:]
k1 = h*g(r,t[i])
k2 = h*g(r +.5*k1,t[i]+.5*h)
k3 = h*g(r +.5*k2,t[i]+ .5*h)
k4 = h*g(r +k3,t[i]+h)
r += (1/6)*(k1 + 2*k2 + 2*k3 + k4)
return t,x1,x2,x3,vx1,vx2,vx3
cla() # Clear axis
clf()
r_i0 = r_2 + array([0,R_e+300e3,0],float) # Satellite is initially in Earth orbit
vr_i0 = array([11000,0,0],float) # Initial velocity of satellite
t,R1,R2,R3,vR1,vR2,vR3 = RungeKutta4o(F,r_1,omega,r_2,omega,r_i0,vr_i0,0,10*T_a )
plot(R3[:,0],R3[:,1]) #Plots the Satellite orbit R3[:,0] is the x position of the satellite and R3[:,1] is the y position of the satellite
show()
anything**(1/2)returns1.0and1/6returns zero, since Python 2 uses integer math on integer types, it does not upcast to floats. Now I can get your spiral, yay! First suggestion, your time stephis 0.01 month, or about 6.5 hours. You need to use much smaller time steps of the order of minutes to handle the orbit around the earth. Second, in the CR3BP you just leave the Earth and Moon fixed, you don't calculate their motion. Only the spacecraft's position should appear in the calculation of acceleration (don't say "Force"). – uhoh Apr 18 '17 at 07:09list. When you are done, convert the list to an array. And like I said, don't let any information about the Earth or Moon even get into the calculation. The whole point of using the CR3BP formalism is to reduce the problem to a one body problem in a rotating frame, with an effective force that includes centrifugal. – uhoh Apr 18 '17 at 07:16scipy.odeintinstead of the manual RK4. I recommend you that to first make sure that your RK4 gives the same answer asodeint. – uhoh Apr 18 '17 at 07:20+1for math! Since you will need time steps of 1 minute or probably a few seconds with RK4, a two-week mission will need 200,000 to 1,000,000 time steps. If you want to maintain a dozen 3 x 1,000,000 arrays in memory and pop them in and out of cache, this program will be extremely slow. – uhoh Apr 18 '17 at 19:32scipy.odeintfirst, then try to see if it is even possible with RK4, which needs much much smaller time steps than the much higher order, variable step-size algorithms inodeint. Again, look at the python here. Fastest way to get a program to work is to get small parts of it to work, one at a time. Much faster than just posting a doesn't work script in stackexchange. – uhoh Apr 18 '17 at 19:37n=1000only, and I did that because it is# written for readability, not speed. Forn=1000000you drop the big arrays and just.append()each new position to a list. Each position could be a small array, or a list. btw I just checked it and it still works. If you download blender and paste that script into the script editor and hit run, it generates an animation that works. – uhoh Apr 18 '17 at 20:24