Как программно рассчитать элементы орбиты, используя векторы положения/скорости?

Даны центрированные на Земле, инерциальные (ECI) положение и векторы скорости$\vec{r}$а также$\vec{v}$, вы можете напрямую решить для классических орбитальных элементов$(a,e,i,\Omega,\omega,\nu)$следующим образом (сначала алгоритмы, а затем псевдокод внизу):

Сначала найдите угловой момент

$$\vec{h}=\vec{r} \times \vec{v}$$

тогда вектор узла

$$\hat{n}=\hat{K} \times \vec{h}$$ который будет использован позже.

Тогда вектор эксцентриситета равен

$$\vec{e} = \frac{(v^2-\mu /r)\vec{r}-(\vec{r} \cdot \vec{v})\vec{v}}{\mu}$$

а также$e=|\vec{e}|$.

Удельная механическая энергия

$$E=\frac{v^2}{2}-\frac{\mu}{r}$$

Если$e\neq 1$, тогда

$$a = -\frac{\mu}{2E}$$ $$p=a(1-e^2)$$ В противном случае, $$p=\frac{h^2}{\mu}$$ $$a=\infty$$

В настоящее время,

$$i=\cos^{-1}{\frac{h_K}{h}}$$ $$\Omega=\cos^{-1}{\frac{n_I}{n}}$$ $$\omega=\cos^{-1}{\frac{\vec{n}\cdot\vec{e}}{ne}}$$ $$\nu=\cos^{-1}{\frac{\vec{e}\cdot\vec{r}}{er}}$$

И вам нужно будет сделать следующие проверки: Если$n_J<0$, тогда$\Omega=360^{\circ}-\Omega$,

Если$e_K<0$, тогда$\omega=360^{\circ}-\omega$, а также

Если$\vec{r}\cdot\vec{v}<0$, тогда$\nu=360^{\circ}-\nu$.

Обратите внимание, что в некоторых случаях вы столкнетесь с проблемами (особенностями): круговые орбиты ($e\approx 0$), и экваториальные орбиты ($i\approx 0$), особенно. В этих случаях вы обычно вводите новую, менее хлопотную переменную, например, среднюю долготу или истинную долготу перигея.

h=cross(r,v)
nhat=cross([0 0 1],h)

evec = ((mag(v)^2-mu/mag(r))*r-dot(r,v)*v)/mu
e = mag(evec)

energy = mag(v)^2/2-mu/mag(r)

if abs(e-1.0)>eps
   a = -mu/(2*energy)
   p = a*(1-e^2)
else
   p = mag(h)^2/mu
   a = inf

i = acos(h(3)/mag(h))

Omega = acos(n(1)/mag(n))

if n(2)<0
   Omega = 360-Omega

argp = acos(dot(n,evec)/(mag(n)*e))

if e(3)<0
   argp = 360-argp

nu = acos(dot(evec,r)/(e*mag(r))

if dot(r,v)<0
   nu = 360 - nu

Примечание : это следует из метода, изложенного в «Основах астродинамики и приложений » Вальядо, 2007.

Первый ключ к выяснению этого — правильная система координат. Для таких вещей обычно используются две системы координат. Это фиксированные кадры с центрированием на земле ( ECEF ) и кадры с центрированием на земле (ECI). В полночь эти две точки точно совпадают, но в другое время они расходятся в зависимости от вращения Земли. ECEF лучше всего работает для вещей на Земле (если вы не двигаетесь, у вас должна быть нулевая скорость. ECEF принимает это во внимание, скорость ECI заставит вас двигаться вместе с вращением Земли), ECI лучше всего работает для вещей на орбите (орбитальное движение). объекты не заботятся о вращении Земли, по крайней мере, физике все равно). Убедитесь, что системы координат правильные!

Итак, у вас есть положение и скорость в координатах ECI, что вы делаете? Есть отличная статья , описывающая весь процесс, конечные формулы которой я скопирую сюда. Есть также несколько хороших источников здесь , здесь и здесь . Я настоятельно рекомендую внимательно их прочитать. Неопределенность гораздо сложнее, так что давайте просто предположим, что вы в совершенстве знаете скорость и положение. В частности, 6 классических элементов Кеплариана — это эксцентриситет (е), наклон (i), прямое восхождение восходящего узла ($\Omega$), аргумент перигея ($\omega$), большая полуось (а) и время прохождения перигея ($T_O$).

Я должен упомянуть, что я в первую очередь следую методу определения орбиты Лапласа , существует конкурирующая методология, известная как метод Гаусса. Но в конце концов дело дошло до расшифровки кода Matlab .

Большая полуось

$W_s = \frac{1}{2}*v^2s - \text{mus}./r;$

$a = -mus/2./W_s$; % большая полуось

Эксцентриситет

L = [rs(2,:).*vs(3,:) - rs(3,:).*vs(2,:);...
     rs(3,:).*vs(1,:) - rs(1,:).*vs(3,:);...
     rs(1,:).*vs(2,:) - rs(2,:).*vs(1,:)]; %angular momentum

$p = \sum{L^2}./mus;$% полуширокая прямая кишка

$e = \sqrt{1 - p/a}; %eccentricity$

наклон

$I = atan(\frac{\sqrt{L(1,:)^2 + L(2,:)^2}}{L(3,:)});$

Аргументы Перицентра

$\omega = atan2(\frac{(vs(1,:).*L(2,:) - vs(2,:).*L(1,:))./mus - rs(3,:)./r)./(e.*sin(I))}{((\sqrt{L2s}.*vs(3,:))./mus - (L(1,:).*rs(2,:) - L(2,:).*rs(1,:))./(\sqrt{L2s}.*r))./(e.*sin(I)))}$

Долгота восходящего узла

$\Omega = atan2(-L(2,:),L(1,:));$

Время прохождения перигея:

$T_0 = -(E - e.*sin(E))./\sqrt{mus.*a.^-3}$

OrbitalPy имеет удобную elements_from_state_vectorфункцию, которая делает именно это:

https://github.com/RazerM/orbital/blob/0.7.0/orbital/utilities.py#L252

Вы можете проверить, что математика такая же, как в ответе пользователя 29.