Вычислительное решение задачи двух тел для аппроксимации заплатанных коник

Вы, вероятно, уже нашли ответ к тому времени, когда разместили этот вопрос, но я все же даю ответ на случай, если это поможет.

Использование орбитальных элементов - "правильный способ" сделать это. Их 6, описанных в этой статье Википедии :

• эксцентриситет$e$
• большая полуось$a$
• склонность$i$
• долгота восходящего узла$\Omega$
• аргумент перицентра$\omega$
• настоящая аномалия$\theta$(что зависит от времени)

Первые пять из них полностью описывают геометрию орбиты в 3D. Истинная аномалия — это угол, относящийся к положению тела/космического корабля на этой орбите. Все эти параметры можно рассчитать по векторам начального состояния (положения и скорости), зная стандартный гравитационный параметр тела-аттрактора. Математика описана в этом PDF-файле Рене Шварцем о том, как преобразовать декартовы орбитальные векторы состояния в орбитальные элементы.

Когда у вас есть элементы орбиты, вы можете использовать их для вычисления векторов состояния (вектор положения/радиуса и вектор скорости) в любой момент времени, следуя этому другому PDF-файлу о том, как преобразовать элементы орбиты в декартовы векторы состояния. Однако некоторые формулы, описанные в документе, работают только для эллиптических орбит. Для гиперболических орбит (когда$e > 1$) вам придется использовать другие формулы для эксцентрической и истинной аномалии (см. эту статью в Википедии и второй ответ на этот пост ), а также адаптировать формулу для вектора скорости из документа (см. этот пост ).

«Сложно выглядящие» матрицы вращения, которые задают положение и скорость космического корабля/тела в трехмерном пространстве, могут быть разработаны за следующие шаги:

• Выберите опорные векторы для «вправо» и «вверх», например, направление x (1, 0, 0) и направление y (0, 1, 0), я называю их$\vec{u_x}$и$\vec{u_y}$.
• Вычислить вектор восходящего узла$\vec{u_{asc}}$: это$\vec{u_x}$вращается вокруг$\vec{u_y}$под углом$\Omega$;
• Установить положение$\vec{r} = \vec{o}$и скорость$\vec{v} = \dot{\vec{o}}$как описано в документе (с измененной формулой в случае$e > 1$);
• Повернуть$\vec{r}$и$\vec{v}$вокруг$\vec{u_y}$угла$\Omega + \omega$;
• Повернуть снова$\vec{r}$и$\vec{v}$вокруг$\vec{u_{asc}}$угла$i$.

(Обратите внимание, что в документах Рене Шварца он считает ось Z осью «вверх», а не осью Y)

Случай параболической орбиты ($e = 1$) можно, как сказал нотовны в комментарии, вежливо проигнорировать (применить алгоритм страуса ).

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

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