Klein-Nishina formula

klein_nishina.jl

using Unitful
using Unitful: MeV, NoUnits, cm
using UnitfulRecipes

const h     = 6.626_070_040e-34*u"J*s"
const h_bar = h / (2pi)
const m_e   = 9.10938356e-31 * u"kg"
const c     = 299_792_458.0 * u"m/s"
const r_c   = h_bar / (c*m_e) # reduced compton wavelength of electron
const alpha = 1/137.04 # fine structure constant
function klein_nishina_cross_section(E,theta)
    P_E_theta = photon_energy_ratio(E, theta)
    ret = alpha^2 * r_c^2 * P_E_theta^2 * (P_E_theta + 1/P_E_theta - (sin(theta)^2)) / 2
    uconvert(cm^2, ret)
end

function photon_energy_ratio(E, theta)
    ret = 1 / (1+ (E/(m_e*c^2))*(1 - cos(theta)))
    uconvert(NoUnits, ret)
end

using Plots
thetas = range(0.0, stop=π, length=100)
plt = plot(xlabel="scatter angle/deg")

for E in [0.6MeV, 1MeV, 50MeV]
    xcses = [klein_nishina_cross_section(E, theta) for theta in thetas]
    ys = uconvert.(NoUnits, xcses .* (cm^-2)) * 10^24
    xs = rad2deg.(thetas)
    plot!(xs,ys, label=string(E))
end
display(plt)

plt = plot(xlabel="scatter angle/deg")
for E in [0.5MeV, 1MeV, 5MeV]
    ys = E * photon_energy_ratio.(E, thetas)
    xs = rad2deg.(thetas)
    plot!(xs,ys, label=string(E))
end
display(plt)