- #1
jbay
- 9
- 0
So I am writing a program in python to do RK4 for the two body problem. I want it to display a sphere moving around another. It currently displays one sphere for a split second and then it goes blank. Any suggestions?
Code:
from __future__ import division
from visual import *
from visual.graph import *
G = 6.67E-11
msat = 1500E21
mearth = 5.98E24
i = 0
h = 100
earth = sphere(pos = (0,0,0), radius = 6378100, color = color.blue)
sat = sphere(pos = (384403000,0,0), radius = 3000000, color = color.red)
pearth = vector(0,-1000000E21,0)
psat = vector(0,1000000E21,0)
trail = curve(color = sat.color)
satvelocity = psat / msat
earthvelocity = pearth / mearth
satacceraltion = G * mearth * norm(sat.pos - earth.pos) / mag2(sat.pos - earth.pos)
earthacceraltion = G * msat * norm(earth.pos - sat.pos) / mag2(earth.pos - sat.pos)
Ag = arrow(pos=sat.pos, axis = 10000000000 * satacceraltion, color = color.red)
vel = arrow(pos=sat.pos, axis = 100000 * (psat / msat), color = color.red)
z = [earth.pos, earthvelocity, sat.pos, satvelocity]
Force = [earthvelocity, earthacceraltion, satvelocity, satacceraltion]
k1 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)]
k2 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)]
k3 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)]
k4 = [vector(0,0,0),vector(0,0,0), vector(0,0,0), vector(0,0,0)]while 1:
(100000)
psat = psat + msat * Force[3] * h
pearth = pearth - msat * Force[3] * h
sat.pos = sat.pos + Force[2] * h
earth.pos = earth.pos + Force[0] * h
trail.append(pos = sat.pos)
Ag.pos = sat.pos
Ag.axis = 10000000000 * Force[3]
vel.pos = sat.pos
vel.axis = 100000 * Force[3]
while (i < 4):
k1[i] = h * Force[i]
k2[i] = h *(z[i] + .5 * k1[i])
k3[i] = h *(z[i] + .5 * k2[i])
k4[i] = h *(z[i] + k3[i])
Force[i] = z[i] + h / 6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i])
i = i + 1
i = 0
Last edited by a moderator: