Alcuni dettagli matematici del metodo di Herget

Click here for English version

Il metodo di Herget è in realtà molto semplice; questa è la ragione principale per cui ho deciso di usarlo per determinare l'orbita iniziale nel mio programma FIND_ORB per la determinazione delle orbite. Se qualcuno mi chiedesse di mettermi alla lavagna e spiegare come funziona, oppure di sedermi al computer e scrivere un programma su questo argomento, lo potrei fare senza aprire un libro e senza mettermi le mani nei capelli quando si arriva alle parti più difficili...  questo non avviene per i metodi di Laplace, Olbers, o Gauss. Eppure il mio programma può utilizzare molte osservazioni e tener conto delle perturbazioni orbitali; gli altri programmi utilizzano soltanto tre osservazioni e ipotizzano un'orbita non perturbata.

Il codice sorgente in C/C++ che ho usato per FIND_ORB è disponibile su questo sito; si può fare clic qui per maggiori dettagli. Se siete buoni programmatori con il linguaggio C, analizzare il codice può essere un modo utile per capire il metodo di Herget (e anche per vedere come si effettua la correzione di un'orbita con il metodo dei minimi quadrati, oltre ad altri concetti fondamentali). Tuttavia, la descrizione matematica che segue è probabilmente di maggiore interesse, anche per i programmatori in linguaggio C.

Quasi tutti gli argomenti descritti più avanti sono tratti dall'ottimo libro Fundamentals of Celestial Mechanics di J.M.A. Danby, pubblicato dalla Willmann-Bell.

Il metodo funziona nel modo seguente. Supponiamo di avere un insieme di osservazioni per un certo oggetto, espresse in ascensione retta e in declinazione. Si inizia ipotizzando di conoscere la distanza dell'oggetto per due di quelle osservazioni...  sono i valori R1 e R2 che il programma FIND_ORB richiede. FIND_ORB sceglie anche la prima e l'ultima osservazione, cosa che sembra aiutare il metodo a evitare strani andamenti divergenti.

Ad ogni modo, poiché si conosce la distanza, l'AR e la declinazione per quei due punti e la loro collocazione temporale, si è in grado di sapere esattamente dove si trova l'oggetto in quelle due date. Chiamiamo le due posizioni P1 = (x1, y1, z1) e  P2 = (x2, y2, z2). Con queste ipotesi, esiste una e una sola orbita che possa portare l'oggetto da P1 a P2 in un intervallo di tempo delta_t = (t2 - t1). Il testo di Danby fornisce vari modi per calcolare quale dovrebbe essere una tale orbita. Personalmente, ho trovato quasi tutti quei modi poco chiari, per cui ho escogitato un mio metodo molto semplice. (Sono sicuro, tuttavia, che non si tratta di un'idea originale. Non ho però trovato nessun riferimento bibliografico in proposito).

(1) Si inizia con una velocità iniziale "ipotizzata" V1 = (vx1, vy1, vz1),

   vx1 = (x2 - x1) / delta_t
   vy1 = (y2 - y1) / delta_t
   vz1 = (z2 - z1) / delta_t

Se non ci fossero effetti gravitazionali da parte del Sole e dei pianeti, e l'oggetto si muovesse semplicemente su una linea retta, questa velocità sarebbe esattamente corretta. (Infatti, il programma TWOOBS di Dave Tholen fa sostanzialmente così. Se il valore di delta_t è abbastanza piccolo, diciamo di pochi giorni, questa è un'approssimazione "abbastanza buona"; l'oggetto non precipita tanto rapidamente verso il sole nei tempi successivi. Potete provare a usare prima questa approssimazione, e poi usare i raffinamenti descritti più avanti, in un secondo tempo).

(1b) (Completamente opzionale) Il successivo raffinamento rispetto all'approssimazione precedente, quella cioè di considerare l'ipotesi di un "moto in linea retta", consiste nell'includere gli effetti gravitazionali del primo ordine. A tale scopo va calcolata l'accelerazione dell'oggetto nel punto approssimativamente a metà dell'orbita: (x1+x2)/2, (y1+y2)/2, (z1+z2)/2. Se si assume che questa accelerazione (ax, ay, az) sia molto prossima all'accelerazione costante tra t1 e t2, si può ricavare un valore molto migliore per la velocità iniziale V1:

   vx1 = (x2 - x1) / delta_t - ax1 * delta_t / 2
   vy1 = (y2 - y1) / delta_t - ay1 * delta_t / 2
   vz1 = (z2 - z1) / delta_t - az1 * delta_t / 2

