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