2 Elementy orbitalne, równanie Keplera i ślad naziemny

0 Prolog: Gauss odnajduje Ceres

1 stycznia 1801 roku Giuseppe Piazzi w Palermo zauważył nową „gwiazdę”, która przesuwała się między gwiazdami. Obserwował ją przez kilka tygodni, aż zginęła w blasku Słońca. Miał pomiary obejmujące niewielki łuk orbity i nikt nie umiał z nich przewidzieć, gdzie szukać obiektu za kilka miesięcy. Dwudziestoczteroletni Carl Friedrich Gauss wymyślił metodę wyznaczania orbity z trzech obserwacji, policzył elementy orbity i przewidział położenie. W ostatnich dniach roku Franz von Zach odnalazł Ceres prawie dokładnie tam, gdzie wskazał Gauss.

Gauss musiał rozwiązać dokładnie te problemy, którymi zajmuje się ta lekcja: opisać orbitę w przestrzeni kilkoma liczbami i z nich policzyć położenie w dowolnej chwili. Dziś robi to każdy odbiornik GPS dla każdego satelity, kilka razy na sekundę.

Pytanie na start. Stacja kosmiczna ma orbitę nachyloną o $51{,}6^\circ$ do równika. Warszawa leży na szerokości $52{,}2^\circ$ . Czy stacja kiedykolwiek przeleci dokładnie nad Warszawą? Odpowiedź w części 4.

1 Sześć elementów orbitalnych

Orbitę opisujemy w układzie inercjalnym związanym z Ziemią (ECI): początek w środku Ziemi, oś $z$ wzdłuż osi obrotu na północ, oś $x$ w płaszczyźnie równika w kierunku punktu Barana (gdzie Słońce przecina równik niebieski w równonocy wiosennej), oś $y$ dopełnia układ prawoskrętny. Układ się nie obraca razem z Ziemią. Orbitę opisuje sześć liczb, klasycznych elementów orbitalnych:

Klasyczne elementy orbitalne

Kształt i rozmiar: wielka półoś $a$ i mimośród $e$ (lekcja 1 Zagadnienie dwóch ciał i równanie orbity).

Płaszczyzna orbity: nachylenie $i$ (kąt między $\mathbf h$ a osią $z$ , od $0^\circ$ do $180^\circ$ ) i rektascensja węzła wstępującego $\Omega$ (kąt od osi $x$ do punktu, w którym satelita przecina równik, lecąc na północ).

Orientacja elipsy w płaszczyźnie: argument perygeum $\omega$ (kąt od węzła wstępującego do perygeum, w kierunku ruchu).

Położenie na orbicie: anomalia prawdziwa $\nu$ albo anomalia średnia $M$ w danej chwili (część 3).

Przejście z układu peryfokalnego (oś $p$ do perygeum, $q$ w płaszczyźnie orbity, $w$ wzdłuż $\mathbf h$ ) do ECI: $\;Q = R_3(\Omega)\,R_1(i)\,R_3(\omega)$ .

Macierz $Q$ to złożenie trzech obrotów, jak kąty Eulera w lekcji o orientacji samolotu (dynamika lotu, lekcja 1 Orientacja w 3D: macierz obrotu, kąty Eulera, oś i kąt, kwaterniony) i obroty z algebry liniowej (matematyka, lekcja 2.2 Przekształcenia liniowe: obroty, zmiana bazy, transformacja Lorentza): najpierw obrót o $\omega$ wokół normalnej do orbity, potem o $i$ wokół linii węzłów, na koniec o $\Omega$ wokół osi $z$ . Ma też te same osobliwości. Dla orbity równikowej ( $i = 0$ ) linia węzłów nie istnieje i $\Omega$ nie jest określone; dla kołowej ( $e = 0$ ) nie ma perygeum i nie ma $\omega$ . W takich przypadkach używa się kombinacji, które pozostają określone: argumentu szerokości $u = \omega + \nu$ (orbita kołowa) albo długości perygeum $\varpi = \Omega + \omega$ (orbita równikowa).

Zadania · część 1

3 zadania · ok. 20 min

Orbita ma nachylenie $98^\circ$ . W którą stronę względem obrotu Ziemi leci satelita i do jakiej szerokości geograficznej sięga jego ślad?

Odpowiedź
Nachylenie powyżej $90^\circ$ oznacza ruch wsteczny: rzut prędkości na płaszczyznę równika jest skierowany na zachód. Ślad sięga szerokości $180^\circ - 98^\circ = 82^\circ$ . Takie nachylenia mają orbity heliosynchroniczne (lekcja 3 Perturbacje: spłaszczenie Ziemi, orbita heliosynchroniczna i opór).
Które elementy orbitalne nie są określone dla idealnej orbity geostacjonarnej? Czym je zastąpić?

