import matplotlib.pyplot as plt
import numpy as np

# Potential energy from Rutherford's 1911 paper.
# V(r) = Ne(1/r - 3/2R + r^2/2R^3)
# Note:
# V(R) = Ne(1/R - 3/2R + 1/2R) = 0
# V(0) -> Ne (1/r) 

# Units of b = 2NeE/mv^2
b = 1
R = 10000*b
rs = np.linspace(R, b, 500)
nrs = np.linspace(-b, -R, 500)
rs = np.concatenate([rs,nrs])
Vs = [(1/np.absolute(r) - 3/(2*R) + r*r/(2*R*R*R)) for r in rs

fig, ax = plt.subplots()
ax.plot(rs, Vs)

t = ax.text(-R/2, 0.5, "Atomic radius",
            ha="center", va="center",  size=18,
            bbox=dict(boxstyle="darrow,pad=0.3",
                      fc="white", ec="steelblue", lw=2))
t2 = ax.text(R/2, 0.5, "Incoming kinetic energy",
            ha="center", va="center",  size=17,rotation=90,
            bbox=dict(boxstyle="darrow,pad=0.3",
                      fc="white", ec="steelblue", lw=2))

ax.set(xlabel='r/b', ylabel='$V(r)/E_{in}$',
       title='Highly concentrated Rutherford atom potential')
fig.savefig("RutherfordConcentrated.png")
plt.show()