import numpy as np
import scipy as sc
import matplotlib
import matplotlib.pyplot as plt
import datetime
import pandas as pd
import plotly.express as px
import time
#Option d'affichage et taille de police sur les figures:
fs=20
plt.style.use('fivethirtyeight')
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
plt.rc('text', usetex=True)
plt.rcParams["font.family"] = "Ubuntu"
## Constantes utilisées:
c0=299792458 # Vitesse de la lumière dans le vide en m/s
# Diametre solaire en m
Dsol=1392680000
# Masse solaire en kg
Msol=1.9884e30
# Constante de Gravitation
G=6.6743015e-11
# Parametre solaire GM0:
mu0=G*Msol
Cette première partie du TN présente les principales caractéristique physique et orbitales des 8 planétes du système solaire.
Le fichier Syst_solaire.csv contient les principales données orbitales et physiques des 8 planètes du système solaire que l'on peut mémoriser dans l'ordre d'éloignement par les première lettre de la phrase "Mercredi Viendra Tu Manger Avec Jean Sur Une Nappe ?" - Le A étant la ceinture d'Astéroïdes.
En regardant le système depuis le pôle nord terrestre, toutes les planètes tounent dans le sens antihoraire autour du soleil. Leur sens de rotation autour de leur axe est également antihoraire sauf pour Venus et Uranus pour lesquelles un observateur à leur surface verrait le soleil se lever à l'Ouest et se coucher à l'Est.
Comme dans le TN1 on utilise pandas pour explorer facilement les données du fichier:
#! wget https://hpp.education/Lessons/MecaSpace/Files/Syst_solaire.csv
SystSol=pd.read_csv('Syst_solaire.csv',delimiter=" ")
SystSol["Size"]=SystSol["Diamètre(km)"]/SystSol["Diamètre(km)"][2]
SystSol
Planete | Masse(kg) | Diamètre(km) | Demi-GA(km) | Trot(s) | exentricity | Axial-Tilt(°) | Orbit-incl | Size | |
---|---|---|---|---|---|---|---|---|---|
0 | Mercure | 3.290000e+23 | 4879.0 | 5.790905e+07 | 5.066928e+06 | 0.20560 | 0.035200 | 7.00000 | 0.382487 |
1 | Venus | 4.870000e+24 | 12104.0 | 1.082095e+08 | -2.099892e+07 | 0.00678 | 177.360000 | 3.39471 | 0.948887 |
2 | Terre | 5.970000e+24 | 12756.0 | 1.495979e+08 | 8.616040e+04 | 0.01671 | 23.436691 | 0.00000 | 1.000000 |
3 | Mars | 6.290000e+23 | 6794.0 | 2.279440e+08 | 8.864266e+04 | 0.09339 | 25.190000 | 1.85000 | 0.532612 |
4 | Jupiter | 1.900000e+27 | 142080.0 | 7.783400e+08 | 3.572730e+04 | 0.04839 | 3.120000 | 1.30400 | 11.138288 |
5 | Saturne | 5.680000e+26 | 120536.0 | 1.426700e+09 | 3.798000e+04 | 0.05390 | 26.730000 | 2.48600 | 9.449357 |
6 | Uranus | 8.680000e+25 | 50800.0 | 2.870700e+09 | -6.206371e+04 | 0.04726 | 97.800000 | 0.77300 | 3.982440 |
7 | Neptune | 1.020000e+26 | 50538.0 | 4.498400e+09 | 5.760660e+04 | 0.00859 | 28.320000 | 1.77000 | 3.961900 |
fig=px.scatter(SystSol, hover_name="Planete",x="Demi-GA(km)", y="Masse(kg)",color="Trot(s)",size="Size",color_continuous_scale='reds_r')
fig.show()
Il est également important de connaitre les ordres de grandeur des données orbitales par rapport à celles de la terre:
# Diamètres par rapport à celui de la terre:
SystSol["Diamètre(km)"]/SystSol["Diamètre(km)"][2]
0 0.382487 1 0.948887 2 1.000000 3 0.532612 4 11.138288 5 9.449357 6 3.982440 7 3.961900 Name: Diamètre(km), dtype: float64
#Distances en UA:
SystSol["Demi-GA(km)"]/SystSol["Demi-GA(km)"][2]
0 0.387098 1 0.723336 2 1.000000 3 1.523711 4 5.202881 5 9.536899 6 19.189442 7 30.069943 Name: Demi-GA(km), dtype: float64
Calculer les demi grand axes des planètes en minutes lumière:
#Distance en minutes lumière:
dml=...
print(dml)
Ellipsis
L'équation d'une conique (Cf Chap 3.) s'écrit en coordonnées polaires:
$$ \rho=\dfrac{p}{1+e\cos(\theta)} $$Avec:
soit donc:
$$ \rho=\dfrac{a(1-e^2)}{1+e\cos(\theta)} $$En prenant un vecteur azimutale $\theta$ entre $0$ et $2\pi$ on peut tracer l'orbite des planètes (on prendra ici les planètes telluriques):
# Vecteur Azimutale:
theta=np.linspace(0,2*np.pi,100)
# Creation de la figure
fig=plt.figure(figsize=(16,12))
plt.axes(projection='polar')
# On représente les orbites des 4 premières planètes
for i,plan in SystSol[0:4].iterrows():
a=SystSol["Demi-GA(km)"][i]/SystSol["Demi-GA(km)"][2]
p=a*(1-SystSol["exentricity"][i]**2)
plt.plot(theta,p/(1+SystSol["exentricity"][i]*np.cos(theta)),label=SystSol["Planete"][i])
plt.legend(fontsize=fs)
<matplotlib.legend.Legend at 0x7f0892601970>
Le soleil est un des foyers des orbites planétaires. On s'interresse maintenant à l'orbite terrestre et à la trajectoire apparente du soleil vu depuis la terre:
terre=SystSol.iloc[2]
Exercice: En utilisant les éphémérides fournis dans le repère héliocentrique ($\lambda$,$\beta$) du plan de l'écliptique, représenter l'état actuel du système solaire et calculer le nouvel état du système dans 1 an. On considérera des orbites circulaires pour chacune des planètes.
Le jour solaire est défini comme la durée que met le soleil pour revenir à la même position dans le ciel d'un jour à l'autre (1->3 sur la figure). La durée du jour solaire est fixée à 86400 secondes (soir 24 heures).
Le jour sidéral est défini comme la durée que met la terre pour faire un tour sur elle même (1->2 sur la figure). En pratique ce n'est pas tout à fait la même chose que le jour solaire car d'un jour à l'autre elle a parcouru un petit arc de cercle dans sa révolution autour du soleil:
On peut donc calculer la durée du jour sidéral en remarquant qu'au bout d'un an la terre aura bien fait un tour de plus sur elle même soit $366.25$ jours. Ainsi on a:
Tsid=86400*(365.25/366.25)
print("Tsid =",Tsid)
Tsid = 86164.09556313994
La durée d'un jour sidéral est donc de 86164 seconde soit 23h 56 min et 4 secondes.
L'orbite terrestre permet de décrire sa position autour du soleil. Cette orbite ne fixe pas de calendrier (on ne peut pas a priori positioner le 1er Janvier ou le 13 Mai sur la figure précédente). Le calendrier est principalement fixé par deux paramètres:
A partir de la troisième loi de Kepler, la periode de révolution d'une planète en fonction de son demi-grand axe $a$ est donnée par:
$$ T^2=\dfrac{4\pi^2 a^3}{GM_{sol}} $$soit:
T=np.sqrt(4*np.pi**2*(1000*terre["Demi-GA(km)"])**3/mu0)
Soit en terme de jour solaire:
T0=T/(86400);print(T0)
365.2578254600518
Ou en terme de nombre de rotation complète:
T1=T/terre["Trot(s)"];print(T1)
366.27355629440524
En pratique la durée réelle de l'année sidérale est de $365,256363004$ jours. Un calendrier consiste donc à fixer un nombre de jour entier dans une année en s'approchant au plus de la valeur réelle:
Les accords entre l'orbite et le calendrier sont alors fixés par la postion de l'axe de rotation de la terre par rapport à l'éclitptique:
Il est important de noter que les positions des solstices ne correspondent pas au périhélie et à l'aphélie qui sont respectivement atteint le 03 janvier et le 03 juillet:
Si l'on considère le système de coordonnée écliptique $(\beta,\lambda)$ défini dans le TN1 dans sa version héliocentrique, les coordonnées de la Terre sont donc $(\beta=0,\lambda)$. Par définition $\lambda=0$ lors de l'équinoxe de printemps. On peut alors définir un calendrier à l'aide d'un entier $N$ variant de $1$ (premier jour du calendrier) à $365$ ou $366$. On considère dans toute la suite une année à 365 jours. La longitude écliptique de la terre est donc donnée en degré par:
$$ \lambda=\dfrac{360}{365} (N-80) $$ou en radian par:
$$ \lambda=\dfrac{2\pi}{365} (N-80) $$En considérant une orbite circulaire de vitesse constante, la hauteur du soleil $h$ dans le système de coordonnée horizontal définis dans le TN1 dépend de l'angle $\delta$ entre les rayons du soleil et le plan de l'équateur (la déclinaison). Cet angle est nulle aux équinoxes et maximum aux solstices. Il s'exprime simplement:
$$ \delta=\epsilon \cos\left(\dfrac{2\pi}{365}(N+10)\right) $$La formule exacte tenant compte de l'execentricité de l'orbite est donnée par:
$$ \delta=\arcsin\left(\sin(-\epsilon)\cos\left(\dfrac{360}{365.2425}(N+10)+\dfrac{360}{\pi} e \sin\left(\dfrac{360}{365.2425}(N-3) \right) \right) \right) $$Avec $\epsilon = 23.44$, l'obliquité de la terre et $e$ l'excentricité de son orbite.
La hauteur $h$ en degrés du soleil dans le ciel au midi solaire s'exprime donc pour une latitude donnée $L$:
$$ h=90-L-\delta $$def Hsol(N,L):
long_ec=2*np.pi*(N+10)/365
delta=23.44*np.cos(long_ec)
return 90-L-delta
N=np.arange(1,366)
fig=plt.figure(figsize=(16,8))
plt.plot(N,Hsol(N,23.44),label="Lat=23.44°")
plt.plot(N,Hsol(N,48.8),label="Lat=48.8° (Paris)")
plt.plot(N,Hsol(N,60),label="Lat=60°")
plt.xlabel("Jour de l'année",fontsize=fs)
plt.ylabel("Hauteur du soleil",fontsize=fs)
plt.legend(fontsize=fs)
<matplotlib.legend.Legend at 0x7f0892503a60>
L'équation du temps désigne l'ensemble des corrections qu'il faut apporter à l'azimute du soleil dans le système de coordonée azimutale pour passer de la position moyenne (basée sur une orbite circulaire uniforme) à la position réelle (basée sur une orbite elliptique).
Ces corrections sont souvent données en minutes et permettent alors de passer de l'heure solaire moyenne à l'heure solaire réelle. Ces corrections sont dues d'une part à l'orbite elliptique de la terre mais également à l'inclinaison de l'axe de rotation de la terre.
L'équation du temps s'écrit alors: $$ \Delta z=C_E+C_O $$
Que l'on exprime souvent en minute en considérant que 24H (24*60 minutes) correspondent à $360°$ on a donc 4 minutes par degré:
$$ \Delta T (')= 4\Delta z $$La version plus précise de l'équation du temps est développé ici.
def correction_excentrique(N):
B=2*np.pi*(N-81)/365
E=7.678*np.sin(B+1.374)
return E
def correction_oblique(N):
B=2*np.pi*(N-81)/365
O=-9.87*np.sin(2*B)
return O
def equation_du_tps(N):
return correction_excentrique(N)+correction_oblique(N)
def Real_time_eq(L,N):
rad=np.pi/180
# Ecart au 1er Janvier de l'époque J2000:
d0=(datetime.date(2022,1,1)-datetime.date(2000,1,1)).days
N=N+d0
M=357.5291+0.98560028*N
#Excentricité de l'orbite Terrestre
e=0.01671
#Obliquité
eps=23.43929
C=(2*e-0.25*e**3)*np.sin(M*rad)+1.25*e**2*np.sin(2*M*rad)+(13/12)*e**3*np.sin(3*M*rad)
l=280.4665+0.98564736*N+C
y=np.tan(eps*rad/2)
R=-y**2*np.sin(2*l*rad)+0.5*y**4*np.sin(4*l*rad)-(y**6/3)*np.sin(6*l*rad)
T=60*(C+R)
#Hauteur solaire
P=Hsol(N,L)
return 4*T,P
#Lattitude:
Lat=48.8667
fig=plt.figure(figsize=(16,8))
plt.plot(N,equation_du_tps(N),'--',label='Version Simplifiée')
plt.plot(N,Real_time_eq(Lat,N)[0],'-',label='Version Précise')
plt.xlabel('Jour de l\'année',fontsize=fs)
plt.ylabel('Décallage en min',fontsize=fs)
plt.grid(True)
plt.legend(fontsize=fs)
plt.show()
On peut également représenter l'équation du temps dans les coordonnées polaires en fonction de la longitude ecliptique (nulle le 21 Mars):
fig=plt.figure(figsize=(14,14))
long_eq=2*np.pi*(N-80)/365
plt.polar(long_eq,equation_du_tps(N),'--',long_eq,Real_time_eq(Lat,N)[0],'-')
[<matplotlib.lines.Line2D at 0x7f0892423d60>, <matplotlib.lines.Line2D at 0x7f0892423f70>]
En représentant sur la même figure la déclinaison du soleil en fonction de l'équation du temps, on obtient une Analemme. Cette figure peut également être obtenue en marquant quotidiennement sur le sol la position de l'ombre d'un gnomon à la même heure ou en faisant un cliché de la position d'un astre tous les jours à la même heure.
En pratique on utilise l'analemme solaire pour passer de l'heure solaire vrai à l'heure légale qui s'obtient en ajoutant à l'équation du temps une correction de longitude pour tenir compte de l'heure au méridien local. En effet, l'heure solaire réelle change avec la longitude. On considérant toujours une correction de 4 minutes par degrés on a donc:
$$ H_{legale}=H_{solaire} + \delta_f + 4l + \Delta T $$Ou $\delta_f$ est la correction due au fuseau horaire (+1 pour l'Europe en hivers, +2 en été), $l$ est la longitude en degrés et $\Delta T$ est l'équation du temps.
# Longitude (CNAM Paris)
long=2.3559
T,P=Real_time_eq(Lat,N)
# On repère le premier jour du calendrier
T0,P0=Real_time_eq(Lat,1)
# On repère le premier jour de chaque mois
N1=np.cumsum([1,31,28,31,30,31,30,31,31,30,31,30]);
T1,P1=Real_time_eq(Lat,N1)
# Liste des mois
month=['Jan','Fev','Mar','Avr','Mai','Juin','Juil','Aout','Sep','Oct','Nov','Dec']
def dtr(d):
return 1/np.tan(d*np.pi/180)
def rtd(r):
return np.arctan2(1,r)
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
fig=plt.figure(figsize=(16,22))
ax=fig.add_subplot(111)
plt.plot(T,P,'.k')
plt.plot(T0,P0,'om',markersize=12)
plt.plot(T1,P1,'or')
plt.title('Analemme solaire 2022 (Lattitude: $'+str(Lat)+'$°) \n Correction Longitude cl: '+str(-4*long)+' min \n \n Heure légale = Heure solaire + 1h(hiver ou 2h été) + cl(min) + $\Delta T$(min)',fontsize=fs)
plt.xlabel('$\Delta T$ (min)',fontsize=fs)
plt.ylabel('Elevation au midi solaire (°)',fontsize=fs)
ax.xaxis.set_minor_locator(AutoMinorLocator())
ax.yaxis.set_minor_locator(MultipleLocator(1))
secaxy = ax.secondary_yaxis('right', functions=(dtr, rtd))
secaxy.set_ylabel(r'$l_o=l/tan(\delta)$',fontsize=fs)
plt.grid(True,which='major',color='k',linestyle='-')
plt.grid(True,which='minor',color='gray',linestyle='--')
for i,n in enumerate(month):
ax.annotate('%s' % n, xy=(T1[i],P1[i]),xytext=(T1[i],P1[i]), xycoords='data',textcoords='data',fontsize=fs)
fig.savefig('Analeme_CNAM_2022.png',bbox_inches='tight',dpi=120)
Applications:
Retrouver les grandeurs associées au calendrier des planètes telluriques du système solaire (durée de l'année, du jour sidéral, du jour solaire, équation du temps...)
Un exemple de mesure du temps sur Mars:
https://fr.wikipedia.org/wiki/Mesure_du_temps_sur_Mars
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)