# Earth and Moon computer simulation

Edit: I made a test model with 2 bodies. The sun and earth or w/e. But I now have the problem of zero divisions every time that one of the coordinates is zero.

import numpy as np
from numpy import exp, sqrt, pi, linspace
from time import clock,time
import pylab as pylab

S1 = 1E3                   #Mass Star
M1 = 1                     #Mass Planet

Rs = np.array([0, 0])      #position
R1 = np.array([0,200])

Vs = np.array([0,0])       #velocity
V1 = np.array([10,0])

As = np.array([0,0])       #acceleration
A1 = np.array([0,0])

dt = 0.1                   #timestep
N = 1000                   #steps
G = 1                       #G const

for i in range(1,N):

deltaR = Rs - R1
Rsq = deltaR**2
F12 = -G*M1*S1/(Rsq)
F21 = -F12

As = F21/S1
A1 = F12/M1

Vs += As*dt
V1 += A1*dt

Rs += Vs*dt + 0.5*As*dt**2
R1 += V1*dt + 0.5*A1*dt**2

x = [ Rs, R1]
y = [ Rs, R1]
pylab.ion()
pylab.scatter( x, y )
pylab.axis([-200, 200, -200, 200])
pylab.plot
pylab.clf()


Any tips?

So I want to simulate the solar system but want to start simple with one orbiting body. However, I never did anything like this before and was wondering if anyone here could give me some hints.

Clearly the first thing we need to use is the law of gravitation and Kepler. Then I'd like to throw everything into an array and just calculate the timesteps. However, I am not sure on how to get the position of the planet/moon as a function of time. (Can I just integrate the acceleration via the law of gravitation?)

• Look at answers to this question as you will need a somewhat specialized time advancement technique for this problem. – Godric Seer May 15 '13 at 15:24
• To summarize: you're dealing with the $n$-body problem, $n=2$. The DEs are in that article. As @Godric hints, you'll need a symplectic integrator for your Hamiltonian system. – J. M. May 15 '13 at 15:36

If you want to start simple, you can just use simple kinematics: \begin{align*} a(t) &= \text{Gravity}[x(t)] \\ v(t+ \Delta t) &= v(t) + a(t) \Delta t \\ x(t+\Delta t) &= x(t) + v(t) \Delta t + a(t) (\Delta t)^2/2 \\ \end{align*} That should give you something that looks OK for small time-steps, $\Delta t$, and short duration. I also think it helps demonstrate the connection between continuous physics and explicit methods.