Отображение орбиты с помощью vpython с использованием уравнения Кеплера, но планета не будет вращаться - PullRequest
Новая
       103

Отображение орбиты с помощью vpython с использованием уравнения Кеплера, но планета не будет вращаться

0 голосов
/ 15 октября 2019

Я пытался сделать своего рода симуляцию солнечной системы, и я наткнулся на уравнение Кеплера, чтобы дать мне положение планеты в любой момент времени. Я использовал метод Ньютона-Рафсона для расчета эксцентрической аномалии. Тем не менее, я не могу заставить планету действительно вращаться должным образом, и я не уверен, что не так. Планета просто движется туда-сюда в очень крошечном месте.

Я заранее ценю помощь. 

from datetime import datetime, timedelta
from math import *
from vpython import *

 scene = canvas(title='Solar system simulation',
 width=1024, height=720,
 center=vector(0,0,0), background = color.black)

 sun = sphere(pos = vector(0,0,0), radius = 5, color = color.yellow)
 earth = sphere(pos=vector(1,1,0), radius = 1, color = color.blue, make_trail = True)
 earth.trail_color = color.white

 d1 = datetime.now()

while True:

    d1 = d1 + timedelta(minutes = 1)
    d2 = datetime(2000, 1, 1, 00, 00)

    SecondsFrom2000 = (d1 - d2).total_seconds()  #seconds between today and 01.01.2000 at midnight
    CenturiesFrom2000 = SecondsFrom2000/(60*60*24*365.25*100) #centuries between today and 2000
    e0Earth = 0.01671123 #eccentricity
    edotEarth = -0.00004392
    eEarth = e0Earth + edotEarth*CenturiesFrom2000

    a0Earth = 1.00000018 #semi-major axis[au], from 3000 BC - 3000 AD
    adotEarth = -0.00000003 #[au/century]
    aEarth = a0Earth + adotEarth*CenturiesFrom2000

    L0Earth = 100.46691572 #mean longitude [deg]
    LdotEarth = 35999.37306329 #[deg/century]
    LEarth = (L0Earth + LdotEarth*CenturiesFrom2000)*(pi/180) #[rad/century]

    Pi0Earth = 102.93005885 #longitude of the perihelion [deg]
    PidotEarth = 0.31795260 #[deg/century]
    PiEarth = (Pi0Earth + PidotEarth*CenturiesFrom2000)*(pi/180)

    W0Earth = -5.11260389 #longitude of the ascending node [deg]
    WdotEarth = -0.24123856
    WEarth = (W0Earth + WdotEarth*CenturiesFrom2000)*(pi/180)

    I0Earth = -0.00054346 #inclination [deg]
    IdotEarth = -0.01337178
    IEarth = (I0Earth + IdotEarth*CenturiesFrom2000)*(pi/180)

    MEarth = LEarth - PiEarth #mean anomaly
    wEarth = PiEarth - WEarth #argument of perihelion

    E0 = 0.1
    E1 = MEarth #E = eccentric anomaly
    error = 1

    while error > 0.000001:
        E2 = E1 - (E1 - eEarth*sin(E1)/(1 - eEarth*cos(E1)))
        E1 = E2
        erreur = abs(E2 - E1)

    E_Earth = E2

    PEarth = aEarth * (cos(E_Earth) - eEarth) #[au]
    QEarth = aEarth * sin(E_Earth)*sqrt(1 - eEarth*eEarth) #[au]

    xEarth = cos(wEarth)*PEarth - sin(wEarth)*QEarth
    yEarth = sin(wEarth)*PEarth + cos(wEarth)*QEarth
    zEarth = sin(IEarth)*xEarth
    xEarth = cos(IEarth)*xEarth
    xtempEarth = xEarth
    xEarth = cos(WEarth)*xtempEarth - sin(WEarth)*yEarth
    yEarth = sin(WEarth)*xtempEarth + cos(WEarth)*yEarth

    earth.pos = vector(xEarth*10,yEarth*10,zEarth*10)
...