bethebloch-phy5601hw2q4

import numpy as np
import matplotlib.pyplot as plt
#everything in SI

n=8.49*10**28 #electron density of copper
e=1.60217663*10**-19
c=299792458
me=9.1093837*10**-31
epsilon=8.85418782*10**-12
K1=(n*e**4)/(me*c**2*4*np.pi*epsilon**2)

I=322*e #https://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html #the 'e' is to convert to SI
K2=2*me*c**2/I

def dEdx(E, z, m): #Sigmund, Peter. Particle Penetration and Radiation Effects.
    gamma=E/(m*c**2)
    beta=np.sqrt(1-gamma**-2)
    return K1*z**2*beta**-2*(np.log(K2*beta**2*gamma**2)-beta**2)

m_mu=1.883531627*10**-28
z_mu=1
m_p=1.67262192*10**-27
z_p=1
m_a=6.6446573357*10**-27
z_a=2

betagamma=3
E_mu=m_mu*c**2*np.sqrt(1+betagamma**2)
E_p=m_p*c**2*np.sqrt(1+betagamma**2)
E_a=m_a*c**2*np.sqrt(1+betagamma**2)

L=0.1 #10cm in m
ndx=1000
#super basic numerical integration haha
def energy_left(E, z, m, L, ndx): #L=length of copper block, ndx=number of steps
    energy_lost=0
    dx=L/ndx
    for i in range(ndx):
        dE=dx*dEdx(E-energy_lost,z,m)
        energy_lost+=dE
    return E-energy_lost

energy_left_mu=energy_left(E_mu, z_mu, m_mu, L, ndx)
energy_left_p=energy_left(E_p, z_p, m_p, L, ndx)
energy_left_a=energy_left(E_a, z_a, m_a, L, ndx)
print('Muon started with ', round(10**-6*E_mu/e, 2), ' MeV. Muon energy left = ', round(10**-6*energy_left_mu/e, 2), ' MeV')
print('Proton started with ', round(10**-6*E_p/e, 2), ' MeV. Proton energy left = ', round(10**-6*energy_left_p/e, 2), ' MeV')
print('Alpha started with ', round(10**-6*E_a/e, 2), ' MeV. Alpha energy left = ', round(10**-6*energy_left_a/e, 2), ' MeV')
print('')

L=1
ndx=9000
def rangelength(E, z, m, L, ndx):
    stoppower=[]
    i=0
    energy_lost=0
    dx=L/ndx
    threshold=m*c**2 #'at rest' condition
    while E-energy_lost>=threshold:
        dE=dx*dEdx(E-energy_lost,z,m)
        stoppower.append(dE/dx)
        energy_lost+=dE
        i+=1
    return i*dx, 10**-8*np.array(stoppower)/e

range_mu, sp_mu=rangelength(E_mu, z_mu, m_mu, L, ndx)
print('Muon range = ', round(100*range_mu, 2), ' cm.')
range_p, sp_p=rangelength(E_p, z_p, m_p, L, ndx)
print('Proton range = ', round(100*range_p, 2), ' cm.')
range_a, sp_a=rangelength(E_a, z_a, m_a, L, ndx)
print('Alpha range = ', round(100*range_a, 2), ' cm.')

#trying to plot bragg peaks but they look weird so never mind :(
'''plt.plot(100*np.array(range(len(sp_mu)))/ndx, sp_mu, label='mu bragg')
plt.plot(100*np.array(range(len(sp_p)))/ndx, sp_p, label='p bragg')
plt.plot(100*np.array(range(len(sp_a)))/ndx, sp_a, label='a bragg')
plt.xlabel('path length (cm)')
plt.ylabel('stopping power (MeV/cm)')
plt.legend()
plt.show()'''