VPython неправильно вращается при использовании фактических значений - PullRequest
0 голосов
/ 25 ноября 2018

Я пытаюсь создать симуляцию из трех тел, состоящую из Солнца, Марса и планеты-изгоя, совершающей облет, который повлияет на орбиту Марса.Пока что то, что я имею, чем-то напоминает симуляцию.

Константа гравитации больше, чем должна, и радиус орбиты неверен.Марс вращается вокруг Солнца на отметке 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
...