Troisième loi de Kepler (Python-Skyfield)

Où l’on calcule le demi grand-axe et la période orbitale des planètes, pour mettre en évidence la troisème loi de Kepler.

Ce code Python utilise les modules suivants :

  • Skyfield (calculs astronomiques)
  • Matplotlib (tracé de graphiques)
  • sklearn (régressions linéaires)
  • numpy (calcul numérique)
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Créé le Sat Feb 26 09:57:44 2022

@auteur: David ALBERTO
(www.astrolabe-science.fr)
"""
import matplotlib.pyplot as plt
import matplotlib.patches as patches# rectangles du tableau
import numpy as np
from sklearn.metrics import r2_score#pour calculer le coeff R2

# éphémérides astronomiques :
from skyfield.api import load
from skyfield.elementslib import osculating_elements_of

plt.rc('text', usetex=True)
# plt.rcParams["font.weight"] = "bold"
# plt.rcParams["font.family"] = "Ubuntu"# ou autre police installée

ts = load.timescale()
t = ts.now()#date de la compilation

planets = load('de421.bsp')# charge un fichier d'éphémérides (valable 1900-2050)
soleil=planets['Sun']
listeplanetes=['mercury','venus','earth','mars']+['jupiter','saturn','uranus','neptune']
listePlanetesFR=['Mercure','Vénus','Terre', 'Mars']+['Jupiter', 'Saturne', 'Uranus', 'Neptune']

# Création des listes de valeurs de a et de T :
liste_T=[]#liste des périodes orbitales
liste_a=[]#liste des demi grands-axes

for nom in listeplanetes:
    planete=planets[nom+' barycenter']
    position=soleil.at(t).observe(planete)#définir la position de la planète depuis le centre du Soleil
    elements = osculating_elements_of(position)#calcule a et T de la planète
    a = elements.semi_major_axis.au
    liste_a.append(a)
    T = elements.period_in_days/365.25 #conversion de T en années terrestres
    liste_T.append(T)
# ----------------------------------------------
# Paramètres du graphique :

fig,axs=plt.subplots(nrows=2, ncols=2,tight_layout=True)
axs=axs.flatten()#pour accéder à  chaque graphe par axs[N]
fig.set_size_inches(29.7/2.54,21/2.54)#format A4 paysage

plt.suptitle(r"$\textbf{Troisième loi de Kepler}$",fontsize=18)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ax0=axs[0]# AXE :  TABLEAU
ax0.set_ylim(-2.5,len(listePlanetesFR))
ax0.invert_yaxis()
ax0.set_xlim(-1,2.5)
ax0.set_xticks([])
ax0.set_yticks([])

proptexte=dict(fontsize=14,ha='right')#format du texte du tableau
proptitre=dict(fontsize=14,ha='right',c='crimson',fontweight='bold')#format des titres
ax0.text(0,-1.5,r'$\textbf{planète}$',**proptitre)
ax0.text(1,-1.5,r'$\textbf{a (u.a.)}$',**proptitre)
ax0.text(2,-1.5,r'$\textbf{T (ans)}$',**proptitre)

ax0.plot([-0.5,2.05],[-1,-1],lw=0.7,c='k')

for N in range(len(listePlanetesFR)):
    ax0.text(0,N,listePlanetesFR[N],**proptexte)
    ax0.text(1,N,'%.3f'%(liste_a[N]),**proptexte)
    ax0.text(2,N,'%.3f'%(liste_T[N]),**proptexte)
    if N%2==0:#une bande bleue en fond (planètes de rang pair)
        ax0.add_patch(patches.Rectangle(
        (-0.5, N+0.25),2.55,-1,
        facecolor = 'lightsteelblue', alpha=0.5,zorder=-1,
     ) )

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ax1=axs[1]# AXE : T en fonction de a

prop_scatter=dict(s=15,edgecolor='k')#format des points
ax1.scatter(liste_a,liste_T,**prop_scatter)
ax1.set_xlabel("demi grand-axe a    (u.a.)")
ax1.set_ylabel("période orbitale T (ans)")

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ax2=axs[2]# AXE :  ln (T) = f (ln(a))

prop_scatter=dict(s=30,edgecolor='k')#format des points
c_regression='green'# couleur régression

ax2.set_xlabel(r"ln (a)")
ax2.set_ylabel(r"ln (T)")
ln_T=np.log(liste_T)
ln_a=np.log(liste_a)
ax2.scatter(ln_a,ln_T,**prop_scatter)#trace les points

# Régression linéaire :
x=ln_a
y=ln_T
fit=np.polyfit(x, y, 1)#polynôme degré 1
fit_fn=np.poly1d(fit)
pente=fit[0]
ordo=fit[1]
R2=r2_score(y, fit_fn(x))
ax2.plot(x, fit_fn(x), '--',c=c_regression, zorder=0)#trace la droite de régression
ax2.text(0.99, 0.01, r'$R^2$={:.8f}'.format(R2), ha='right', transform=ax2.transAxes)
ax2.text(0.99, 0.2, 'pente : {:.5f}'.format(pente), ha='right', transform=ax2.transAxes)
ax2.text(0.99, 0.1, 'ordonnée : {:.5f}'.format(ordo), ha='right', transform=ax2.transAxes)

ax2.scatter(0,0,s=10, c='r')#point origine du graphe
ax2.text(0,0,'  (0,0)',ha='left',va='top', c='r')
ax2.text(0.5,0.90,r"$\ln (T)=\frac{3}{2}\cdot \ln (a)$",c=c_regression,
         transform=ax2.transAxes, ha='center', fontsize=16)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ax3=axs[3]# AXE :  T^2 = a^3

c_regression='crimson'# couleur régression

# Pour calculer les puissances, il faut convertir les listes en séries Numpy :
liste_T=np.array(liste_T)
liste_a=np.array(liste_a)
ax3.set_xticks([])
ax3.set_yticks([])
ax3.set_xlabel(r"$a^3$")
ax3.set_ylabel(r"$T^2$")
ax3.scatter(liste_a**3,liste_T**2,**prop_scatter)

# Régression linéaire :
x=liste_a**3
y=liste_T**2
fit=np.polyfit(x, y, 1)
fit_fn=np.poly1d(fit)
pente=fit[0]
ordo=fit[1]
R2=r2_score(y, fit_fn(x))
ax3.plot(x, fit_fn(x), '--',c=c_regression, zorder=0)#trace la droite de régression
ax3.text(0.99, 0.01, r'$R^2$={:.8f}'.format(R2), ha='right', transform=ax3.transAxes)
ax3.text(0.99, 0.2, 'pente : {:.5f}'.format(pente), ha='right', transform=ax3.transAxes)
ax3.text(0.99, 0.1, 'ordonnée : {:.5f}'.format(ordo), ha='right', transform=ax3.transAxes)
ax3.text(0.5,0.90,r"$T^2=a^3$",c=c_regression, transform=ax3.transAxes,
         ha='center',fontsize=16)

fig.text(0.01,0.01,r"D. \textsc{Alberto} (www.astrolabe-science.fr)",fontsize=8)

fig.savefig("Kepler3.pdf",dpi=300)
fig.savefig("Kepler3.png",dpi=300)

Code Python (compressé) :

Soyez le premier à commenter

Laisser un commentaire