Odpowiedź
Orbita jest kołowa i równikowa: nie istnieje ani perygeum ( $\omega$ ), ani linia węzłów ( $\Omega$ ). Zostają $a$ (i $e = 0$ , $i = 0$ ), a położenie podaje się jedną liczbą, długością prawdziwą $\ell = \Omega + \omega + \nu$ , mierzoną od osi $x$ . Satelity geostacjonarne opisuje się w praktyce elementami „ekwinoktialnymi”, które nie mają tych osobliwości.
Satelita ma $\mathbf r = [-6045;\ -3490;\ 2500]$ km i $\mathbf v = [-3{,}457;\ 6{,}618;\ 2{,}533]$ km/s. Wektor $\mathbf h = \mathbf r\times\mathbf v$ ma składowe $[-25\,385;\ 6670;\ -52\,071]$ km²/s (sprawdź jedną z nich). Policz nachylenie i rektascensję węzła wstępującego.

Odpowiedź
$h = 58\,311$ km²/s; $\cos i = h_z/h = -0{,}893$ , $i = 153{,}2^\circ$ (orbita wsteczna). Wektor węzłów $\mathbf n = \hat{\mathbf z}\times\mathbf h = [-h_y;\ h_x;\ 0] = [-6670;\ -25\,385;\ 0]$ ; $\Omega = \operatorname{atan2}(-25\,385,\ -6670) = 255{,}3^\circ$ (trzecia ćwiartka: sam $\arccos$ dałby $104{,}7^\circ$ ).

2 Od wektorów do elementów i z powrotem

Algorytm z lekcji 1 rozszerzamy o trzy kąty. Oprócz $\mathbf h$ i $\mathbf e$ potrzebny jest wektor węzłów $\mathbf n = \hat{\mathbf z}\times\mathbf h$ , skierowany do węzła wstępującego.

Algorytm $(\mathbf r, \mathbf v)\to(a, e, i, \Omega, \omega, \nu)$

1. $\mathbf h = \mathbf r\times\mathbf v$ , $\;\mathbf n = \hat{\mathbf z}\times\mathbf h$ , $\;\mathbf e = \mathbf v\times\mathbf h/\mu - \mathbf r/r$ , $\;a = -\mu/(2\varepsilon)$ .

2. $\cos i = h_z/h$ .

3. $\cos\Omega = n_x/n$ ; jeśli $n_y < 0$ , to $\Omega = 360^\circ - \Omega$ .

4. $\cos\omega = \mathbf n\cdot\mathbf e/(ne)$ ; jeśli $e_z < 0$ , to $\omega = 360^\circ - \omega$ .

5. $\cos\nu = \mathbf e\cdot\mathbf r/(er)$ ; jeśli $\mathbf r\cdot\mathbf v < 0$ , to $\nu = 360^\circ - \nu$ .

Z powrotem: $\mathbf r_{pf} = \dfrac{p}{1 + e\cos\nu}[\cos\nu;\ \sin\nu;\ 0]$ , $\;\mathbf v_{pf} = \sqrt{\mu/p}\,[-\sin\nu;\ e + \cos\nu;\ 0]$ , $\;\mathbf r = Q\mathbf r_{pf}$ , $\mathbf v = Q\mathbf v_{pf}$ .

Każdy z trzech kroków z warunkiem rozstrzyga ćwiartkę, bo $\arccos$ daje kąty tylko od $0$ do $180^\circ$ . Pominięcie tych warunków to najczęstszy błąd w programach studentów. Dobrym testem jest obieg: elementy $\to(\mathbf r, \mathbf v)\to$ elementy musi zwrócić te same liczby dla orbit we wszystkich ćwiartkach.

Zadania · część 2

3 zadania · ok. 30 min

Dla satelity na orbicie kołowej $e = 0$ , więc $\mathbf e = 0$ . Który krok algorytmu się psuje i co zrobić?

Odpowiedź
Krok 4 i 5 dzielą przez $e$ . Dla orbity kołowej przyjmuje się $\omega = 0$ i liczy argument szerokości $u$ : $\cos u = \mathbf n\cdot\mathbf r/(nr)$ , z warunkiem $r_z < 0 \Rightarrow u = 360^\circ - u$ . W praktyce przy bardzo małym $e$ (np. $10^{-4}$ ) $\omega$ i $\nu$ są osobno źle uwarunkowane, ale ich suma $u$ jest dokładna.
Dla tego samego satelity co w zadaniu 1.3 (dane jak w przykładzie u Curtisa) policz $\varepsilon$ i $a$ .

Odpowiedź
$r = 7414$ km, $v^2 = 11{,}95 + 43{,}80 + 6{,}42 = 62{,}17$ km²/s². $\varepsilon = 31{,}08 - 53{,}76 = -22{,}68$ km²/s², $a = -\mu/(2\varepsilon) = 8788$ km.
Dla tych samych danych $\mathbf e = [-0{,}0916;\ -0{,}1422;\ 0{,}0264]$ (sprawdź, jeśli chcesz). Policz $e$ , $\omega$ i $\nu$ , pilnując ćwiartek.

