MODULE shared
  USE, INTRINSIC :: ISO_FORTRAN_ENV, dp=>REAL64 !modern DOUBLE PRECISION
  INTEGER :: j, k, n, atomic, atomic_read, count
  INTEGER, DIMENSION(:), ALLOCATABLE :: choice, amount
  REAL(dp) :: Z, A, wZA, wZA_tot, omega_p, rho, rho_read, sum, Cbar, X_0, X_1, aa, delta_0, eta, m, I, E, f
  CHARACTER(*), PARAMETER :: element_format =                                                                                      &
    '(i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2, i4, a1, a2)',    &
    my_format = '(a8, 2f8.3, f6.1, 2e11.4, f8.4, f7.4, f8.5, f7.4, f5.2, f4.1, e11.4)'
  CHARACTER(LEN=2), DIMENSION(100) :: element
  DATA element /"H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",     &
    "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb",    &
    "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm",   &
    "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi",   &
    "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm"/
  REAL(dp), DIMENSION(100) :: mass, potential, total
  CHARACTER(LEN=1) :: gas, gas_read
END MODULE shared

PROGRAM Sternheimer !Int J Appl Radiat Isot 33, 1189 (1982)
  USE shared
  IMPLICIT NONE
  mass = [1.00794, 4.002602, 6.941, 9.012182, 10.811, 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797, 22.98976928, 24.3050,        &
    26.9815386, 28.0855, 30.973762, 32.065, 35.453, 39.948, 39.0983, 40.078, 44.955912, 47.867, 50.9415, 51.9961, 54.938045,       &
    55.845, 58.933195, 58.6934, 63.546, 65.38, 69.723, 72.63, 74.92160, 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585, 91.224,   &
    92.90638, 95.96, 98.0, 101.07, 102.90550, 106.42, 107.8682, 112.411, 114.818, 118.710, 121.760, 127.60, 126.90447, 131.293,    &
    132.9054519, 137.327, 138.90547, 140.116, 140.90765, 144.242, 145.0, 150.36, 151.964, 157.25, 158.92535, 162.500, 164.93032,   &
    167.259, 168.93421, 173.054, 174.9668, 178.49, 180.94788, 183.84, 186.207, 190.23, 192.217, 195.084, 196.966569, 200.59,       &
    204.3833, 207.2, 208.98040, 210.0, 210.0, 222.0, 223.0, 226.0, 227.0, 232.03806, 231.03588, 238.02891, 237.0, 244.0, 243.0,    &
    247.0, 247.0, 251.0, 252.0, 257.0] !amu
  potential = [19.2, 41.8, 40.0, 63.7, 76.0, 78.0, 82.0, 95.0, 115.0, 137.0, 149.0, 156.0, 166.0, 173.0, 173.0, 180.0, 174.0,      &
    188.0, 190.0, 191.0, 216.0, 233.0, 245.0, 257.0, 272.0, 286.0, 297.0, 311.0, 322.0, 330.0, 334.0, 350.0, 347.0, 348.0, 343.0,  &
    352.0, 363.0, 366.0, 379.0, 393.0, 417.0, 424.0, 428.0, 441.0, 449.0, 470.0, 470.0, 469.0, 488.0, 488.0, 487.0, 485.0, 474.0,  &
    482.0, 488.0, 491.0, 501.0, 523.0, 535.0, 546.0, 560.0, 574.0, 580.0, 591.0, 614.0, 628.0, 650.0, 658.0, 674.0, 684.0, 694.0,  &
    705.0, 718.0, 727.0, 736.0, 746.0, 757.0, 790.0, 790.0, 800.0, 810.0, 823.0, 823.0, 830.0, 825.0, 794.0, 827.0, 826.0, 841.0,  &
    847.0, 878.0, 890.0, 902.0, 921.0, 934.0, 939.0, 952.0, 966.0, 980.0, 994.0] !eV
  total = [0.01359844, 0.07900519, 0.20348619, 0.39914912, 0.67098452, 1.03010566, 1.48606618, 2.04380319, 2.71586454, 3.51164078, &
    4.41997248, 5.45113072, 6.60498998, 7.88829137, 9.30580079, 10.85952461, 12.55480964, 14.39786659, 16.37945006, 18.50896698,   &
    20.78588723, 23.2208425, 25.8192349, 28.5853251, 31.51804801, 34.6045282, 37.915054, 41.38009864, 45.03842478, 8.0465816,      &
    0.1212235, 0.19727122, 0.2971326, 0.41010778, 0.56101381, 25.22016146, 0.94176213, 1.17041503, 1.4317741, 0.1574419,           &
    0.36202585, 18.44316403, 0.05208, 0.0525905, 0.0565989, 0.0606969, 0.0638962, 0.06338212, 0.10668616, 0.16549377, 0.25863891,  &
    0.3594296, 0.06258256, 0.06546259, 0.02705135, 0.0152156, 0.1473642, 0.2164947, 0.134157, 0.078765, 0.079882, 0.0815136,       &
    0.0845314, 0.0828701, 0.0790838, 0.0818789, 0.0831615, 0.0834777, 0.08461431, 0.08704026, 0.1523353, 0.07835507, 0.0075496,    &
    0.007864, 0.0078335, 0.0084382, 0.008967, 0.0275217, 0.0297255, 0.0633935, 0.0563662, 0.16550616, 0.2391356, 0.008417,         &
    0.009534973, 0.0107485, 0.0040727, 0.01542556, 0.01727, 0.0666067, 0.00589, 0.00619405, 0.0062657, 0.0060262, 0.0059738,       &
    0.0059915, 0.0061979, 0.0062817, 0.00642, 0.0065] !keV
  PRINT *, " "
  WRITE(6, '(a39)') "  Sternheimer Density Effect Parameters"
  PRINT *, " "
  WRITE(6, '(a64)', advance='no') "  Number of atomic species in the chemical (empirical) formula: "
  READ *, atomic_read
  atomic = atomic_read !number of chemical species that are not empty
  ALLOCATE(choice(atomic))
  ALLOCATE(amount(atomic))
  count = 1
  menu: DO
    IF (count > atomic) EXIT
    PRINT *, " "
    WRITE(6, '(a13, i1)') "  Selection #", count
    PRINT element_format, (j, "-", element(j), j = 1, 100) !print menu
    WRITE(6, '(a48)', advance='no') "  Select element and amount: [element] [amount] "
    READ *, choice(count), amount(count)
    count = count + 1
  END DO menu
  PRINT *, " "
  WRITE(6, '(a24)', advance='no') "  Mass density [g/cm3]: "
  READ *, rho_read
  rho = rho_read
  WRITE(6, '(a31)', advance='no') "  Is the material a gas? (y/n) "
  READ *, gas_read
  gas = gas_read
  Z = 0.0_dp
  A = 0.0_dp
  sum_Z_and_A: DO k = 1, atomic
    Z = Z + FLOAT(amount(k)) * FLOAT(choice(k))
    A = A + FLOAT(amount(k)) * mass(choice(k))
  END DO sum_Z_and_A
  wZA_tot = 0.0_dp
  sum = 0.0_dp
  Table_6: DO n = 1, atomic
    wZA = ((amount(n) * mass(choice(n))) / A) * (FLOAT(choice(n)) / mass(choice(n))) !weighted quantity of Z and A
    wZA_tot = wZA_tot + wZA !total weighted quantity of Z and A
    IF (gas == 'n') THEN !use adopted values from Table 6
      potential(1) = 19.2_dp !H
      potential(6) = 81.0_dp !C
      potential(7) = 82.0_dp !N
      potential(8) = 106.0_dp !O
      potential(9) = 112.0_dp !F
      potential(17) = 180.0_dp !Cl
    ELSE
      potential(1) = 19.2_dp !H
      potential(6) = 70.0_dp !C
      potential(7) = 82.0_dp !N
      potential(8) = 97.0_dp !O
    END IF
    IF (choice(n) == 1 .OR. choice(n) == 6 .OR. choice(n) == 7 .OR. choice(n) == 8 .OR. choice(n) == 9 .OR. choice(n) == 17) THEN !Table 6
      f = 1.00_dp
    ELSE
      f = 1.13_dp
    END IF
    sum = sum + wZA * LOG(f * potential(choice(n)))
  END DO Table_6
  I = EXP(sum / wZA_tot)
  omega_p = 28.816_dp * SQRT(rho * wZA_tot) !plasma frequency
  Cbar = 1.0_dp + 2.0_dp * LOG(I / omega_p)
  IF (gas == 'n') THEN !solids and liquids
    IF (I < 100.0_dp) THEN
      IF (Cbar < 3.681_dp) THEN
        X_0 = 0.2_dp
        X_1 = 2.0_dp
      ELSE IF (Cbar > 3.681_dp) THEN
        X_0 = 0.326_dp * Cbar - 1.0_dp
        X_1 = 2.0_dp
      END IF
    ELSE IF (I > 100.0_dp) THEN
      IF (Cbar < 5.215_dp) THEN
        X_0 = 0.2_dp
        X_1 = 3.0_dp
      ELSE IF (Cbar > 5.215_dp) THEN
        X_0 = 0.326_dp * Cbar - 1.5_dp
        X_1 = 3.0_dp
      END IF
    END IF
  ELSE !gases
    IF (Cbar < 10.0_dp) THEN
      X_0 = 1.6_dp
      X_1 = 4.0_dp
    ELSE IF (Cbar < 10.5_dp .AND. Cbar >= 10.0_dp) THEN
      X_0 = 1.7_dp
      X_1 = 4.0_dp
    ELSE IF (Cbar < 11.0_dp .AND. Cbar >= 10.5_dp) THEN
      X_0 = 1.8_dp
      X_1 = 4.0_dp
    ELSE IF (Cbar < 11.5_dp .AND. Cbar >= 11.0_dp) THEN
      X_0 = 1.9_dp
      X_1 = 4.0_dp
    ELSE IF (Cbar < 12.25_dp .AND. Cbar >= 11.5_dp) THEN
      X_0 = 2.0_dp
      X_1 = 4.0_dp
    ELSE IF (Cbar < 13.804_dp .AND. Cbar >= 12.25_dp) THEN
      X_0 = 2.0_dp
      X_1 = 5.0_dp
    ELSE IF (Cbar >= 13.804_dp) THEN
      X_0 = 0.326_dp * Cbar - 2.5_dp
      X_1 = 5.0_dp
    END IF
  END IF
  aa = (Cbar - 2.0_dp * LOG(10.0_dp) * X_0) / (X_1 - X_0)**3
  m = 3.0_dp
  delta_0 = 0.0_dp !for compounds, mixtures, and some elemental substances
  IF (gas == 'y') THEN
    eta = 1.0_dp !gas density ratio
  ELSE
    eta = 0.0_dp
  END IF
  IF (atomic == 1) THEN !elemental substances
    E = total(choice(atomic))
  ELSE !compounds and mixtures
    E = 0.0_dp
  END IF
  IF (atomic > 1) THEN !compounds and mixtures
    Z = Z / 3.0_dp !weird output
    A = A / 3.0_dp !weird output
  END IF
  PRINT *, " "
  !output precision here matches Target.dat START
  WRITE(6, '(a22, f7.3)') "  Total atomic charge:", Z
  WRITE(6, '(a20, f8.3)') "  Total atomic mass:", A
  WRITE(6, '(a28, f6.1, a3)') "  Mean ionization potential:", I, " eV"
  WRITE(6, '(a15, e11.4, a6)') "  Mass density:", rho, " g/cm3"
  WRITE(6, '(a19, e11.4, a3)') "  Plasma frequency:", omega_p, " eV"
  WRITE(6, '(a7, f8.4)') "  X_0 =", X_0
  WRITE(6, '(a7, f7.4)') "  X_1 =", X_1
  WRITE(6, '(a5, f8.5)') "  a =", aa
  WRITE(6, '(a5, f7.4)') "  m =", m
  WRITE(6, '(a11, f5.2)') "  delta_0 =", delta_0
  WRITE(6, '(a20, f4.1)') "  Gas density ratio:", eta
  WRITE(6, '(a34, e11.4, a4)') "  Total electronic binding energy:", E, " keV"
  !output precision here matches Target.dat END
  PRINT *, " "
  DEALLOCATE(choice)
  DEALLOCATE(amount)
END PROGRAM Sternheimer