Autorem tego opracowania jest dr Piotr A. Dybczyński z Instytutu Obserwatorium Astronomiczne UAM w Poznaniu.
Najpierw prezent
Zanim przejdziemy do właściwego zagadnienia musimy stworzyć podstawowy warsztat rachunkowy dla ruchu keplerowskiego. Skorzystamy tu z formuł dla heliocentrycznego ruchu ciała o zaniedbywalnej masie (np. komety czy planetoidy). Oznacza to, że jednostką odległości będzie tzw. jednostka astronomiczna (AU), jednostką czasu doba a jednostką masy masa Słońca. W tych jednostkach stała grawitacji jest równa kwadratowi tzw. stałej Gaussa: G = k2 = (0.01720209895)2 . Zagubionym w obszarze mechaniki nieba polecam wykład Profesora Sławomira Breitera.
Musimy umieć wyliczyć pozycję ciała na zadany moment czasu z danych elementów orbity. Najtrudniejszym miejscem jest tu rozwiązanie równania Keplera. Sposobów praktycznego rozwiązania tego zagadnienia jest wiele, ja zaproponuję użycie jednej procedury, która wylicza anomalię prawdziwą po podaniu mimośrodu, odległości peryhelium i interwału czasu w dniach (może być ujemny) od momentu przejścia praz peryhelium. Pomysł i algorytm tej procedury pochodzi z pracy Profesora Grzegorza Sitarskiego (1968, Acta Astronomica,vol.18,no.2,strony 197-205).
Poniżej prezent, czyli gotowa i przetestowana procedura realizująca ten rachunek:
/****************************************************************************
******** ANOMALIAT.C *********
*****************************************************************************
******** **********
******** **********
*****************************************************************************
** **
** Pomysł i algorytm: **
** G.Sitarski (1968), Acta Astronomica,vol.18,no.2,pp 197-205, **
** **
*****************************************************************************
** Programmed by: Piotr A. Dybczyński, Last modified: 3.03.2012 **
****************************************************************************/
#define GAUSS 0.01720209895
#define E_SMALL (1.e-12)
#include<math.h>
double anomaliat(double dt,double q,double e)
{
double x,M,M1,E,E1,w,F,v,u,F1,z,n,mi;
x = (1.-e)/(1.+e);
mi=GAUSS*GAUSS; /* tylko dla ciał o zaniedbywalnej masie */
if(e<(1.-E_SMALL)) /* elipsa */
{
n = sqrt(mi*(1.-e)*(1.-e)*(1.-e)/(q*q*q));
M = n*dt;
M = fmod(M,2.0*M_PI);
M1 = fabs(M);
E1 = M_PI;
do {
E = E1;
E1 = E + (v=M1-E+e*sin(E))/(1.-e*cos(E));
} while( fabs(v) > fabs(M1-E1+e*sin(E1) ));
if(1.!=1.+cos(E1/2.)) w = sqrt(1./x)*tan(E1/2.);
else w = 1.0E+70; /* theta = 180 deg*/
}
else if(e>(1.+E_SMALL)) /*hiperbola */
{
M = sqrt(mi*(e-1.)*(e-1.)*(e-1.)/(q*q*q))*dt;
M1 = fabs(M);
F1 = 1;
while( (z = M1+F1-e*sinh(F1)) > 0.0 ) F1+=1.;
do {
F = F1;
F1 = F + (v=z)/(e*cosh(F)-1.);
} while( fabs(v) > fabs(z=M1+F1-e*sinh(F1)));
w = sqrt(-1./x)*tanh(F1/2.);
}
else /* bliskoparabola: e = 1 +/- E_SMALL */
{
M = sqrt(mi/(2.*q*q*q))*dt;
u = fabs(M)*1.5;
v = sqrt(1.+u*u);
w = pow((v+u),1./3.) - pow((v-u),1./3.);
}
return( M<0? -2.*atan(w): 2.*atan(w) );
}
Kilka słów wyjaśnienia: procedura definiuje stałą E_SMALL i jeśli mimośród różni się od jedynki o mniej niż E_SMALL korzysta ze wzorów dla paraboli, w pozostałych przypadkach rozwiązywane jest eliptyczne lub hiperboliczne równanie Keplera. Daje to optymalną dokładność numeryczną.
Pozostaje napisać prosty program testowy, który poprosi o podanie q, e oraz dt i wylicza anomalię prawdziwą f. W cytowanej pracy Profesora Sitarskiego jest tabelka z danymi testowymi ale prościej jest chyba porównać działanie takiego programu z moim, zamieszczonym na pracowni w katalogu /home/COMMON/varia programem anomalia.
Położenie ciała na orbicie keplerowskiej - metoda klasyczna
Orbitę keplerowską opisuje się przy pomocy tzw. elementów orbity. Zależnie od zagadnienia czy typu orbity stosuje się różne, równoważne zestawy sześciu elementów. Dla komet najwygodniejszym zestawem jest:
- odległość peryhelium w AU q
- mimośród e
- moment przejścia peryhelium (data lub dzień juliański) T
- długość węzła wstępującego Ω
- argument peryhelium ω
- nachylenie (najczęściej do płaszczyzny ekliptyki) i
W innych zagadnieniach (ale i w ruchu komet) używane bywają również inne elementy:
- parametr orbity p = q(1 + e)
- wielka półoś orbity a, dla elipsy: a = q / (1 - e), dla hiperboli: a = q / (e - 1)
- odległość aphelium (dla elipsy) Q = q (1 + e) / (1 - e)
- ruch średni n = k / a3
- anomalia średnia Mo na jakiś ustalony moment czasu
Położenie ciała na orbicie keplerowskiej opisujemy tradycyjnie przy pomocy jednej z anomalii: prawdziwej f, mimośrodowej E lub średniej M. Ponieważ będziemy korzystać z załączonej wyżej procedury użyjemy anomalii prawdziwej.
W płaskim (tzw. orbitalnym) układzie odniesienia Oξη związanym z płaszczyzną orbity, ze środkiem w masie centralnej (Słońcu), z osią Oξ skierowaną do peryhelium a osią Oη prostopadle do niej i zgodnie z kierunkiem ruchu ciała, składowe położenia i prędkości pozycję prostymi wzorami:
` r = p/(1+e cos f )`
Aby uzyskać prostokątne składowe położenia w układzie równikowym przechodzimy po drodze przez układ ekliptyczny (względem którego mamy określone elementy) a cała transformacja wymaga w sumie czterech obrotów:
r(α, δ) = Rx(-ε) Rz(-Ω) Rx(-i) Rz(-ω) r(ξ,η)
gdzie ε jest nachyleniem ekliptyki do równika niebieskiego. Proszę używać wartości dla epoki J2000: ε = 23° 26' 21.448" .
Oczywiście składowe prędkości uzyskujemy z takiej samej transformacji.
Położenie ciała na orbicie keplerowskiej - użycie zmiennych uniwersalnych
Zupełnie inaczej wylicza się położenie i prędkość na dowolnej orbicie keplerowskiej przy użyciu zmiennych uniwersalnych. Szczegóły można znaleźć w tym materiale i cytowanej tam pracy.
Teraz zadanie
Napisać program, który pobiera jako dane elementy dowolnej (niezdegenerowanej) orbity keplerowskiej oraz dwa dowolne momenty czasu. Następnie wylicza metodą klasyczną położenia i prędkości na moment pierwszy, wylicza przy użyciu zmiennych uniwersalnych położenia i prędkości na moment drugi i przy użyciu metody konwencjonalnej je sprawdza.
Dane testowe można sobie wygenerować moim programem, do znalezienia na pracowni: /home/COMMON/varia/unitest