Численно интегрируйте и визуализируйте траекторию движения спутника под действием двух больших масс.
В орбитальной механике задача трех тел (ЗТТ) — это изучение того, как незначительная масса движется под действием двух больших масс (планеты, луны, звезды и т. д.). Результаты используются для расчета траекторий космических аппаратов, когда необходимо учитывать гравитационное притяжение двух больших масс. Например, траектория и текущая орбита телескопа Джеймса Уэбба были разработаны с использованием 3BP. Обычно, когда люди думают об орбитах, они думают о незначительной массе с эллиптической орбитой вокруг одной планеты или звезды. В задаче трех тел к системе добавляется дополнительное массивное тело, что усложняет вывод уравнений, описывающих движение меньшего тела.
Чтобы создать и визуализировать орбиты в 3BP, мы должны вывести уравнения движения для ничтожно малой массы. Я опускаю здесь детали вывода, но если вам интересно, в этой статье описано, как вывести уравнения для задачи трех тел. Я рекомендую вам изучить эту статью, чтобы полностью понять переменные и уравнения, которые вы собираетесь решать с помощью Python. Мы будем изучать задачу трех тел с круговыми ограничениями (CR3BP). В этой ограниченной версии более широкой задачи трех тел две большие массы (или основные массы) вращаются вокруг своего центра масс по круговым орбитам. Это упрощает уравнения движения космического корабля (или другой незначительной массы). На следующей диаграмме показана установка CR3BP.
Уравнения движения CR3BP выводятся для m₃ (космического корабля, астероида и т. д.) во вращающейся системе отсчета (x, y, кадр z-hat):
Эти уравнения немного пугают, но как только вы узнаете, что такое каждая переменная, они не так уж плохи. Здесь x, y и z (вектор ρ) — координаты положения m₃относительно центра масс основных цветов в вращающийся кадр (x, y, z-hat кадр), который перемещается вместе с основными цветами. Точки над этими членами обозначают производную по времени, поэтому одна точка обозначает скорость, а две точки обозначают ускорение. µ – безразмерное отношение масс основных цветов. r₁ и r₂ – расстояния м₃ от соответствующей первичной.
Приведенные выше уравнения движения получены с использованием безразмерных переменных. В упомянутой выше статье подробно рассказывается о том, как использовать безразмерные параметры, что важно, если вы действительно пытаетесь изучить проблему трех тел. Не вдаваясь в подробности, нижеприведенные параметры используются для удаления и добавления размеров (км, кг, с и т. д.) по мере необходимости. Мы можем использовать M* для управления единицами массы (кг), L* для единиц длины (км), и T* для единиц времени (с). Это будет объяснено далее в разделе кодирования.
Чтобы решить CR3BP, мы будем использовать решатель ОДУ для численного интегрирования уравнений движения для m₃. Чтобы сделать это, нам нужно определить вектор состояния и производный по времени вектор состояния. Для нашей задачи вектор состояния включает векторы положения и скорости в системе отсчета вращения. Производный по времени вектор состояния — это просто производный по времени вектор состояния. Они будут использоваться с помощью числового интегратора Python odeint. Здесь две точки x, y и z взяты из уравнений движения, определенных выше.
Мы можем преобразовать вектор во вращающуюся рамку (x, y и z ) до единицы в инерциальной системе отсчета (X, Y и Z), используя следующее уравнение:
Я знаю, что это был краткий обзор проблемы, но вам не нужно полностью понимать CR3BP, чтобы иметь возможность кодировать решение. Я предполагаю, что вы немного разбираетесь в численном интегрировании и в том, как реализовать его в Python. Если вам нужно какое-то руководство по этому поводу, я написал еще одну статью, которая может помочь, и в Интернете есть множество других ресурсов. Без лишних слов, давайте перейдем к коду!
Импорт пакетов
Для этого кода нам нужно будет импортировать несколько пакетов. В следующем списке кратко описаны пакеты, используемые в этом коде:
- NumPy используется для создания массивов и управления ими (определяется как
npдля простоты вызова). odeintиз библиотеки SciPy используется для численного интегрированияpyplotиз matplotlib используется для визуализации результатов численного интегрирования (определяется какpltдля простоты вызова)
# Importing Packages import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt
Создание функции для числового интегрирования
Следующим шагом является создание определяемой пользователем функции Python, model_CR3BP, которая будет использоваться в odeint для численного интегрирования нашего вектора состояния. Функция odeint использует текущий вектор состояния для создания производного по времени вектора состояния. Обратите внимание, что мы извлекаем компоненты положения и скорости из вектора состояния, чтобы создать производную по времени вектора состояния.
# CR3BP Model
def model_CR3BP(state, t):
x = state[0]
y = state[1]
z = state[2]
x_dot = state[3]
y_dot = state[4]
z_dot = state[5]
x_ddot = x+2*y_dot-((1-mu)*(x+mu))/((x+mu)**2+y**2+z**2)**(3/2)\
-(mu*(x-(1-mu)))/((x-(1-mu))**2+y**2+z**2)**(3/2)
y_ddot = y-2*x_dot-((1-mu)*y)/((x+mu)**2+y**2+z**2)**(3/2)\
-(mu*y)/((x-(1-mu))**2+y**2+z**2)**(3/2)
z_ddot = -((1-mu)*z)/((x+mu)**2+y**2+z**2)**(3/2)\
-(mu*z)/((x-(1-mu))**2+y**2+z**2)**(3/2)
dstate_dt = [x_dot, y_dot, z_dot, x_ddot, y_ddot, z_ddot]
return dstate_dt
Определение безразмерных параметров
Здесь мы определяем несколько параметров для нашей выбранной системы. Для этой статьи я выбрал систему Земля-Луна, но если вы хотите выбрать другую систему (Плутон-Харон, Солнце-Юпитер и т. д.), вы можете соответствующим образом настроить массы и переменные большой полуоси. Здесь m₁ и m₂ — массы Земли и Луны соответственно. a — большая полуось орбиты Луны (если бы у них была идеально круговая орбита). Я использовал среднее расстояние по орбите Луны для a, поскольку Луна вращается вокруг Земли практически по круговой орбите. Далее безразмерные параметры определяются с помощью уравнений из начала статьи.
# Defining ND Parameters
G = 6.67408E-20 # Univ. Gravitational Constant [km3 kg-1 s-2]
mEarth = 5.97219E+24 # Mass of the Earth [kg]
mMoon = 7.34767E+22 # Mass of the Moon [kg]
a = 3.844E+5 # Semi-major axis of Earth and Moon [km]
m1 = mEarth
m2 = mMoon
Mstar = m1+m2 # ND Mass Parameter
Lstar = a # ND Length Parameter
Tstar = (Lstar**3/(G*Mstar))**(1/2) # ND Time Parameter
mu = m2/Mstar
print('\u03BC = ' + str(mu))
Определение входных данных решателя ODE
Теперь, наряду с определяемой пользователем моделью, odeint нужны еще два входа: начальные условия и временной интервал для интегрирования. Я определил произвольный набор начальных условий, чтобы создать орбиту, которую вы видели в начале. Вы можете попробовать свои собственные начальные условия или использовать данные об эфемеридах из этого инструмента NASA JPL. Затем мы определяем вектор начального состояния state_0 и интересующий временной массив t. Обратите внимание, что размеры удалены с использованием параметров ND, указанных выше. Километры и секунды удаляются путем деления или умножения на T* и L*.
# Initial Conditions [Initial State Vector] X_0 = 50000/Lstar # ND x Y_0 = 0 # ND y Z_0 = 0 # ND z VX_0 = 1.08*Tstar/Lstar # ND x_dot VY_0 = 3.18*Tstar/Lstar # ND y_dot VZ_0 = 0.68*Tstar/Lstar # ND z_dot state_0 = [X_0, Y_0, Z_0, VX_0, VY_0, VZ_0] # ND ICs # Time Array t = np.linspace(0, 15, 1000) # ND Time
Модель численного интегрирования
Следующий шаг — объединение того, что мы уже закодировали, и использование всего этого в функции odeint. Здесь мы используем model_CR3BP, state_0 и t. Результатом численного интегрирования является вектор состояния для каждого шага, который мы определили в массиве времени. Из результатов мы можем извлечь x, y и zкоординаты (в системе вращения), определяющие положение m₃. Далее, используя уравнения из начала статьи, мы можем преобразовать вращательное состояние в инерционное состояние.
# Numerically Integrating sol = odeint(model_CR3BP, state_0, t) # Rotational Frame Position Time History X_rot = sol[:, 0] Y_rot = sol[:, 1] Z_rot = sol[:, 2] # Inertial Frame Position Time History X_Iner = sol[:, 0]*np.cos(t) - sol[:, 1]*np.sin(t) Y_Iner = sol[:, 0]*np.sin(t) + sol[:, 1]*np.cos(t) Z_Iner = sol[:, 2]
Добавление истории времени основной позиции
Этот шаг необязателен, но учет истории двух основных цветов может быть важен для понимания движения m₃. Положения первичных элементов являются стационарными во вращающейся раме (для CR3BP) и могут быть получены просто путем использования центра масс первичных элементов и преобразования в единицы ND. Затем, используя то же преобразование вращающейся системы отсчета в инерциальную систему координат, что и раньше, мы можем определить временную историю инерциального движения основных цветов.
# Constant m1 and m2 Rotational Frame Locations for CR3BP Primaries m1_loc = [-mu, 0, 0] m2_loc = [(1-mu), 0, 0] # Moving m1 and m2 Inertial Locations for CR3BP Primaries X_m1 = m1_loc[0]*np.cos(t) - m1_loc[1]*np.sin(t) Y_m1 = m1_loc[0]*np.sin(t) + m1_loc[1]*np.cos(t) Z_m1 = m1_loc[2]*np.ones(len(t)) X_m2 = m2_loc[0]*np.cos(t) - m2_loc[1]*np.sin(t) Y_m2 = m2_loc[0]*np.sin(t) + m2_loc[1]*np.cos(t) Z_m2 = m2_loc[2]*np.ones(len(t))
Визуализация данных
Последним шагом является построение вращательного и инерционного движения каждой массы. Чтобы различать орбиты, мы будем использовать зеленый для m₃, черный для Луны и синий для Земли. Мы также можем установить оси так, чтобы они отображали одинаковую длину осей, чтобы дать более точное представление о системе.
# Rotating Frame Plot
fig = plt.figure()
ax = plt.axes(projection='3d')
# Adding Figure Title and Labels
ax.set_title('Rotating Frame CR3BP Orbit (\u03BC = ' + str(round(mu, 6)) + ')')
ax.set_xlabel('x [ND]')
ax.set_ylabel('y [ND]')
ax.set_zlabel('z [ND]')
# Plotting Rotating Frame Positions
ax.plot3D(X_rot, Y_rot, Z_rot, c='green')
ax.plot3D(m1_loc[0], m1_loc[1], m1_loc[2], c='blue', marker='o')
ax.plot3D(m2_loc[0], m2_loc[1], m2_loc[2], c='black', marker='o')
# Setting Axis Limits
xyzlim = np.array([ax.get_xlim3d(), ax.get_ylim3d(), ax.get_zlim3d()]).T
XYZlim = np.asarray([min(xyzlim[0]), max(xyzlim[1])])
ax.set_xlim3d(XYZlim)
ax.set_ylim3d(XYZlim)
ax.set_zlim3d(XYZlim * 3 / 4)
# Inertial Frame Plot
fig = plt.figure()
ax = plt.axes(projection='3d')
# Adding Figure Title and Labels
ax.set_title('Inertial Frame CR3BP Orbit (\u03BC = ' + str(round(mu, 6)) + ')')
ax.set_xlabel('X [ND]')
ax.set_ylabel('Y [ND]')
ax.set_zlabel('Z [ND]')
# Plotting Inertial Frame Positions
ax.plot3D(X_Iner, Y_Iner, Z_Iner, c='green')
ax.plot3D(X_m1, Y_m1, Z_m1, c='blue')
ax.plot3D(X_m2, Y_m2, Z_m2, c='black')
# Setting Axis Limits
ax.set_xlim3d(XYZlim)
ax.set_ylim3d(XYZlim)
ax.set_zlim3d(XYZlim * 3 / 4)
plt.show()
Вывод кода:
Графики выше показывают результаты инерциальной и вращательной системы отсчета для орбиты m₃и первичных масс. Как видите, космический корабль имеет хаотичную орбиту для этого набора начальных условий. Если вы отрегулируете начальные условия, вы можете создать различные различные орбиты (некоторые даже периодические). Я призываю вас возиться с начальными условиями и экспериментировать самостоятельно. Чтобы статья и код были как можно короче, я не упомянул, как анимировать сюжет. Если вы заинтересованы в его добавлении, ознакомьтесь с пошаговым руководством ниже:
Спасибо за прочтение статьи! Это был краткий обзор задачи трех тел с ограничением по кругу и того, как численно интегрировать уравнения движения. Это продвинутая тема орбитальной механики, поэтому ее изучение может быть сложным. Если у вас есть какие-либо вопросы, обращайтесь сюда или в LinkedIn. Подпишитесь на еженедельные статьи об орбитальной механике, программировании и машинном обучении!
Если эта статья показалась вам интересной, вам могут быть интересны и следующие:
