Calcolo della traiettoria geocentrica di un
METEOROIDE
tramite il metodo vettoriale di DUBYAGO
e il teorema di LAMBERT

(Giuseppe Matarazzo - Agosto 2000)

Scopo di questa trattazione è quello di determinare tutti i parametri necessari a individuare l'orbita geocentrica di un meteoroide fotografato da due stazioni terrestri con l'applicazione di un noto metodo vettoriale e l'utilizzo di una moderna tecnica astrodinamica.
Quando in epoche passate gli elaboratori elettronici erano una chimera per i comuni mortali, ogni tentativo di risolvere il problema per via vettoriale appariva velleitario, ma adesso, con i personal computer a disposizione di chiunque, le cose si sono radicalmente invertite facendo sembrare sorpassati i metodi trigonometrici di un tempo e le necessarie mappe stellari a proiezione centrografica (gnomonica).
Il metodo di calcolo qui esposto è quello descritto nel libro [1] dell'autore russo A.D.Dubyago, e da qui in avanti verrα chiamato, per comoditα lessicale, METODO DUBYAGO. Esso determina le coordinate geocentriche di due punti qualsiasi della scia luminosa della meteora, da cui, con facili passaggi algebrici, si arriva a conoscere l'altezza rispetto alla superficie terrestre, la distanza e l'azimut dalle due stazioni osservative e il radiante di provenienza del corpo. Alcune apparecchiature di ripresa (fotocamere a otturatore rotante) consentono anche di conoscere con molta precisione la durata della traccia luminosa (o di parte di essa), per cui l'orbita geocentrica del meteorite è completamente determinata sfruttando il teorema di Lambert, applicato in astronautica negli omonimi 'piani di guida'. Questo ci dice che gli elementi orbitali di una conica sono univocamente determinati quando si conoscono i due raggi vettori r1 ed r2 di un corpo, l'angolo tra essi compreso e il tempo di percorrenza dell'arco tra i due punti; si calcolano anche, con i noti algoritmi della meccanica celeste, i vettori velocitα agli estremi dell'arco considerato, cosicchè il panorama completo delle informazioni sul meteoroide è perfettamente delineato. Inoltre, poichè sappiamo che dinamicamente la sua traiettoria è, nell'ambito della sfera di influenza terrestre, di tipo iperbolico, i dati vettoriali consentono anche di ricavare rapidamente la direzione dell'asintoto e di conoscere il radiante geocentrico VERO del meteorite. Quando però le apparecchiature fotografiche non consentono di determinare la durata della traccia luminosa, nel programma al computer [2] si introduce un suo valore arbitrario, per esempio 1 secondo, e il radiante vero, insieme alle velocitα iniziale e finale del tratto di traccia considerata, viene ugualmente calcolato anche se, com'è ovvio, in maniera imprecisa.
Ma la grande VERSATILITA' del metodo DUBYAGO, che lo fa sicuramente preferire agli altri in circolazione, è la possibilitα di determinare l'orbita geocentrica del meteoroide anche quando i quattro punti scelti sulle due tracce fotografiche NON sono a due a due OMOLOGHI. Ciò è molto importante in quanto elimina una delle INCERTEZZE maggiori che si presentano all'atto dell'impostazione del problema: quella di non poter individuare chiaramente i punti corrispondenti a causa della frequente labilitα delle tracce fotografiche. Con il metodo esposto, che è ovviamente applicabile a punti omologhi (inizio e fine scia per esempio), si arriva a determinare SEMPRE gli elementi orbitali dell'iperbole percorsa dal meteoroide, consentendo di avere i parametri necessari per il calcolo eventuale della sua orbita eliocentrica; procedura quest'ultima che esula dalla presente trattazione.

Il metodo DUBYAGO

Nella figura riportata accanto A₁ e A₂ sono i punti scelti sulla prima traccia del corpo vista dalla stazione S₁ con d₁₁ e d₁₂ le due distanze da S₁, mentre analogamente B₁ e B₂ sono i punti della seconda traccia visti dalla stazione S₂ con d₂₁ e d₂₂ le rispettive distanze da S₂. Dalle tracce fotografiche, simili a quelle dell'esempio riportato pi· avanti, si ricavano con il procedimento ivi descritto le coordinate equatoriali topocentriche dei 4 punti:
$\begin{array} {rlcrl} 1^a traccia: & A_1(\alpha_{11},\delta_{11}) & & 2^a traccia: & B_1(\alpha_{21},\delta_{21}) \\ & A_2(\alpha_{12},\delta_{12}) & & & B_2(\alpha_{22},\delta_{22}) \end{array}$ Siano inoltre (A,D) le coordinate equatoriali del radiante del meteorite, S il tempo siderale a Greenwich all'istante dello scatto delle foto e (x₁,y₁,z₁),(x₂,y₂,z₂) le coordinate rettangolari geocentriche delle due stazioni di rispresa S₁ ed S₂.

