Scattering in Central Potentials¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
La funzione simula il fenomeno di scattering di una particella in un potenziale centrale $V(r)$ utilizzando una formulazione hamiltoniana e un integratore simplettico.
Si considera una particella di massa (m) con momento angolare conservato $$\large L = p_\theta, $$ e si lavora in coordinate polari $(r,\theta)$, riducendo il problema a un moto radiale accoppiato all’evoluzione angolare.
L’Hamiltoniana del sistema è $$\large H(r,p_r) = \frac{p_r^2}{2m}\frac{L^2}{2mr^2}V(r), $$
da cui derivano le equazioni di Hamilton $$\large \dot r = \frac{p_r}{m}, \qquad \dot p_r = \frac{L^2}{m r^3} - \frac{dV}{dr}, \qquad \dot\theta = \frac{L}{m r^2}. $$
La derivata del potenziale (\frac{dV}{dr}) viene calcolata numericamente mediante differenze finite centrate.
L’integrazione temporale è effettuata con uno schema simplettico di tipo ``kick--drift--kick'' (leapfrog):
- aggiornamento di mezzo passo del momento radiale $p_r$;
- aggiornamento di un passo completo della posizione radiale $r$ e dell’angolo $\theta$;
- aggiornamento finale di mezzo passo del momento radiale.
Questo metodo preserva la struttura hamiltoniana del sistema e garantisce una buona conservazione dell’energia nel lungo tempo, caratteristica essenziale nello studio dello scattering.
La funzione restituisce le traiettorie temporali $$\large r(t), \quad p_r(t), \quad \theta(t), $$ che permettono di ricostruire la traiettoria spaziale della particella e di analizzare quantitativamente il processo di diffusione, come l’angolo di scattering in funzione del parametro d’impatto o del momento angolare $L$.
import numpy as np
import matplotlib.pyplot as plt
def central_potential_scattering(mass, V_func, initial_conditions, L, delta_t, T):
"""
Function to simulate scattering in a central potential using
symplectic integration.
:param mass: Mass of the particle
:param V_func: Function defining the potential V(r)
:param initial_conditions: Initial radial position and momentum (r0, p_r0)
:param L: Angular momentum
:param delta_t: Time step for integration
:param T: Total time of simulation
:return:
r(t), p_r(t), theta(t)
"""
r0, p_r0, theta0 = initial_conditions
r, p_r, theta = r0, p_r0, theta0
trajectory_r = [r]
trajectory_p_r = [p_r]
trajectory_theta = [theta]
def dV_dr(r):
epsilon = 1e-5
return (V_func(r + epsilon) - V_func(r - epsilon)) / (2 * epsilon)
time = 0.0
while time < T:
# half step momento radiale
p_r += (delta_t / 2) * ((L**2) / (mass * r**3) - dV_dr(r))
# step posizione
r += delta_t * (p_r / mass)
theta += delta_t * (L / (mass * r**2))
# half step momento radiale
p_r += (delta_t / 2) * ((L**2) / (mass * r**3) - dV_dr(r))
trajectory_r.append(r)
trajectory_p_r.append(p_r)
trajectory_theta.append(theta)
time += delta_t
return (np.array(trajectory_r),
np.array(trajectory_p_r),
np.array(trajectory_theta))
def V_inverse_square(r, k=1.0):
"""
Inverse-square law potential function.
:param r: Radial distance
:param k: Coupling constant
:return:
Potential value
"""
return -k / r
# =======================
# Parametri simulazione
# =======================
mass = 1.0
L = 1.0
delta_t = 0.01
T = 20.0
# condizioni iniziali (r, p_r, theta)
initial_conditions = (5.0, -0.5, 0.0)
# integrazione
r, p_r, theta = central_potential_scattering(
mass, V_inverse_square, initial_conditions, L, delta_t, T
)
# =======================
# Traiettoria cartesiana
# =======================
x = r * np.cos(theta)
y = r * np.sin(theta)
# =======================
# Grafico
# =======================
plt.figure(figsize=(4, 4))
plt.plot(x, y, label="Traiettoria")
plt.scatter([0], [0], color="red", s=50, label="Centro")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Scattering in un campo centrale (1/r)")
plt.axis("equal")
plt.grid(True)
plt.legend()
plt.show()
