Untitled
unknown
plain_text
6 months ago
2.3 kB
7
Indexable
import numpy as np import matplotlib.pyplot as plt # Constants m_e = 9.11 * 10**-31 # Electron mass in kg e = 1.60 * 10**-19 # Electron charge in C k = 8.99 * 10**9 # Coulomb's constant in N·m²/C² pi = 3.1416 # Value of π c = 299792458 # Speed of light in m/s # Alpha particle SP calculations for different mediums def alpha_B(T): m_p = 6.644657230 * 10**-27 # Alpha particle mass in kg c_mev = 1.60 * 10**-13 # MeV to kg conversion constant lorentz = ((T * c_mev) + (m_p * c**2)) / (m_p * c**2) B = np.sqrt(1 - (1 / lorentz**2)) B2 = np.round(B**2, 7) return B2 def alpha_SP_H2O(T): # Constants for H2O n_H2O = 4.96 * 10**28 # Number density for H2O in m^-3 I_H2O = 74.59 * e # Mean excitation potential for H2O in J z_a = 2 # Charge number for alpha particle B = alpha_B(T) if B < 1: sp_h2o = ((4 * pi * k**2 * z_a**2 * e**4 * n_H2O) / (m_e * c**2 * B)) * \ (np.log((2 * m_e * c**2 * B) / (I_H2O * (1 - B))) - B) else: sp_h2o = 0 # Avoid division by zero in extreme cases return sp_h2o def alpha_SP_H2(T): # Constants for H2 n_H2 = 3.34 * 10**29 # Number density for H2 in m^-3 I_H2 = 4.36 * e # Mean excitation potential for H2 in J z_a = 2 # Charge number for alpha particle B = alpha_B(T) if B < 1: sp_h2 = ((4 * pi * k**2 * z_a**2 * e**4 * n_H2) / (m_e * c**2 * B)) * \ (np.log((2 * m_e * c**2 * B) / (I_H2 * (1 - B))) - B) else: sp_h2 = 0 # Avoid division by zero in extreme cases return sp_h2 # Kinetic energy values (0.2 to 400 MeV in 100 steps) T_values = np.linspace(0.2, 400, 100) # Calculate stopping power for each T value SP_H2O_values = [alpha_SP_H2O(T) for T in T_values] SP_H2_values = [alpha_SP_H2(T) for T in T_values] # Plotting plt.figure(figsize=(12, 6)) plt.plot(T_values, np.log10(SP_H2O_values), label='alpha_SP_H2O', color='blue') plt.plot(T_values, np.log10(SP_H2_values), label='alpha_SP_H2', color='red') plt.xlabel('Kinetic Energy (MeV)') plt.ylabel('Log10(Stopping Power)') plt.title('Stopping Power vs. Kinetic Energy (Logarithmic Scale)') plt.legend() plt.grid(True) plt.show()
Editor is loading...
Leave a Comment