Я пытаюсь создать симуляцию из трех тел, состоящую из Солнца, Марса и планеты-изгоя, совершающей облет, который повлияет на орбиту Марса.Пока что то, что я имею, чем-то напоминает симуляцию.
Константа гравитации больше, чем должна, и радиус орбиты неверен.Марс вращается вокруг Солнца на отметке 1,52 а.е., но в этой программе они разделены на 9 а.е.Если я использую правильные значения для всех, орбиты не работают, особенно когда я изменяю импульс в строке 23 с (0,2,0) на (0,24100,0), что составляет орбитальную скорость 24,1 км./s.
Что я делаю не так?
from visual import *
from visual.controls import *
import math
scene = display(width=1200, height=1200) #set animation window size
scene.autoscale = 0 #disable zoom if orbit size increases
#constants
m1 = 1.989e30 #mass Sun
m2 = 6.39e23 #mass Mars
m3 = 1.898e27 #mass intruding object
G = 6.67e-7 #gravitational constant, supposed to be 6.67e-11 but doesn't work
AU = 1.496e11 #astronomical unit
R = 0.1 #object render size
T = 1000.00 #max run time
dt = 0.001 #run rate
#momenta
p1 = m1*vector(0,0,0) #Sun
p2 = m2*vector(0,2,0) #Mars
p3 = m3*vector(0,10,0) #Intruder
p0 = m1*vector(0,0,0) #zero vector to make something stationary
#object positions must be scaled down, but use correct for calculations
obj1 = sphere(pos=vector(-4.50,0), radius = R, color = color.orange, make_trail = True) #Sun
obj2 = sphere(pos=vector(4.5,0), radius = R/3, color = color.red, make_trail = True) #Mars
obj3 = sphere(pos=vector(10,-10), radius = R, color = color.white, make_trail = True) #Intruder
t = 0
while t < T:
rate(1000) #change this to make animation faster
#seperation vectors
r12 = (obj2.pos - obj1.pos)*AU
r21 = -r12
r13 = (obj3.pos - obj1.pos)*AU
r31 = -r13
r23 = (obj3.pos - obj2.pos)*AU
r32 = -r23
#forces, norm = unit vector x/|x|, mag = magnitude |x|
F12 = -G*m1*m2*norm(r12)/(mag(r12)**2)
F21 = -F12
F13 = -G*m1*m3*norm(r13)/(mag(r13)**2)
F31 = -F13
F23 = -G*m2*m3*norm(r23)/(mag(r23)**2)
F32 = -F23
#updating object positions
obj1.pos = obj1.pos + p0*dt/m1 #Sun is stationary
obj2.pos = obj2.pos + p2*dt/m2
obj3.pos = obj3.pos + p3*dt/m3
#updating momenta
p1 = p1 + (F21 + F31)*dt
p2 = p2 + (F12 + F32)*dt
p3 = p3 + (F13 + F23)*dt
t = t + dt