Эллиптическая орбита в vpython - PullRequest
1 голос
/ 05 апреля 2019

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

import math
from vpython import *
lamp = local_light(pos=vector(0,0,0), color=color.yellow)
# Data in units according to the International System of Units
G = 6.67 * math.pow(10,-11)

# Mass of the Earth
ME = 5.973 * math.pow(10,24)
# Mass of the Moon
MM = 7.347 * math.pow(10,22)
# Mass of the Mars
MMa = 6.39 * math.pow(10,23)
# Mass of the Sun
MS = 1.989 * math.pow(10,30)
# Radius Earth-Moon
REM = 384400000
# Radius Sun-Earth
RSE = 149600000000
RMS = 227900000000
# Force Earth-Moon
FEM = G*(ME*MM)/math.pow(REM,2)
# Force Earth-Sun
FES = G*(MS*ME)/math.pow(RSE,2)
# Force Mars-Sun
FEMa = G*(MMa*MS)/math.pow(RMS,2)

# Angular velocity of the Moon with respect to the Earth (rad/s)
wM = math.sqrt(FEM/(MM * REM))
# Velocity v of the Moon (m/s)
vM = wM * REM
print("Angular velocity of the Moon with respect to the Earth: ",wM," rad/s")
print("Velocity v of the Moon: ",vM/1000," km/s")

# Angular velocity of the Earth with respect to the Sun(rad/s)
wE = math.sqrt(FES/(ME * RSE))
# Angular velocity of the Mars with respect to the Sun(rad/s)
wMa = math.sqrt(FEMa/(MMa * RMS))

# Velocity v of the Earth (m/s)
vE = wE * RSE
# Velocity v of the Earth (m/s)
vMa = wMa * RMS
print("Angular velocity of the Earth with respect to the Sun: ",wE," rad/s")
print("Velocity v of the Earth: ",vE/1000," km/s")


# Initial angular position
theta0 = 0

# Position at each time
def positionMoon(t):                                     
    theta = theta0 + wM * t
    return theta

def positionMars(t):                                     
    theta = theta0 + wMa * t
    return theta

def positionEarth(t):
    theta = theta0 + wE * t
    return theta


def fromDaysToS(d):
    s = d*24*60*60
    return s

def fromStoDays(s):
    d = s/60/60/24
    return d

def fromDaysToh(d):
    h = d * 24
    return h

# Graphical parameters
print("\nSimulation Earth-Moon-Sun motion\n")
days = 365
seconds = fromDaysToS(days)
print("Days: ",days)
print("Seconds: ",seconds)

v = vector(384,0,0)
E = sphere(pos = vector(1500,0,0), color = color.blue, radius = 60, make_trail=True)
Ma = sphere(pos = vector(2300,0,0), color = color.orange, radius = 30, make_trail=True)
M = sphere(pos = E.pos + v, color = color.white,radius = 10, make_trail=True)
S = sphere(pos = vector(0,0,0), color = color.yellow, radius=700)

t = 0
thetaTerra1 = 0
dt = 5000
dthetaE = positionEarth(t+dt)- positionEarth(t)
dthetaM = positionMoon(t+dt) - positionMoon(t)
dthetaMa = positionMars(t+dt) - positionMars(t)
print("delta t:",dt,"seconds. Days:",fromStoDays(dt),"hours:",fromDaysToh(fromStoDays(dt)),sep=" ")
print("Variation angular position of the Earth:",dthetaE,"rad/s that's to say",degrees(dthetaE),"degrees",sep=" ")
print("Variation angular position of the Moon:",dthetaM,"rad/s that's to say",degrees(dthetaM),"degrees",sep=" ")

while t < seconds:
    rate(500)
    thetaEarth = positionEarth(t+dt)- positionEarth(t)
    thetaMoon = positionMoon(t+dt) - positionMoon(t)
    thetaMars = positionMars(t+dt) - positionMars(t)
    # Rotation only around z axis (0,0,1)
    E.pos = rotate(E.pos,angle=thetaEarth,axis=vector(0,1,0))
    Ma.pos = rotate(Ma.pos,angle=thetaMars,axis=vector(0,1,0))
    v = rotate(v,angle=thetaMoon,axis=vector(0,1,0))
    M.pos = E.pos + v
t += dt

Мне интересно Как изменить путь орбиты на эллиптическую?Я пробовал несколько способов, но мне не удалось найти решение.

Спасибо.Спасибо

1 Ответ

2 голосов
/ 25 июля 2019

Это похоже больше на физическую проблему, чем на проблему программирования.Проблема в том, что вы предполагаете, что каждая из орбит круговая при расчете скорости и линейном интегрировании позиции (например, v * dt).Это не то, как вы бы рассчитывали траекторию орбитального тела.

Для простоты мы будем предполагать, что все массы являются точечными, поэтому нет никаких странных градиентов гравитации или динамики ориентации для учета.

Оттуда вы можете обратиться к этой странице MIT.(http://web.mit.edu/12.004/TheLastHandout/PastHandouts/Chap03.Orbital.Dynamics.pdf) по динамике орбиты. На 7-й странице приведено уравнение, связывающее радиальное положение вашего центрального тела как функцию множества параметров орбиты. Кажется, что у вас есть все параметры, кроме эксцентриситета орбитыВы можете посмотреть это онлайн или рассчитать, если у вас есть подробные эфемерные данные или информация об апоапсисе / периапсисе.

Из этого уравнения вы увидите термин phi-phi_0 в знаменателе. Это в разговорной речи называетсяистинная аномалия спутника. Вместо времени, вы бы итерировали этот параметр истинной аномалии от 0 до 360, чтобы найти ваше радиальное расстояние, и от истинной аномалии, наклона, прямого угла к восходящему узлу и аргумента периапсисов выможно найти трехмерные декартовы координаты для конкретной истинной аномалии.

Переход от истинной аномалии немного менее тривиален. Вам нужно будет найти эксцентрическую аномалию, а затем среднюю аномалию на каждом этапе эксцентричной аномалии. Теперь у вас естьзначит аномалия как афпринцип времени.Вы можете линейно интерполировать между «узлами», в которых вы рассчитываете положение с помощью v * dt.Вы можете рассчитать скорость, используя уравнение vis-viva, и dt будет разницей между рассчитанными временными шагами.

На каждом временном шаге вы можете обновлять положение спутника в вашей программе на Python, и он будет правильно рисовать ваши траектории.

Для получения дополнительной информации об истинной аномалии, Wikipedia имеет хорошее описание этого: https://en.wikipedia.org/wiki/True_anomaly

Для получения дополнительной информации об орбитальных элементах (которые необходимы для преобразования из радиальной позиции в декартовы координаты): https://en.wikipedia.org/wiki/Orbital_elements

...