наклонение планеты, ось вращения трассы по орбите с питоном - PullRequest
Новая
       29

наклонение планеты, ось вращения трассы по орбите с питоном

0 голосов
/ 04 апреля 2019

Я пытаюсь сделать простое моделирование планеты, которая вращается вокруг Луны. Пока у меня есть проблема с двумя телами, которая решает планету и орбиту Луны. Теперь я хотел бы добавить фиксированную ось вращения к планете и посмотреть, как на нее влияет луна. Любая идея, как это можно сделать с помощью Python?

Проблема двух тел может быть запущена с кодом ниже: 

import pylab
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# Set Constants
G = 6.67e-11
AU = 1.5e11
daysec = 24.0*60*60

Ma =5.972e24   # Planet mass in Kg
Mb = 7.348e22   # Moon mass in Kg


gravconst = G*Ma*Mb

# Set up starting conditions

# Planet
xa = 0.0
ya = 0.0
za = 0.0

xva = 0.0
yva = 0.0
zva = 0.0

# Moon
xb = 384400000
yb = 0.0
zb = 0.0

xvb = 0.0
yvb = 1000.0
zvb = 0.0

# Time steps
t = 0.0
dt = 0.01*daysec

# Coordinate lists
xalist = []
yalist = []

xblist = []
yblist = []

zalist = []
zblist = []

# Loop
while t < 100.0*daysec:
    # Compute Force
    rx = xb-xa
    ry = yb-ya
    rz = zb-za

    modr3 = (rx**2+ry**2+rz**2)**1.5

    fx = -gravconst*rx/modr3
    fy = -gravconst*ry/modr3
    fz = -gravconst*rz/modr3

    # Update quantities
    # Moon
    xvb += fx*dt/Mb
    yvb += fy*dt/Mb
    zvb += fz*dt/Mb

    xb += xvb*dt
    yb += yvb*dt
    zb += zvb*dt

    # Planet
    xva += -fx*dt/Ma
    yva += -fy*dt/Ma
    zva += -fz*dt/Ma

    xa += xva*dt
    ya += yva*dt
    za += zva*dt

    t += dt

    # Saving them in lists
    xalist.append(xa)
    yalist.append(ya)
    zalist.append(za)

    xblist.append(xb)
    yblist.append(yb)
    zblist.append(zb)


xalist[:] = [x / 1e6 for x in xalist]
yalist[:] = [x / 1e6 for x in yalist]
zalist[:] = [x / 1e6 for x in zalist]

xblist[:] = [x / 1e6 for x in xblist]
yblist[:] = [x / 1e6 for x in yblist]
zblist[:] = [x / 1e6 for x in zblist]

#Creating the point to represent the planet at the origin (not to       scale),
plt.scatter(0,0,s=200,color='blue')
plt.annotate('Planet', xy=(-45,-50))
plt.scatter(xblist[0],0,s=100,color='grey')
plt.annotate('Mond', xy=(xblist[0]-45,-50))

# Plotting
pylab.plot(xalist, yalist, "-g")
pylab.plot(xblist, yblist, "-r")
plt.axhline(0, color='black')
plt.axvline(0, color='black')
pylab.axis("equal")
pylab.xlabel("X (Mio. Meter)")
pylab.ylabel("Y (Mio. Meter)")
pylab.show()

1 Ответ

0 голосов
/ 05 апреля 2019

Не ответ , так как я не эксперт в этом вопросе, а только некоторые подсказки (было невозможно прочитать в форме комментария)

Материал, который вы хотите добавить, вполнесложный , так как вам нужно было бы принять во внимание:

  1. движущиеся массы обоих тел

    , поэтому вам нужно "превышение контактных поверхностей"тело, которое имеет какие-либо движущиеся массы (такие как океаны, магма, вращающееся ядро ​​и т. д.), так что вы можете вычислять реальный центр тяжести каждый раз.Также вам нужно приложить силы и к самой движущейся массе (которая движется не только гравитацией и вращением, но также резонансом и в основном инерцией), не забывайте, что у Земли есть также ядро ​​и магма вокруг нее, а не только океаны, поэтому вам нужно принять во вниманиеучитывать как минимум 3 поверхности ...

  2. неоднородность распределения массы обоих тел

    , чтобы вы могли вычислить центр тяжести, иквадратичная инерция массы для вращений относительно фактической оси вращения

  3. планета / луна обычно составляет не менее 3 проблем тела, а не только 2

    , поскольку существуетЯвляется ли местная звезда слишком сильно влияющей на орбиту Луны ...

В зависимости от чисел некоторые эффекты будут настолько малы, что их можно будет отбросить, но некоторые - нет (особенно с инерцией +резонанс на месте).

Уравнения вращения похожи на положение / скорость / ускорение, как вы уже получили.Она называется интеграция Ньютона Даламбера / физика , но для этого вам потребуется реализовать матрицы преобразования .

См. Несколько связанных QAs :

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

Как вы можете видеть, это очень много вещей для обработки, и большинство программ моделирования, которые я видел, не делают этого (включая мою) ... онивместо этого имитируйте его константами нутации и прецессии.

...