Odpowiedź
$e = 0{,}171$ . $\cos\omega = \mathbf n\cdot\mathbf e/(ne) = (611 + 3610)/(26\,247\cdot 0{,}1712) = 0{,}939$ , $\omega = 20{,}1^\circ$ ( $e_z > 0$ , więc bez poprawki). $\cos\nu = \mathbf e\cdot\mathbf r/(er) = (554 + 496 + 66)/(0{,}1712\cdot 7414) = 0{,}879$ , $\nu = 28{,}4^\circ$ ( $\mathbf r\cdot\mathbf v = 4133 > 0$ ). Razem z zadaniem 1.3: $a = 8788$ km, $e = 0{,}171$ , $i = 153{,}2^\circ$ , $\Omega = 255{,}3^\circ$ , $\omega = 20{,}1^\circ$ , $\nu = 28{,}4^\circ$ .

3 Czas na orbicie: równanie Keplera

Satelita porusza się po elipsie nierównomiernie: szybko przy perygeum, wolno przy apogeum. Jak znaleźć $\nu$ po czasie $t$ od przelotu przez perygeum? Kepler rozwiązał to geometrycznie. Wokół elipsy opisujemy okrąg o promieniu $a$ i rzutujemy satelitę pionowo na ten okrąg; kąt widziany ze środka elipsy to anomalia mimośrodowa $E$ . Drugie prawo Keplera (pole $\propto$ czas) przekłada się na:

Równanie Keplera

M = E - e\sin E,\qquad M = n(t - t_p),\qquad n = \sqrt{\mu/a^3}.

\tan\frac{\nu}{2} = \sqrt{\frac{1 + e}{1 - e}}\,\tan\frac{E}{2},\qquad r = a(1 - e\cos E).

$M$ (anomalia średnia) rośnie jednostajnie z czasem; $n$ to ruch średni. Dla hiperboli: $M = e\sinh H - H$ .

Z czasu na położenie trzeba rozwiązać równanie Keplera dla $E$ : równanie przestępne, bez rozwiązania w postaci jawnej. Metoda Newtona (matematyka, lekcja 6.2 Metody numeryczne: Python i C++) zbiega tu kwadratowo:

E_{k+1} = E_k - \frac{E_k - e\sin E_k - M}{1 - e\cos E_k} .

Dla $e = 0{,}3$ i $M = 1$ rad: $E_0 = 1$ , $E_1 = 1{,}3013$ , $E_2 = 1{,}28812$ , $E_3 = 1{,}2880913$ , a $E_4$ zgadza się z $E_3$ do 10 cyfr. Start $E_0 = M$ jest dobry dla małych mimośrodów; dla dużych bezpieczniejszy jest $E_0 = \pi$ , bo przy perygeum krzywa $M(E)$ jest prawie płaska i pierwszy krok Newtona potrafi skoczyć daleko.

Zadania · część 3

3 zadania · ok. 30 min

Dla $e = 0{,}3$ i $M = 1$ rad rozwiązanie to $E = 1{,}28809$ rad. Policz $\nu$ i $r/a$ .

Odpowiedź
$\tan(\nu/2) = \sqrt{1{,}3/0{,}7}\tan(0{,}64405) = 1{,}3628\cdot 0{,}7509 = 1{,}0233$ , $\nu = 91{,}3^\circ$ . $r/a = 1 - 0{,}3\cos 1{,}28809 = 0{,}916$ . W 16% okresu ( $M = 1$ rad) satelita przebył ponad ćwierć kąta pełnego: przy perygeum leci szybciej.
Ile trwa lot po orbicie GTO (perygeum 250 km, $a = 24\,396$ km, $e = 0{,}728$ ) od perygeum do $\nu = 90^\circ$ ? Okres wynosi 10,53 h.

Odpowiedź
$\tan(E/2) = \sqrt{0{,}2717/1{,}7283}\tan 45^\circ = 0{,}3965$ , $E = 0{,}7549$ rad; $M = E - e\sin E = 0{,}7549 - 0{,}7283\cdot 0{,}6852 = 0{,}2559$ rad. $t = M/n = 0{,}2559/(1{,}657\cdot10^{-4}) = 1544$ s $= 25{,}7$ min. Ćwierć kąta pełnego w 4% okresu.
Orbita Mołnia: $e = 0{,}74$ , okres 11,96 h. Jaką część okresu satelita spędza w „dalszej połowie” orbity, $90^\circ < \nu < 270^\circ$ ?

Odpowiedź
Czas od perygeum do $\nu = 90^\circ$ : $\tan(E/2) = \sqrt{0{,}26/1{,}74} = 0{,}3866$ , $E = 0{,}7378$ , $M = 0{,}7378 - 0{,}74\cdot 0{,}6727 = 0{,}2400$ rad, czyli $0{,}2400/2\pi = 3{,}8\%$ okresu. Bliższa połowa trwa $2\cdot 3{,}8\% = 7{,}6\%$ , dalsza $92{,}4\%$ : 11,05 h z 11,96 h. Dlatego Mołnia „wisi” nad wysokimi szerokościami północnymi.