Tutto questo non è assolutamente necessario. Fisicamente, tutto questo somiglia al caso di un tiratore scelto con il suo fucile che vuole sparare un colpo verso un bersaglio lontano e pensa: "Il proiettile si abbasserà di dieci centimetri, a causa della gravità, su una distanza così lunga; mirerò allora dieci centimetri più in alto per compensare tale abbassamento". In questo caso, il tiratore ignora il fatto che la gravità (per questo speciale tipo di "proiettile") cambia in intensità e direzione mentre il proiettile viaggia nell'aria; tuttavia i cambiamenti sono abbastanza piccoli, per cui si ha una buona approssimazione, che è comunque enormemente migliore dell'ipotesi di trascurare completamente la gravità. E' questa una comoda analogia, e ci tornerò sopra più avanti.

Se si include questa grossolana approssimazione per tenere conto della gravità, si otterranno valori accettabili per un intervallo di tempo delta_t di alcune settimane.

(2) Va poi eseguita una integrazione numerica dell'orbita dall'istante t1 all'istante "finale" t2. (Se si include nel calcolo soltanto l'effetto del sole, si può anche usare una soluzione analitica. Occorre però essere in grado di integrare numericamente l'orbita per tenere conto delle perturbazioni. FIND_ORB usa la versione modificata da Fehlberg del metodo di Runge-Kutta, come è descritta nelle pagg. 297-300 del testo Fundamentals of Celestial Mechanics; comunque il tipo di algoritmo che si usa per integrare l'orbita è del tutto irrilevante, purché dia risultati accurati).

(3) Il risultato del calcolo non sarà naturalmente (x2, y2, z2), a causa degli effetti gravitazionali durante l'intervallo di tempo...  l'oggetto in realtà non si muoverà così semplicemente come è stato suggerito fin qui. Supponendo di trovare uno scostamento di una certa entità (delta_x, delta_y, delta_z), occorre impostare una nuova "congettura" per la velocità,

   vx1 = vx1 - delta_x / delta_t
   vy1 = vy1 - delta_y / delta_t
   vz1 = vz1 - delta_z / delta_t

(4) Occorre infine ripetere i passi (2) e (3) finché questi valori delta non diventino sufficientemente piccoli.

Questo metodo ha il vantaggio di essere semplice e fisicamente ovvio. Il passo (3) corrisponde a dire "Il mio primo tentativo è stato 0.13 UA 'troppo alto' e 0.08 UA 'troppo a sinistra'; pertanto il mio prossimo tentativo dovrà essere 'più basso' di 0.13 UA e 'più a destra' di 0.08 UA", anche qui allo stesso modo con cui il tiratore spara un colpo, e poi usa il risultato per aggiustare il tiro. Inoltre l'integrazione numerica può includere gli effetti perturbatori dei pianeti; i metodi più sofisticati descritti da Danby non ne tengono conto. Si potrebbe dire che questo metodo è più "orientato all'aspetto fisico" del problema, mentre gli altri metodi sono più "orientati all'aspetto matematico".

Il metodo ha però lo svantaggio che, se l'arco di orbita considerata è maggiore di 1/8 dell'orbita totale, i valori delta_xyz non convergono, e non vi è soluzione! Recentemente ho modificato il mio programma per tenere conto di questo fatto. A questo punto, gli algoritmi matematici cominciano a diventare un po' più difficili e la connessione con l'aspetto fisico del problema molto meno ovvia. La ragione è che, su lunghi periodi di tempo, la relazione tra i cambiamenti della velocità iniziale e la posizione in cui va a finire l'oggetto non è molto semplice. Nell'approssimazione precedente, avevamo ipotizzato che se si fa andare l'oggetto più veloce, esso arriverà oltre il punto voluto; se si mira al di sopra del bersaglio, il proiettile passerà al di sopra del bersaglio stesso, e così via. Questo cessa di essere vero per orbite più lunghe.

Potrò descrivere in futuro le modifiche necessarie, se qualcuno me lo chiederà. Nel frattempo, la parte di programma che fa questi calcoli si trova nella funzione find_transfer_orbit, nel file del codice sorgente ORB_FUNC.CPP. La parte di codice interessata è lunga una trentina di righe, ma è piuttosto dura da comprendere.

In ogni caso...  sia se vi accontentiate di un metodo semplice ma di modesta accuratezza, oppure vogliate un metodo più complesso ma molto accurato, seguendo i passi descritti sopra, avete determinato l'orbita di un oggetto, partendo dall'istante t1, dalla posizione (x1, y1, z1), e dalla velocità (vx1, vy1, vz1), che si accorda alle vostre osservazioni negli istanti t1 e t2. Nei tempi intermedi vi saranno degli errori in AR e in declinazione, per le osservazioni 1, 2, 3, ... n.  

Chiamiamo questi errori

   delta_ra(i)  = observed_ra(i)  - computed_ra(i)
   delta_dec(i) = observed_dec(i) - computed_dec(i),  i = 1, 2,  ... n.

L'osservazione 1 si accorda con l'istante t1, e l'osservazione n si accorda con l'istante t2, e in questi due istanti l'errore sarà nullo. Pertanto: delta_ra(1) = delta_dec(1) = delta_ra(n) = delta_dec(n) = 0. Viceversa, per tutti gli altri valori di i, si può assumere che le posizioni osservate differiscano da quelle calcolate. 

