Задачи моделирования орбиты двух тел

Question

Перейдите к Обновлению 2 ниже, если вы не хотите читать слишком много фона.

Я пытаюсь реализовать модель для простого орбитального моделирования (два тела).

Однако, когда я пытаюсь использовать код, который я написал, графики, сгенерированные из результата, выглядят довольно странно.

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

Кажется, это работает нормально и само по себе строит графики правильно, пока я держу график на орбитальной плоскости. Но я хотел бы повернуть график к системе отсчета (родительскому телу), чтобы я мог видеть крутой трехмерный вид того, как выглядят орбиты (obvs).

Прямо сейчас я подозреваю, что ошибка заключается в том, что я преобразую два вектора состояния в орбитальной плоскости в поворот их в систему отсчета. Я использую уравнения из шага 6 этого документа для создания следующего кода из (но применяя отдельные матрицы вращения [ скопировано здесь ]):

from numpy import sin, cos, matrix, newaxis, asarray, squeeze, dot

def Rx(theta):
    """
    Return a rotation matrix for the X axis and angle *theta*
    """
    return matrix([
        [1, 0,           0           ],
        [0, cos(theta), -sin(theta)  ],
        [0, sin(theta), cos(theta)   ],
    ], dtype="float64")

def Rz(theta):
    """
    Return a rotation matrix for the Z axis and angle *theta*
    """
    return matrix([
        [cos(theta), -sin(theta),   0],
        [sin(theta), cos(theta),    0],
        [0,          0,             1],
    ], dtype="float64")

def rotate1(vector, O, i, w):
    # The starting value of *vector* is just a 1-dimensional numpy
    # array.
    # Transform into a column vector.
    vector = vector[:, newaxis]
    # Perform the rotation
    R = Rz(-O) * Rx(-i) * Rz(-w)
    res2 = dot(R, vector)
    # Transform back into a row vector (because that's what
    # the rest of the program uses)
    return squeeze(asarray(res2))

(Для контекста это полный класс, который я использую для модели орбиты.)

Когда я строю координаты X и Y из результата, я получаю это:

Но когда я меняю матрицу вращения на R = Rz(-O) * Rx(-i) Я получаю этот более правдоподобный сюжет (хотя, очевидно, отсутствует один поворот и немного смещен от центра):

И когда я уменьшу его дальше R = Rx(-i), как и следовало ожидать, я получаю это:

Итак, как я уже сказал, я совершенно уверен, что это не код орбитального вычисления, который ведет себя странно, а скорее некоторая ошибка в коде вращения. Но я не уверен, где сузить это, поскольку я довольно плохо знаком с математикой и математикой в ​​целом.

---

Обновление: основываясь на ответе Стохастика, я переставил матрицы (R = Rz(-O).T * Rx(-i).T * Rz(-w).T), но потом получил этот сюжет:

что заставило меня задуматься, было ли мое преобразование в координаты экрана каким-то образом неправильным, но для меня это выглядит правильно (и это тот же код, что и для более правильных графиков с меньшим поворотом), а именно:

def recenter(v_position, viewport_width, viewport_height):
    x, y, z = v_position
    # the size of the viewport in meters
    bounds = 20000000
    # viewport_width is the screen pixels (800)
    scale = viewport_width/bounds
    # Perform the scaling operation
    x *= scale
    y *= scale
    # recenter to screen X and Y measured from the top-left corner
    # of the viewport
    x += viewport_width/2
    y = viewport_height/2 - y
    # Cast to int, because we don't care about pixel fractions
    return int(x), int(y)

---

Обновление 2

Хотя я трижды проверил свою реализацию уравнений, а также вращения с помощью Стохастика, я все еще не могу заставить орбиты выйти правильно. Они по-прежнему выглядят в основном так же, как на графиках выше.

Используя данные из системы NASA Horizon, я настроил орбиту с определенными векторами состояний от МКС (2457380.183935185 = AD 2015-Dec-23 16:24:52.0000 (TDB)) и проверил их по элементам орбиты Кеплера для того же момент времени, который дает такой результат:

inclination :
  0.900246137041
  0.900246137041
true_anomaly :
  0.11497063007
  0.0982485984565
