# -*- coding: utf-8 -*-
# Chapitre 15 — Tester les deuxième et troisième lois de Kepler (CORRIGÉ)
# Mercure : aires balayées pendant des durées égales.
# Six planètes : régression affine de T^2 en fonction de a^3.

import math
import matplotlib.pyplot as plt
from scipy.stats import linregress

# --- deuxième loi : Mercure, une position tous les 5 jours ---
r = [0.3075, 0.315, 0.336, 0.363, 0.392, 0.418, 0.440, 0.455,
     0.464, 0.467, 0.462, 0.450, 0.432, 0.408, 0.381, 0.352,
     0.326, 0.310, 0.309]                       # u.a.
theta = [0, 31, 60, 85, 106, 124, 140, 155, 169, 183, 197, 211,
         227, 244, 263, 286, 312, 342, 13]       # degrés

aires = []
for i in range(len(r) - 1):
    dtheta = theta[i + 1] - theta[i]
    if dtheta < 0:
        dtheta = dtheta + 360
    S = 0.5 * r[i] * r[i + 1] * math.sin(math.radians(dtheta))
    aires.append(S)

somme_aires = 0.0
for S in aires:
    somme_aires = somme_aires + S
aire_moyenne = somme_aires / len(aires)
somme_ecarts = 0.0
for S in aires:
    somme_ecarts = somme_ecarts + (S - aire_moyenne) ** 2
ecart_type = math.sqrt(somme_ecarts / len(aires))

print("DEUXIEME LOI — aires balayees en 5 jours")
for i in range(len(aires)):
    print("S", i + 1, "=", aires[i], "ua2")
print("moyenne =", aire_moyenne, "ua2")
print("ecart-type relatif =", 100 * ecart_type / aire_moyenne, "%")

# --- troisième loi : six planètes ---
noms = ["Mercure", "Vénus", "Terre", "Mars", "Jupiter", "Saturne"]
a = [0.39, 0.72, 1.00, 1.52, 5.2, 9.52]      # u.a.
T = [87.9, 224.7, 365.25, 687.0, 4331.0, 10751.0]  # jours

G  = 6.674e-11
UA = 1.496e11    # m
JOUR = 86400.0   # s

# --- les deux listes de la régression ---
x = []
y = []
for i in range(len(a)):
    x.append(a[i] ** 3)   # a au cube  # À COMPLÉTER
    y.append(T[i] ** 2)   # T au carré  # À COMPLÉTER

# --- régression linéaire avec la méthode de la fiche M2 ---
reg = linregress(x, y)
pente = reg.slope          # en j2/ua3
ordonnee = reg.intercept   # en j2
n = len(x)
sx = 0.0
for valeur in x:
    sx = sx + valeur

# incertitude-type estimée de l'ordonnée à l'origine
x_moyen = sx / n
Sxx = 0.0
residus = []
for i in range(n):
    Sxx = Sxx + (x[i] - x_moyen) ** 2
    residus.append(y[i] - (pente * x[i] + ordonnee))
somme_residus = 0.0
for e_residu in residus:
    somme_residus = somme_residus + e_residu ** 2
ecart_residuel = math.sqrt(somme_residus / (n - 2))
u_ordonnee = ecart_residuel * math.sqrt(1 / n + x_moyen ** 2 / Sxx)

print("\nTROISIEME LOI — regression affine")
print("T2 =", pente, "a3 +", ordonnee)
print("b =", ordonnee, "+/-", u_ordonnee, "j2")
print("zero appartient a [b-u(b) ; b+u(b)] :",
      ordonnee - u_ordonnee <= 0 <= ordonnee + u_ordonnee)
for i in range(len(a)):
    print(noms[i], ": T2/a3 =", y[i] / x[i])

# la pente pèse le Soleil : k = 4 pi^2 / (G M), en unités SI
k_SI = pente * JOUR**2 / UA**3
M = 4 * math.pi**2 / (G * k_SI)
print("M_Soleil =", M, "kg   (valeur admise : 1,989e30 kg)")

plt.figure()
plt.plot(range(1, len(aires) + 1), aires, "o-", color="teal")
plt.axhline(aire_moyenne, color="orange", label="aire moyenne")
plt.xlabel("intervalle de 5 jours")
plt.ylabel("aire du triangle (u.a.2)")
plt.title("Deuxième loi de Kepler")
plt.grid(True)
plt.legend()
plt.show()

plt.figure()
plt.scatter(x, y, color="orange", zorder=3, label="planètes")
xs = [0, x[len(x) - 1] * 1.05]
ys = []
for valeur in xs:
    ys.append(pente * valeur + ordonnee)
plt.plot(xs, ys, color="teal", label="modèle T2 = k a3 + b")
plt.xlabel("a3 (u.a.3)")
plt.ylabel("T2 (jours2)")
plt.title("Troisième loi de Kepler")
plt.grid(True)
plt.legend()
plt.show()