PROGRAM Gamma_profile
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  REAL(dp) :: E_i, depth, PE, E_over_MEC2, Compton, b, PP, sigma, mu_m
  REAL(dp), PARAMETER :: PI = ATAN2(0.0_dp, -1.0_dp) !mathematical constants
  REAL(dp), PARAMETER :: FINE = 0.0072973525664D+00, MEC2 = 0.5109989461D+00, R_E = 2.8179403227D-13, N_A = 6.022140857D+23,       &
    EXPONENT = 3.5D+00 !physical constants
  REAL(dp), PARAMETER :: DENSITY = 2.648D+00 !mass density of SiO2 in g/cm3 (α-quartz)
  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
  PRINT *, " "
  WRITE(6, '(a25)', advance='no') "  Energy of gamma [MeV]: "
  READ *, E_i
  PRINT *, " "
  WRITE(6, '(a23)') "  Target material: SiO2"
  PRINT *, " "
  IF (E_i > 0.1_dp .AND. E_i < 0.35_dp) THEN !photon energy in MeV above the K-absorption edge
    PE = 4.0_dp * FINE**4 * SQRT(2.0_dp) * Z**5 * (8.0_dp * PI * R_E**2 / 3.0_dp) * (MEC2 / E_i)**EXPONENT !cross section for photoelectric effect (PE)
    PRINT *, " Photoelectric effect:", PE, "cm2"
  ELSE
    PE = 0.0_dp
  END IF
  E_over_MEC2 = E_i / MEC2 !dimensionless
  Compton = (PI * R_E**2 / E_over_MEC2) * ((1.0_dp - ((2.0_dp * E_over_MEC2 + 2.0_dp) / E_over_MEC2**2)) * LOG(2.0_dp              &
    * E_over_MEC2 + 1.0_dp) + 0.5_dp + (4.0_dp / E_over_MEC2) - 1.0_dp / (2.0_dp * (2.0_dp * E_over_MEC2 + 1.0_dp)**2)) !cross section for Compton scattering
  PRINT *, "   Compton scattering:", Compton, "cm2"
  IF (E_i > 2.0_dp * MEC2) THEN !photon energy in MeV below which pair production cannot occur
    b = SQRT(1.0_dp - (MEC2 / (E_i + MEC2))**2) !beta from special relativity
    PP = (PI * R_E**2 / 2.0_dp) * (1.0_dp - b**2) * ((3.0_dp - b**4) * LOG((1.0_dp + b)/(1.0_dp - b)) - 2.0_dp * b                 &
      * (2.0_dp - b**2)) !Gould and Schréder (1967), cross section for pair production (PP)
    PRINT *, "      Pair production:", PP, "cm2"
  ELSE
    PP = 0.0_dp
  END IF
  sigma = PE + Z * Compton + PP !cm2
  PRINT *, "  Total cross section:", sigma, "cm2"
  mu_m = N_A * sigma / A !cm2/g
  depth = 1.0_dp / mu_m !probability of non-absorption has dropped to 1/e
  PRINT *, " "
  PRINT *, " Attenuation length [cm]:", depth / DENSITY
  PRINT *, " "
END PROGRAM Gamma_profile