Untitled

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

G = 6.7e-11
mu = 2e29

Sun = sphere(pos = vec(-1.1e11, 0, 0), mass = 2e24, radius = 10*7e8, color = color.orange, make_trail = True)

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

#Pre loop
t = 0
dt = 2 * 30 * 24 * 60 * 60

Sun.v = vec(0,0,0)
Sun.p = -Earth.p

Earth.v = vec(0, 10, 0)
Earth.p = Earth.mass*Earth.v
Sun.p = -Earth.p

Ki_earth = Earth.p.mag**2/(2*Earth.mass)
Ki_sun = Sun.p.mag**2/(2*Sun.mass)

R = Earth.pos - Sun.pos

w = 0

#Work Energy
w_graph = graph(title = 'Work-Energy', xtitle = 'time(s)', ytitle = 'Joules')
eswx = gcurve(graph= w_graph, color=color.red)
eswy = gcurve(graph= w_graph, color=color.blue)
eswz = gcurve(graph= w_graph, color=color.green)

#Total Kinetic Energy
k_graph = graph(title = 'Total Change in Kinetic Energy', xtitle = 'time(s)', ytitle = 'Joules')
eskex = gcurve(graph= k_graph, color=color.red)
eskey = gcurve(graph= k_graph, color=color.blue)
eskez = gcurve(graph= k_graph, color=color.green)

#Potential Energy
pe_graph = graph(title = 'Total Energy vs (t)', xtitle = 'time(s)', ytitle = 'Joules')
estex = gcurve(graph= pe_graph, color=color.red)
estey = gcurve(graph= pe_graph, color=color.blue)
estez = gcurve(graph= pe_graph, color=color.green)

myrate = 2000

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

#Loop
while R.mag > Earth.radius + Sun.radius:
    rate(myrate)
    R = Earth.pos - Sun.pos   
    F_g = -G * Sun.mass * Earth.mass / R.mag2 * R.hat
    
    Fe = mu * (F_g.mag/Earth.mass)**2 * (-Earth.p.hat)
    Fs = mu * (F_g.mag/Sun.mass)**2 * (-Sun.p.hat)
    
    Earth.p = Earth.p + Fe * dt + F_g * dt
    Earth.v = Earth.p / Earth.mass
    Earth.pos = Earth.pos + Earth.v * dt
    
    Sun.p = Sun.p + Fs * dt - F_g * dt
    Sun.v = Sun.p / Sun.mass
    Sun.pos = Sun.pos + Sun.v * dt
    
    w = w + dot(F_g, Earth.v) * dt + dot(-F_g, Sun.v) * dt
    
    Ket_sun = Earth.p.mag**2/(2*Earth.mass)
    Ket_earth = Sun.p.mag**2/(2*Sun.mass)
    
    Ke_sun = Ket_sun - Ki_sun
    Ke_earth = Ket_earth - Ki_earth
    
    Ke = Ke_sun + Ke_earth
    
    Pe = -G * Sun.mass * Earth.mass / R.mag
    
    E = Pe + Ke
    
    
    eswx.plot(pos = (t, w))
    eswy.plot(pos = (t, Ke))
    eswz.plot(pos = (t, Pe))
    
    eskex.plot(t, Ke)
    eskey.plot(t, Pe)
    eskez.plot(t, Ke + Pe)
    
    estex.plot(R.mag, Ke)
    estey.plot(R.mag, Pe)
    estez.plot(R.mag, Ke + Pe)
    
    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