4 Ślad naziemny

Ślad naziemny to punkt na powierzchni Ziemi bezpośrednio pod satelitą. Dla orbity kołowej, z argumentem szerokości $u = n(t - t_0)$ liczonym od węzła wstępującego (satelita jest w węźle w chwili $t_0$ ):

\sin\varphi = \sin i\,\sin u,\qquad \lambda = \Omega + \operatorname{atan2}(\cos i\sin u,\ \cos u) - \omega_\oplus(t - t_0) - \lambda_0 ,

gdzie $\varphi$ to szerokość, $\lambda$ długość geograficzna, $\lambda_0$ to kąt między osią $x$ układu ECI a południkiem Greenwich w chwili $t_0$ (czas gwiazdowy Greenwich), a $\omega_\oplus = 7{,}292\cdot10^{-5}$ rad/s to prędkość obrotu Ziemi (względem gwiazd). Wynikają z tego trzy reguły. Ślad sięga szerokości $\pm i$ (albo $\pm(180^\circ - i)$ dla orbit wstecznych). W czasie jednego obiegu Ziemia obraca się o $\omega_\oplus T$ , więc kolejne przejścia przez równik przesuwają się na zachód: dla stacji kosmicznej na wysokości 420 km ( $T = 93{,}0$ min) o $23{,}3^\circ$ . Jeśli w ciągu całkowitej liczby dób gwiazdowych satelita robi całkowitą liczbę obiegów, ślad się powtarza: 15 obiegów na dobę gwiazdową to wysokość 554 km (bez poprawki na spłaszczenie Ziemi z lekcji 3).

Odpowiedź na pytanie z prologu: nie. Ślad stacji sięga szerokości $51{,}6^\circ$ , a Warszawa leży $0{,}6^\circ$ (ok. 66 km) dalej na północ. W najlepszym przelocie stacja jest jednak widoczna prawie w zenicie, ok. $80^\circ$ nad horyzontem (dla kulistej Ziemi; różnica szerokości geodezyjnej i geocentrycznej przesuwa te liczby o ok. 20 km, przykład 6.3).

Ślady prawdziwych satelitów liczy się z elementów publikowanych w formacie TLE (ang. two-line element set). To nie są elementy keplerowskie, tylko „średnie” elementy modelu SGP4, który uwzględnia spłaszczenie Ziemi i opór. Wstawienie ich do równania Keplera daje błędy rzędu kilometrów do dziesiątek kilometrów; do TLE trzeba używać SGP4 (biblioteki są dostępne dla Pythona i C++).

Zadania · część 4

3 zadania · ok. 20 min

Ile obiegów na dobę robi stacja kosmiczna na wysokości 420 km i o ile stopni przesuwa się jej ślad na obieg?

Odpowiedź
$T = 2\pi\sqrt{6798^3/\mu} = 5578$ s $= 92{,}97$ min, $86\,400/5578 = 15{,}49$ obiegu na dobę. Przesunięcie: $\omega_\oplus T = 7{,}292\cdot10^{-5}\cdot 5578 = 0{,}4068$ rad $= 23{,}3^\circ$ na zachód.
Jaka wysokość orbity kołowej daje ślad powtarzający się po 16 obiegach w ciągu doby gwiazdowej? Dlaczego takich orbit nie stosuje się w praktyce?

Odpowiedź
$T = 86\,164/16 = 5385$ s, $a = (\mu T^2/4\pi^2)^{1/3} = 6640$ km, wysokość 262 km. Na takiej wysokości opór atmosfery skraca życie orbity do dni lub tygodni (lekcja 3), więc satelita musiałby stale pracować silnikiem.
Satelita geosynchroniczny ( $T$ równe dobie gwiazdowej) ma nachylenie $10^\circ$ i $e = 0$ . Opisz jego ślad.

Odpowiedź
Szerokość oscyluje między $\pm 10^\circ$ z okresem doby. Długość też oscyluje, bo rzut ruchu na równik nie jest jednostajny: kąt w płaszczyźnie równika spełnia $\tan\lambda' = \cos i\tan u$ , więc $\lambda' - u$ zmienia się z amplitudą ok. $i^2/4 = 0{,}0076$ rad $= 0{,}44^\circ$ , dwa razy na dobę. Ślad to wąska ósemka z przecięciem na równiku, $\pm 10^\circ$ wysokości i ok. $\pm 0{,}44^\circ$ szerokości.

5 Propagacja numeryczna: RK4 kontra metoda symplektyczna

Gdy na satelitę działają siły inne niż grawitacja punktowa (lekcja 3), równania Keplera nie wystarczą i orbitę trzeba całkować numerycznie. Do tego potrzebny jest wybór metody. Porównajmy RK4 (rząd 4, cztery obliczenia przyspieszenia na krok) z metodą skokową (ang. leapfrog, Störmer–Verlet; rząd 2, jedno obliczenie przyspieszenia na krok), która jest symplektyczna: zachowuje strukturę równań Hamiltona, a w praktyce trzyma błąd energii ograniczony na zawsze.