I coseni direttori delle 4 direzioni (S1-A1, S1-A2, S2-B1, S2-B2) e del radiante (a,b,c) sono (cfr.[3]):
$\begin{array} {lcl} a_{ik}=\cos\delta_{ik}\cdot \cos\alpha_{ik} & & a=\cos D\cdot \cos A \\ b_{ik}=\cos\delta_{ik}\cdot \sin\alpha_{ik} & & b=\cos D\cdot \sin A \\ c_{ik}=\sin\delta_{ik} & & c=\sin D \\ con \; i=1,2; \quad k=1,2 & & \end{array}$ Note le coordinate geografiche delle due stazioni S₁(λ₁,φ₁) e S₂(λ₂,φ₂), si calcolano quelle rettangolari geocentriche mediante le formule, sempre riprese da [3]: $\begin{array} {ll} x_i=RT\cdot \cos\phi_i\cdot \cos(\lambda_i+S) & essendo\; RT\; il\; raggio\; LOCALE \\ y_i=RT\cdot \cos\phi_i\cdot \sin(\lambda_i+S) & della\; sfera\; terrestre\; ed\; i=1,2 \\ z_i=RT\cdot \sin\phi_i & \end{array}$ Il fatto di considerare sferica la Terra è una semplificazione accettabile ai fini di questo studio, così come supporre ininfluente l'attrito dell'atmosfera per quote non inferiori a una settantina di km; l'orbita che ne verrα fuori sarα pertanto quella classica di tipo kepleriano.
Come raggio terrestre si è assunto quello della sfera locale tangente, nel punto di osservazione di latitudine í, all'ellissoide di Bessel, ma con i parametri dell'ellisse del meridiano corretti secondo le indicazioni date nel 1976 dall'Unione Astronomica Internazionale (IAU); esso si ricava con l'espressione: $\begin{array} {ll} RT= a \cdot \sqrt{\frac{1-e^2}{1-e^2\cdot\sin^2\phi}} & con\; a=6378.140\; km;\quad e=0.08181922 \end{array}$ Il valore di RT così calcolato corrisponde ad 1 Unitα Canonica di lunghezza, mentre in base al noto parametro gravitazionale μ=398600.5 km³/sec² l'Unitα Canonica di tempo(TU) è uguale a 806.81153 sec e l'Unitα di velocitα, derivata dalle precedenti, è pari a 1 RT/TU= 7.9053644 km/sec.
Imponiamo adesso le condizioni di appartenenza allo stesso piano sia delle direzioni (A1-A2-R, S1-A1, S1-A2) relative alla stazione S₁ che a quelle (B1-B2-R, S2-B1, S2-B2) riguardanti la S₂. Si ottengono le 2 equazioni: $\left | \begin{array} {ccc} a & b & c \\ a_{11} & b_{11} & c_{11} \\ a_{12} & b_{12} & c_{12} \\ \end{array} \right | = 0 \qquad \left | \begin{array} {ccc} a & b & c \\ a_{21} & b_{21} & c_{21} \\ a_{22} & b_{22} & c_{22} \\ \end{array} \right | = 0 $$ \begin{center} che, sviluppate, formano il seguente sistema: \end{center} $$ \left\{ \begin{array} {cc} l_1\cdot a + m_1\cdot b + n_1\cdot c & = 0 \\ l_2\cdot a + m_2\cdot b + n_2\cdot c & = 0 \\ \end{array} \right. \qquad [1] $$ \begin{center} essendo $(l_1,m_1,n_1)$ i minori della $1^a$ matrice e \\ $(l_2,m_2,n_2)$ i minori della $2^a$, cio\`e \end{center} $$ l_1= \left | \begin{array} {cc} b_{11} & c_{11} \\ b_{12} & c_{12} \\ \end{array} \right | \quad m_1= \left | \begin{array} {cc} c_{11} & a_{11} \\ c_{12} & a_{12} \\ \end{array} \right | \quad n_1= \left | \begin{array} {cc} a_{11} & b_{11} \\ a_{12} & b_{12} \\ \end{array} \right | $$ $$ l_2= \left | \begin{array} {cc} b_{21} & c_{21} \\ b_{22} & c_{22} \\ \end{array} \right | \quad m_2= \left | \begin{array} {cc} c_{21} & a_{21} \\ c_{22} & a_{22} \\ \end{array} \right | \quad n_2= \left | \begin{array} {cc} a_{21} & b_{21} \\ a_{22} & b_{22} \\ \end{array} \right |$ Per risolvere il sistema [1] dividiamo tutto per c e calcoliamo i rapporti a'=a/c e b'=b/c, dopo aver aggiunto una terza equazione, quella riguardante una proprietα dei coseni direttori: a²+b²+c²=1. Con alcuni passaggi algebrici si ottengono i coseni direttori (a,b,c) del radiante: $\begin{array} {ccc} a=a'/\sqrt{a'^2+b'^2+1} & b=b'/\sqrt{a'^2+b'^2+1} & c=c'/\sqrt{a'^2+b'^2+1} \\ \end{array} $$ \begin{center} con \end{center} $$ a'= \frac{ \left | \begin{array} {cc} -n_1 & m_1 \\ -n_2 & m_2 \\ \end{array} \right |} { \left | \begin{array} {cc} l_1 & m_1 \\ l_2 & m_2 \\ \end{array} \right |} \qquad b'= \frac{ \left | \begin{array} {cc} l_1 & -n_1 \\ l_2 & -n_2 \\ \end{array} \right |} { \left | \begin{array} {cc} l_1 & m_1 \\ l_2 & m_2 \\ \end{array} \right |}$ Quindi l'equazione della traccia rettilinea della meteora può essere scritta in forma parametrica: $$$ \left\{ \begin{array} {c} x = x_0 - a\cdot s \\ y = y_0 - b\cdot s \\ z = z_0 - c\cdot s \\ \end{array} \right. $$ \begin{center} dove $(x_0,y_0,z_0)$ sono le coordinate geocentriche \\ della postazione ed $s$ il parametro che esprime la \\ distanza lungo la traccia \end{center}$ Si scrive quindi il seguente sistema di 12 equazioni, raggruppabili a 3 a 3, di facile soluzione per confronto, come sarα mostrato appresso. Le coordinate dei 4 punti sono: A₁(x₁₁,y₁₁,z₁₁) B₁(x₂₁,y₂₁,z₂₁)
A₂(x₁₂,y₁₂,z₁₂) B₂(x₂₂,y₂₂,z₂₂) $\{1\} \left\{ \begin{array} {c} x_{11} = x_1 + a_{11}\cdot d_{11} \\ y_{11} = y_1 + b_{11}\cdot d_{11} \\ z_{11} = z_1 + c_{11}\cdot d_{11} \\ \end{array} \right. \qquad \{2\} \left\{ \begin{array} {c} x_{21} = x_2 + a_{21}\cdot d_{21} = x_{11} - a\cdot s_{21} \\ y_{21} = y_2 + b_{21}\cdot d_{21} = y_{11} - b\cdot s_{21} \\ z_{21} = z_2 + c_{21}\cdot d_{21} = z_{11} - c\cdot s_{21} \\ \end{array} \right. $$ $$ \{3\} \left\{ \begin{array} {c} x_{12} = x_1 + a_{12}\cdot d_{12} = x_{11} - a\cdot s_{12} \\ y_{12} = y_1 + b_{12}\cdot d_{12} = y_{11} - b\cdot s_{12} \\ z_{12} = z_1 + c_{12}\cdot d_{12} = z_{11} - c\cdot s_{12} \\ \end{array} \right. \qquad \{4\} \left\{ \begin{array} {c} x_{22} = x_2 + a_{22}\cdot d_{22} = x_{11} - a\cdot s_{22} \\ y_{22} = y_2 + b_{22}\cdot d_{22} = y_{11} - b\cdot s_{22} \\ z_{22} = z_2 + c_{22}\cdot d_{22} = z_{11} - c\cdot s_{22} \\ \end{array} \right.$ Confrontando le {1} con le {2} si ottiene il sistema lineare di tre equazioni nelle tre incognite (d₁₁,d₂₁,s₂₁): $\left\{ \begin{array} {c} a_{11}\cdot d_{11} - a_{21}\cdot d_{21} - a\cdot s_{21}= x_2 - x_1 \\ b_{11}\cdot d_{11} - b_{21}\cdot d_{21} - b\cdot s_{21}= y_2 - y_1 \\ c_{11}\cdot d_{11} - c_{21}\cdot d_{21} - c\cdot s_{21}= z_2 - z_1 \\ \end{array} \right.$ e quindi, determinato d₂₁, dalle {2} si ricavano le coordinate di B₁ e, noto d₁₁, dalle {1} quelle di A1. Dal confronto delle {1} con le {4} si ottiene un altro sistema 3x3 nelle incognite (d₁₁,d₂₂,s₂₂): $\left\{ \begin{array} {c} a_{11}\cdot d_{11} - a_{22}\cdot d_{22} - a\cdot s_{22}= x_2 - x_1 \\ b_{11}\cdot d_{11} - b_{22}\cdot d_{22} - b\cdot s_{22}= y_2 - y_1 \\ c_{11}\cdot d_{11} - c_{22}\cdot d_{22} - c\cdot s_{22}= z_2 - z_1 \\ \end{array} \right.$ e quindi, determinato d₂₂, dalle{4} si ricavano le coordinate di B₂. Infine dal confronto delle {1} con le {3} si calcolano le rimanenti incognite (d₁₂,s₁₂) per cui bastano solo 2 equazioni, mentre la terza si utilizza come controllo: $\left\{ \begin{array} {c} a_{12}\cdot d_{12} + a\cdot s_{12} = a_{11}\cdot d_{11} \\ b_{12}\cdot d_{12} + b\cdot s_{12} = b_{11}\cdot d_{11} \\ \end{array} \right.$ L'equazione di controllo è: c₁₂·d₁₂=c₁₁· d₁₁-c·s₁₂. Sostituendo d₁₂ nelle {3} si ottengono le coordinate geocentriche di A₂.
A questo punto, note le coordinate rettangolari geocentriche delle 4 posizioni, si passa immediatamente al calcolo del radiante geocentrico del meteoroide individuato dalla direzione del vettore (1°punto-2°punto): si può scegliere indifferentemente la coppia (A1,A2) o quella (B1,B2).

Radiante Apparente

TAN α_app = (y₁₁-y₁₂)/(x₁₁-x₁₂) con α_app nel corretto quadrante
SIN δ_app= (z₁₁-z₁₂)/dist.(A1-A2) essendo dist.(A1-A2) =
[(x₁₁-x₁₂)²+ (y₁₁-y₁₂)²+ (z₁₁-z₁₂)²]^½

Come si vede dalla figura a lato, le altezze dei due punti della meteora sulla superficie terrestre sono immediatamente determinate come differenza tra le distanze geocentriche dei punti A1 e A2 ed il raggio terrestre; si ha:
hA1= [x₁₁²+ y₁₁²+ z₁₁²]^½ - RT
hA2= [x₁₂²+ y₁₂²+ z₁₂²]^½ - RT
e analogamente per le quote di B1 e B2:
hB1= [x₂₁²+ y₂₁²+ z₂₁²]^½ - RT
hB2= [x₂₂²+ y₂₂²+ z₂₂²]^½ - RT
Le lunghezze delle scie luminose viste dalle due stazioni sono anch' esse di facile calcolo e valgono:
lA=[(x₁₁-x₁₂)²+ (y₁₁-y₁₂)²+ (z₁₁-z₁₂)²]^½
lB=[(x₂₁-x₂₂)²+ (y₂₁-y₂₂)²+ (z₂₁-z₂₂)²]^½

E' evidente che, se i punti sono stati scelti in maniera omologa, cioè B1 corrispondente ad A1(inizio scia per es.) e B2 ad A2 (fine scia), le lunghezze dei due tratti devono, nell'ambito delle tolleranze sperimentali, risultare UGUALI, così come la quota di A1 uguale a quella di B1 e la quota di A2 pari a quella di B2.
I dati vettoriali finora acquisiti permettono di calcolare facilmente altre grandezze complementari, quali le elevazioni angolari dei punti, le distanze sulla superficie terrestre e gli azimut. Esse non sono strettamente necessarie per l'obiettivo di questo lavoro, ma servono per avere un esauriente quadro d'insieme del fenomeno, per cui procediamo a stimarle. In figura è rappresentato un generico punto P della traccia distante d dalla stazione S ed r dal centro C della Terra. I 3 vettori (RT, d, r) sono noti ed abbiamo giα calcolato l'altezza h di P sulla superficie terrestre (sferica) come differenza dei moduli dei vettori r e RT. Il coseno dell'angolo al centro β si determina tramite il prodotto scalare di questi due vettori e vale:
COS β = (x·X+y·Y+z·Z) / (RT·r) da cui, noto β, si ricava la distanza S-P0 = RT·β(radianti)
Introdotto un sistema di coordinate orizzontali (Sud-Est-Zenit) con centro nella stazione S, le componenti (dS,dE,dZ) del vettore topocentrico d si ricavano dalle seguenti equazioni, tratte da [3]: $\left | \begin{array} {c} dS \\ dE \\ dZ \\ \end{array} \right | = \left | \begin{array} {ccc} \sin\phi \cos\theta & \sin\phi \sin\theta & -\cos\phi \\ -\sin\theta & \cos\theta & 0 \\ \cos\phi \cos\theta & \cos\phi \sin\theta & \sin\phi \\ \end{array} \right | \cdot \left | \begin{array} {c} x \\ y \\ z \\ \end{array} \right |$ essendo (x,y,z) le coordinate di P e θ il tempo siderale locale, che è uguale a quello di Greenwich + la longitudine est di S. Dalla terza equazione si determina dZ e quindi l'elevazione angolare (k) di P sull'orizzonte di S, mentre con le altre due si calcola l'azimut Nord-Est del punto P: $$$ \begin{array} {cc} \sin k=\frac{dZ}{d} & \tan A_z = \frac{dE}{-dS} = \frac{dE}{dN} \\ \end{array} $$ \begin{center} con dN=-dS componente nord di d \end{center}$

Il teorema di LAMBERT

Questo teorema, noto fin dalla metα del XVIII secolo, trova diffusissima applicazione nella definizione e correzione di orbite di sonde interplanetarie e in anni recenti molti sono stati i ricercatori di ogni parte del mondo che vi si sono dedicati, a testimonianza del FASCINO che esso esercita sugli studiosi contemporanei di astrodinamica. La generalitα applicativa dell'algoritmo creato da Battin (cfr.[4]) e la sua straordinaria convergenza anche in casi limite come θ=π (laddove fallì il grande Friedrich GAUSS), lo fanno preferire ai metodi degli altri autori e pertanto lo adottiamo per sviluppare i calcoli orbitali della meteora in esame;l'algoritmo di Battin è riportato nell'Appendice 1 e pronto per essere facilmente trascritto in un programma per computer.

Per trovare gli elementi orbitali di una conica basta allora conoscere i due raggi vettori r1 e r2 dei due estremi dell'arco P1 e P2,l'angolo tra essi compreso e il tempo di percorrenza τ=t2-t1. Nel caso della meteora considerata noi conosciamo le componenti vettoriali dei due estremi A1 e A2 (si possono scegliere anche B1 e B2, non c'è differenza) che abbiamo calcolato e sono: A₁(x₁₁,y₁₁,z₁₁) A₂(x₁₂,y₁₂,z₁₂) Il coseno dell'angolo θ compreso tra i raggi vettori C-A1 e C-A2 è noto dal calcolo vettoriale e si ottiene dividendo il prodotto scalare dei due vettori per il prodotto dei due raggi, cioè in formula COS θ = (x₁₁·x₁₂+ y₁₁·y₁₂+ z₁₁·z₁₂) / (r₁₁·r₁₂)
essendo r₁₁= (x₁₁²+ y₁₁²+ z₁₁²)^½ e r₁₂= (x₁₂²+ y₁₂²+ z₁₂²)^½ La durata τ è un dato sperimentale che si ricava da una traccia fotografica quando è spezzettata da un otturatore rotante; se non è conosciuto si ricorre a un valore arbitrario e ci si accontenta, a malincuore, di dati imprecisi sulle velocitα del corpo e sul suo radiante vero.

Gli elementi di forma (a,e) della conica, che deve essere necessariamente un'IPERBOLE in quanto la costante energetica del moto (-μ/a) è sempre positiva, si calcolano allora mediante l'algoritmo [4] e insieme ad essi si trovano i parametri DINAMICI della meteora, vale a dire le due velocitα V₁ e V₂ e gli angoli di giacitura orbitale (i,ω,Ω). Una rappresentazione grafica, naturalmente fuori scala, di un'orbita di un meteoroide è riportata nella figura a lato e in essa è evidenziato il versore n˜ dell'asintoto del ramo superiore dell'iperbole che ci accingiamo a calcolare.

Chiamati (p˜,q˜) i versori perifocali della conica, si determinano i coseni direttori (p,q) del versore n˜ dell'asintoto, sapendo che il coseno della sua inclinazione vale: COSβ=-1/e. $\left | \begin{array} {c} p \\ q \\ 0 \\ \end{array} \right | = \left | \begin{array} {c} \cos\beta \\ -\sin\beta \\ 0 \\ \end{array} \right | = \left | \begin{array} {c} -\/e \\ -\sqrt{1-1/e^2} \\ 0 \\ \end{array} \right |$ Si prosegue con il calcolo dei versori (x˜,y˜,z˜) dell'asintoto rispetto al piano equatoriale mediante la seguente matrice di trasformazione: $\left | \begin{array} {c} \tilde{x} \\ \tilde{y} \\ \tilde{z} \\ \end{array} \right | = \left | \begin{array} {ccc} R_{11} & R_{12} & ... \\ R_{21} & R_{22} & ... \\ R_{31} & R_{32} & ... \\ \end{array} \right | \cdot \left | \begin{array} {c} p \\ p \\ 0 \\ \end{array} \right | $$ \begin{center} che in forma esplicita diventano: \end{center} $$ \left\{ \begin{array} {c} \tilde{x} = R_{11}\cdot p + R_{12}\cdot q \\ \tilde{y} = R_{21}\cdot p + R_{22}\cdot q \\ \tilde{z} = R_{31}\cdot p + R_{32}\cdot q \\ \end{array} \right. $$ \begin{center} con: \end{center} $$ \left\{ \begin{array} {l} R_{11} = \cos\Omega\cdot \cos\omega - \sin\Omega\cdot \sin\omega\cdot \cos i \\ R_{12} =-\cos\Omega\cdot \sin\omega - \sin\Omega\cdot \cos\omega\cdot \cos i \\ R_{21} = \sin\Omega\cdot \cos\omega + \cos\Omega\cdot \sin\omega\cdot \cos i \\ R_{22} =-\sin\Omega\cdot \sin\omega + \cos\Omega\cdot \cos\omega\cdot \cos i \\ R_{31} = \sin\omega\cdot \sin i \\ R_{32} = \cos\omega\cdot \sin i \\ \end{array} \right.$ Le coordinate polari (α_vero,δ_vero) del radiante vero sono quindi determinate con le formule:

Radiante Vero

TAN α_vero = y˜/x˜ con α_vero nel corretto quadrante SIN δ_vero= z˜

1° ESEMPIO

Le foto della Perseide 3, gentilmente fornite da [5], sono state riprese da due localitα situate nell'Italia centrale e aventi le seguenti coordinate geografiche:

   A≡S1 (44°07'35" N; 10°47'05" E)    B≡S1 (44°12'20" N; 10°44'10" E)
            Data di osservazione: 12 08 1991 ore 22:58:15 TU

1^a FOTO

2^a FOTO

Entrambe sono state processate con lo scanner del PC in modo da avere dei file grafici sui quali poter marcare i punti di inizio e fine traccia e numerare le stelle di riferimento.
Per determinare le coordinate equatoriali dei 4 punti OMOLOGHI, cioè Ini e Fin della 1^a traccia corrispondenti ad A1 e A2 e i rispettivi B1 e B2 della 2^a traccia, si fa ricorso al metodo di confronto noto come RIDUZIONE ASTROMETRICA che trasforma le posizioni degli oggetti della lastra in coordinate standard; un buon algoritmo di calcolo per la riduzione stellare è riportato in [6]. Per la nostra meteora abbiamo preferito adottare un metodo pi· dettagliato ricorrendo al software [7], che utilizza automaticamente il precisissimo catalogo PPM. Esso può apparire esagerato ai fini di una valutazione abbastanza grezza delle misure in gioco, ma dα sicuramente la certezza, attraverso l'analisi ponderata dei residui in AR e declinazione, che non si è incorsi in errori grossolani.

  Riduzione  dati  lastre  fotografiche  con  calcolo  posizione  topocentrica
     del  corpo  celeste   ( min.3   max.15  stelle ) ,  di  Augusto  TESTA
                         * Osservatorio  di  Sormano *
+------------------------------------------------------------------------------+
| Titolo :  Perseide 3 - Staz.A - Inizio Traccia         Data 1991.08 12.95712 |
| Misur. VM      cat. PPM  n stelle abil.   4   ARC  18.1800     DC  46.2400   |
+------------------------------------------------------------------------------+
|st.n  S/N  Catal.n    Ascen.Retta  m.p. Declinazione  m.p.   X (mm)    Y (mm) |
|   1   1     57038    18.2132611   -33    49.071764    52     4.580     1.960 |
|   2   1     57761    19.0126424    27    46.560562   -79     4.790     4.860 |
|   3   1     57635    18.5520111    21    43.564599    82     6.090     4.800 |
|   4   0     36267    17.5636367    -8    51.292021   -19     3.820     0.100 |
|   5   1     57119    18.3033901    12    48.175239   -17     4.800     2.630 |
|       1               Misure del Pianetino o Cometa          4.760     2.640 |
+------------------------------------------------------------------------------+
                              h m s.dd                     ' ".d
                Ascen.Retta  18.304983    Declinazione  48.22421
                RMS secondi       1.537 tempo                4.27 arco
                                 23.06  arco
                   Residui stelle (catalogo PPM )  in < mm >
                       n   1 ->    X  -0.003   Y   0.001
                       n   2 ->    X  -0.000   Y   0.000
                       n   3 ->    X  -0.001   Y   0.000
                       n   5 ->    X   0.004   Y  -0.001
                       Deviazione      0.003       0.000

+------------------------------------------------------------------------------+
| Titolo :  Perseide 3 - Staz.A - Fine Traccia           Data 1991.08 12.95712 |
| Misur. VM      cat. PPM  n stelle abil.   5   ARC  17.5800     DC  32.3900   |
+------------------------------------------------------------------------------+
|st.n  S/N  Catal.n    Ascen.Retta  m.p. Declinazione  m.p.   X (mm)    Y (mm) |
|   1'  1     80778    17.5724269     0    33.240337   -12    11.800     0.770 |
|   2'  1     80788    17.5804995   -32    32.385314   -27    12.170     0.860 |
|   3'  1     80831    18.0036319   -13    33.125010   -11    11.830     1.090 |
|   4'  1     80918    18.0549554     7    32.135049   -17    12.180     1.650 |
|   5'  1     80796    17.5830135    -7    30.112124    -3    13.370     0.990 |
|       1               Misure del Pianetino o Cometa         12.330     0.520 |
+------------------------------------------------------------------------------+
                              h m s.dd                     ' ".d
                Ascen.Retta  17.543595    Declinazione  32.28276
                RMS secondi       1.887 tempo                6.56 arco
                                 28.31  arco
                   Residui stelle (catalogo PPM )  in < mm >
                       n   1 ->    X  -0.008   Y  -0.002
                       n   2 ->    X   0.006   Y   0.001
                       n   3 ->    X   0.005   Y   0.001
                       n   4 ->    X  -0.002   Y  -0.000
                       n   5 ->    X  -0.001   Y  -0.000
                       Deviazione      0.005       0.001

+------------------------------------------------------------------------------+
| Titolo :  Perseide 3 - Staz.B - Inizio Traccia         Data 1991.08 12.95712 |
| Misur. VM      cat. PPM  n stelle abil.   4   ARC  18.4800     DC  45.1500   |
+------------------------------------------------------------------------------+
|st.n  S/N  Catal.n    Ascen.Retta  m.p. Declinazione  m.p.   X (mm)    Y (mm) |
|   1   1     57038    18.2132611   -33    49.071764    52     3.960    13.280 |
|   2   1     57761    19.0126424    27    46.560562   -79     4.490    16.000 |
|   3   1     57635    18.5520111    21    43.564599    82     5.720    15.800 |
|   4   1     57641    18.5533377   -50    47.262653   -60     4.370    15.580 |
|   5   0     57665    18.5632119     5    46.455848   -42     4.620    15.700 |
|       1               Misure del Pianetino o Cometa          5.200    15.270 |
+------------------------------------------------------------------------------+
                              h m s.dd                     ' ".d
                Ascen.Retta  18.490394    Declinazione  45.27546
                RMS secondi       1.332 tempo               12.16 arco
                                 19.98  arco

                   Residui stelle (catalogo PPM )  in < mm >
                       n   1 ->    X  -0.001   Y   0.000
                       n   2 ->    X  -0.003   Y   0.002
                       n   3 ->    X   0.000   Y  -0.000
                       n   4 ->    X   0.003   Y  -0.002
                       Deviazione      0.002       0.001

+------------------------------------------------------------------------------+
| Titolo :  Perseide 3 - Staz.B - Fine Traccia           Data 1991.08 12.95712 |
| Misur. VM      cat. PPM  n stelle abil.   4   ARC  18.0700     DC  28.4500   |
+------------------------------------------------------------------------------+
|st.n  S/N  Catal.n    Ascen.Retta  m.p. Declinazione  m.p.   X (mm)    Y (mm) |
|   1'  1     81056    18.1154108   -17    31.241907    22    11.210    12.960 |
|   2'  1     80945    18.0701566   -70    30.334314    74    11.610    12.520 |
|   3'  1    106802    18.0732555     1    28.454501     9    12.390    12.580 |
|   4'  0     80796    17.5830135    -7    30.112124    -3    11.830    11.740 |
|   5'  1    106582    17.5745895    64    29.145246   -17    12.240    11.660 |
|       1               Misure del Pianetino o Cometa         12.010    12.960 |
+------------------------------------------------------------------------------+
                              h m s.dd                     ' ".d
                Ascen.Retta  18.114047    Declinazione  29.33553
                RMS secondi       0.631 tempo                4.33 arco
                                  9.46  arco
                   Residui stelle (catalogo PPM )  in < mm >
                       n   1 ->    X   0.001   Y   0.001
                       n   2 ->    X  -0.002   Y  -0.001
                       n   3 ->    X   0.000   Y   0.000
                       n   5 ->    X   0.001   Y   0.000
                       Deviazione      0.001       0.001

                              ---------------------
                                R I E P I L O G O
                              ---------------------

             Staz.A -> Inizio Traccia  AR  18 30 49.83  =  277.7076°
                                       De  48 22 42.1   =   48.3784°

             Staz.A ->   Fine Traccia  AR  17 54 35.95  =  268.6498°
                                       De  32 28 27.6   =   32.4743°
             -------------------------------------------------------
             Staz.B -> Inizio Traccia  AR  18 49  3.94  =  282.2664°
                                       De  45 27 54.6   =   45.4652°

             Staz.B ->   Fine Traccia  AR  18 11 40.47  =  272.9186°
                                       De  29 33 55.3   =   29.5654°

Le schermate della riduzione astrometrica sono riportate sopra unitamente ad un riepilogo delle coordinate equatoriali dei 4 punti. L'utilizzo dei file grafici ha permesso inoltre una precisione ottimale della misure delle stelle e delle posizioni delle tracce grazie alle letture micrometriche effettuate in PIXEL; successivamente i valori letti sono stati divisi per una costante arbitraria (100) ed inseriti come millimetri in [7]. Poichè infine non si conosce la durata della traccia della perseide esaminata si è supposto che τ fosse verosimilmente di 0.63 sec, quindi il quadro dei dati è completo e trascritto sotto nello stesso standard del file dati DUBY.dat di [2]. I valori di α, δ degli estremi delle tracce sono riferiti all'equinozio medio della data

   1991  8  12   22  58  15
   44.1264   10.7847  44.2055   10.7361
   0  0.63
   277.7076  48.3784  282.2664  45.4652
   268.6498  32.4743  272.9186  29.5654

   -------------  Sono lette le prime 5 RIGHE ----------------------------------
   1^ riga: Data dell'osservaz.: Anno, mese, giorno; Ore, Min, Sec (in TU)
   2^ riga: Latitudine, Longit.Est in GRADI 1^ Staz.(A)
            Latitudine, Longit.Est in GRADI 2^ Staz.(B)
   3^       Tempi in SEC  [t1=0; t2=...]
   4^ riga: (α,δ) in GRADI del primo Punto visto da A e da B; Punti NON necessa-
   5^ riga: (α,δ) in GRADI del secon.Punto visto da A e da B  riamente omologhi
   -----------------------------------------------------------------------------

  Vettore geocentrico di una METEORA fotografata da 2 stazioni terrestri, note
  le coord. (α,δ) di 2+2 punti NON necessariamente OMOLOGHI  -  Metodo DUBYAGO
  ----------------------------------------------------------------------------
                              RT= 6367.109 km

 QUOTE met. ->   A1= 112.1 km     A2= 90.0 km    [Dist_Stazioni= 9.6 km]
                 B1= 112.1 km     B2= 90.3 km
 LUNGH.scie ->   lA=  37.6 km     lB= 37.0 km

 ... Per punti OMOLOGHI Q_A1=Q_B1, Q_A2=Q_B2 + L_scie Uguali

     RADIANTE Apparente: α= 46.7°   δ= 58.6°

         Coordinate geocentriche di A1 e A2 (punti visti dalla Staz_1)
         -------------------------------------------------------------
           x11= 0.520134 RT    y11=-0.509452 RT    z11= 0.710940 RT
           x12= 0.518025 RT    y12=-0.511687 RT    z12= 0.705904 RT
           Dist_Topocentr.   ro_A1= 125.2 km     ro_A2= 114.6 km
           Dist_Terrestre     d_A1=  55.3 km      d_A2=  70.4 km
           Elevazione        El_A1=  63.31°      El_A2=  51.47°
           Azimut (N-E)      Az_A1= 292.91°      Az_A2= 269.01°


        Coordinate geocentriche di B1 e B2 (punti visti dalla Staz_2)
        -------------------------------------------------------------
           x21= 0.520131 RT    y21=-0.509456 RT    z21= 0.710932 RT
           x22= 0.518051 RT    y22=-0.511659 RT    z22= 0.705967 RT
           Dist_Topocentr.   ro_B1= 122.4 km     ro_B2= 112.7 km
           Dist_Terrestre     d_B1=  48.7 km      d_B2=  67.0 km
           Elevazione        El_B1=  66.09°      El_B2=  52.94°
           Azimut (N-E)      Az_B1= 285.05°      Az_B2= 261.61°


  PROBLEMA di LAMBERT per traiettorie di METEORE Terrestri (Metodo BATTIN)
  ------------------------------------------------------------------------
                        Tempo di volo   τ [sec] =  0.630
                        Angolo Trasfer. θ [ gradi] =   0.26924
                        r1=  1.01760372 RT  r2=  1.01413739 RT
        N_iter.  2)  yit=  1.00000002       Residuo = -6.86539713967704E-0012
       i= 117.4684   ω= 164.9049  Ω= 105.0879
                               RISULTATI FINALI
                      Semiasse magg.  a = -0.01814490 RT  =  -116 km
                      Eccentricitα    e = 46.11572481
                      Sistema EQUATORIALE Geocentrico
 POS/VEL_1       r1= 1.017604 RT =  6479 km ♦ V1 = 7.55495 RT/TU = 59.621 km/s
     x1=  0.520134 RT =   3312 km   ♦   V1x= -2.70073 RT/TU = -21.313 km/s
     y1= -0.509452 RT =  -3244 km   ♦   V1y= -2.86222 RT/TU = -22.588 km/s
     z1=  0.710940 RT =   4527 km   ♦   V1z= -6.44911 RT/TU = -50.894 km/s

 POS/VEL_2       r2= 1.014137 RT =  6457 km ♦ V2 = 7.55540 RT/TU = 59.625 km/s
     x2=  0.518025 RT =   3298 km   ♦   V2x= -2.70112 RT/TU = -21.316 km/s
     y2= -0.511687 RT =  -3258 km   ♦   V2y= -2.86184 RT/TU = -22.585 km/s
     z2=  0.705904 RT =   4495 km   ♦   V2z= -6.44964 RT/TU = -50.899 km/s

     RADIANTE Apparente: α= 46.7°   δ= 58.6°     Vero: α= 47.5°   δ= 58.4°

I risultati finali dell'altezza della Perseide 3, dei due radianti, delle distanze e degli azimut, oltre agli elementi orbitali della traiettoria iperbolica e delle velocitα agli estremi della traccia luminosa sono stati ripresi dal file DUBY.ris di [2] e copiati sopra.

2° ESEMPIO

Viene preso in considerazione il bolide avvistato e fotografato da diverse localitα dell'Italia settentrionale nell'agosto del '93, i cui dati sperimentali sono stati cortesemente forniti da [8].
Poichè il corpo è stato visto da 5 posti diversi, si sono scelte le due stazioni pi· vicine ad esso e le coordinate di 2 punti scia pi· riconoscibili. I dati del problema sono i seguenti:

       A≡S1 (45°54'15.6" N; 09°29'18" E)    B≡S1 (46°03'15" N; 11°18'48" E)
                Data di osservazione: 11 08 1993 ore 23:13:20 TU

            Punto Scia 5        A1  α=00h34m11s      B1  α=19h27m10s
              (in alto)             δ=20°47'00"          δ=18°53'00"

            Punto Scia 1        A2  α=00h18m27s      B2  α=19h01m15s
             (in basso)             δ=12°17'28"          δ=07°12'00"

Con l'ausilio del programma computerizzato si ottengono le quote dei 2 punti (visti dalla stazione A, per es.) pari rispettivamente a 88.7 e 71.8 km, gli azimut nord-est uguali a 104.19° e 115.42° e le distanze terrestri di 96.1 e 88.6 km. Con questi elementi (e due note formule di trigonometria sferica) si determinano le coordinate geografiche delle località che si trovano sulle verticali dei due punti scia considerati, cioè:

                              Latit.        Long.Est
                          ----------------------------
            Punto Scia 5    45°41'09"       10°41'18"
                          [ 45 41 23 ]    [ 10 40 51 ]

            Punto Scia 1    45°33'27"       10°31'23"
                          [ 45 33 35 ]    [ 10 30 06 ]

I valori entro parentesi quadre sono quelli pubblicati nella Tabella I di [9], che riporta anche le altezze dei due punti, pari rispettivamente a 90.0 km e 73.6 km, tutti in OTTIMO accordo con quelli calcolati con questo metodo.

C O N C L U S I O N E

Il calcolo vettoriale possiede indubbiamente delle enormi potenzialità e come tale riesce a semplificare moltissimo la soluzione di problemi che con altri procedimenti, quali quelli trigonometrici, risulta laboriosa. La Meccanica Celeste offre poi gli algoritmi appropriati per definire e completare esaustivamente lo studio di un'orbita.

ing. Giuseppe MATARAZZO
03.08.2000 ♦ A corredo l'Appendice 1

Riferimenti Bibliografici

[1] A.D.Dubyago - The Determination Of Orbits - Mc Millan Company, New
    York, 1961; in inglese tradotto dal russo. Cap.12: Determination of
    Meteor Orbits

[2] Giuseppe Matarazzo - Programma MY_DUBY in linguaggio Pascal (Giu.00),
    riadattato dal Basic BATTBEST.BAS (Ago.95) per la parte riguardante il
    teorema di Lambert risolto con l'algoritmo di Battin

[3] R.R.Bate, D.D.Mueller, J.P.White - Fundamentals of Astrodynamics -
    Dover, New York 1971; p.117: coseni direttori vettore topocentrico;
    p.98-99: coord. vettore geocentrico stazione; p.79: matrice di rota=
    zione per passaggio da sist. coord. geocentriche a topocentriche

[4] Richard H. Battin, Robin M. Vaughan
    An Elegant Lambert Algorithm - Vol.7,n.6 del Journal of Guidance and
    Control - Nov/Dic.1984

[5] Mirco Villa - comunicazioni personali

[6] O.Montenbruck, T.Pfleger - Astronomy on the Personal Conputer -
    Springer-Verlag editore: cap.10 Astrometry p.197-210

[7] Augusto Testa (Osservatorio di Sormano) - Programma MAPPA2, vers.6.7-
    2000, Database del catalogo PPM fino alla 12^ magn. e programma di
    riduzione astrometrica RID2 (1992)

[8] Roberto Serpilli - comunicazioni personali

[9] Roberto Serpilli - Osservazione di un bolide l'11-08-1993 - Astronomia
    UAI n.4 - 1994: p.15

Appendice 1

         ALGORITMO di BATTIN per RISOLVERE il PROBLEMA di LAMBERT
         --------------------------------------------------------
Si procede per successive approssimazioni della grandezza y, uguale
al rapporto tra il settore C-P1-P2 e il triangolo sotteso C-P1-P2. Ecco
nell'ordine ed in sintesi, i vari steps del procedimento iterativo.
Sono noti: r1(x1,y1,z1), r2(x2,y2,z2), θ, τ=(t2-t1).

     1) Calcolo dei parametri ausiliari: (c,s,λ,í,L,rop,m)
        c= √(r1² + r2² - 2 * r1 * r2 * COS(θ))  <- corda dell'arco
        s= (r1 + r2 + c) / 2                    <- semi-perimetro
        λ= √(r1 * r2) / s * COS(θ/2)
        COS(í)= √(r1·r2)·COS(θ/2) / [«·(r1+r2)]    da cui -> í
        L= TAN²(í/2)
        rop=¬·(r1+r2+2·√(r1·r2)·COS(θ/2))
        m=μ·τ²/(8·rop^3)  con μ=parametro gravitaz. =1 in unita' canoniche

     2) Scelta del valore iniziale x0 dell'iterazione (per es. x0=0,
        oppure x0=L), essendo il parametro x compreso tra [-1 e +ì[
        {saltare al punto 4}

     3) Calcolo di x dopo la prima iterazione
        x= √[(1 - L) / 2)² + m / y²] - (1 + L) / 2

     4) Calcolo della Funzione Ipergeometrica (csi)
            Posto:  η= x/[1 + √(1 + x)]²
     gs3= 169/25/27*η/(1 + 196/27/29*η/(1 + 225/29/31*η / (1 + 0)))
     gs2= 100/19/21*η/(1 + 121/21/23*η/(1 + 144/23/25*η / (1 + gs3)))
     gs1= 49/13/15*η/(1 + 64/15/17*η/(1 + 81/17/19*η/(1 + gs2)))
     gs0= 16/7/9*η/(1 + 25/9/11*η/(1 + 36/11/13*η/(1 + gs1)))
                                si ottiene
     csi= 8*(1 + sqrt(1 + x))/(3 + 1/(5 + η + 9/7*η/(1+gs0)))

     5) Calcolo dei parametri (h1,h2)
        Posto: hden= (1 + 2 * x + L) * (4 * x + csi * (3 + x)) si ha:
     h1= (L + x)² * (1 + 3 * x + csi)/hden;  h2= m * (x - L + csi)/hden

     6) Calcolo della Funzione Ipergeometrica (ku)
        Posto: B= 27 * h2 / (4 * (1 + h1)*(1 + h1)²)
               u= -B / 2 / (1 + √(1 + B))
        fs3= 928/3591*u/(1 - 1054/4347*u/(1 - 266/1035*u/(1 - 0)))
        fs2= 418/1755*u/(1 - 598/2295*u/(1 - 700/2907*u/(1 - fs3)))
        fs1= 22/81*u/(1 - 208/891*u/(1 - 340/1287*u/(1 - fs2)))
        fs0= 4/27*u/(1 - 8/27*u/(1 - 2/9*u/(1 - fs1)))
                                si ottiene
         ku= 1/3/(1 - fs0)

     7) Radice positiva y del polinomio di 3° grado di Gauss
        y_iterativo= (1 + h1) / 3 * (2 + √(1 + B) / (1 - 2 * u * ku²)
        Se | y   -  y | < î  (=1·E-6 per es.)  fine iterazione -> punto 8)
              i+1    i       altrimenti saltare al punto 3)

     8) Elementi dell'orbita (a,p,e)
        Detto yf il valore di y dell'ultima iterazione, si ottiene:
        a= m * s * (1 + λ)² / (8 * x * yf²)
        p= 2 * r1 * r2 * (yf * (1 + x) * SIN(θ/2))² / ((1 + λ)² * m * s)
        e= √(1 - p / a)

     9) Coefficienti lagrangiani (f,g,f',g')
        f= 1 - r2 * 2 * SIN²(θ/2) / p
        g= r1 * r2 * SIN(θ) / √(μ * p)
        f'= (2 / p * SIN²(θ/2) - 1/r1 - 1/r2) * TAN(θ/2) * √(μ/p)
        g'= 1 - r1 * 2 * SIN²(θ/2) / p

    10) Componenti delle VELOCITA' ai 2 estremi dell'arco P1 e P2
        Vx1=(x2-f*x1)/g;   Vy1=(y2-f*y1)/g;   Vz1=(z2-f*z1)/g;
        Vx2=f'*x1+g'*Vx1; Vy2:=f'*y1+g'*Vy1; Vz2:=f'*z1+g'*Vz1;
        Moduli delle VELOCITA':  V1=√(Vx1²+Vy1²+Vz1²); V2=√(Vx2²+Vy2²+Vz2²)

(.... fine)

Calcolo della traiettoria geocentrica di un METEOROIDE tramite il metodo vettoriale di DUBYAGO e il teorema di LAMBERT