Определение орбитальной позиции в будущем моменте времени

Постановка задачи

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

Как вы обнаружили, единственный способ распространить декартовы элементы вперед во времени — это численное интегрирование. Это работает нормально, но может быть медленным, если вам нужна высокая точность, и есть некоторые проблемы с числами (ошибки, вызванные округлением [медленно накапливаются, и многие интеграторы вызывают дрейф энергии ). Вы можете обойти некоторые из этих проблем, используя интегратор более высокого порядка ( популярным является Рунге-Кутта ), который позволяет делать большие шаги для того же уровня точности или получать лучшую точность для того же размера шага. Тем не менее, это излишество для простого моделирования.

Если вашу симуляцию можно рассматривать как задачу двух тел , то все значительно упрощается. Задача двух тел является хорошим упрощением, если на объекты симуляции в основном влияет один большой объект. Например, Земля, движущаяся вокруг Солнца, или космический корабль, движущийся по низкой околоземной орбите, хорошо моделируются как задача двух тел; однако космический корабль, летящий с Земли на Луну, - нет (об этом позже).

Поскольку вы пытаетесь смоделировать положение планет со средней точностью, сведение к задаче двух тел должно вам помочь.

Значение терминов

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

1. большая полуось,$a$: половина максимального диаметра эллиптической орбиты (= радиус цикла, если орбита круговая). Энергия и период обращения зависят только от$a$. полуширокая прямая кишка $\ell$, «ширина» орбиты, может быть лучшим выбором для орбит, близких к параболическим (например, для астероидов) или меняющихся с эллиптических на гиперболические (например, межпланетные космические корабли). Эти два связаны$\ell=a(1-e^2)$.

2. Эксцентриситет,$e$: "Заостренность" орбиты. Диапазоны от$e=0$для идеально круговой орбиты, чтобы$e=1$для параболической орбиты, чтобы$e>1$для гиперболических орбит. Меркурий — самая эксцентричная планета с$e \approx 0.2$. Космические аппараты на околоземной орбите обычно имеют$e<0.01$.

Помимо$e$а также$a$мы можем определить самую дальнюю и самую близкую точки орбиты, апоцентр и перицентр (вместе апсиды ):

$$r_a = a(1+e) \\ r_p = a(1-e)$$ _{Названия этих точек немного забавны: апоцентр и перицентр — общие термины, но орбиты вокруг конкретных тел имеют специфические термины: космический корабль вокруг Земли имеет апогей и перигей , а Земля (на орбите вокруг Солнца) имеет афелий и перигелий .}

Два параметра$a$а также$e$достаточно, чтобы определить форму орбиты. Следующие три параметра определяют ориентацию орбиты относительно системы координат, состоящей из базовой плоскости и базового направления (параллельного плоскости).

Почти для всех орбит в Солнечной системе используемой системой координат является эклиптическая система координат . Базовой плоскостью является плоскость эклиптики , плоскость обращения Земли вокруг Солнца. Опорным направлением является точка весеннего равноденствия , направление от Земли к Солнцу в момент весеннего равноденствия. Поскольку обе эти ссылки медленно дрейфуют во времени, мы должны указать конкретное время, в которое определяются эти ссылки, называемое эпохой . Наиболее распространенным является J2000 , полдень 1 января 2000 года (UTC).

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

Следующие параметры определяют орбиту относительно орбиты Земли:

3. Наклон,$i$: угол между плоскостью орбиты и базовой плоскостью. Наклон от 90 до 180 градусов относится к ретроградной орбите, которая вращается «назад» от обычного направления.

4. Долгота восходящего узла,$\Omega$: восходящий узел находится там, где орбита пересекает опорную плоскость снизу вверх. (Находится на пересечении плоскости орбиты и плоскости отсчета)$\Omega$угол между этой точкой и опорным направлением, измеренный против часовой стрелки.

5. Аргумент перицентра,$\omega$: угол между восходящим узлом и перицентром (самая низкая точка орбиты). Для орбит с очень малым наклонением, где трудно определить положение восходящего узла (поскольку это пересечение двух почти параллельных плоскостей), вместо этого мы используем долготу перицентра. $\varpi = \Omega + \omega$.

Шестой параметр определяет положение объекта на его орбите. Есть несколько вариантов, но самый распространенный:

• Средняя аномалия,$M$: «мнимый» угол, который равен нулю в перицентре и увеличивается с постоянной скоростью 360 градусов на орбиту.

Скорость, с которой$M$изменения называется средним движением ,$n$, равно$2\pi/T$. Обычно у вас есть измерение$M$в определенную эпоху$t_0$, называемая (что неудивительно) средней аномалией в эпоху ,$M_0$.

Точно так же, как аргумент перицентра, для орбит с малым наклонением мы используем связанное значение, среднюю долготу ,$L=\varpi + M$.

Фактический угол между вращающимся телом и перицентром называется истинной аномалией .$\nu$. Это угол, который нам нужен для вычисления положения тела. К сожалению, нет возможности напрямую вычислить$\nu$из$M$. Вместо этого мы сначала находим эксцентрическую аномалию . $E$:

$$M = E - e \sin E$$

Это называется уравнением Кеплера , и его нельзя решить аналитически. Как только у нас есть$E$хотя существует относительно простое выражение для$\nu$.

Вычисление положения по элементам орбиты

Мы выполним это вычисление в три этапа: сначала решим уравнение Кеплера. Во-вторых, мы вычислим 2-е положение тела в плоскости орбиты. Наконец, мы повернем нашу 2D-позицию в 3D-координаты. Я приведу некоторый «псевдокод» на Javascript для большинства этих задач.

Я предполагаю, что вы используете набор подобных элементов с веб-сайта JPL . Они используют$L$а также$\varpi$вместо$M$а также$\omega$. Таблица дает два значения для каждого из элементов; вторая производная по времени. Если вы используете значения из этой таблицы, вы также должны использовать производные.

Рассчитать время$t$в веках от J2000:

// month is zero-indexed, so 0 is January
var tMillisFromJ2000 = Date.now() - Date.UTC(2000, 0, 1, 12, 0, 0);
var tCenturiesFromJ2000 = tMillisFromJ2000 / (1000*60*60*24*365.25*100);

Теперь вычисляем текущие значения каждого из параметров орбиты. Например, большая полуось Земли с использованием значений из таблицы 1 (действительна с 1800 по 2500 год):

// a0 = 1.00000261; adot = 0.00000562
var a = a0 + adot * tCenturiesFromJ2000;

(Обратите внимание, что значения на самом деле даны для «EM Barycenter », центра масс системы Земля-Луна. Земля находится примерно в 4600 километрах от барицентра в направлении, противоположном Луне. Если вы хотите исправить это неточность, вам также нужно будет имитировать движение Луны, но это, вероятно, излишне.)

В таблице 2а приведены элементы с точностью от 3000 г. до н.э. до 3000 г. н.э.; однако, если вы используете элементы из таблицы 2а, вы должны дополнить их поправками на$L$из таблицы 2b! Например, вот вычисление долготы Сатурна:

// L0 = 34.33479152; Ldot = 3034.90371757
// b = -0.00012452          
// c =  0.06064060
// s = -0.35635438
// f = 38.35125000
var L = L0 + Ldot * tCenturiesFromJ2000
           + b * Math.pow(tCenturiesFromJ2000, 2)
           + c * Math.cos(f * tCenturiesFromJ2000)
           + s * Math.sin(f * tCenturiesFromJ2000);

Нам не нужно явно вычислять среднее движение и добавлять его в$L$, так как обе таблицы включают его в$\dot L$.

Теперь мы готовы вычислить$M$а также$\omega$( w):

var M = L - p \\ p is the longitude of periapsis
var w = p - W \\ W is the longitude of the ascending node

Переходим к шагу 2: нам нужно решить уравнение Кеплера:

$$M = E - e \sin E$$

Мы можем решить это численно, используя метод Ньютона . Решение уравнения Кеплера эквивалентно нахождению корней$f(E) = E - e \sin E - M$. Данный$E_i$, оценка$E$, мы можем использовать метод Ньютона, чтобы найти лучшую оценку:

$$E_{i+1} = E_i - f(E_i) / f'(E_i) \\ f'(E) = 1 - e \cos E$$

Поскольку нелинейная часть$e \sin E$очень мало, мы можем начать с оценки$E=M$. Наш код выглядит примерно так:

E = M;
while(true) {
  var dE = (E - e * Math.sin(E) - M)/(1 - e * Math.cos(E));
  E -= dE;
  if( Math.abs(dE) < 1e-6 ) break;
}

Теперь есть два способа вычислить положение по эксцентричной аномалии. Мы можем сначала вычислить истинную аномалию и радиус (положение объекта в полярных координатах), а затем преобразовать в прямоугольные координаты; однако, если мы применим немного геометрии, мы можем вместо этого вычислить координаты непосредственно из$E$:

var P = a * (Math.cos(E) - e);
var Q = a * Math.sin(E) * Math.sqrt(1 - Math.pow(e, 2));

( Pи Qсформируйте двумерную систему координат в плоскости орбиты с +Pуказанием на перицентр.)

Наконец, мы можем повернуть эти координаты в полную трехмерную систему координат:

// rotate by argument of periapsis
var x = Math.cos(w) * P - Math.sin(w) * Q;
var y = Math.sin(w) * P + Math.cos(w) * Q;
// rotate by inclination
var z = Math.sin(i) * y;
    y = Math.cos(i) * y;
// rotate by longitude of ascending node
var xtemp = x;
x = Math.cos(W) * xtemp - Math.sin(W) * y;
y = Math.sin(W) * xtemp + Math.cos(W) * y;

( x, y, и zбудет в единицах AU.)

И вы сделали!

---

Несколько советов:

• Если вы хотите также вычислить скорость, вы можете сделать это одновременно с вычислением.$P$а также$Q$, затем поверните его таким же образом.

  $$\dot M = n = \dot L \\ \dot M = \dot E - e (\cos E) \dot E \\ \dot E = \dot M / (1 - e \cos E) \\ \dot P = -a (\sin E) \dot E \qquad \dot Q = a (\cos E) \dot E \sqrt{1 - e^2}$$ Примечание. Я не включаю какие-либо производные (кроме$\dot L$) в этом расчете, так как они не сильно влияют на результат. Вы можете закодировать это как:

   var vP = - a * Math.sin(E) * Ldot / (1 - e * Math.cos(E));
   var vQ = a * Math.cos(E) * Math.sqrt(1 - e*e) * Ldot / (1 - e * Math.cos(E));
  

  Обратите внимание, что скорости будут в а.е. за столетие.

• Если вы очень часто обновляете позиции, вы можете использовать предыдущее значение$E$засеять метод Ньютона и выполнить фиксированное количество итераций (вероятно, будет достаточно только одной). Однако обратите внимание, что вам нужно сохранить это значение$E$локально для каждого объекта!

• Вы также можете просто использовать фиксированное количество итераций для исходного решения. Даже для$e=0.2$, после трех итераций ошибка в$E$только о$10^{-13}$, и после четырех итераций ошибка меньше, чем ошибка округления IEEE в два раза до$e=0.42$.

---

Если вам нужна дополнительная информация, вы можете поискать в Интернете, но если вы действительно заинтересованы, вам следует прочитать вводный текст по орбитальной механике. Лично я рекомендую «Основы астродинамики » Бейта, Мюллера и Уайта (pdf) . Мой отец пользовался этой книгой, когда учился в колледже, и я нашел ее более читабельной, чем мой учебник для колледжа. Вам будет интересна глава 4 «Положение и скорость как функция времени».

Поскольку это всего лишь игра, были бы вы довольны круговыми орбитами и орбитами планет, на которые влияет только центральное тело? В этом случае распространение довольно просто. В плоскости орбиты с центральным телом в точке (0,0) положение как функция времени:

$$x(t)=a\cos\left({2\pi \left(t-t_0\right)\over T}\right)$$

$$y(t)=a\sin\left({2\pi \left(t-t_0\right)\over T}\right)$$

куда$a$это большая полуось или в данном случае просто радиус орбиты,$T$- период обращения, а$t_0$определяет фазировку орбиты, где при$t=t_0$, планета находится на оси x с положительной стороны.

Чтобы орбиты разных планет согласовывались друг с другом, вам просто нужно определить$GM$центрального тела, которое мы будем называть$\mu$. Тогда для любого радиуса орбиты$a$, период обращения связан с$a$по:

$$T=2\pi\sqrt{a^3\over\mu}$$

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

---

Если вы не занимаетесь физикой N тел и планеты не взаимодействуют, вы можете использовать аналитические решения проблемы Кеплера. В конце концов вы поймете, что в какой-то момент вам нужно решить и гиперболические орбиты. Это приведет вас к формулировкам универсальных переменных для решения проблемы Кеплера.

Лучшим решением для этого, вероятно, будет метод Goodyear:

В. Гудиер, «Полностью общее решение в замкнутой форме для координат и частных производных задачи двух тел», The Astronomical Journal, Vol. 70, № 3, 1965, стр. 189–192 (или документ NASA NTRS TD по тому же материалу )

Метод Шепперда:

Шепперд, SW Небесная механика (1985) 35: 129. https://doi.org/10.1007/BF01227666

Или Дэнби-Стумпф:

Дэнби, JMA Celestial Mechanics (1987) 40: 303. https://doi.org/10.1007/BF01235847

Здесь есть некоторый код MATLAB, который может быть полезен (и гораздо более доступен), хотя случайные фрагменты кода на matlabcentral далеко не гарантированно свободны от ошибок, и похоже, что этому коду может не хватать полезной нормализации его входных данных (обычно вы собираетесь хотеть нормализовать масштаб вашей проблемы, чтобы вы выполняли математику в единицах, где r0-bar = 1,0 и mu-bar = 1,0, а где v-bar = 1 — это скорость на круговой орбите при r0 или что-то в этом роде) .

Если вы собираетесь выполнять интеграцию движения планет по N телам, то, я думаю, вам придется использовать численное интегрирование. Рунге-Кутта нарушит закон сохранения энергии, поэтому вы, вероятно, захотите использовать симплектическую интеграцию . Симплектический интегратор 4-го порядка в этой статье не так уж сложен в программировании, хотя это оставляет вам трудности с угадыванием правильного временного шага (опять же, нормализация помогает, потому что круговая планетарная орбита и круговая НОО — это одна и та же проблема, только с разными масштабами расстояний). ) и с интерполяцией внутренних точек (и вам нужно остерегаться феномена Рунге , но я не боролся с этим, поэтому не знаю, какой подход там выбрать).

Если вы собираетесь использовать Рунге-Кутту, то Дорманд-Принц с динамической стороной шага и его интерполянтом 3-го порядка будет очень удобным, и это то, что Matlab использует в своем решателе ode45.

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