Mając na Ziemi, bezwładnościowe (ECI) położenie i wektory prędkości $ \ vec {r} $ i $ \ vec {v} $, możesz bezpośrednio obliczyć klasyczne elementy orbity $ (a, e, i, \ Omega, \ omega, \ nu) $ w następujący sposób (najpierw algorytmy, a na dole pseudokod):

Najpierw znajdź moment pędu $$ \ vec {h} = \ vec {r} \ times \ vec {v} $$

to wektor węzłowy $$ \ hat {n} = \ hat {K} \ times \ vec {h} $$, który zostanie użyty później.

Wektorem mimośrodu jest zatem $$ \ vec {e} = \ frac {(v ^ 2- \ mu / r) \ vec {r} - (\ vec {r} \ cdot \ vec {v} ) \ vec {v}} {\ mu} $$

i $ e = | \ vec {e} | $.

Specyficzna energia mechaniczna wynosi $$ E = \ frac {v ^ 2} {2} - \ frac {\ mu} {r} $$

Jeśli $ e \ neq 1 $, to $$ a = - \ frac {\ mu} {2E} $$$$ p = a (1-e ^ 2) $$ W przeciwnym razie $$ p = \ frac {h ^ 2} {\ mu} $$ $$ a = \ infty $$

Teraz $$ 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}} $$

I będziesz musiał dokonać następujących kontroli: Jeśli $ n_J<0 $, to $ \ Omega = 360 ^ {\ circ} - \ Omega $,

Jeśli $ e_K<0 $, to $ \ omega = 360 ^ {\ circ} - \ omega $ i

Jeśli $ \ vec {r} \ cdot \ vec {v} <0 $, a następnie $ \ nu = 360 ^ {\ circ} - \ nu $.

Zauważ, że napotkasz problemy (osobliwości) z pewnymi przypadki: szczególnie orbity kołowe ($ e \ approx 0 $) i orbity równikowe ($ i \ ok. 0 $). W takich przypadkach zwykle wprowadzasz nową, mniej kłopotliwą zmienną, taką jak średnia długość lub rzeczywista długość geograficzna perygeum.

  h = cross (r, v) nhat = cross ([0 0 1], h) evec = ((mag (v) ^ 2-mu / mag (r)) * r-dot (r, v) * v) / mue = mag (evec) energia = mag (v) ^ 2 / 2- mu / mag (r) if abs (e-1.0) >eps a = -mu / (2 * energia) p = a * (1-e ^ 2) else p = mag (h) ^ 2 / mu a = infi = acos (h (3) / mag (h)) Omega = acos (n (1) / mag (n)) if n (2) <0 Omega = 360-Omegaargp = acos (dot (n, evec) / (mag ( n) * e)) if e (3) <0 argp = 360-argpnu = acos (dot (evec, r) / (e * mag (r)) if dot (r, v) <0 nu = 360 - nu  kod> 

Uwaga : wynika to z metody przedstawionej w Fundamentals of Astrodynamics and Applications , Vallado, 2007.

Pierwszym kluczem do zrozumienia tego jest poprawienie układu współrzędnych. W takich przypadkach istnieją dwa powszechnie używane układy współrzędnych. Są to ramy zwierciadlane (ECEF) i Earth Centered Earth Fixed (ECI). O północy te dwie dokładnie ustawiają się w jednej linii, ale w innych momentach różnią się w zależności od obrotu Ziemi. ECEF działa najlepiej w przypadku obiektów na Ziemi (jeśli się nie poruszasz, powinieneś mieć prędkość 0. ECEF bierze to pod uwagę, prędkość ECI sprawi, że będziesz się poruszać wraz z obrotem Ziemi), ECI działa najlepiej w przypadku obiektów na orbicie (Orbitowanie obiekty nie przejmują się obrotem Ziemi, a przynajmniej fizyka to nie obchodzi). Upewnij się, że układy współrzędnych są poprawne!

OK, więc masz pozycję i prędkość we współrzędnych ECI, co robisz? Istnieje doskonały artykuł, który opisuje cały proces, którego końcowe formuły powielę tutaj. Jest też kilka dobrych źródeł tutaj, tutaj i tutaj. Gorąco polecam uważne ich przeczytanie. Niepewność jest znacznie trudniejsza, więc załóżmy, że masz doskonałą wiedzę o prędkości i pozycji. Konkretnie, 6 klasycznych elementów keplariańskich to mimośrodowość (e), nachylenie (i), rektascensja węzła wstępującego ($ \ Omega $), argument perygeum ($ \ omega $), półoś wielka (a) i czas przejścia perygeum ($ T_O $).

Powinienem wspomnieć, że przede wszystkim stosuję metodę Laplace'a określania orbity, istnieje konkurencyjna metodologia znana jako metoda Gaussa. Ale ostatecznie sprowadzało się to do rozszyfrowania kodu Matlab.

Półoś główna oś

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

$ a = -mus / 2. / W_s $; % półoś wielkiej

Mimośrodowość

  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, :)]; % momentu pędu  

$ p = \ sum {L ^ 2} ./ mus; $% semi-latus rectum

$ e = \ sqrt {1 - p / a}; % ekscentryczności $

Nachylenie

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

Argumenty Pericenter

$ \ 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)))} $

Długość węzła wstępującego

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

Czas przejścia Perygeum:

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

OrbitalPy ma przydatną funkcję elements_from_state_vector , która właśnie to robi:

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

Możesz sprawdzić, czy matematyka jest taka sama jak w odpowiedzi user29.