PROGRAM Alpha_dEdx
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  REAL(dp) E_i, second, a1, b, b2, capA, capB, cbar, delt, f1, f2, f4, f6, g, X, X0, X1, z0, z1, z23, dEdx !f5 not used
  REAL(dp), PARAMETER :: ALPHA = 7.297352569824D-03, 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
  CHARACTER(LEN=1) FD, FE
  PRINT *, " "
  WRITE(6, '(a42)', advance='no') "  Kinetic energy of alpha particle [MeV]: "
  READ *, E_i
  PRINT *, " "
  WRITE(6, '(a23)') "  Target material: SiO2"
  PRINT *, " "
  FD = 'Y' !density effect correction
  FE = 'Y' !electron capture
  second = E_i / 4.0_dp !MeV/n
  g = 1.0_dp + (second / 931.49406121_dp) !gamma from special relativity
  IF (FD == 'Y') THEN !density effect correction, DELT
    PRINT *, " FD"
    X0 = X_0
    X1 = X_1
    cbar = 2.0_dp * LOG(I / plasma_freq) + 1.0_dp !Seltzer, eqn. (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
  END IF !end FD
  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))
  IF (FE == 'Y') THEN !electron capture, Z1
    PRINT *, " FE"
    IF (INT(ATOMIC_NUMBER) == 4) THEN !beryllium target
      z1 = z0 * (1.0_dp - (2.045_dp + 2.000_dp * EXP(-0.04369_dp * z0)) * EXP(-7.000_dp * EXP(0.2643_dp * LOG(second))             &
        * (-0.4171_dp * LOG(z0))))
    ELSE IF (INT(ATOMIC_NUMBER) == 6) THEN !carbon target
      z1 = z0 * (1.0_dp - (2.584_dp + 1.910_dp * EXP(-0.03958_dp * z0)) * EXP(-6.933_dp * EXP(0.2433_dp * LOG(second))             &
        * EXP(-0.3969_dp * LOG(z0))))
    ELSE !everything else: Weaver, eqn. (59)
      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))))
    END IF
  ELSE !empirical formula by Anthony and Landford
    capA = 1.16D+00 - ATOMIC_NUMBER * (1.91D-03 - 1.26D-05 * ATOMIC_NUMBER)
    capB = (1.18D+00 - ATOMIC_NUMBER * (7.5D-03 - 4.53D-05 * ATOMIC_NUMBER)) / ALPHA
    z1 = z0 * (1.0_dp - capA * EXP(-capB * b / z23))
  END IF !end FE
  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 *, " "
  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