Détection et prévision d’une situation d’éclipse

Dans le cadre de mon projet de logiciel de simulation de mouvements planétaires, j’ai du trouver un moyen de prévoir à quelle date certains astres se trouveraient en situation de syzygie (éclipse, transit).

On considère trois sphères (nos trois astres), S1, S2, S3, de rayons respectifs R1, R2, R3. On connait les vecteurs position des centres \vec{p}_1, \vec{p}_2, \vec{p}_3 de chacune des sphères. On note S1 et S2 les astres depuis lesquels le phénomène d’éclipse peut être observé, en appelant S1 l’astre de plus petit rayon parmi S1 et S2 (R_1 \leq R_2). On appelle S3 l’astre dont on cherche à savoir s’il se trouve « entre » S1 et S2.

 Première méthode : solution « approchée » géométrique

La première méthode de résolution proposée ici donne une bonne approximation du résultat exact, dans les conditions étudiées.
La condition d’occultation est l’intersection entre le volume de la sphère 3 et le volume du cône tangent aux sphères 1 et 2.

On se place dans le plan (S1, S2, S3). Le problème consiste désormais à montrer que le cercle se situe en partie entre les deux tangentes communes aux deux sphères.
On construis la droite (O1, O2), puis les points T1 et T2, tels que (T_1O_1) \perp (O_1O_2), T_1 \in S_1 et (T_2O_2) \perp (O_1O_2), T_2 \in S_2.
Ainsi, la droite (T1, T2) est une bonne approximation d’une tangente commune aux deux droites. Cela est vrai tant que la différence des rayons des astres 1 et 2 est petite devant la distance qui la sépare.

Schéma sphères 1

