Untitled
unknown
plain_text
a year ago
2.7 kB
14
Indexable
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
Editor is loading...
Leave a Comment