Déterminer la masse du Soleil en Python

Objectif de la fiche : utiliser un langage de programmation pour étudier la relation approchée entre la variation du vecteur vitesse d'un système modélisé par un point matériel entre deux instants voisins et la somme des forces appliquées sur celui-ci — et en déduire une estimation de la masse du Soleil.

Cette fiche prolonge directement la méthode de la fiche n°3 (vecteurs vitesse et accélération par différences centrées), appliquée cette fois à un mouvement réel à l'échelle du Système solaire. Elle est autonome : tu n'as besoin d'aucun autre document. Chaque étape contient un script incomplet à compléter, puis à copier-coller dans Thonny (ou dans la console Basthon tout en bas de la page si tu n'as pas Thonny sous la main).

Les positions de Mercure

À savoir

Dans l'activité originale, les positions de Mercure sont relevées sur le site des éphémérides de l'IMCCE, pour 11 dates démarrant le 1er janvier 2017 et espacées de 8 jours (soit 80 jours au total, presque un tour complet, puisque Mercure met environ 88 jours à boucler son orbite).
Pour que cette fiche reste autonome sans avoir besoin d'aller chercher ces données en ligne, on les remplace ici par un modèle simplifié : une orbite parfaitement circulaire, de rayon égal au demi-grand axe réel de Mercure (0,387 UA) et de période réelle (87,969 jours). Cela reproduit la même physique (une trajectoire quasi circulaire autour du Soleil) sans passer par le pointage.
Les distances sont exprimées en unités astronomiques (UA), les temps en jours : ce sont des unités pratiques pour manipuler ces grandeurs, mais il faudra les convertir en unités SI (mètres, secondes) avant de pouvoir en tirer une masse en kilogrammes.
Le script ci-dessous est déjà complet : copie-le tel quel dans Thonny (ou dans Basthon) pour commencer.

etape0_donnees.py

import numpy as np
import matplotlib.pyplot as plt

# Modèle d'orbite circulaire pour Mercure (mêmes rayon et période que l'orbite réelle)
r = 0.387        # demi-grand axe, en UA
Torb = 87.969    # période de révolution, en jours

# 11 dates espacées de 8 jours, comme sur le site des éphémérides
t = 8 * np.arange(11)         # jours écoulés depuis le 1er janvier 2017
theta = 2 * np.pi * t / Torb   # angle parcouru sur l'orbite (rad)

x = r * np.cos(theta)          # abscisse (UA)
y = r * np.sin(theta)          # ordonnée (UA)

print('x :', x)
print('y :', y)

Résultat attendu

Deux tableaux de 11 valeurs chacun, oscillant entre environ -0,39 et 0,39 : les coordonnées successives de Mercure sur son orbite, en unités astronomiques.

Représenter les vecteurs vitesse

Consigne

Comme dans la fiche n°3, la vitesse au point i+1 s'estime par différence centrée : (x[i+2]-x[i])/(t[i+2]-t[i]) pour la composante horizontale, et de même pour la composante verticale.
Un facteur d'échelle de 5 est nécessaire ici pour que les flèches soient visibles : sans lui, les vecteurs vitesse (de l'ordre de 0,03 UA/jour) seraient minuscules comparés à la taille de l'orbite (0,39 UA).
Trace ces vecteurs avec plt.arrow, superposés à la trajectoire (plt.plot(x,y,'bo')). Utilise plt.axis('equal') pour que l'orbite ne soit pas déformée par des échelles différentes en x et en y.

etape1_vecteurs_vitesse.py

m = np.arange(len(t) - _)
for i in m :
    plt.arrow(x[i+_], y[i+_], _*(x[i+_]-x[_])/(t[i+_]-t[_]), _*(y[i+_]-y[_])/(t[i+_]-t[_]), head_width=0.05)

plt.plot(x, y, '__')
plt.axis('_____')
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.title("__________________")
plt.grid()
plt.show()

Indice