Questi valori 'delta_ra' e 'delta_dec' sono anche noti come 'errori residui' o anche valori 'O-C' (osservati meno calcolati).

A questo punto, abbiamo un modo per misurare l'errore per ciascun oggetto, e possiamo quindi calcolare un "errore totale" per la nostra congettura dell'orbita dell'oggetto: 

   total_error = sum( delta_ra(i) ^ 2 + delta_dec(i) ^ 2),  i = 1 to n.

Questo valore è generalmente mostrato da FIND_ORB come "errore quadratico medio", o errore RMS:

   rms_error = sqrt( total_error / n)

Il nostro scopo è quello di minimizzare l'errore totale (o, come dire, l'errore RMS). Allora potremo dire di avere un'orbita che "si accorda al meglio" con le osservazioni. 

La procedura usuale per ottenere questo risultato è il metodo dei minimi quadrati. Nella maggior parte dei casi, il metodo dei minimi quadrati è considerato tremendo dai programmatori...  peccato! poiché si tratta di un algoritmo matematico di estrema utilità. (Personalmente ho "incapsulato" questo metodo all'interno di molti algoritmi di calcolo usati nei miei programmi di astrometria, nel calcolo delle orbite, e nell'analisi degli errori delle montature dei telescopi. Tutto questo si trova nel file LSQUARE.CPP che fa parte del codice sorgente di Find_Orb. Penso che vi occorrerà un buon testo di matematica per comprendere i dettagli di tali algoritmi).

In questo caso particolare, tuttavia, si tratta soltanto di ricavare due variabili incognite, R1 e R2. Questo semplifica un po' il problema, anche se non tanto quanto farebbe piacere. La descrizione completa dei dettagli si trova nel testo di Danby.

Dave Tholen e Rob McNaught hanno escogitato ciascuno un diverso metodo per ricavare R1 e R2. Il programma TWOOBS di Dave Tholen (il nome significa "due osservazioni", ed è costituito dalle due sillabe "two" e "obs" delle due parole inglesi "two","due" e "observations","osservazioni") funziona prendendo in considerazione inizialmente un gran numero di valori di R1 e R2. Gli attuali PC sono abbastanza veloci, per cui questo non è un grosso onere. TWOOBS riesce immediatamente a scartare le soluzioni assurdamente iperboliche e (anche se è un po' più rischioso) le soluzioni con moto retrogrado. 

TWOOBS non dà un giudizio sulle restanti orbite. Si ottiene una grande quantità di possibili soluzioni, che, tra l'altro, dà un'idea dell'incertezza delle effemeridi. Rob McNaught ha realizzato un procedimento che ha chiamato "Pangloss", poiché trova la "migliore di tutte le orbite possibili". Anche in questo metodo si ottiene una grande quantità di valori di R1 e R2, ma va considerato quanto sia "ragionevole" un'orbita. Gli asteroidi si suddividono in gruppi e famiglie, per cui si può dire (per esempio) "questo asteroide sembra che sia di un tipo con un'orbita nella fascia principale", oppure "praticamente nessuna orbita conosciuta ha un semiasse maggiore di questo genere...  probabilmente siamo molto lontani da una soluzione corretta".

Attualmente FIND_ORB aggiusta i valori di R1 e R2 per minimizzare la somma dei residui. Spero in futuro di rendere il programma un po' più simile al "Pangloss", minimizzando la somma dei residui a meno di un "ragionevole fattore". Questo fattore sarà zero per un'orbita realmente strana, ma di una certa entità per un'orbita che appaia ragionevole. Se sono disponibili soltanto due osservazioni (nel qual caso, qualunque sia il valore di R1 e R2, la somma dei residui sarà zero), il risultato sarà puramente del tipo "Pangloss". Quando poi saranno disponibili altri dati osservativi, il procedimento sarà una via di mezzo tra il metodo di Herget e quello del "Pangloss". Inoltre, con un arco di orbita abbastanza lungo, il giudizio sulla attendibilità di un'orbita sarà basato sui dati osservativi, e non verranno prese in considerazione idee circa la "ragionevolezza". 

Più in generale, una volta che si abbiano dei "buoni" valori di R1 e R2, e quindi si sia ottenuta un'orbita che corrisponda bene ai dati osservativi, si può passare al metodo generale dei minimi quadrati, quello che in FIND_ORB si attiva con il bottone "Passo completo" (nella versione inglese "Full step"). Con questo metodo, FIND_ORB riaggiusta i valori (x1, y1, z1) e (vx1, vy1, vz1), cercando di nuovo di trovare il minimo per il valore dell'errore totale ("total_error"). Poiché ora non si ipotizza più che la prima e l'ultima osservazione abbiano un errore nullo (cioè, poiché ora si aggiustano sei valori anziché soltanto due, R1 e R2), normalmente si ottiene un errore totale più piccolo rispetto al caso precedente.