# Problems with Sun-Earth-Moon simulations in vPython

I'm trying to simulate the Sun-Earth-Moon system using the Verlet algorithm with vpython. I'm using natural units and a cartesian coordinate system, and am having trouble working out the physics.

I tried calculating the initial velocity of the moon by adding Earth's velocity with respect to the sun to the moon's velocity with respect to Earth, but what I'm getting is just an elliptical orbit of the moon around the sun.

I will point out that the simulation seems ok for just Earth, but I don't know how to make the moon orbit it.
I guess I need help with determining the correct initial velocity for the moon and any other remarks on my calculations and explanations on the physics would be greatly appreciated.

The code:

from vpython import *
from math import *

earth_mass = (5.97237 / 1.9855) * 10 ** (-6)  # Earth's proportional mass to the sun in natural units
moon_mass = (7.24767306 / 1.9855) * 10 ** (-8)  # Moon's proportional mass to the sun in natural units
moon_period = 27.32152778 / 365.25  # Moon's period around Earth in natural units
earth_moon_dist = 0.00256955529  # AU, Moon is ~384400 km from Earth

scene = canvas(fov=0.01, background=color.black)
dt = 0.001  # chosen dt
half_dt = dt / 2
half_dt_squared = (dt ** 2) / 2
M = 1  # Natural unit - 1 Sun mass
G = 4 * (pi ** 2)  # Gravitational constant in natural units

def acceleration(body_1, body_2):
# Calculates acceleration of body_1
sun_acc = (G * M / body_1.pos.mag2) * (-body_1.pos.hat)  # Force sun exerts on body_1 divided by body_1's mass
body_1_acc = (G * body_2.m / (body_1.pos - body_2.pos).mag2) * (
-(body_1.pos - body_2.pos).hat)  # Force body_2 exerts on body_1 divided by body_1's mass
return sun_acc + body_1_acc

# Bodies initialization

# Sun
sun = sphere(color=color.yellow, radius=.05, m=M, texture=r"https://i.imgur.com/lrI62Ot.jpg")

# Earth
earth_r = 1  # Natural unit - 1 AU
earth = sphere(pos=vector(earth_r, 0, 0), radius=.005, texture=textures.earth, make_trail=True, interval=dt,
trail_color=color.blue, m=earth_mass, v=vector(0, 0, 0), a=vector(0, 0, 0))
earth_initial_velocity = sqrt((G * M) * (1 / earth_r))
earth.v = vector(0, earth_initial_velocity, 0)

# Moon
moon_r = earth.pos.x + earth.radius + earth_moon_dist
moon = sphere(pos=vector(moon_r, 0, 0), radius=.001, texture=r"https://i.imgur.com/0lAj5pJ.jpg", make_trail=True,
interval=dt,
trail_color=color.white, m=moon_mass, v=vector(0, 0, 0), a=vector(0, 0, 0))

moon_initial_velocity = earth_initial_velocity + sqrt(G * earth.m / earth_moon_dist)
moon.v = vector(0, moon_initial_velocity, 0)

# Calculate accelerations
earth.a = acceleration(earth, moon)
moon.a = acceleration(moon, earth)
# Simulation
t = 0
while t <= 31557600:  # 1 Year in seconds
rate(100)
# Verlet algorithm - update positions
earth.pos += (earth.v * dt + earth.a * half_dt_squared)
moon.pos += (moon.v * dt + moon.a * half_dt_squared)
# Verlet algorithm - partially update velocities
earth.v += (earth.a * half_dt)
moon.v += (moon.a * half_dt)
# Verlet algorithm - update accelerations according to new positions
earth.a = acceleration(earth, moon)
moon.a = acceleration(moon, earth)
# Verlet algorithm - partially update velocities according to new accelerations
earth.v += (earth.a * half_dt)
moon.v += (moon.a * half_dt)
t += dt


EDIT: To emphasize, each time I'm calculating the acceleration of the moon, I use:
$$\vec{a_{M}}=\vec{a_{EM}}+\vec{a_{SM}}=\frac{\vec{F_{EM}}}{m_{M}}+\frac{\vec{F_{SM}}}{m_{M}}=\frac{G*M_{S}}{|\vec{r_{SM}}|^2}\hat{r_{SM}}+\frac{G*M_{E}}{|\vec{r_{EM}}|^2}\hat{r_{EM}}$$,

where $${F_{SM}}$$ is the force the sun exerts on the moon, $$M_{S}$$ is the mass of the sun, $$m_{M}$$ is the mass of the moon, and $$r_{SM}=\vec{S}-\vec{M}$$ with $$\vec{S}$$ the position vector of the sun (in my case the origin) and $$\vec{M}$$ the position vector of the moon.
The same for the force the earth exerts on the moon.

Are there any other forces I should be taking into consideration?

• Related. You may need a more careful visualization to see the Earth’s effect on the Moon’s orbit around the Sun, which is convex everywhere. – rob Jun 14 at 11:08
• I haven't looked at this in any detail, but you might try expanding some "+=". I can't remember the details, but I once hit a case where += did something unexpected. It took a long time to figure out what was going on. Interestingly, it was also in a Vpython program. – garyp Jun 14 at 11:51
• @rob I am pretty certain there is something wrong with my simulation as the moon is supposed to orbit around the earth, I tried looking in the link you sent but I didn't understand how it could help me. Could you please elaborate? – orav94 Jun 14 at 12:02
• @garyp Tried that just now, the code seems to function propely, it's the physics that are the problem. – orav94 Jun 14 at 12:03
• @rob I have added a clarification as to how I'm calculating the acceleration. They do not stay together and act like independent planets orbiting around the sun, the orbits don't even cross. This is why I think there is something wrong with my calculations. – orav94 Jun 14 at 12:51