Le facteur d'échelle 5 se place directement devant chacune des deux longueurs de flèche, comme le facteur 0.1 dans la fiche n°3 — seule sa valeur change, la logique reste identique.

Résultat attendu

Une orbite quasi circulaire, ponctuée de flèches tangentes à la trajectoire, toutes à peu près de la même longueur : une conséquence directe de l'orbite quasi circulaire du modèle utilisé, où la vitesse change surtout de direction, très peu de norme.

Construire les listes de vitesse et d'accélération

Consigne

Pour la suite, il faut cette fois stocker les valeurs de vitesse et d'accélération dans des listes, exactement comme dans les étapes 2 et 3 de la fiche n°3 : d'abord vxl, vyl à partir de x, y, t, puis axl, ayl à partir de vxl, vyl, t.
Le script ci-dessous est déjà complet : c'est une reprise directe de la méthode déjà maîtrisée, tu peux le copier-coller tel quel.

etape2_listes_vitesse_acceleration.py

vxl = []
m = np.arange(len(t) - 2)
for i in m :
    vx = (x[i+2] - x[i]) / (t[i+2] - t[i])
    vxl.append(vx)

vyl = []
m = np.arange(len(t) - 2)
for i in m :
    vy = (y[i+2] - y[i]) / (t[i+2] - t[i])
    vyl.append(vy)

axl = []
m = np.arange(len(t) - 4)
for i in m :
    ax = (vxl[i+2] - vxl[i]) / (t[i+2] - t[i])
    axl.append(ax)

ayl = []
m = np.arange(len(t) - 4)
for i in m :
    ay = (vyl[i+2] - vyl[i]) / (t[i+2] - t[i])
    ayl.append(ay)

Résultat attendu

Aucun affichage, mais à la fin de cette étape, vxl et vyl contiennent 9 valeurs (en UA/jour), et axl, ayl en contiennent 5 (en UA/jour²) — les listes dont on va avoir besoin pour la suite.

Passer aux unités du Système International

Consigne

axl et ayl sont pour l'instant de simples listes Python : convertis-les d'abord en tableaux numpy avec np.array(...), puis passe des UA/jour² aux m/s² en multipliant par 1.5*10**11 (mètres par UA) et en divisant par (24*3600)**2 (secondes par jour, au carré).
Calcule ensuite la norme du vecteur accélération, norma, à partir de ses deux composantes SI : (ayln**2 + axln**2)**(1/2).
Calcule aussi la distance Mercure-Soleil, en mètres, à partir de x et y : dx = x[:-2]*1.5*10**11 (le slicing [:-2] retire les 2 derniers termes, pour que la liste ait la même longueur que les vitesses), puis d = (dx**2+dy**2)**(1/2).
Enfin, redimensionne encore une fois : dr = d[2:] (retire les 2 premiers termes cette fois) pour que les distances correspondent bien aux mêmes instants que les accélérations. Calcule alors GM = dr**2 * norma : d'après la deuxième loi de Newton, ce produit devrait être à peu près constant, égal à G × M_soleil.

etape3_unites_si.py

# TODO : conversion de l'accélération en m/s²
ayln = ___*10**__ * np.array(___) / (24*3600)**_
axln = ___*10**__ * np.array(___) / (24*3600)**_

# TODO : norme du vecteur accélération
norma = (ayln**_ + axln**_)**(1/_)

# TODO : distance Mercure-Soleil, en mètres
dx = x[:-_] * 1.5*10**11
dy = y[:-_] * 1.5*10**11
d = (dx**2 + dy**2)**(1/2)

# TODO : redimensionne et calcule le produit GM
dr = d[_:]
GM = dr**_ * _____

Indice

Deux redimensionnements différents interviennent ici, dans des sens opposés : x[:-2] retire les 2 derniers éléments, d[2:] retire les 2 premiers. Relis bien la consigne pour comprendre pourquoi ce ne sont pas les mêmes bornes qu'on retire à chaque fois.

