Calcul
d'éphémérides dans le systéme solaire
par
Ghislain de FROMENT
Cet exposé, sans prétendre
à l'exhaustivité, est destiné à présenter
diverses méthodes permettant d'obtenir la position d'un astre mobile
dans le systéme solaire.
A) Rappels et définitions
1) Repéres
Pour l'établissement d'éphémérides
dans le systéme solaire, deux repéres sont d'usage courant
:
-
le repére écliptique, fixe
ou quasiment fixe, dont le plan fondamental est le plan moyen du mouvement
de la Terre autour du Soleil, dans lequel sont définis les éléments
osculateurs;
-
le repére équatorial, dont
le plan fondamental est le plan équatorial terrestre, et auquel
sont rattachées les observations;
Ces plans, non paralléles, se coupent
selon une droite : celle ci sert à leurs orientations respectives.
Sur cette droite se trouve le point vernal g,
défini comme suit : c'est la direction du Soleil, vu de la Terre,
lorsque, dans son mouvement apparent, il passe de l'hémisphére
Sud à l'hémisphére Nord du repére équatorial.
Sous l'effet de diverses perturbations
(essentiellement liées à la forme de la Terre donnant prise
aux attractions de la Lune et du Soleil), ces repéres sont mobiles
l'un par rapport à l'autre, ce qui donne lieu aux phénoménes
de précession et de nutation.
Les perturbations d'origine planétaires
induisent également une variation des éléments osculateurs
des orbites. Il est, de ce fait, indispensable de préciser la date
de référence des éléments et repéres
utilisés.
D'autre part, la vitesse finie de la
lumiére induit les phénoménes connus sous le nom d'aberration
:
-
aberration du Soleil et des étoiles
(aberration dite "des fixes"), liée à la vitesse de la Terre
sur son orbite (la constante de l'aberration vaut 20"49552 pour J2000.0).
Il importe de noter des "composantes secondaires" dues à l'excentricité
de l'orbite terrestre (voisine de 0"342) et à l'influence des autres
planétes. Cette aberration n'intervient que pour la réduction
d'observations astrométriques. L'aberration solaire est à
inclure dans le calcul de la position de la Terre au moment de l'observation.
-
aberration planétaire, due au temps
de trajet de la lumiére entre le corps observé et la Terre
(499.0047815 secondes/unité astronomique), pendant lequel ce corps
a bougé. L'expression "temps de lumiére" est également
utilisée pour désigner ce phénoméne.
Ces divers phénoménes ont
amené à définir plusieurs qualificatifs pour les coordonnées,
afin de préciser lesquels sont pris en compte :
-
coordonnées moyennes : liées
à une date, elles ne prennent en compte que la précession;
-
coordonnées vraies : se déduisent
des précédentes par la prise en compte de la nutation;
-
coordonnées apparentes : se déduisent
des précédentes par la prise en compte de l'aberration "des
fixes".
Les coordonnées astrométriques
publiées sont les coordonnées moyennes corrigées du
temps de lumiére (aberration planétaire).
Le repére horizontal peut accessoirement
être utilisé pour s'assurer de la visibilité locale
de l'objet.
2) Eléments osculateurs
Ce sont les éléments
orbitaux transmis dans les circulaires UAI :
-
T instant du périhélie;
-
q distance du périhélie;
-
e excentricité;
-
w argument
du périhélie : distance angulaire dans l'orbite entre le
noeud ascendant (cf ci après) et le périhélie;
-
W longitude
du noeud ascendant : point où le corps entre dans l'hémisphére
écliptique Nord;
-
i inclinaison de l'orbite sur l'écliptique;
Il est possible de trouver le demi grand
axe a à la place de la distance du périhélie q [q=a.(1-e)].
En toute rigueur, on devrait trouver également la date pour laquelle
sont donnés les éléments osculateurs, ce qui permettrait
de donner également l'anomalie moyenne pour cette date à
la place de l'instant du périhélie. En l'absence de cette
information, il me semble logique de prendre l'instant du périhélie
comme date pour les éléments osculateurs, en considérant
qu'ils doivent être proches, par convention sinon dans les faits.
3) Loi de magnitude
Pour les cométes :
m1 = H + 5 log(D)
+ 2.5 G log(R)
avec :
-
H : magnitude à 1 UA de la Terre
(D) et du Soleil (R)
-
G : facteur d'activité (souvent
égal à 4, mais reste à déterminer par régression
linéaire chaque fois que possible).
Pour les astéroïdes :
m = H + 5 log(R) + 5 log(D)
+ 2.5 log[(1-G).F1+G.F2]
avec
-
H : magnitude à 1 UA de la Terre
(D) et du Soleil (R)
-
G : "paramétre de pente" (slope
parameter), pris égal à 0.15 par défaut.
-
F1 : exp(-3.33 tg(b/2)^0.63)
-
F2 : exp(-1.87 tg(b/2)^1.22)
-
^ : puissance
-
b : angle de phase
B) Cas non perturbé
1) Orbites elliptiques
On pose : M = k.(t-T)/a.sqrt(a)
-
a : demi grand axe => a=q/(1-e).
-
q : distance au périhélie.
-
t : instant de l'observation.
-
T : instant du périhélie.
-
k = 0.01720209895 (constante de Gauss,
à la base du systéme UAI).
-
M : anomalie moyenne.
-
u : anomalie excentrique.
-
v : anomalie vraie.
-
sqrt() = racine carrée.
On résoud ensuite l'équation
de Képler : u-e.sin(u)=M par itérations successives, en posant
d'abord :
u(0) = M
puis :
u(n+1) = u(n) + [(M - u(n)
+ e.sin(u(n)))/(1 - e.cos(u(n)))]
jusqu'à ce que la valeur absolue
de la différence u(n+1)-u(n) soit inférieure à une
limite fixée à l'avance.
On passe ensuite aux quantités
:
X = r.cos(v) = a.(cos(u)
- e)
Y = r.sin(v) = a.sqrt(1-e*e).sin(u)
2) Orbites paraboliques
On pose l'équation : s^3 + 3.s
= 3.M = 3.k.(t-T)/(q.sqrt(2.q)), que l'on résout par les quantités
annexes :
-
A = 3.M/2
-
B = (A+sqrt(1+A.A))^(1/3)
-
s = B-(1/B) = tan(v/2)
On obtient ensuite :
3) Orbites hyperboliques
On pose :
M = k.(t-T)/a.sqrt(a)
avec a=q/(e-1)
On résoud ensuite l'équation
de Képler : u-e.sin(u)=M par itérations successives, en posant
d'abord :
u(0) = M
puis :
u(n+1) = u(n) + [(M - u(n)
+ e.sh(u(n)))/(e.ch(u(n)) - 1)]
jusqu'à ce que la valeur absolue
de la différence u(n+1)-u(n) soit inférieure à une
limite fixée à l'avance.
On passe ensuite aux quantités
:
-
X = r.cos(v) = a.(e - ch(u))
-
Y = r.sin(v) = a.sqrt(e*e - 1).sh(u)
4) Passage à la position
apparente
Ayant déterminé la position
dans l'orbite, on passe à la position apparente par la suite de
calculs :
On définit les quantités
:
-
Px' = cos(W).cos(w)-sin(W).sin(w).cos(i)
-
Py' = sin(W).cos(w)+cos(W).sin(w).cos(i)
-
Pz' = sin(w).sin(i)
-
Qx' = -cos(W).sin(w)-sin(W).cos(w).cos(i)
-
Qy' = -sin(W).sin(w)+cos(W).cos(w).cos(i)
-
Qz' = -cos(w).sin(i)
-
Wx' = sin(W).sin(i)
-
Wy' = -cos(W).sin(i)
-
Wz' = cos(i)
Puis on passe de la position dans l'orbite
à la position dans le repére écliptique par les relations
:
-
x' = Px'.X + Qx'.Y
-
y' = Py'.X + Qy'.Y
-
z' = Pz'.X + Qz'.Y
On calcule ensuite la position de la Terre
dans le même repére au même instant. On prend la différence
(le vecteur étant orienté de la Terre vers la cométe)
en tenant compte du temps de lumière entre le corps et la Terre,
puis l'on transpose dans le repére équatorial.
Note : Les quantités Px',
Py', Pz' et Qx', Qy', Qz' sont équivalentes aux éléments
de Thiele - Innes - van den Bos utilisés dans le calcul des éphémérides
d'étoiles doubles. Dans ce contexte, elles sont multipliées
par le demi grand axe.
C) Cas perturbé
1) Conditions initiales
On prend comme position de départ
celle correspondant à l'instant pour lequel sont donnés les
éléments osculateurs. Dans le cas où cette date n'est
pas donnée, on prendra, par défaut, l'instant du périhélie
comme instant de référence. Outre la position initiale, il
faut connaître la vitesse à l'instant origine. On obtient
celle dans l'orbite par les relations :
-
Vx = k.(e.X-p).Y/(sqrt(p).R^2)
-
Vy = k.(e.Y^2+p.X)/(sqrt(p).R^2)
où p est le paramètre de
l'orbite, défini par : p = q.(1+e) ou p = a.(1-e^2).
On transpose ensuite ce vecteur vitesse
dans le repére écliptique (ou équatorial) par des
relations identiques à celles utilisées pour les positions.
2) Intégration numérique
Connaissant les conditions initiales,
il est possible de passer à l'intégration numérique
des équations différentielles du mouvement. Celles ci étant
du second degré, sans termes dépendant de la vitesse (sauf
peut être pour l'orientation instantanée de l'orbite, utile
dans le cas des forces non gravitationnelles), j'utilise la méthode
décrite par Runge Kutta et Nyström :
-
k1 = 0.5*h*f(t,f,df)
-
k2 = 0.5*h*f(t+0.5*h,f+0.5*h*(df+0.5*k1),df+k1)
-
k3 = 0.5*h*f(t+0.5*h,f+0.5*h*(df+0.5*k1),df+k2)
-
k4 = 0.5*h*f(t+h,f+h*(df+k3),df+2*k3)
avec :
-
f = f + h*(df+(k1+k2+k3)/3)
-
df = df + (k1+2*k2+2*k3+k4)/6
où f représente le vecteur
position, df le vecteur vitesse et h le pas d'intégration. k1, k2,
k3 et k4 sont des vecteurs intermédiaires de travail. La valeur
exacte de ce pas est fixée, au moment du calcul, de telle maniére
qu'il y ait un nombre entier de pas.
Note sur les forces non gravitationnelles
:
Si j'ai bien compris les articles de
B.G. MARSDEN et al. dans "Astronomical Journal", les forces non gravitationnelles
s'expriment sous la forme :
Ai.exp(-Bi.dt).exp(-0.5*R^2)/R^2
avec :
-
i = 1 pour la composante radiale;
-
i = 2 pour la composante tangentielle
dans le plan de l'orbite;
-
i = 3 pour la composante normale à
l'orbite;
Je n'ai pas relevé de valeurs pour
cette dernière composante dans les articles de MARSDEN : c'est pour
cette raison je ne demande et n'utilise pas dans mon programme les quantités
A3 et B3. Les paramètres A1, A2, B1 et B2 ne sont pris en compte
que
pour le principe, n'ayant pas relevé les valeurs correspondantes
dans les circulaires ou messages en ma possession.