PROGRAM Beta !Beta_dEdx USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION IMPLICIT NONE REAL(dp) :: E_i, b, dEdx, dEdx_c, dEdx_b, Coulomb, projectile, absorber REAL(dp), PARAMETER :: PI = ATAN2(0.0_dp, -1.0_dp) !mathematical constants REAL(dp), PARAMETER :: MEC2 = 0.5109989461D+00, R_E = 2.8179403227D-13, N_A = 6.022140857D+23 !physical constants REAL(dp), PARAMETER :: A = 21.648886845D+00 !effective atomic mass of SiO2 in g/mol REAL(dp), PARAMETER :: Z = 10.80461D+00 !effective atomic number of SiO2 REAL(dp), PARAMETER :: I = 139.2D-06 !mean ionization potential of SiO2 in MeV PRINT *, " " WRITE(6, '(a41)', advance='no') " Kinetic energy of beta particle [MeV]: " READ *, E_i PRINT *, " " WRITE(6, '(a23)') " Target material: SiO2" PRINT *, " " b = SQRT(1.0_dp - (MEC2 / (E_i + MEC2))**2) !beta from special relativity Coulomb = MEC2 * R_E**2 !MeV-cm2 projectile = (2.0_dp * PI * Z) / b**2 !dimensionless absorber = N_A / A !leave out rho to get units of ENERGY/DEPTH ... g-1 dEdx_c = Coulomb * projectile * absorber * (LOG((E_i * (E_i + MEC2)**2 * b**2) / (2.0_dp * I**2 * MEC2)) + 1.0_dp - b**2 & - (2.0_dp * SQRT(1.0_dp - b**2) - 1.0_dp + b**2) * LOG(2.0_dp) + 0.125_dp * (1.0_dp - SQRT(1.0_dp - b**2))**2) !collision stopping power PRINT *, " dEdx_c =", dEdx_c, "MeV-cm2/g" Coulomb = R_E**2 !cm2 projectile = ((E_i + MEC2) * Z**2) / 137.0_dp !MeV dEdx_b = Coulomb * projectile * absorber * (4.0_dp * LOG((2.0_dp * (E_i + MEC2)) / MEC2) - 4.0_dp / 3.0_dp) !bremsstrahlung (radiative) stopping power PRINT *, " dEdx_b =", dEdx_b, "MeV-cm2/g" dEdx = dEdx_c + dEdx_b !total stopping power PRINT *, " dEdx = ", dEdx, "MeV-cm2/g" PRINT *, " " END PROGRAM Beta