Dati i vettori $ \ vec {r} $ e $ \ vec {v} $ di posizione e velocità inerziali (ECI) centrati sulla Terra, è possibile risolvere direttamente gli elementi orbitali classici $ (a, e, i, \ Omega, \ omega, \ nu) $ come segue (algoritmi prima, seguito da pseudocodice in basso):

Per prima cosa, risolvi il momento angolare $$ \ vec {h} = \ vec {r} \ times \ vec {v} $$

quindi il vettore del nodo $$ \ hat {n} = \ hat {K} \ times \ vec {h} $$ che verrà utilizzato in seguito.

Il vettore di eccentricità è quindi $$ \ vec {e} = \ frac {(v ^ 2- \ mu / r) \ vec {r} - (\ vec {r} \ cdot \ vec {v} ) \ vec {v}} {\ mu} $$

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

L'energia meccanica specifica è $$ E = \ frac {v ^ 2} {2} - \ frac {\ mu} {r} $$

Se $ e \ neq 1 $, allora $$ a = - \ frac {\ mu} {2E} $$$$ p = a (1-e ^ 2) $$ Altrimenti, $$ p = \ frac {h ^ 2} {\ mu} $$ $$ a = \ infty $$

Ora, $$ 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}} $$

E dovrai effettuare i seguenti controlli: If $ n_J<0 $, allora $ \ Omega = 360 ^ {\ circ} - \ Omega $,

Se $ e_K<0 $, allora $ \ omega = 360 ^ {\ circ} - \ omega $ e

Se $ \ vec {r} \ cdot \ vec {v} <0 $, quindi $ \ nu = 360 ^ {\ circ} - \ nu $.

Nota che incontrerai problemi (singolarità) per alcuni casi: orbite circolari ($ e \ circa 0 $) e orbite equatoriali ($ i \ circa 0 $), in particolare. In questi casi normalmente si introduce una nuova variabile meno problematica, come longitudine media o longitudine reale del perigeo.

  h = cross (r, v) nhat = cross ([0 0 1], h) evec = ((mag (v) ^ 2-mu / mag (r)) * r-punto (r, v) * v) / mue = 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 = infi = acos (h (3) / mag (h)) Omega = acos (n (1) / mag (n)) if n (2) <0 Omega = 360-Omegaargp = acos (punto (n, evec) / (mag ( n) * e)) if e (3) <0 argp = 360-argpnu = acos (punto (evec, r) / (e * mag (r)) if punto (r, v) <0 nu = 360 - nu  codice> 

Nota : questo deriva dal metodo esposto in Fundamentals of Astrodynamics and Applications di Vallado, 2007.

La prima chiave per capirlo è ottenere il corretto sistema di coordinate. Ci sono due sistemi di coordinate comunemente usati per queste cose. Sono i frame Earth Centered Earth Fixed (ECEF) e Earth Centered Interial (ECI). A mezzanotte, questi due si allineano esattamente, ma divergono in altri tempi, in base alla rotazione della Terra. ECEF funziona meglio per le cose sulla Terra (se non ti muovi, dovresti avere una velocità 0. ECEF ne tiene conto, la velocità ECI ti farà muovere con la rotazione della Terra), ECI funziona meglio per le cose in Orbita (Orbita agli oggetti non interessa la rotazione della Terra, almeno, alla fisica non interessa). Assicurati che i sistemi di coordinate siano corretti!

Ok, quindi hai una posizione e una velocità in coordinate ECI, cosa fai? C'è un documento eccellente che descrive l'intero processo, di cui copierò qui le formule finali. Ci sono anche alcune buone fonti qui, qui e qui. Consiglio vivamente di leggerli attentamente. L'incertezza è molto più difficile, quindi supponiamo che tu abbia una perfetta conoscenza della velocità e della posizione. Nello specifico, i 6 elementi keplarian classici sono l'eccentricità (e), l'inclinazione (i), l'ascensione retta del nodo ascendente ($ \ Omega $), l'argomento del perigeo ($ \ omega $), il semiasse maggiore (a) e tempo di passaggio perigeo ($ T_O $).

Devo dire che sto principalmente seguendo il Metodo Laplace di determinazione orbitale, esiste una metodologia concorrente nota come Metodo Gauss. Ma alla fine si è trattato di decifrare il codice Matlab.

Semiasse maggiore

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

$ a = -mus / 2. / W_s $; % semiasse maggiore

Eccentricità

  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, :)]; % momento angolare  

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

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

Inclinazione

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

Argomenti di 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)))} $

Longitudine del nodo ascendente

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

Tempo di passaggio del Perigeo:

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

OrbitalPy ha una pratica funzione elements_from_state_vector che fa proprio questo:

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

Puoi verificare che la matematica sia la stessa della risposta user29.