Kepler Problem¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
Nel problema di Kepler si considera una particella di massa $m$ soggetta a un potenziale centrale gravitazionale $$\large V(r) = -\frac{k}{r}, $$ e si introduce l’Hamiltoniana del sistema in coordinate polari piane $$\large H(r,\theta,p_r,p_\theta) = \frac{p_r^2}{2m}\frac{p_\theta^2}{2mr^2}\frac{k}{r}. $$
Poiché l’angolo $\theta$ è una coordinata ciclica, il momento coniugato $$\large p_\theta = L $$ è una costante del moto e coincide con il momento angolare.
Il problema si riduce quindi a un moto unidimensionale nella variabile radiale, governato dall’Hamiltoniana efficace $$\large H_{\text{eff}}(r,p_r) = \frac{p_r^2}{2m} + V_{\text{eff}}(r), \qquad V_{\text{eff}}(r) = \frac{L^2}{2mr^2} - \frac{k}{r}. $$
Le equazioni di Hamilton per il moto radiale sono $$\large \dot r = \frac{\partial H_{\text{eff}}}{\partial p_r} = \frac{p_r}{m}, \qquad \dot p_r = -\frac{\partial H_{\text{eff}}}{\partial r}. $$
L’evoluzione temporale si ottiene integrando queste equazioni. In particolare, fissata l’energia totale $E$, il moto radiale è determinato dalla relazione $$\large \frac{1}{2} m \dot r^2 + V_{\text{eff}}(r) = E, $$ che consente di ricavare il tempo come funzione del raggio tramite $$\large t - t_0 = \int \frac{dr}{\sqrt{\frac{2}{m}\left(E - V_{\text{eff}}(r)\right)}}. $$
Una volta nota $r(t)$, l’evoluzione angolare segue da $$\large \dot\theta = \frac{\partial H}{\partial p_\theta} = \frac{L}{mr^2}, \qquad \theta(t) = \theta_0 + \int \frac{L}{m r(t)^2}\, dt. $$
In questo modo l’evoluzione temporale completa del sistema è ottenuta risolvendo separatamente il moto radiale e quello angolare, sfruttando la conservazione dell’energia e del momento angolare caratteristica del problema di Kepler.
import numpy as np
def reduced_mass(m1, m2):
"""
Calculate the reduced mass of two interacting bodies.
:param ml: Mass of the first body.
:param m2: Mass of the second body.
:return:
Reduced mass.
"""
return (m1 * m2) / (m1 + m2)
def lagrangian(mu, r, theta_dot, phi_dot):
"""
Compute the Lagrangian of the system.
:param mu: Reduced mass.
:param r: Radial distance.
:param theta_dot: Rate of change of theta.
:param phi_dot: Rate of change of phi.
:return:
Lagrangian value.
"""
T = 0.5 * mu * (r**2 * theta_dot**2 + r**2 * np.sin(theta)**2 * phi_dot**2)
U = -G * m1 * m2 / r
return T + U
def hamiltonian(p_r, p_theta, p_phi, mu, r, theta):
"""
Compute the Hamiltonian of the system.
:param p_r: Radial momentum.
:param p_theta: Theta momentum.
:param p_phi: Phi momentum.
:param mu: Reduced mass.
:param r: Radial distance.
:param theta: Theta angle.
:return:
Hamiltonian value.
"""
T = (p_r**2 / (2 * mu)) + (p_theta**2 / (2 * mu * r**2)) + (p_phi**2 / (2 * mu * r**2 * np. sin(theta)**2))
U = -G * m1 * m2 / r
return T + U
def kepler_orbit_simulation(r0, r_dot0, phi0, L, dt, T, mu, G, m1, m2):
t = 0.0
r = r0
r_dot = r_dot0
phi = phi0
r_list = [r]
phi_list = [phi]
while t < T:
# forza radiale efficace
F_r = (-G * m1 * m2 / r**2) + (L**2 / (mu * r**3))
# integrazione (Euler, coerente con il tuo codice)
r_dot += F_r * dt / mu
r += r_dot * dt
phi_dot = L / (mu * r**2)
phi += phi_dot * dt
r_list.append(r)
phi_list.append(phi)
t += dt
return np.array(r_list), np.array(phi_list)
# Parameters for a two-body system
m1 = 5.972e24 # Mass of the Earth (kg)
m2 = 1.989e30 # Mass of the Sun (kg)
G = 6.67430e-11 # Gravitational constant (m''3 kg'-l s'-2)
mu = reduced_mass(m1, m2)
r0 = 1.47e11 # Approximate distance at perihelion (m)
r_dot0 = 0.0 # Initial radial velocity (m/s)
L = 2.66e40 # Assumed angular momentum (kg m^/s)
dt = 86400 # Time step of one day (s)
T = 365.25 * dt # One year in seconds
phi0 = 0.0
r, phi = kepler_orbit_simulation(
r0, r_dot0, phi0, L, dt, T, mu, G, m1, m2
)
x = r * np.cos(phi)
y = r * np.sin(phi)
import matplotlib.pyplot as plt
plt.figure(figsize=(4, 4))
plt.plot(x, y, label="Orbita")
plt.scatter(0, 0, color="orange", s=100, label="Sole")
plt.xlabel("x [m]")
plt.ylabel("y [m]")
plt.axis("equal")
plt.grid()
plt.legend(loc="center left", bbox_to_anchor=(0.6, 0.5))
plt.tight_layout()
plt.title("Orbita kepleriana Terra–Sole")
plt.show()
# Simulating the orbit
#positions, velocities = kepler_orbit_simulation(r0, r_dot0, L, dt, T, mu, G, m1, m2)
# Print results for demonstration
#for i in range(len(positions)):
# print(f"Day {i}: Position (r) = {positions[i]} m. Velocity (r_dot) = {velocities [i]} m/s")
