Untitled

Web VPython 3.2
scene = canvas(background = color.black)

G = 6.7e-11
Sun = sphere(pos = vec(0, 0, 0), mass = 1.989e30, radius = 10*7e8, color = color.orange, make_trail = True)

Earth = sphere(pos = vec(1.49e11, 0, 0), mass = 5.972e24, radius = 6.64e8, color = color.blue, make_trail = True) 

#Sun momentum
Sun.p = Sun.mass * vec(0,0,0)

#Earth momentum
Earth.p = Earth.mass * vec(0,2.978e4, 0)

#Pre loop
t = 0
dt = 60 * 60

#Graph
gr_p = graph(title = "Earth & Sun momentum", xtitle = "time(s)", ytitle = "p(t)")

#curve of Earth
p_plot_E_x = gcurve(graph = gr_p, color = color.red)
p_plot_E_y = gcurve(graph = gr_p, color = color.blue)
p_plot_E_z = gcurve(graph = gr_p, color = color.green)

#curve of Sun
p_plot_S_x = gcurve(graph = gr_p, color = color.red)
p_plot_S_y = gcurve(graph = gr_p, color = color.blue)
p_plot_S_z = gcurve(graph = gr_p, color = color.green)

#momentum plot total
gr_pt = graph(title = "Total p(t)", xtitle = "time(s)", ytitle = "p(t)")

p_plot_t_x = gcurve(graph = gr_pt, color = color.red)
p_plot_t_y = gcurve(graph = gr_pt, color = color.blue)
p_plot_t_z = gcurve(graph = gr_pt, color = color.green)

#graph and curves of velocity
gr_v = graph(title = "Earth & Sun velocity", xtitle = "time(s)", ytitle = "p(t)")
v_plot_E_x = gcurve(graph = gr_v, color = color.red)
v_plot_E_y = gcurve(graph = gr_v, color = color.blue)
v_plot_E_z = gcurve(graph = gr_v, color = color.green)

v_plot_S_x = gcurve(graph = gr_v, color = color.red)
v_plot_S_y = gcurve(graph = gr_v, color = color.blue)
v_plot_S_z = gcurve(graph = gr_v, color = color.green)

myrate = 1500

Earth.old_p = Earth.p
Earth.T = 0

#Loop
while t < 5 * 365 * 24 * 60 * 60: 
    rate(myrate)
    R = Earth.pos - Sun.pos 
    F_g = -G * Sun.mass * Earth.mass / R.mag2 * R.hat
    
    Earth.p = Earth.p + F_g * dt
    Earth.v = Earth.p / Earth.mass
    Earth.pos = Earth.pos + Earth.v * dt
    
    Sun.p = Sun.p - F_g * dt
    Sun.v = Sun.p / Sun.mass
    Sun.pos = Sun.pos + Sun.v * dt

    p_plot_E_x.plot(pos = (t,Earth.p.x))
    p_plot_E_y.plot(pos = (t,Earth.p.y))
    p_plot_E_z.plot(pos = (t,Earth.p.z))
    
    p_plot_S_x.plot(pos = (t,Sun.p.x))
    p_plot_S_y.plot(pos = (t,Sun.p.y))
    p_plot_S_z.plot(pos = (t,Sun.p.z))

    p_plot_t_x.plot(pos = (t,Earth.p.x+Sun.p.x))
    p_plot_t_y.plot(pos = (t,Earth.p.x+Sun.p.y))
    p_plot_t_z.plot(pos = (t,Earth.p.x+Sun.p.z))
    
    v_plot_E_x.plot(pos = (t,Earth.v.x))
    v_plot_E_y.plot(pos = (t,Earth.v.y))
    v_plot_E_z.plot(pos = (t,Earth.v.z))
  
    v_plot_S_x.plot(pos = (t,Sun.v.x))
    v_plot_S_y.plot(pos = (t,Sun.v.y))
    v_plot_S_z.plot(pos = (t,Sun.v.z))
    
    t = t + dt
    
    if Earth.old_p.x > 0 and Earth.p.x < 0:
        print('year: ', Earth.T/(60*60*24), ' days')
        Earth.T=0 
        
    Earth.old_p = Earth.p 
    Earth.T = Earth.T + dt