Orbite de la comète de Halley (Python)

Le code Python proposé ici utilise les modules Python suivants :

• matplotlib (graphiques)
• datetime (création et gestion des dates)
• ephem (éphémérides astronomiques)
• numpy (gestion de listes)

Le résultat :

Quelques commentaires sur le code :

import matplotlib.pyplot as plt
import datetime as dt
import numpy as np
import ephem

#coordonnées de la comète de Halley (code du logiciel XEphem,avec les paramètres orbitaux)
# Source : http://astro.vanbuitenen.nl/cometelements?format=xephem&mag=obs

line="1P/Halley,e,162.2239,58.9763,111.9047,17.872265,0.0130447,0.966321,171.4710,02/10.0/2022,2000,g 4.0,6.0"
halley = ephem.readdb(line)#attention : pas de tirets dans le nom de la comète

angle=0.93#trouvé par tâtonnements, pour rendre horizontal le demi grand-axe de l'ellipse
#cet angle est à adapter à l'époque choisie pour les dates.

Pour les petits corps célestes, les commandes du module Ephem sont moins évidentes que pour les planètes et la Lune. Je vous renvoie à la notice du module. Pour la comète souhaitée, il faut rentrer une obscure suite de nombres, qui obéissent à un format spécifique du logiciel XEphem (sur lequel se base le module Ephem).

Heureusement, j’ai trouvé sur le web les codes de plusieurs comètes avec ce format.

Pour l’exemple de Halley, je crois reconnaître quelques valeurs comme le demi grand-axe (17.87 ua), l’excentricité (0.966), peut-être la distance du périhélie (0.013 ua ?).

On pourra remplacer par une autre comète, mais toute la mise en page du graphique sera à adapter.

Avec le module datetime, si vous avez besoin d’afficher le nom du mois, il sera par défaut en anglais. Pour afficher les noms des jours et des mois en français, il faut ajouter ces deux lignes en début de code :

import locale
locale.setlocale(locale.LC_ALL, 'fr_FR.UTF-8')#dates en français

# -------------------Création de la liste de dates :

#Liste de dates au 1er janvier de l'année :
listeDates=[dt.datetime(1985+N,1,1) for N in range(76)]

Voici une possibilité pour créer la liste de dates. Si on prend plus d’années que la période, les points se superposent et rendent peu lisible le document final. De plus, on remarque que la comète ne repasse pas exactement au même point.

Les coordonnées de positions utilisées ici sont la longitude écliptique et la distance au Soleil. Pour une fois, vue l’orbite très allongée, je transforme ces coordonnées en coordonnées rectangulaires.

prop1=dict(s=taillePoint,c=couleur_comete,edgecolor='k',lw=0.2,zorder=0)#format du point année
prop2=dict(s=2*taillePoint,c='crimson',edgecolor='k',lw=0.5,zorder=1)#format du point 5 ans

for date in listeDates:
    halley.compute(date)
    R=halley.sun_distance
    theta=halley.hlon#longitude écliptique
    # Passage aux coordonnées cartésiennes :
    X=R*np.cos(theta+angle)
    Y=R*np.sin(theta+angle)
    plt.scatter(X,Y,**prop1)
    if int(date.year)%5==0:#pour les années en multiple de 5 : étiquette de l'année
    # calcul du décalage des étiquettes par rapport aux points
        X_offset=(X+18)/abs(X+18)#décalage +1 si X>-18, sinon -1
        Y_offset=Y/abs(Y)#décalage +1 si Y positif, sinon -1
        plt.annotate(date.year,xy=(X,Y),xytext=(X+X_offset,Y+Y_offset),fontsize=6,arrowprops=dict(arrowstyle="-", linewidth=0.5,color='gray'))
        plt.scatter(X,Y,**prop2)

Sur le site suivant, on peut accéder à des paramètres orbitaux plus détaillés (mais peu expliqués…), et simuler la trajectoire au fil du temps.

Rendez-vous en 2061 pour le prochain passage !