Orbit of Earth: Plotting Trajectory

6 replies · 2K views · last reply
Graham87
Messages
72
Reaction score
16
Homework Statement
Write a program to calculate the orbit of the Earth using leapfrog and plot it
Relevant Equations
-GM r/r^3
Printed homework problem text asking to use the leapfrog method to compute Earth's orbit around the Sun, with the equation d2r/dt2 = -GM r/r^3 and its constants.


I am attempting this homework exercise but 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.
Matplotlib line plot of a green curve starting at the origin and bending steadily downward and to the left, ending near y = -5 on the scaled y-axis. The y-axis offset reads 1e6 + 1.4710000000e11.

Cheers!This is my code:
Python:
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()
 
drmalawi said:
vx = 3.0287*20**4
Ah!

Should be vx=3.0287*10**4 and
t_stop = 24 * 365 * 0.99 * h

Now it's correct.

Thanks!
 
drmalawi said:
Plot works now?
Works neatly !
Plot of a closed green elliptical orbit around a yellow dot at the origin, axes in units of 1e11 m ranging from -1.5 to 1.5; a small gap is visible at the top of the ellipse.
 
  • Like
  • Love
Likes   Reactions: PeroK and malawi_glenn
Great job!

Could you fix the plot so that the non-cicrular shape is manifest?
Perhaps also you can plot a circular orbit centered around the sun with the Earth's average orbital distance as well
 
  • Like
Likes   Reactions: Graham87
drmalawi said:
Great job!

Could you fix the plot so that the non-cicrular shape is manifest?
Perhaps also you can plot a circular orbit centered around the sun with the Earth's average orbital distance as well
Great idea.
This is much clearer.
Matplotlib plot with axes scaled by 1e11 showing a green Earth 5-year trajectory nearly circular around a yellow dot for the Sun at the origin, with a thin red perfect circle nearly overlapping it.
 
  • Like
Likes   Reactions: malawi_glenn