Python · jeden krok RK4 i jeden krok metody skokowej

import numpy as np
mu = 398600.4418
acc = lambda r: -mu * r / np.linalg.norm(r)**3

def krok_rk4(r, v, dt):
    k1r, k1v = v, acc(r)
    k2r, k2v = v + dt/2*k1v, acc(r + dt/2*k1r)
    k3r, k3v = v + dt/2*k2v, acc(r + dt/2*k2r)
    k4r, k4v = v + dt*k3v, acc(r + dt*k3r)
    return (r + dt/6*(k1r + 2*k2r + 2*k3r + k4r),
            v + dt/6*(k1v + 2*k2v + 2*k3v + k4v))

def krok_skokowy(r, v, dt):              # Störmer–Verlet, symplektyczny
    v = v + dt/2*acc(r)
    r = r + dt*v
    v = v + dt/2*acc(r)
    return r, v

C++ · ten sam krok skokowy (g++ -O2 -std=c++17)

#include <array>
#include <cmath>
using V3 = std::array<double, 3>;
const double mu = 398600.4418;

V3 acc(const V3& r) {
    double R = std::sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]), k = -mu / (R*R*R);
    return {k*r[0], k*r[1], k*r[2]};
}
void krok_skokowy(V3& r, V3& v, double dt) {
    V3 a = acc(r);
    for (int j = 0; j < 3; ++j) { v[j] += dt/2*a[j]; r[j] += dt*v[j]; }
    a = acc(r);
    for (int j = 0; j < 3; ++j) v[j] += dt/2*a[j];
}

Test: orbita z $a = 8000$ km i $e = 0{,}2$ , 200 kroków na okres, 1000 okresów (ok. 82 doby). Wyniki:

	RK4	metoda skokowa
największy błąd energii, okresy 1–10	$1{,}8\cdot10^{-7}$	$2{,}6\cdot10^{-4}$
największy błąd energii, okresy 991–1000	$1{,}6\cdot10^{-5}$ (rośnie)	$2{,}6\cdot10^{-4}$ (nie rośnie)
błąd położenia na orbicie po 1000 okresach	$6{,}7^\circ$ (przed)	$251^\circ$ (za)

Metoda symplektyczna rzeczywiście nie pozwala energii dryfować: satelita nie spada i nie ucieka, nawet po bardzo długim czasie (teoria gwarantuje czasy wykładniczo długie względem $1/\Delta t$ ). Ale faza, czyli położenie satelity na orbicie, rozjeżdża się w niej szybciej niż w RK4, bo jej błąd jest drugiego rzędu. Do śledzenia satelity liczy się faza, więc w praktyce używa się metod wysokiego rzędu z adaptacyjnym krokiem (np. Dormanda–Prince’a) albo rozwiązania Keplera z dodanymi poprawkami od zaburzeń. Metody symplektyczne są w swoim żywiole przy badaniu długoterminowej stabilności, np. Układu Słonecznego przez miliony lat.

Zadania · część 5

2 zadania · ok. 25 min

Metoda Eulera $\mathbf r_{k+1} = \mathbf r_k + \Delta t\,\mathbf v_k$ , $\mathbf v_{k+1} = \mathbf v_k + \Delta t\,\mathbf a(\mathbf r_k)$ dla orbity kołowej. Pokaż, że prędkość rośnie w każdym kroku, i powiedz, jak będzie wyglądał numeryczny tor.

Odpowiedź
W pierwszym kroku (start na orbicie kołowej) $\mathbf a\perp\mathbf v$ , więc $|\mathbf v_1|^2 = v^2 + \Delta t^2a^2 > v^2$ i $|\mathbf r_1|^2 = r^2 + \Delta t^2v^2 > r^2$ . Ogólniej: ruch po okręgu to dwa sprzężone oscylatory, a metoda Eulera mnoży amplitudę oscylatora w każdym kroku przez $\sqrt{1 + (n\Delta t)^2} > 1$ . Energia rośnie systematycznie i satelita odlatuje po spirali na zewnątrz. Zamiana kolejności (najpierw prędkość, potem położenie z nową prędkością, czyli półjawny Euler) daje metodę symplektyczną z ograniczonym błędem energii.
Przy równej liczbie obliczeń przyspieszenia metoda skokowa może zrobić cztery razy więcej kroków niż RK4. Czy przy 800 krokach na okres dogoni RK4 z 200 krokami w błędzie fazy ( $251^\circ$ wobec $6{,}7^\circ$ przy 200 krokach)? Oszacuj z rzędów metod.

