Блог

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

 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 motionn")
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
  

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

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

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

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

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

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

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

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

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

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