long_of_asc_node :
  3.80727461492
  3.80727461492
eccentricity :
  0.000429082122137
  0.000501850615905
semi_major_axis :
  6778560.7037
  6779057.01374
mean_anomaly :
  0.114872215066
  0.0981501816537
argument_of_periapsis :
  0.843226618347
  0.85994864996

Верхние значения - мои (рассчитанные) значения, а нижние значения - значения НАСА. Очевидно, что следует ожидать некоторой ошибки точности с плавающей запятой, но вариации mean_anomaly а также true_anomaly действительно ударил меня как больше, чем я ожидал. (В настоящее время я выполняю все свои вычисления с использованием float128 номера в 64-битной системе).

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

Answer 1

Я полагаю, у вас есть как минимум две проблемы.

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

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

О транспонированных матрицах вращения

В цитируемом вами документе не указывается, должны ли R_x и R_z быть правосторонними поворотами осей или вектора, на который они будут умножаться, хотя вы можете понять это из уравнения 9 (или 10). Оказывается, это должны быть правосторонние повороты осей, а не вектора. Это означает, что они должны быть определены следующим образом:

    return matrix([
        [1, 0,           0           ],
        [0, cos(theta), sin(theta)  ],
        [0,-sin(theta), cos(theta)   ],
    ], dtype="float64")

вместо этого:

    return matrix([
        [1, 0,           0           ],
        [0, cos(theta),-sin(theta)  ],
        [0, sin(theta), cos(theta)   ],
    ], dtype="float64")

Я выяснил это, воспроизведя уравнение 9 от руки на бумаге.

• В этом уравнении посмотрите на первый компонент вектора r(t).
• Есть два термина: один с o_x в нем и один с o_y.
• Посмотрите на вещь мультипликационный o_y. Это: -(sin(omega)*cos(Omega)+cos(omega)*cos(i)*sin(Omega)),
• Этот ведущий знак минус является ключом. Это происходит от знака минус в первом ряду вашей матрицы Rz.
• Так как все омега, i и омега в уравнении 9 являются отрицательными, это означает, что знак минуса должен находиться во второй строке R_z, что будет означать, что R_z представляет правостороннее вращение осей, а не вектор.
• Точно так же мы можем посмотреть на компонент o_y последнего члена и увидеть, что знак минус должен находиться во втором ряду R_x, что означает (слава богу, здравомыслие) как R_z, так и R_x правостороннее вращение осей.

Ваши функции Rx и Rz в настоящее время определяют правостороннее вращение вектора, а не осей.

Вы можете исправить это любым (все три эквивалентны):

• Удаление знаков минус на углах Эйлера: Rz(O) * Rx(i) * Rz(w)

• транспонирование ваших матриц вращения: Rz(-O).T * Rx(-i).T * Rz(-w).T

• перемещая - войдите в определение Rx и Rz к синусу второго ряда, как показано выше

Answer 2

Я собираюсь отметить ответ Стохастика как правильный, потому что а) он заслуживает баллов за то, что он был таким полезным, и б) его совет был в основном верным.

Однако источником странного сюжета на самом деле оказались эти строки в связанном Orbit учебный класс:

    self.v_position = self.rotate(v_position, self.long_of_asc_node, self.inclination, self.argument_of_periapsis)
    self.v_velocity = self.rotate(v_velocity, self.long_of_asc_node, self.inclination, self.argument_of_periapsis)

Обратите внимание, что self.v_position свойство обновляется до вызова поворота вектора скорости; при чтении кода можно также заметить, что я в своем уме решил сделать все методы значений орбитальных элементов обернутыми в @property декораторы, чтобы сделать расчеты более понятными.

Но, конечно, это также означает, что методы вызываются - и значения пересчитываются - каждый раз при обращении к свойству. Итак, второй звонок self.rotate() происходит с немного отличающимися значениями орбитальных элементов от первого вызова и, что более важно, со значениями, которые не совпадают на 100% правильно с векторами "текущего" положения и скорости!

Так что после нескольких дней удара головой об эту ошибку, я понял это из-за небольшого бритья, которое я делал в форме рефакторинга, и теперь все работает отлично.