Odpowiedź
Błąd fazy metody skokowej $\propto\Delta t^2$ : czterokrotnie mniejszy krok zmniejsza go 16 razy, z $251^\circ$ do ok. $16^\circ$ (symulacja: $19^\circ$ ). To wciąż ok. trzy razy więcej niż RK4 z 200 krokami ( $6{,}7^\circ$ ), przy tym samym koszcie obliczeń. Przy wymaganiu wysokiej dokładności metoda wyższego rzędu wygrywa: podwojenie kroków RK4 zmniejsza błąd 16 razy (przy 400 krokach $0{,}2^\circ$ ), metody skokowej tylko 4 razy.

6 Metoda i przykłady rozwiązane

Metoda: gdzie jest satelita w chwili $t$ ?

Elementy $a, e, i, \Omega, \omega$ i chwila przelotu przez perygeum $t_p$ (albo $M_0$ w epoce $t_0$ ).
Anomalia średnia: $M = M_0 + n(t - t_0)$ , sprowadzona do $[0, 2\pi)$ .
Równanie Keplera metodą Newtona: $E$ ; potem $\nu$ i $r = a(1 - e\cos E)$ .
Wektory: $\mathbf r_{pf}$ , $\mathbf v_{pf}$ , a potem $\mathbf r = Q\mathbf r_{pf}$ ; dla śladu jeszcze obrót o $-\omega_\oplus t$ wokół osi $z$ i przejście na szerokość i długość.
Sprawdzenie: obieg elementy $\to$ wektory $\to$ elementy; $|\mathbf r|$ z vis-viva.

Przykład 6.1 GTO godzinę po perygeum

Satelita na GTO ( $r_p = 6628$ km, $r_a = 42\,164$ km) minął perygeum godzinę temu. Gdzie jest i jak szybko leci?

$a = 24\,396$ km, $e = 0{,}7283$ , $n = \sqrt{\mu/a^3} = 1{,}657\cdot10^{-4}$ rad/s; $M = n\cdot 3600 = 0{,}5965$ rad.
Newton od $E_0 = M$ : po 6 krokach $E = 1{,}2978$ rad.
$\nu = 2\arctan\left(\sqrt{1{,}7283/0{,}2717}\tan 0{,}6489\right) = 124{,}8^\circ$ ; $r = a(1 - e\cos E) = 19\,606$ km.

Wynik: po godzinie (9,5% okresu) satelita jest 13 228 km nad Ziemią, $125^\circ$ za perygeum, z prędkością $v = \sqrt{\mu(2/r - 1/a)} = 4{,}93$ km/s.

Sprawdzenie: $r = p/(1 + e\cos\nu) = 11\,455/(1 - 0{,}4157) = 19\,604$ km ✓ (zgodnie z zaokrągleniami).

Przykład 6.2 Wektory na elementy

Zbierz wyniki zadań 1.3, 2.2 i 2.3 w pełny zestaw elementów i sprawdź obiegiem.

$a = 8788$ km, $e = 0{,}1712$ , $i = 153{,}25^\circ$ , $\Omega = 255{,}28^\circ$ , $\omega = 20{,}07^\circ$ , $\nu = 28{,}45^\circ$ .
Z elementów: $p = a(1 - e^2) = 8530$ km, $r = p/(1 + e\cos\nu) = 7414$ km, a $Q\mathbf r_{pf}$ daje z powrotem $[-6045;\ -3490;\ 2500]$ km.

Wynik: orbita wsteczna, lekko eliptyczna, z perygeum 905 km nad Ziemią.

Sprawdzenie: te same liczby podaje Curtis (rozdz. 4) dla tych danych; skrypt a2.py wykonuje obieg elementy $\to$ wektory $\to$ elementy z dokładnością do $10^{-4}$ ✓.

Przykład 6.3 Stacja nad Polską

Stacja kosmiczna: 420 km, $i = 51{,}64^\circ$ . Pod jakim kątem nad horyzontem widać ją z Warszawy ( $52{,}23^\circ$ N) w najlepszym przelocie?

Najbliższy punkt śladu jest $0{,}59^\circ$ na południe, czyli $s = 0{,}59^\circ\cdot\frac{\pi}{180}\cdot 6378 = 66$ km po powierzchni.
Kąt środkowy $\theta = 0{,}59^\circ$ ; satelita jest na promieniu $r_s = 6798$ km. Względem obserwatora: w poziomie $r_s\sin\theta = 70$ km, w pionie $r_s\cos\theta - R = 420$ km.
Elewacja $= \arctan(420/70) = 80{,}5^\circ$ .

Wynik: ok. $80^\circ$ : prawie w zenicie, ale nigdy dokładnie.

Sprawdzenie: rachunek dla kulistej Ziemi; szerokość geograficzna Warszawy jest geodezyjna, a ślad liczony geocentrycznie, co na tej szerokości zmienia odległość o ok. 20 km. Wniosek się nie zmienia ✓.

7 Błędne intuicje, których trzeba się oduczyć