On notera d = [O_1,O_2] = ||\vec{p}_2-\vec{p}_1|| et \varepsilon = \frac{R_2-R_1}{d} (d’après notre choix de notation, \varepsilon > 0. La condition se réécrit, mathématiquement : |\varepsilon| \ll 1. En pratique, dans notre cas, elle est vérifiée : les distances qui séparent deux astres sont d’un ordre de grandeur bien plus élevé que celui de leur taille.

Dès lors, la condition d’occultation devient alors (approximativement) : O_3B \leq HB+R_3
Le calcul de O_3B est simple. C’est la distance de O_3 à la droite (O1, O2). D’où :

O_3B = \dfrac{||(\vec{p}_3-\vec{p}_1) \wedge (\vec{p}_2-\vec{p}_1)||}{d}

Le calcul de HB se décompose en HA+AB. AB n’est que le rayon de la sphère 1 (car c’est la plus petite, d’après notre choix de notation), donc AB = R_1.

D’autre part, le théorème de Thalès donne : HA = (R_2-R_1)\dfrac{O_1B}{d} = \varepsilon.O_1B

Par projection scalaire : O_1B = |\dfrac{(\vec{p}_3-\vec{p}_1).(\vec{p}_2-\vec{p}_1)}{d}|

La condition devient donc, si \varepsilon \ll 1 :

\dfrac{||(\vec{p}_3-\vec{p}_1) \wedge (\vec{p}_2-\vec{p}_1)||}{d} \leq R_1+R_3+\varepsilon. |\dfrac{(\vec{p}_3-\vec{p}_1).(\vec{p}_2-\vec{p}_1)}{d}|, soit :

||(\vec{p}_3-\vec{p}_1) \wedge (\vec{p}_2-\vec{p}_1)|| \leq d(R_1+R_3)+\varepsilon. |(\vec{p}_3-\vec{p}_1).(\vec{p}_2-\vec{p}_1)|

Cette condition est celle que j’utilise actuellement dans mon programme

Deuxième méthode : solution « exacte » (analytique)

Bien que la première solution donne des résultats satisfaisants et permette de prévoir des éclipses et transits à quelques mois/ans avec une précision de 2-3 heures, il est intéressant d’étudier le cas général. J’ai pour cela demandé de l’aide à mon prof de maths qui m’a proposé d’utiliser les équations des tangentes dans un système de coordonnés simplifié. Cette démonstration s’inspire de sa méthode mais s’applique à n’importe quel système de coordonnées.

On se place toujours dans le même plan que précédemment. Mais cette fois, on le muni d’un repère (O_1, \vec{i}, \vec{j}) orthonormé et direct. Volontairement, on se place dans une situation ou la condition de validité de l’approximation précédente n’est plus respectée.

Schéma sphères

On appelle A=(x_A, 0) l’intersection des tangentes avec l’axe des abscisses. Puisque (O_1T_1)\parallel(O_2T_2) alors d’après le théorème de Thalès :

\dfrac{x_A}{x_A+d} = \dfrac{R_1}{R_2} d’où x_A = -\dfrac{R_1}{\varepsilon}. Lorsque \varepsilon \to 0, A est logiquement rejeté à l’infini.

On souhaite déterminer l’équation des tangentes communes aux deux cercles, de la forme xx_0 \pm yy_0 = R_1^2, sachant que T_1 = (x_0,y_0) \in S_1. Puisque A appartient à ces tangentes on peut écrire : -x_0 . \dfrac{R_1}{\varepsilon} = R_1^2 d’où x_0 = -R_1.\varepsilon.

Mais T_1 \in S_1 \Rightarrow R_1^2\varepsilon^2 + y_0^2 = R_1^2 d’où y_0 = \pm R_1 \sqrt{1-\varepsilon^2}. Ainsi les équations des deux tangentes commune sont :

\varepsilon.x \pm \sqrt{1-\varepsilon^2}.y = R_1

On rappelle la condition d’occultation : O_3B \leq HB+R_3

Soit \mathcal{D} la perpendiculaire à une de ces tangentes (prenons par exemple celle de pente positive) passant par O_3 = (x_3,y_3). Elle coupe la tangente en H et l’axe des abscisses en B = (x_B, 0).

Cette droite a donc pour équation : \sqrt{1-\varepsilon^2}x-\varepsilon.y = \sqrt{1-\varepsilon^2}x_3-\varepsilon. y_3.

On en déduit les coordonnées de B\in\mathcal{D} \mbox{ : } B = (x_3-\dfrac{\varepsilon}{\sqrt{1-\varepsilon^2}}y_3,0). On en déduit :

O_3B = \sqrt{y_3^2 + y_3^2.\dfrac{\varepsilon^2}{1-\varepsilon^2}} = \dfrac{y_3}{\sqrt{1-\varepsilon^2}}

Reste à calculer HB. Cela se fait par application du théorème de Thalès :

\dfrac{HB}{R_2} = \dfrac{x_B+AO_1}{d+AO_1} soit

HB = R_2.\dfrac{x_3-\dfrac{\varepsilon}{\sqrt{1-\varepsilon^2}}y_3+R_1/|\varepsilon|}{d+R_1/\varepsilon} = R_1 + |\varepsilon| ( x_3-\dfrac{\varepsilon}{\sqrt{1-\varepsilon^2}}y_3 )

Nous devons calculer x_3 et y_3, les coordonnées de O_3 dans le repère (O_1, \vec{i}, \vec{j}), sachant ses coordonnées dans une base (O, \vec{u}, \vec{v}, \vec{\omega}) quelconque.
L’abscisse x_3 est en fait la projection de \vec{p}_3-\vec{p}_1 sur la droite (O_1O_2). On  a alors :

x_3 = \dfrac{(\vec{p}_3-\vec{p}_1).(\vec{p}_2-\vec{p}_1)}{d}

Quant à y_3, son signe n’a pas d’importance (par symétrie par rapport à Ox), aussi on peut le calculer comme étant la distance de O_3 à (O_1O_2) :

y_3 = \dfrac{||(\vec{p}_3-\vec{p}_1) \wedge (\vec{p}_2-\vec{p}_1)||}{d}

La condition d’occultation s’écrit alors \dfrac{y_3}{\sqrt{1-\varepsilon^2}} \leq R_1 + R_3 + |\varepsilon| ( x_3-\dfrac{\varepsilon}{\sqrt{1-\varepsilon^2}}y_3 ) soit :

\dfrac{||(\vec{p}_3-\vec{p}_1) \wedge (\vec{p}_2-\vec{p}_1)||}{\sqrt{1-\varepsilon^2}} \leq d(R_1 + R_3) + \varepsilon \left( (\vec{p}_3-\vec{p}_1).(\vec{p}_2-\vec{p}_1)-\dfrac{||\varepsilon(\vec{p}_3-\vec{p}_1) \wedge (\vec{p}_2-\vec{p}_1)||}{\sqrt{1-\varepsilon^2}} \right)

Voilà le résultat sans approximation que nous recherchions. Si nous négligeons les termes en epsilon d’ordre 2 (\varepsilon^2) nous retrouvons notre première approximation. L’erreur est alors faible et le calcul plus rapide (moins de racine(s) carrée(s) à calculer). Dans le cas d’une planète passant entre la Terre et le Soleil, \varepsilon \approx 0,0046.

Application

Je me suis servi de la formule approchée dans mon programme, ce qui a permis de prévoir la date du transit de vénus de juin 2012 à partir de la disposition du système solaire en mars de la même année. Le programme utilise des données issues de kernels SPICE. L’état initial du système solaire a été pris le 2012 MAR 1 12:00 puis l’évolution en a été déduite par résolution du problème à N corps (méthode d’Euler, pas de 1 seconde). J’ai de plus tenu compte des effets dus au temps de propagation de la lumière. (Pour cela, il suffit d’effectuer les calculs à partir de \vec{p}_2 = \vec{p}_2(t - ||\vec{p}_2 - \vec{p}_1||/c) et idem pour \vec{p}_3).

Valeur calculée Observé
Début transit (UTC) 05/06/2012 22:04 05/06/2012 22:09
Fin transit (UTC) 06/06/2012 04:57 06/06/2012 04:49
Diamètre apparent Vénus 58″ 58″

Les résultats sont donc satisfaisants. L’heure de décalage s’explique par l’erreur accumulée sur les positions des planètes.

J’ai aussi effectué un test pour le transit de mercure de 2016 à l’aide de l’état initial du système solaire au premier mai de cette année. Ce transit devant avoir lieu le 9 mai, l’interpolation par calcul des effets gravitationnels ne se fait que sur 9 jours, ce qui permet de réduire les erreurs dues à celle-ci.

Valeur calculée
vitesse lumière infinie
Valeur calculée Prédiction NASA
Début transit (UTC) 09/05/2016 11:03 09/05/2016 11:10 09/05/2016 11:12
Fin transit (UTC) 09/05/2016 18:37 09/05/2012 18:43 06/06/2012 18:42

Les résultats très bons (2 minutes de décalage seulement). Cet exemple montre bien que pour une telle précision il est nécessaire de tenir compte du temps de propagation de la lumière.