# Code to create example hyperbolic trajectories for Rutherford scattering.
%matplotlib notebook

import numpy as np
import math
import matplotlib.pyplot as plot
import matplotlib.patches as mpatches
from IPython.display import HTML
# r(phi) = 2(p/b)^2/(e cos(phi) - 1)
# Where p is the impact parameter
# b is the closest approach radius
# and e = Sqrt((2p/b)^2 + 1)

# Units of "b", so r/b, p/b.
# b = 34fm
# Rutherford's table uses 
# p/b = 10, 5, 2, 1, 0.5, 0.25, 0.125

def trajectory(p2b):
    """Array of r, theta describing hyperbolic trajectory fo p/b"""
    p2bSq = p2b*p2b
    e = np.sqrt(4* p2bSq + 1)
    # Rutherford says "the eccentricity is sec theta"
    # Reduce a bit to avoid very large r
    phi_in = np.arccos(1/e) - 0.001
    phi_out = -phi_in
    phis = np.linspace(phi_in, phi_out, 50)
    r2bs = [2 * p2bSq / (e * np.cos(phi) - 1) for phi in phis]
    # Rotate to make the alpha particle come in from left side.
    phi_lefts = [phi - phi_in + math.pi for phi in phis]
    return (r2bs, phi_lefts)

def nucleus(r):
    """Array of r, theta describing a circle of radius r"""
    phis = np.linspace(0, 2*math.pi, 50)
    rs = [r for phi in phis]
    return (rs, phis)

fig, ax = plot.subplots(subplot_kw={'projection': 'polar'})

r_max = 10
p2bs = [5, 2, 1, .5, .25, .125]
for p2b in p2bs:
    r, t = trajectory(p2b)
    display(f'p/b {p2b} angle {round(math.degrees(t[-1]))}')
    ax.plot(t, r)
    
# 7.3fm Modern Au radius / 34fm Rutherford r_min = b    
r, t = nucleus(7.3/34) 
ax.plot(t,r,color='gray', fillstyle='full')
        
ax.set_rmax(r_max)
ax.set_axis_off()

plot.show()
plot.savefig('RutherfordHyperbolasRaw.png')