„Satelita krąży nad równikiem, a nachylenie mówi, jak bardzo się od niego oddala.”: Płaszczyzna każdej orbity (bez zaburzeń) przechodzi przez środek Ziemi, więc rzut toru na sferę to koło wielkie: orbita przecina płaszczyznę równika dwa razy na obieg. Nachylenie mówi, jak daleko na północ i południe sięga ślad.
„Ślad naziemny przesuwa się, bo orbita się obraca.”: Orbita (bez zaburzeń) stoi w przestrzeni; obraca się Ziemia pod nią. Spłaszczenie Ziemi rzeczywiście obraca orbitę, ale dużo wolniej (lekcja 3).
„Anomalia średnia to średni kąt.”: To miara czasu w radianach: kąt, jaki zakreśliłby satelita poruszający się jednostajnie z tym samym okresem. Zgadza się z $\nu$ tylko w perygeum, apogeum i na orbicie kołowej.
„Metoda symplektyczna jest dokładniejsza.”: Jest wierniejsza jakościowo (energia nie dryfuje), ale jej błąd fazy może być dużo większy niż w metodzie wyższego rzędu.
„Elementy z TLE wstawiam do równania Keplera.”: To elementy średnie modelu SGP4; bez SGP4 wynik jest błędny o kilometry do dziesiątek kilometrów.

8 Zadania do przemyślenia

Te zadania wymagają więcej myślenia niż liczenia.

Zadania do przemyślenia

5 zadań · ok. 1,5 godziny

Zbieżność Newtona. Dlaczego metoda Newtona dla równania Keplera zawsze zbiega przy starcie $E_0 = \pi$ , a przy $E_0 = M$ i $e$ bliskim 1 potrafi skakać?
Rozwiązanie
1. $f(E) = E - e\sin E - M$ ma $f' = 1 - e\cos E > 0$ (rosnąca) i $f'' = e\sin E$ , które na $(0, \pi)$ ma stały znak. Funkcja rosnąca i wypukła na przedziale z pierwiastkiem: start z tej strony, gdzie $f$ i $f''$ mają ten sam znak, daje ciąg monotoniczny, bez przeskoków. Dla $M \in (0, \pi)$ to $E_0 = \pi$ .
2. Przy $E_0 = M$ blisko zera i $e \to 1$ : $f'(M) = 1 - e\cos M$ jest prawie zerem, więc pierwszy krok $-f/f'$ jest olbrzymi i może wylądować daleko za pierwiastkiem.
Czas do tego samego miejsca. Satelita na orbicie o okresie 92,97 min przelatuje nad punktem na równiku. Kiedy następny raz przeleci dokładnie nad nim, pomijając spłaszczenie Ziemi?
Rozwiązanie
1. Po $k$ obiegach ślad przesunął się o $k\cdot 23{,}31^\circ$ . Dokładny powrót nad ten sam punkt (w tym samym węźle) wymaga, żeby $k\cdot 23{,}31^\circ$ było wielokrotnością $360^\circ$ . Przy okresie niewspółmiernym z dobą to się nigdy nie zdarza: ślad z czasem gęsto pokrywa cały pas szerokości $\pm i$ .
2. Prawie-powtórzenia: $360^\circ/23{,}31^\circ = 15{,}45$ ; najlepsze bliskie całkowitym wielokrotności to $k = 139$ obiegów (ok. 9 dób, błąd $0{,}43^\circ$ , czyli ok. 50 km na równiku) i $k = 278$ (18 dób, $0{,}87^\circ$ ). Satelity obserwacyjne dobiera się tak, żeby powtarzalność była dokładna (np. 15 obiegów w dobie gwiazdowej), z poprawką na spłaszczenie Ziemi z lekcji 3.
Dlaczego atan2. W kroku 3 algorytmu zamiast $\arccos$ z warunkiem można użyć $\Omega = \operatorname{atan2}(n_y, n_x)$ . Czy podobnie da się policzyć $\nu$ i $\omega$ ? Co jest zaletą atan2 numerycznie?
Rozwiązanie
1. Tak: $\nu = \operatorname{atan2}\big((\mathbf e\times\mathbf r)\cdot\hat{\mathbf h},\ \mathbf e\cdot\mathbf r\big)$ i analogicznie $\omega$ z $\mathbf n$ i $\mathbf e$ . Sinus ze znakiem bierze się z iloczynu wektorowego rzutowanego na normalną.
2. $\arccos$ ma nieskończoną pochodną przy $\pm 1$ : dla kątów bliskich $0^\circ$ i $180^\circ$ mały błąd zaokrąglenia argumentu daje duży błąd kąta (a argument nieco powyżej 1 daje NaN). atan2 jest dobrze uwarunkowany dla wszystkich kątów.
Mołnia a GEO. Rosja ma większość terytorium na szerokościach powyżej $55^\circ$ . Wyjaśnij, czemu do łączności używa orbit Mołnia zamiast geostacjonarnej, i co jest ceną.
Rozwiązanie
1. Satelita GEO jest nad równikiem: z szerokości $70^\circ$ widać go tylko kilkanaście stopni nad horyzontem albo wcale. Mołnia z apogeum nad półkulą północną spędza 92% okresu w dalszej połowie orbity, wysoko na niebie nad Rosją.
2. Cena: satelita się porusza, więc anteny muszą go śledzić, a ciągłość łączności wymaga 3–4 satelitów na przesuniętych orbitach. Orbita przecina też pasy radiacyjne, a perygeum jest niskie. Argument perygeum $270^\circ$ wybiera się nie bez powodu: przy $i = 63{,}4^\circ$ spłaszczenie Ziemi go nie obraca (lekcja 3).
Trzy obserwacje Gaussa. Obserwacja z Ziemi daje tylko kierunek (dwa kąty), nie odległość. Dlaczego do wyznaczenia orbity wystarczą trzy obserwacje?
Rozwiązanie
1. Orbita ma sześć niewiadomych (elementy, w tym czas przelotu przez perygeum). Każda obserwacja daje dwa kąty i znany czas: trzy obserwacje to sześć równań.
2. Równania są silnie nieliniowe, dlatego Gauss sprowadził je do jednego równania na odległość w środkowej obserwacji, rozwiązywanego iteracyjnie. Przy krótkim łuku orbity układ jest źle uwarunkowany: małe błędy pomiaru dają duże błędy elementów, co Gauss łagodził metodą najmniejszych kwadratów, kiedy przybywało obserwacji.

