Untitled

Web VPython 3.2
G = 6.8e-11 #universal gravitational constant
work_e = 0
work_s = 0
mu = 2e29

sun = sphere(pos = vec(-1.1e11,0,0), radius = 7e9, color = color.yellow, make_trail = True,
mass = 2e24 )

earth = sphere(pos = vec(1.5e11, 0, 100), radius = 6.4e9, color = color.blue,
make_trail = True, mass = 1e24)

 #initial momentum
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

sun.v = vec(0,0,0)
sun.p = -earth.p

Ki = 0.5 * earth.mass * (earth.v)**2 +1/2*sun.mass*sun.v**2
R = earth.pos - sun.pos
Ui = -G * (earth.mass * sun.mass)/(R.mag**2)
#defining graphing environment

w = 0
delta_k = 0
delta_u = 0
#Momentum Graph
w_graph = graph(title = 'Work-Energy', xtitle = 'time [s]', ytitle = 'J')
epwx = gcurve(graph= w_graph, color=color.red)
epwy = gcurve(graph= w_graph, color=color.blue)
epwz = gcurve(graph= w_graph, color=color.green)

spwx = gcurve(graph= w_graph, color=color.cyan)
spwy = gcurve(graph= w_graph, color=color.yellow)
spwz = gcurve(graph= w_graph, color=color.purple)

k_graph = graph(title = 'Total Energy', xtitle = 'time [s]', ytitle = 'J')
epkx = gcurve(graph= k_graph, color=color.red)
epky = gcurve(graph= k_graph, color=color.blue)
epkz = gcurve(graph= k_graph, color=color.green)

spkx = gcurve(graph= k_graph, color=color.cyan)
spky = gcurve(graph= k_graph, color=color.yellow)
spkz = gcurve(graph= k_graph, color=color.purple)

kr_graph = graph(title = 'Energy vs Seperation', xtitle = 'time [s]', ytitle = 'J')
epkrx = gcurve(graph= kr_graph, color=color.red)
epkry = gcurve(graph= kr_graph, color=color.blue)
epkrz = gcurve(graph= kr_graph, color=color.green)

spkrx = gcurve(graph= kr_graph, color=color.cyan)
spkry = gcurve(graph= kr_graph, color=color.yellow)
spkrz = gcurve(graph= kr_graph, color=color.purple)


p_graph = graph(title = 'Momentum', xtitle = 'Time', ytitle = 'Momentum')#name and title the axis accordingly
#earth
epx = gcurve(graph= p_graph, color=color.red) #earth momentum in x direction
epy = gcurve(graph= p_graph, color=color.cyan)#earth momentum in y direction
epz = gcurve(graph= p_graph, color=color.purple)#earth momentum in z direction

#sun
spx = gcurve(graph= p_graph, color=color.yellow)#sun momentum in x direction
spy = gcurve(graph= p_graph, color=color.green)#sun momentum in y direction
spz = gcurve(graph= p_graph, color=color.blue)#sun momentum in z direction

#total Momentum graph
p_sum = graph(title = 'Total Momentum Graph', xtitle = 'Total Time', ytitle = 'Total Momentum')#name and title the axis accordingly

px = gcurve(graph= p_sum, color=color.red) #sum of momentum in x direction
py = gcurve(graph= p_sum, color=color.cyan)
pz = gcurve(graph= p_sum, color=color.purple)

v_graph = graph(title = 'Velocity', xtitle = 'Time', ytitle = 'Velocity')#name and title the axis accordingly

evx = gcurve(graph= v_graph, color=color.red) #earth velocity in x direction
evy = gcurve(graph= v_graph, color=color.cyan)#earth velocity in y direction
evz = gcurve(graph= v_graph, color=color.purple)#earth velocity in z direction

svx = gcurve(graph= v_graph, color=color.yellow)
svy = gcurve(graph= v_graph, color=color.green)
svz = gcurve(graph= v_graph, color=color.blue)

old_value = earth.v #old value to compare to updated value
T =0 #period variable

t = 0
dt = 24*60*60*30 # why did we chose a time step of 3600?

while t < 365*24*60*60*2000: #fufill the conditon so the simulation runs for 1 year in seconds
    rate(1600)
    
#define a force
#what is the vector between the two bodies?
    R = earth.pos - sun.pos
    F = (-R.hat) * (G * earth.mass * sun.mass)/(R.mag**2) # -R.hat gives direction radial inward
    FG = -G * sun.mass * earth.mass / R.mag**2 * R.hat
    Frad_b1 = -mu * (FG.mag/earth.mass)**2 * (earth.p.hat)
    Frad_b2 = -mu * (FG.mag/sun.mass)**2 * (sun.p.hat)
    
    F1 = (-FG + Frad_b1)
    F2 = (FG + Frad_b2)
    K = 0.5 * (earth.mass + sun.mass) * (sun.v + earth.v)**2
    U = -G * (earth.mass * sun.mass)/(R.mag**2)

#update momentum, velocity, and position for both objects (look at last weeks code for a hint)
    earth.p = earth.p + -F1 * dt
    earth.v = earth.p/earth.mass
    earth.pos = earth.pos + earth.v * dt
#for the earth and the sun momentum. which direction should the forces be?
    sun.p = sun.p + -F2 * dt
    sun.v = sun.p/sun.mass
    sun.pos = sun.pos + sun.v * dt

#use velocity graph to set if statment condition
    if old_value.x > 0 and earth.v.x <0:
        print('year is', T/(60*60*24), 'days') #convert to days
        T=0 #reset period variable
    
    epx.plot(pos = (t, earth.p.x ) )
    epy.plot(pos = (t, earth.p.y ) )
    epz.plot(pos = (t, earth.p.z ) )

    spx.plot(pos = (t, sun.p.x ) )
    spy.plot(pos = (t, sun.p.y ) )
    spz.plot(pos = (t, sun.p.z ) )

    px.plot(pos = (t, earth.p.x + sun.p.x) ) # how can we plot the total P??
    py.plot(pos = (t, earth.p.y + sun.p.y) )
    pz.plot(pos = (t, earth.p.z + sun.p.z) )

#Plot the each X Y & Z component of the velocity for the earth and the W.R.T
    evx.plot(pos = (t,earth.v.x ) )
    evy.plot(pos = (t,earth.v.y ) )
    evz.plot(pos = (t,earth.v.z ) )

    svx.plot(pos = (t, sun.v.x) )
    svy.plot(pos = (t, sun.v.y) )
    svz.plot(pos = (t, sun.v.z) )  
    
    epwx.plot(pos = (t, w))
    epwy.plot(pos = (t, K - Ki))
    epwz.plot(pos = (t, U - Ui))
    
    epkx.plot(t, K)
    epky.plot(t, U)
    epkz.plot(t, K + U)
    
    epkrx.plot(R.mag, K)
    epkry.plot(R.mag,U)
    epkrz.plot(R.mag, K + U)
    
    t = t + dt
    T = t
    old_value = earth.v