Untitled

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()