9 Sprawdź się

Odpowiedz w myślach, zanim rozwiniesz odpowiedź.

1 Wymień sześć klasycznych elementów orbitalnych.

$a$ , $e$ , $i$ , $\Omega$ , $\omega$ i $\nu$ (albo $M$ , albo czas przelotu przez perygeum).

2 Kiedy

\Omega

nie jest określone? A

\omega

$\Omega$ dla orbit równikowych ( $i = 0$ lub $180^\circ$ ); $\omega$ dla kołowych ( $e = 0$ ).

3 Napisz równanie Keplera.

$M = E - e\sin E$ , $M = n(t - t_p)$ .

4 Do jakiej szerokości sięga ślad orbity o nachyleniu

i

$\pm i$ (dla $i > 90^\circ$ : $\pm(180^\circ - i)$ ).

5 Dlaczego ślad przesuwa się na zachód?

Ziemia obraca się na wschód o $\omega_\oplus T$ w czasie jednego obiegu.

6 Co zachowuje metoda symplektyczna, a czego nie?

Trzyma ograniczony błąd energii przez bardzo długi czas; nie gwarantuje małego błędu fazy.

10 Karta wzorów i co dalej

Elementy i wektory

$\mathbf n = \hat{\mathbf z}\times\mathbf h$ , $\cos i = h_z/h$

$\cos\Omega = n_x/n$ ( $n_y < 0 \Rightarrow 360^\circ - \Omega$ )

$\cos\omega = \mathbf n\cdot\mathbf e/(ne)$ ( $e_z < 0 \Rightarrow 360^\circ - \omega$ )

$\cos\nu = \mathbf e\cdot\mathbf r/(er)$ ( $\mathbf r\cdot\mathbf v < 0 \Rightarrow 360^\circ - \nu$ )

$\mathbf r = Q\mathbf r_{pf}$ , $Q = R_3(\Omega)R_1(i)R_3(\omega)$

Czas i ślad

$M = n(t - t_p) = E - e\sin E$ , $n = \sqrt{\mu/a^3}$

$\tan\frac\nu2 = \sqrt{\frac{1 + e}{1 - e}}\tan\frac E2$ , $r = a(1 - e\cos E)$

Newton: $E_{k+1} = E_k - \frac{E_k - e\sin E_k - M}{1 - e\cos E_k}$

$\sin\varphi = \sin i\sin u$ , przesunięcie śladu $\omega_\oplus T$ na obieg

$\omega_\oplus = 7{,}292\cdot10^{-5}$ rad/s, doba gwiazdowa 86 164 s

Do przeczytania

Obowiązkowo: H. D. Curtis, Orbital Mechanics for Engineering Students, rozdz. 3 (położenie w funkcji czasu, równanie Keplera) i 4 (orbity w przestrzeni, elementy, ślad naziemny).
Dla chętnych: E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration (Springer), rozdz. I (eksperymenty z problemem Keplera) i IX (wyjaśnienie, dlaczego metody symplektyczne zachowują energię); D. A. Vallado, Fundamentals of Astrodynamics and Applications: algorytmy z wszystkimi przypadkami szczególnymi i SGP4.

Elementy orbitalne są stałe tylko w idealnym świecie dwóch punktów. Ziemia jest spłaszczona, ma atmosferę, a Słońce i Księżyc też ciągną. Kolejna lekcja pokazuje, jak te zaburzenia powoli zmieniają elementy i jak inżynierowie zamieniają je w narzędzie: orbita heliosynchroniczna istnieje tylko dlatego, że Ziemia nie jest kulą.