Résultat attendu

Aucun affichage direct à cette étape, mais GM doit contenir 5 valeurs, toutes du même ordre de grandeur (autour de 1,3×10²⁰) : c'est ce nuage de valeurs presque constantes que la prochaine étape va exploiter.

En déduire la masse du Soleil

Consigne

Trace GM*10**(-20) (mis à une échelle plus pratique à manipuler) en fonction du numéro de point n = np.arange(len(t)-4), en points bleus avec la légende "points expérimentaux".
Modélise ce nuage de points par une constante, avec np.polyfit(n, GM, 0) (même technique qu'à la fiche n°3) : la valeur trouvée, mod[0], est une estimation de G × M_soleil.
Divise cette valeur par G = 6,67×10⁻¹¹ pour obtenir la masse du Soleil, et affiche le résultat avec "%.1e"%M pour ne garder qu'un chiffre après la virgule (notation scientifique).
Superpose la droite horizontale du modèle (GMmod = 0*n + mod[0]) au nuage de points, avec les étiquettes d'axes appropriées.

etape4_masse_soleil.py

n = np.arange(len(t) - _)
plt.plot(n, GM*10**(-20), '__', label="__________________")
plt.legend()
plt.grid()

# TODO : modélise GM par une constante
mod = np.polyfit(_, __, _)
M = float(mod[_] / (6.67*10**(-11)))
print('la masse du soleil est de l ordre de', "%.1e"%M, 'kg')

GMmod = 0*n + mod[_]
plt.plot(n, (GMmod)*10**(-20), '__', label="_______")
plt.legend()
plt.ylim(0, 5)
plt.xlabel("______________")
plt.ylabel("________________________")
plt.show()

Résultat attendu

La console affiche une masse de l'ordre de 2,0e+30 kg : du même ordre de grandeur que la valeur tabulée de la masse du Soleil (1,989×10³⁰ kg), ce qui valide la démarche — même en partant d'un modèle d'orbite simplifié.

Vérifier la loi en 1/d²

Consigne

Une autre façon de vérifier la deuxième loi de Newton : si norma = GM/d², alors norma tracée en fonction de 1/d² devrait donner une droite passant par l'origine, de pente GM.
Calcule drci = 1/(dr)**2, puis trace norma en fonction de drci, en points bleus.
Superpose le modèle théorique : une droite model = G × M_soleil × drci, en utilisant la valeur de masse trouvée à l'étape précédente (par exemple 1.7*10**30 si tu as obtenu cet ordre de grandeur).

etape5_verification.py

drci = 1 / (dr)**_
plt.plot(drci, norma, '__', label="__________________")
plt.legend()
plt.grid()

# TODO : modèle théorique, pente GM avec la masse trouvée à l'étape 4
model = 6.67*10**(-11) * ___*10**__ * _____
plt.plot(drci, model, '__', label="______")
plt.legend()
plt.xlabel("______")
plt.ylabel("____________")
plt.show()

Résultat attendu

Les points expérimentaux s'alignent raisonnablement sur la droite modèle : une confirmation visuelle supplémentaire que l'accélération de Mercure varie bien en 1/d², la signature caractéristique d'une force en 1/d² — la loi de gravitation universelle.

Pense-bête — unités et redimensionnements

1 UA	≈ 1,5×10¹¹ m
1 jour	= 24×3600 s
"%.1e"%M	affiche M en notation scientifique, 1 décimale

x[:-2]	retire les 2 derniers éléments
d[2:]	retire les 2 premiers éléments
GM = a·d²	constant si la force varie en 1/d²

Pas de Thonny sous la main ? Utilise Basthon

Ouvrir dans un nouvel onglet ↗

Basthon est une console Python qui tourne directement dans le navigateur, sans rien installer. Colle ici (avec Ctrl+V ou Cmd+V) le code copié avec les boutons « Copier le code » ci-dessus, puis appuie sur Entrée pour l'exécuter.