PROGRAM Alpha_dEdx USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION IMPLICIT NONE REAL(dp) :: E_i, second, a1, b, b2, cbar, delt, f1, f2, f4, f6, g, X, X0, X1, z0, z1, z23, dEdx !f5 not used REAL(dp), PARAMETER :: EMASS = 0.51099891013D+06 !physical constants REAL(dp), PARAMETER :: plasma_freq = 3.1014D+01, X_0 = 0.1385D+00, X_1 = 3.0025D+00, a = 0.08408D+00, m = 3.5064D+00, & delta_0 = 0.00D+00 !Sternheimer density effect parameters REAL(dp), PARAMETER :: DENSITY = 2.648D+00 !mass density of SiO2 in g/cm3 (α-quartz) REAL(dp), PARAMETER :: ATOMIC_MASS = 21.648886845D+00 !effective atomic mass of SiO2 in g/mol REAL(dp), PARAMETER :: ATOMIC_NUMBER = 10.80461D+00 !effective atomic number of SiO2 REAL(dp), PARAMETER :: I = 139.2D+00 !mean ionization potential of SiO2 in eV PRINT *, " " WRITE(6, '(a42)', advance='no') " Kinetic energy of alpha particle [MeV]: " READ *, E_i PRINT *, " " WRITE(6, '(a23)') " Target material: SiO2" PRINT *, " " second = E_i / 4.0_dp !MeV/n g = 1.0_dp + (second / 931.49406121_dp) !gamma from special relativity X0 = X_0 X1 = X_1 cbar = 2.0_dp * LOG(I / plasma_freq) + 1.0_dp !Seltzer, eq. (14) b = SQRT(1.0_dp - 1.0_dp / (g * g)) !beta from special relativity X = LOG10(b * g) IF (X < X0) THEN delt = delta_0 * EXP(4.6052_dp * (X - X0)) ELSE IF (X >= X0 .AND. X < X1) THEN delt = 4.6052_dp * X - cbar + EXP(LOG(a) + m * LOG(X1 - X)) ELSE delt = 4.6052_dp * X - cbar END IF b2 = 1.0_dp - 1.0_dp / (g * g) !beta squared from special relativity b = SQRT(b2) z0 = 2.0_dp !bare projectile charge a1 = 4.001506466_dp !projectile atomic mass z23 = EXP((2.0_dp / 3.0_dp) * LOG(z0)) z1 = z0 * (1.0_dp - ((1.164_dp + 0.2319_dp * EXP(-0.004302_dp * ATOMIC_NUMBER)) + 1.658_dp * EXP(-0.05170_dp * z0)) & * EXP(-(8.114_dp + 0.9876_dp * LOG(ATOMIC_NUMBER)) * EXP((0.3140_dp + 0.01072_dp * LOG(ATOMIC_NUMBER)) * LOG(second)) & * EXP(-(0.5218_dp + 0.02521_dp * LOG(ATOMIC_NUMBER)) * LOG(z0)))) !everything else: Weaver, eq. (59) f1 = (0.3070722_dp * z1 * z1 * ATOMIC_NUMBER) / (b2 * a1 * ATOMIC_MASS) !constant out front f2 = LOG((2.0_dp * EMASS * b2) / I) !first term in stopping logarithm f4 = 1.0_dp f6 = 2.0_dp * LOG(g) - b2 dEdx = f1 * (f2 * f4 + f6 - (delt / 2.0_dp)) !implement corrections PRINT *, " LET:", dEdx * 4.0_dp, "MeV-cm2/g" !from (MeV/n)(cm2/g) to MeV-cm2/g dEdx = (dEdx * 4.0_dp * DENSITY) / 10.0_dp !from (MeV/n)(cm2/g) to keV/um PRINT *, " ", dEdx, "keV/um" PRINT *, " " END PROGRAM Alpha_dEdx