Domanda:
Come calcolare programmaticamente gli elementi orbitali utilizzando i vettori di posizione / velocità?
Stu
2013-09-11 16:49:49 UTC
view on stackexchange narkive permalink

Vorrei creare da zero un software meccanico orbitale. Ritengo che questo sarebbe un ottimo modo per apprendere i passaggi necessari per calcolare diversi elementi orbitali di Keplero di un oggetto, tracciare le orbite e prevedere dove sarà l'oggetto in un momento futuro.

Nello specifico, voglio iniziare con il calcolo degli elementi keplarian. Gli input che darei al programma sarebbero i vettori di posizione e velocità, insieme a un tempo. Questi vettori di input saranno relativi al centro della Terra, quindi potrei anche dover eseguire un trasferimento di coordinate se voglio utilizzare una posizione specifica sulla superficie come punto di riferimento.

Ho visto il matematica per il calcolo degli elementi orbitali di Keplero da questo libro, e so che molti software sono stati sviluppati nel corso degli anni per calcolarli, ma ho difficoltà a colmare i due. La matematica nel libro è leggermente confusa e penso che sarebbe più facile per me capire se vedessi i passaggi "scritti" in un linguaggio di programmazione.

Posso usare R o Python. Codice vettorizzato dalla maggior parte delle lingue dovrei essere in grado di tradurre in uno di questi due. Grazie! @Chris Aggiornerò ancora una volta la domanda con qualche informazione in più sui miei problemi con la traduzione della matematica in codice.
Controlla http://orsa.sourceforge.net/ per le loro soluzioni / metodi. Lunga discussione qui http://www.physicsforums.com/showthread.php?t=232778. E per la simulazione F = ma di base Python: https://fiftyexamples.readthedocs.org/en/latest/gravity.html
FWIW, se il tuo obiettivo è calcolare la posizione futura, normalmente ci sono poche ragioni per convertire il vettore di posizione e velocità in elementi kepleriani. Calcola e applica la resistenza e la forza di gravità nella posizione e velocità correnti e integrali in avanti a piccoli passi.
Tre risposte:
#1
+40
user29
2013-09-12 19:02:44 UTC
view on stackexchange narkive permalink

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.

Infine ricreato in Python. Grazie per l'aiuto! È stato abbastanza facile e ora capisco molto meglio la matematica.
Mi piacerebbe vedere la prova di queste equazioni così come la loro applicazione in un problema reale derivante da un insieme di elementi orbitali. Gli elementi orbitali sono necessari anche se hai questi vettori di posizione e velocità per cominciare? Sembra che il calcolo vettoriale possa gestirlo abbastanza semplicemente. Inoltre, sono confuso dal fatto che non hai definito piccolo mu. Anche il vettore v è la prima derivata del vettore r. L'n-hat è un vettore unitario ma non hai mostrato i calcoli per renderlo unitario. Quali sono i valori con indice neh? Inoltre, non sei riuscito a mostrare come derivare l'anomalia media e il movimento medio
E per calcolare l'anomalia media?
C'è qualche differenza tra nhat e n? Non sono sicuro di cosa sia n (supponendo che nhat sia il vettore del nodo) ma è stato usato per calcolare l'argomento del periapsis. Presumo che n sia uguale a nhat en è stato utilizzato per errore?
Per la programmazione, preferirei relazioni equivalenti usando la funzione atan2, perché ciò eseguirà automaticamente le risoluzioni dei quadranti e ti protegge dalla necessità di controllare l'intervallo prima della maggior parte degli usi di acos sopra. (A volte, la precisione numerica farà sì che l'argomento di una funzione acos sia leggermente oltre l'intervallo valido da -1.00000 a +1.00000. Quando ciò accade, il programma come mostrato andrà in crash.) È anche possibile che h sia zero, consentendo una divisione per zero errori.
Cosa khat nella seconda riga? Vedo nel codice di esempio che è [0 0 1] è il vettore normale del piano equatoriale terrestre?
Un problema è che si tratta di prendere un punto r e v e trasformarlo in elementi in quell'istante, che non saranno i soliti valori medi in un insieme di elementi a due linee. Quello che le persone fanno per generare una TLE da vettori è propagare il vettore in avanti usando l'integrazione numerica, quindi adattare i termini TLE per abbinarlo. Se sei ragionevolmente alto, probabilmente puoi ignorare il termine di trascinamento. Se non hai bisogno di un'altissima precisione, puoi modellare la terra come una massa puntiforme.
Cos'è "a"? Sto guardando [questa immagine di riferimento] (https://user-images.githubusercontent.com/1646875/47543388-f8cf1f00-d8af-11e8-8775-a97e44df246f.png).
Qual è il vettore del nodo "n ^"? Cos'è "K ^"? Cos'è "μ"? Che cos'è "r" e in che modo è diverso da "r"? Che dire di "v" e "v"? Cos'è "p"?
Cosa sono "h_K" e "n_I"?
@MattJessick Come lo faresti usando Atan2?
#2
+6
PearsonArtPhoto
2013-09-11 20:34:20 UTC
view on stackexchange narkive permalink

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

Il metodo di Laplace è quello della determinazione iniziale dell'orbita dalle misurazioni degli angoli e (credo) ben oltre lo scopo di ciò che OP sta cercando. Se hai posizione e velocità ECI, ottenere elementi kepleriani è solo una semplice [trasformazione di coordinate] (http://ccar.colorado.edu/ASEN5070/primers/cart2kep/cart2kep.htm).
@Chris: Sapevo che doveva esserci un modo più semplice per fare la trasformazione ... Sigh.
E il sistema di coordinate in termini di manualità e orientamento? Da quello che posso dire, la maggior parte delle guide utilizza Z-is-up. Come sarebbero implementate queste formule in un sistema per mancini con Y in alto, come Unity, o in un sistema con la mano destra in Y, come Godot?
#3
+5
alexamici
2016-04-08 11:50:18 UTC
view on stackexchange narkive permalink

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.

Ehi, è fantastico! Non vedo l'ora di provarlo. Nel secondo esempio nella documentazione, intitolato "* Create molniya orbit *", OrbitalPy può implementare la precessione? C'è un posto per aggiungere un J2? (e penso che Molniya dovrebbe essere in maiuscolo - penso che si qualifichi come un nome proprio).
Non ho familiarità con OrbitalPy, ma dall'analisi limitata che ho fatto del codice sorgente sembra che faccia pura propagazione kepleriana a due corpi, senza alcuna perturbazione, quindi nessuna correzione ellissoide.
OK, grazie per le informazioni, farò un giro.
È scritto con un sistema di coordinate "Z-up" in mente? Se ho un sistema di coordinate "Y-up", dovrei sostituire "h.z" con "h.y" e "[0, 0, 1]" con "[0, 1, 0]"?
"L'argomento del periapsis è l'angolo tra il vettore di eccentricità e la sua componente x." Perché?


Questa domanda e risposta è stata tradotta automaticamente dalla lingua inglese. Il contenuto originale è disponibile su stackexchange, che ringraziamo per la licenza cc by-sa 3.0 con cui è distribuito.
Loading...