Orbit of the Earth – leapfrog Python

I am trying to do this homework exercise: Orbit of the Earth

My plot does not show the whole trajectory. I don’t know if it is something wrong with my equations or if it is a plotting matter.

Cheers!

Earth trajectory Plot

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as spi

G = 6.6738*10**-11
M = 1.9891*10**30
h = 3600
y = 1.4710*10**11
vx = 3.0287*20**4

def LeapFrog(f, t_start, t_stop, z0, h):

    t_vec = np.arange(t_start, t_stop, h) 
    n = len(t_vec)
    d = len(z0)  
    z_vec = np.zeros((n, d))
    z_vec[0,:] = z0
    z_half_step=z_vec[0 , :] + (1/2)*h*f(z0,t_vec[0]) 
    
    
    for i in range(0, n - 1):
        z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])
        z_half_step += h*f(z_vec[i+1,:], t_vec[i])   

    return t_vec, z_vec, 


def f(z,t):   
    
    x=z[0]
    y=z[1] 
    vx=z[2] 
    vy=z[3] 
    r=np.sqrt(x**2+y**2)

    dz=np.zeros(4)
    
    dz[0]=vx 
    dz[1]=vy
    dz[2]=-G*M*x/r**3
    dz[3]=-G*M*y/r**3

    return dz

t_start = 0
t_stop = 24*365*5
z0 = np.array([0,y,vx,0])
t_vec, z_vec = LeapFrog(f, t_start, t_stop, z0, h)

plt.plot(z_vec[:,0],z_vec[:,1], 'g', markersize=1, label='Earth trajectory')
plt.plot(0,0,'yo', label = 'Sun positon')  
plt.show()

Answer

Okay, so there’s a pretty funny typo:

vx = 3.0287*20**4

That’s 16 times larger than what you most likely wanted to write and just 20% short of the solar escape velocity.

t_stop is also written as if in hours, but the rest of the code assumes normal SI seconds (as can be verified by using t_stop = 24 * 365 * 0.99 * h, which results in almost a complete orbit). It needs multiplying by h. Sidenote: I like using astropy.units to keep track of physical quantities!

The code itself works nicely, though!