MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER, PARAMETER :: MAXE = 200, Target_dat = 10 !manage UNIT numbers as a resource
  REAL(dp), PARAMETER :: PI = ATAN2(0.0_dp, -1.0_dp), ALPHA = 7.297352569824D-03, EMASS = 0.51099891013D+06,                       &
    TAU = 938.27204621_dp / 931.49406121_dp, COMPTON = 3.0557335D-03
  CHARACTER(*), PARAMETER :: my_format = '(a8, 2f8.3, f6.1, 2e11.4, f8.4, f7.4, f8.5, f7.4, f5.2, f4.1, e11.4)'
  TYPE Projectile
    REAL(dp) z1, a1
  END TYPE Projectile
  TYPE Absorber
    REAL(dp), DIMENSION(:), ALLOCATABLE :: z2, a2, I_pot, rho, pla, X0, X1, a, m, d0, eta, bind
    CHARACTER(LEN=9), DIMENSION(:), ALLOCATABLE :: tname
  END TYPE Absorber
END MODULE Shared

PROGRAM Dreamweaver
  USE Shared, ONLY : dp, MAXE, Target_dat, my_format, Projectile, Absorber
  IMPLICIT NONE
  INTEGER i, j, k, abs, count, line, MAXAB
  REAL(dp) second, second_read, kinetic, third, third_read, depth, rel, ei, energy1, energy2, dE, dx, resultat, en, e1, Benton,    &
    dEdx, dR
  TYPE(Projectile) projektil
  TYPE(Absorber) ziel
  REAL(dp), DIMENSION(MAXE) :: energy_table, range_table, range
  CHARACTER(LEN=1) task, task_read, answer, blank !answer_read not used ... -Wmaybe-uninitialized
  CHARACTER(LEN=1) FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR
  EXTERNAL dEdx, Fbrems, Lindhard, Carg, Chyperg, Clngamma, Benton
  PRINT *, " "
  WRITE(6, '(a51)') "  Dreamweaver: The Berkeley Range-Energy Calculator"
  PRINT *, " "
  WRITE(6, '(a18)') "  Tasks: r = range"
  WRITE(6, '(a19)') "         e = energy"
  WRITE(6, '(a16)') "         d = LET"
  PRINT *, " "
  WRITE(6, '(a25)', advance='no') "  Enter the task letter: "
  READ *, task_read
  task = task_read
  OPEN(Target_dat, file='Target.dat', access='sequential', status='old', action='read')
  line = 0
  file_size: DO !portable method of determining file size
    IF (.false.) EXIT
    READ(Target_dat, *, end=999) blank
    line = line + 1
  END DO file_size
  999 CONTINUE
  MAXAB = line
  REWIND(Target_dat)
  ALLOCATE(ziel%tname(MAXAB), ziel%z2(MAXAB), ziel%a2(MAXAB), ziel%I_pot(MAXAB), ziel%rho(MAXAB), ziel%pla(MAXAB), ziel%X0(MAXAB), &
    ziel%X1(MAXAB), ziel%a(MAXAB), ziel%m(MAXAB), ziel%d0(MAXAB), ziel%eta(MAXAB), ziel%bind(MAXAB))
  material: DO k = 1, MAXAB
    READ(Target_dat, my_format) ziel%tname(k), ziel%z2(k), ziel%a2(k), ziel%I_pot(k), ziel%rho(k), ziel%pla(k), ziel%X0(k),        &
    ziel%X1(k), ziel%a(k), ziel%m(k), ziel%d0(k), ziel%eta(k), ziel%bind(k)
  END DO material
  CLOSE(Target_dat)
  FD = 'Y' !density effect correction
  FE = 'Y' !electron capture
  FSH = 'Y' !shell correction
  FLE = 'N' !Leung relativistic shell correction
  FKIN = 'N' !kinematic correction
  FRAD = 'N' !radiative correction
  FLS = 'Y' !LS replaces BMA
  FNS = 'Y' !finite nuclear size correction accounts for the finite size of atomic nuclei
  FBA = 'Y' !Barkas correction
  FBR = 'N' !projectile bremsstrahlung correction
  IF (task == 'r') THEN
    PRINT *, " "
    WRITE(6, '(a12)') "  RANGE TASK"
  ELSE IF (task == 'e') THEN
    PRINT *, " "
    WRITE(6, '(a13)') "  ENERGY TASK"
  ELSE IF (task == 'd') THEN
    PRINT *, " "
    WRITE(6, '(a29)') "  LINEAR ENERGY TRANSFER TASK"
  END IF
  IF (task == 'r') THEN !RANGE task
    PRINT *, " "
    WRITE(6, '(a39)', advance='no') "  Initial kinetic energy, E_i [MeV/n]: "
    READ *, second_read
    second = second_read
    kinetic = second
    CALL Stuff(MAXAB, projektil, ziel, abs)
  ELSE IF (task == 'e') THEN !ENERGY task
    PRINT *, " "
    WRITE(6, '(a32)', advance='no') "  Initial kinetic energy? (y/n) "
    READ *, answer
    IF (answer == 'y') THEN
      CALL Stuff(MAXAB, projektil, ziel, abs)
      PRINT *, " "
      WRITE(6, '(a34)', advance='no') "  Range in the target, R [g/cm2]: "
      READ *, third_read
      third = third_read
      depth = third
    ELSE
      PRINT *, " "
      WRITE(6, '(a39)', advance='no') "  Initial kinetic energy, E_i [MeV/n]: "
      READ *, second_read
      second = second_read
      kinetic = second
      CALL Stuff(MAXAB, projektil, ziel, abs)
      PRINT *, " "
      WRITE(6, '(a34)', advance='no') "  Depth of the target, x [g/cm2]: "
      READ *, third_read
      third = third_read
      depth = third
    END IF
  ELSE IF (task == 'd') THEN !LET task
    PRINT *, " "
    WRITE(6, '(a39)', advance='no') "  Initial kinetic energy, E_i [MeV/n]: "
    READ *, second_read
    second = second_read
    kinetic = second
    third = 0.0_dp !REL is activated by setting this to a positive number
    CALL Stuff(MAXAB, projektil, ziel, abs)
  END IF
  PRINT *, " "
  WRITE(6, '(a18, a9)') "  Target material:", ziel%tname(abs)
  rel = 0.0_dp !REL is activated by setting this to a positive number
  !passing target parameters were previously here ... now called "ziel"
  energy_table = 0.0_dp !zero-out array
  make_energy_table: DO j = 1, MAXE
    energy_table(j) = EXP(LOG(10.0_dp) * (FLOAT(j - 1) * 6.0_dp / FLOAT(MAXE - 1)))
  END DO make_energy_table
  range_table = 0.0_dp !zero-out array
  IF (task == 'r' .OR. task == 'e' .AND. answer == 'y') THEN
    i = 1
    make_range_table: DO !8 MeV/n upper limit
      IF (energy_table(i) > 8.0_dp) EXIT
      e1 = energy_table(i) !passing variable, energy_table(1) = 1.0 MeV/n
      range_table(i) = Benton(e1, projektil, ziel, abs) !1 is projectile, 2 is target
      i = i + 1
    END DO make_range_table
    range_dR: DO
      IF (i > MAXE) EXIT
      energy2 = energy_table(i) !passing variable
      energy1 = energy_table(i - 1) !passing variable
      range_table(i) = range_table(i - 1)                                                                                &
        + dR(energy2, energy1, rel, projektil, ziel, abs, FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
      i = i + 1
    END DO range_dR
  END IF
  range = 0.0_dp !zero-out array
  IF (task == 'r') THEN !RANGE task
    IF (kinetic > energy_table(1)) THEN !greater than 1.0 MeV/n, use table
      i = 2
      find_energy: DO
        IF (kinetic < energy_table(i)) EXIT
        i = i + 1
      END DO find_energy
      resultat = range_table(i - 1) + (kinetic - energy_table(i - 1)) * (range_table(i) - range_table(i - 1)) / (energy_table(i)   &
        - energy_table(i - 1)) !interpolated range
    ELSE !less than ~1 MeV/n, NIM A 544, 502 (2005)
      resultat = kinetic * range_table(1) / energy_table(1) !energy_table(1) = 1.0 MeV/n
    END IF
    IF (resultat < 0.0_dp) THEN !interpolation failed, now try to integrate
      en = 1.0_dp !1 MeV/n
      ei = 8.0_dp !8 MeV/n
      IF (kinetic > ei) THEN !greater than 8 MeV/n, integrate
        range(1) = Benton(ei, projektil, ziel, abs) !get range for 8 MeV/n
        i = 2
        differential_range: DO !integrate up to input energy
          energy1 = en !save
          en = EXP(LOG(10.0_dp) * (LOG10(ei) + FLOAT(i - 1) * (LOG10(kinetic) - LOG10(ei)) / FLOAT(MAXE - 1)))
          IF (en >= kinetic) EXIT
          energy2 = en !pass
          range(i) = range(i - 1) + dR(energy2, energy1, rel, projektil, ziel, abs, FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA,FBR)
          i = i + 1
        END DO differential_range
        resultat = range(i - 1)
      ELSE IF (kinetic > 1.0_dp .AND. kinetic <= ei) THEN !less than 8 MeV/n, parse table anyway
        resultat = Benton(kinetic, projektil, ziel, abs)
      ELSE
        resultat = kinetic * Benton(en, projektil, ziel, abs) / en
      END IF
    END IF
    CALL Corrections(FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
    PRINT *, " Range:", resultat, "g/cm2"
    resultat = resultat / ziel%rho(abs)
    PRINT *, "       ", resultat, "cm"
    PRINT *, " "
    DEALLOCATE(ziel%tname, ziel%z2, ziel%a2, ziel%I_pot, ziel%rho, ziel%pla, ziel%X0, ziel%X1, ziel%a, ziel%m, ziel%d0, ziel%eta,  &
      ziel%bind)
    STOP
  ELSE IF (task == 'e') THEN !ENERGY task
    IF (answer == 'y') THEN !get initial kinetic energy
      PRINT *, " "
      PRINT *, " INITIAL KINETIC ENERGY"
      depth = third
      IF (depth > range_table(1)) THEN !check if range is within table
        i = 2
        find_range: DO
          IF (depth < range_table(i)) EXIT
          i = i + 1
        END DO find_range
        resultat = energy_table(i - 1) + (depth - range_table(i - 1)) * (energy_table(i) - energy_table(i - 1)) / (range_table(i)  &
          - range_table(i - 1)) !interpolated energy
      ELSE !less than 1 MeV/n, NIM A 544, 502 (2005)
        resultat = depth * energy_table(1) / range_table(1) !energy_table(1) = 1.0 MeV/n
      END IF
      CALL Corrections(FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
      PRINT *, " Energy:", resultat, "MeV/n"
      PRINT *, " "
      DEALLOCATE(ziel%tname, ziel%z2, ziel%a2, ziel%I_pot, ziel%rho, ziel%pla, ziel%X0, ziel%X1, ziel%a, ziel%m, ziel%d0, ziel%eta,&
        ziel%bind)
      STOP
    ELSE !get final kinetic energy
      PRINT *, " "
      PRINT *, " FINAL KINETIC ENERGY"
      en = kinetic
      depth = third
      dx = depth / 20.0_dp
      KE_f: DO count = 1, 20
        dE = dEdx(en, rel, projektil, ziel, abs, FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
        IF (en < 1.0_dp .OR. dE < 0.0_dp) EXIT !less than 1 MeV/n
        en = en - (dx * dE)
      END DO KE_f
      IF (en < 1.0_dp .OR. dE < 0.0_dp) en = 0.0_dp !less than 1 MeV/n
      resultat = en
      CALL Corrections(FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
      PRINT *, " Energy:", resultat, "MeV/n"
      PRINT *, " "
      DEALLOCATE(ziel%tname, ziel%z2, ziel%a2, ziel%I_pot, ziel%rho, ziel%pla, ziel%X0, ziel%X1, ziel%a, ziel%m, ziel%d0, ziel%eta,&
        ziel%bind)
      STOP
    END IF
  END IF
  IF (task == 'd') THEN !LET task
    resultat = dEdx(kinetic, third, projektil, ziel, abs, FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
    CALL Corrections(FD, FE, FSH, FLE, FKIN, FRAD, FLS, FNS, FBA, FBR)
    PRINT *, " LET:", resultat * projektil%a1, "MeV-cm2/g" !from (MeV/n)(cm2/g) to MeV-cm2/g
    resultat = (resultat * projektil%a1 * ziel%rho(abs)) / 10.0_dp !from (MeV/n)(cm2/g) to keV/um
    PRINT *, "     ", resultat, "keV/um"
    PRINT *, " "
    DEALLOCATE(ziel%tname, ziel%z2, ziel%a2, ziel%I_pot, ziel%rho, ziel%pla, ziel%X0, ziel%X1, ziel%a, ziel%m, ziel%d0, ziel%eta,  &
      ziel%bind)
    STOP
  END IF
END PROGRAM Dreamweaver

SUBROUTINE Stuff(MAXABs, projektil, ziel, abss)
  USE Shared, ONLY : dp, Projectile, Absorber
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: MAXABs
  TYPE(Absorber), INTENT(IN) :: ziel !didn't do abs, abs_read because derived-type
  INTEGER, INTENT(OUT) :: abss
  TYPE(Projectile), INTENT(OUT) :: projektil !didn't do z1, z1_read, a1, a1_read because derived-type
  INTEGER i, beam, beam_read
  PRINT *, " "
  WRITE(6, '(a26)') "  Beams: 0. User specified"
  WRITE(6, '(a17)') "         1. 56-Fe"
  WRITE(6, '(a17)') "         2. 28-Si"
  WRITE(6, '(a17)') "         3. 16-O "
  WRITE(6, '(a17)') "         4. 12-C "
  WRITE(6, '(a17)') "         5.  4-He"
  WRITE(6, '(a17)') "         6.  1-H "
  PRINT *, " "
  WRITE(6, '(a23)', advance='no') "  Enter the beam type: "
  READ *, beam_read
  beam = beam_read
  IF (beam == 0) THEN
    PRINT *, " "
    WRITE(6, '(a6)', advance='no') "  Z = "
    READ *, projektil%z1
    WRITE(6, '(a6)', advance='no') "  A = "
    READ *, projektil%a1
  ELSE IF (beam == 1) THEN !Fe
    projektil%z1 = 26.0_dp
    projektil%a1 = 55.845_dp
  ELSE IF (beam == 2) THEN !Si
    projektil%z1 = 14.0_dp
    projektil%a1 = 28.0855_dp
  ELSE IF (beam == 3) THEN !O
    projektil%z1 = 8.0_dp
    projektil%a1 = 15.9994_dp
  ELSE IF (beam == 4) THEN !C
    projektil%z1 = 6.0_dp
    projektil%a1 = 12.0107_dp
  ELSE IF (beam == 5) THEN !He
    projektil%z1 = 2.0_dp
    projektil%a1 = 4.001506466_dp
  ELSE IF (beam == 6) THEN !H
    projektil%z1 = 1.0_dp
    projektil%a1 = 1.00794_dp
  END IF
  PRINT *, " "
  WRITE(6, '(a11, "1.", a9)') "  Targets: ", ziel%tname(1)
  available: DO i = 2, MAXABs
    PRINT '(i12, ".", a9)', i, ziel%tname(i)
  END DO available
  PRINT *, " "
  WRITE(6, '(a25)', advance='no') "  Enter the target type: "
  READ *, abss
END SUBROUTINE Stuff

SUBROUTINE Corrections(FDs, FEs, FSHs, FLEs, FKINs, FRADs, FLSs, FNSs, FBAs, FBRs)
  IMPLICIT NONE
  CHARACTER(LEN=1), INTENT(IN) :: FDs, FEs, FSHs, FLEs, FKINs, FRADs, FLSs, FNSs, FBAs, FBRs
  PRINT *, " "
  WRITE(6, '(a15)', advance='no') "  Corrections: "
  IF (FDs == 'Y') WRITE(6, '(a3)', advance='no') "FD "
  IF (FEs == 'Y') WRITE(6, '(a3)', advance='no') "FE "
  IF (FSHs == 'Y') WRITE(6, '(a4)', advance='no') "FSH "
  IF (FLEs == 'Y') WRITE(6, '(a4)', advance='no') "FLE "
  IF (FKINs == 'Y') WRITE(6, '(a5)', advance='no') "FKIN "
  IF (FRADs == 'Y') WRITE(6, '(a5)', advance='no') "FRAD "
  IF (FLSs == 'Y') WRITE(6, '(a4)', advance='no') "FLS "
  IF (FNSs == 'Y') WRITE(6, '(a4)', advance='no') "FNS "
  IF (FBAs == 'Y') WRITE(6, '(a4)', advance='no') "FBA "
  IF (FBRs == 'Y') WRITE(6, '(a3)', advance='no') "FBR"
  PRINT *, " " !not advancing
  PRINT *, " "
END SUBROUTINE Corrections

REAL(dp) FUNCTION dEdx(secondf, thirdf, projektil, ziel, absf, FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf)
  USE Shared, ONLY : dp, PI, ALPHA, EMASS, Projectile, Absorber !EMASS in eV/c2
  IMPLICIT NONE
  INTEGER absf
  REAL(dp) secondf, thirdf
  TYPE(Projectile) projektil
  TYPE(Absorber) ziel
  CHARACTER(LEN=1) FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf
  REAL(dp) a1, b, b2, Bbr, Blim, cadj, capA, capB, cbar, delt, etam2, f1, f2, f3, f4, f6, f7, f8, f9, fv, g, REL, S, Sbr, v, X, X0,&
    X1, z0, z1, z23, Fbrems, Lindhard !f5 not used
  EXTERNAL Fbrems, Lindhard
  g = 1.0_dp + (secondf / 931.49406121_dp) !gamma from special relativity
  FDf = 'Y' !-Wmaybe-uninitialized
  IF (FDf == 'Y') THEN !density effect correction, DELT
    X0 = ziel%X0(absf)
    X1 = ziel%X1(absf)
    cbar = 2.0_dp * LOG(ziel%I_pot(absf) / ziel%pla(absf)) + 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 (ziel%eta(absf) > 0.0_dp) THEN
      cbar = cbar - 2.303_dp * LOG10(ziel%eta(absf))
      X0 = X0 - 0.5_dp * LOG10(ziel%eta(absf))
      X1 = X1 - 0.5_dp * LOG10(ziel%eta(absf))
    END IF
    IF (X < X0) THEN
      delt = ziel%d0(absf) * EXP(4.6052_dp * (X - X0))
    ELSE IF (X >= X0 .AND. X < X1) THEN
      delt = 4.6052_dp * X - cbar + EXP(LOG(ziel%a(absf)) + ziel%m(absf) * 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 = projektil%z1 !bare projectile charge
  a1 = projektil%a1 !projectile atomic mass
  z23 = EXP((2.0_dp / 3.0_dp) * LOG(z0))
  IF (FEf == 'Y') THEN !electron capture, Z1
    IF (INT(ziel%z2(absf)) == 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(secondf))            &
        * (-0.4171_dp * LOG(z0))))
    ELSE IF (INT(ziel%z2(absf)) == 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(secondf))            &
        * EXP(-0.3969_dp * LOG(z0))))
    ELSE !Weaver, eqn. (59)
      z1 = z0 * (1.0_dp - ((1.164_dp + 0.2319_dp * EXP(-0.004302_dp * ziel%z2(absf))) + 1.658_dp * EXP(-0.05170_dp * z0))          &
        * EXP(-(8.114_dp + 0.9876_dp * LOG(ziel%z2(absf))) * EXP((0.3140_dp + 0.01072_dp * LOG(ziel%z2(absf))) * LOG(secondf))     &
        * EXP(-(0.5218_dp + 0.02521_dp * LOG(ziel%z2(absf))) * LOG(z0))))
    END IF
  ELSE !empirical formula by Anthony and Landford
    capA = 1.16D+00 - ziel%z2(absf) * (1.91D-03 - 1.26D-05 * ziel%z2(absf))
    capB = (1.18D+00 - ziel%z2(absf) * (7.5D-03 - 4.53D-05 * ziel%z2(absf))) / ALPHA
    z1 = z0 * (1.0_dp - capA * EXP(-capB * b / z23))
  END IF !end FE
  f1 = (0.3070722_dp * z1 * z1 * ziel%z2(absf)) / (b2 * a1 * ziel%a2(absf)) !constant out front
  f2 = LOG((2.0_dp * EMASS * b2) / ziel%I_pot(absf)) !first term in stopping logarithm
  IF (FSHf == 'Y') THEN !shell correction, F2
    etam2 = 1.0_dp / (b * b * g * g)
    cadj = 1.0D-06 * ziel%I_pot(absf) * ziel%I_pot(absf) * etam2 * (0.422377D+00 + etam2 * (0.0304043D+00 - etam2                  &
      * 0.00038106D+00)) + 1.0D-09 * ziel%I_pot(absf) * ziel%I_pot(absf) * ziel%I_pot(absf) * etam2 * (3.858019D+00 + etam2        &
      * (-0.1667989D+00 + etam2 * 0.00157955D+00)) !Weaver, eqn. (64)
    f2 = f2 - (cadj / ziel%z2(absf)) !Weaver, eqn. (62)
    IF (FLEf == 'Y') THEN !Leung relativistic shell correction, F2
      f2 = f2 - (5.0D+00 / 3.0D+00) * LOG((2.0D+00 * EMASS * b2) / ziel%I_pot(absf)) * ((1.0D+03 * ziel%bind(absf))                &
        / (ziel%z2(absf) * EMASS)) - ((ziel%I_pot(absf) * ziel%I_pot(absf)) / (4.0D+00 * EMASS * EMASS * b2)) !Weaver, eqn. (65)
    END IF !end FLE
  END IF !end FSH
  f3 = 0.0_dp
  f4 = 1.0_dp
  f6 = 2.0_dp * LOG(g) - b2
  f8 = 0.0_dp
  f9 = 0.0_dp
  Sbr = 0.0_dp
  IF (FRADf == 'Y') THEN !radiative correction
    f9 = (ALPHA / PI) * b2 * (6.0822_dp + LOG(2.0_dp * g) * (LOG(2.0_dp * g) * (2.4167_dp + 0.3333_dp * LOG(2.0_dp * g))           &
      - 8.0314_dp)) !Weaver, eqn. (46)
  END IF !end FRAD, not able to account for "8.0314"
  IF (FKINf == 'Y') THEN !kinematic correction
    f8 = -0.5D+00 * (LOG(1.0D+00 + 2.0D+00 * (g / a1) * 5.4858D-04) + (g / a1) * 5.4858D-04 * b2 / (g * g)) !Weaver, eqn. (48)
  END IF !end FKIN
  IF (FBRf == 'Y') THEN !projectile bremsstrahlung correction
    Blim = 4.0_dp * LOG(183.0_dp / EXP(LOG(ziel%z2(absf)) / 3.0_dp)) + (2.0_dp / 9.0_dp)
    Bbr = (((12.0_dp * g * g + 4.0_dp) / (3.0_dp * b * g * g)) - ((6.0_dp * b + 8.0_dp) / (3.0_dp * b * b * g * g)) * LOG(g + b    &
      * g)) * LOG(g + b * g) - (4.0_dp / 3.0_dp) + (2.0_dp / (b * g * g)) * Fbrems(2.0_dp * b * g * g * (1.0_dp + b))              &
      - 4.828_dp * ALPHA * ALPHA * ziel%z2(absf) * ziel%z2(absf) !Weaver, eqns. (50) and (54)
    IF (Blim < Bbr) Bbr = Blim
    Sbr = 9.7822719D-08 * ((z1 * z1 * z1 * z1 * ziel%z2(absf) * ziel%z2(absf)) / (a1 * a1 * ziel%a2(absf))) * g * Bbr !Weaver, eqn. (49)
  END IF !end FBR
  IF (FLSf == 'Y') f3 = Lindhard(projektil, b, FNSf) !LS replaces BMA
  IF (FBAf == 'Y') THEN !Barkas correction
    v = (b * g) / (ALPHA * SQRT(ziel%z2(absf))) !Weaver, eqn. (61)
    IF (v > 1.0_dp) THEN !Weaver, Table 2
      fv = 0.0019_dp * EXP(-2.0_dp * LOG(v / 10.0_dp)) !replaced fva(9) with 0.0019
      f4 = 1.0_dp + (2.0_dp * z1 * fv) / SQRT(ziel%z2(absf)) !Weaver, eqn. (60)
    ELSE !factor of 2
      f4 = 1.0_dp !no correction
    END IF
  END IF !end FBA
  S = f1 * (f2 * f4 + f3 + f6 - (delt / 2.0_dp) + f8 + f9) + Sbr !implement corrections
  IF (thirdf > 0.0_dp) THEN !compute Restricted Energy Loss (REL) if switched on
    f7 = LOG((2.0_dp * EMASS * b2 * g * g) / thirdf) + b2 * (thirdf / (2.0_dp * EMASS * b2 * g * g) - 1.0_dp) !implement REL plus corrections
    REL = f1 * (f2 * f4 + f3 + f6 - (delt / 2.0_dp) - 0.5_dp * f7 + f8)
    dEdx = REL !Restricted Energy Loss
  ELSE
    dEdx = S
  END IF
END FUNCTION dEdx

REAL(dp) FUNCTION dR(energy2f, energy1f, relf, projektil, ziel, absf, FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf)
  USE Shared, ONLY : dp, Projectile, Absorber
  IMPLICIT NONE
  INTEGER absf
  REAL(dp) energy2f, energy1f, relf
  TYPE(Projectile) projektil
  TYPE(Absorber) ziel
  CHARACTER(LEN=1) FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf
  REAL(dp) e1, e2, e3, e4, de2, dedx1, dedx2, dedx3, dedx4, dEdx
  EXTERNAL dEdx
  de2 = (energy2f - energy1f) / 2.0_dp
  e1 = energy1f + 1.33998104_dp * de2
  dedx1 = dEdx(e1, relf, projektil, ziel, absf, FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf)
  e2 = energy1f + 1.86113631_dp * de2
  dedx2 = dEdx(e2, relf, projektil, ziel, absf, FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf)
  e3 = energy1f + 0.13886369_dp * de2
  dedx3 = dEdx(e3, relf, projektil, ziel, absf, FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf)
  e4 = energy1f + 0.66001869_dp * de2
  dedx4 = dEdx(e4, relf, projektil, ziel, absf, FDf, FEf, FSHf, FLEf, FKINf, FRADf, FLSf, FNSf, FBAf, FBRf)
  dR = de2 * (0.65214515_dp / dedx1 + 0.34785485_dp / dedx2 + 0.34785485_dp / dedx3 + 0.65214515_dp / dedx4)
END FUNCTION dR

REAL(dp) FUNCTION Lindhard(projektil, bf, FNSff) !Lindhard-Sorenson correction replaces Bloch-Mott-Ahlen
  USE Shared, ONLY : dp, PI, ALPHA, COMPTON, Projectile
  IMPLICIT NONE
  TYPE(Projectile) projektil
  REAL(dp) bf
  CHARACTER(LEN=1) FNSff
  INTEGER i, n, max
  REAL(dp) prh, gz, figi, k, fsgs, L, sk, dkm1, grgs, H, term1, term2, term3, sd2, sdm2, eta, rho, dmk, nn, k0, total, b0, a0, a1, &
    bn, anp1, an, anm1, bnm1, asum, bsum, signk, frgr, a3, gg, sdd, z1, Carg
  REAL(dp), DIMENSION(3) :: dk
  COMPLEX(dp) Cexir, Cexis, Cedr, Ceds, Cske, Cmske, Clamr, Clams, Caar, Caas, Cbbr, Cbbs, Czzr, Chyperg, Clngamma
  EXTERNAL Chyperg, Clngamma
  dk = 0.0_dp !zero-out array
  a3 = EXP(LOG(projektil%a1) / 3.0_dp)
  eta = projektil%z1 * ALPHA / bf
  gg = 1.0_dp / SQRT(1.0_dp - bf * bf) !gamma from special relativity
  rho = a3 * COMPTON
  prh = bf * gg * rho
  total = 0.0_dp !the objective of this function is to get this sum
  term1 = 0.0_dp !terms in Weaver, eqn. (34)
  term2 = 0.0_dp
  term3 = 1.0_dp
  IF (gg < 10.0_dp / rho .OR. FNSff == 'N') THEN
    dkm1 = 0.0_dp !JMD
    dmk = 0.0_dp !JMD
    n = 1
    main: DO
      IF (n > 100) EXIT
      k0 = FLOAT(n)
      max = 2
      IF (n == 1) max = 3
      maximum: DO i = 1, max
        IF (i == 1) THEN
          k = k0
        ELSE IF (i == 2) THEN
          k = -k0 - 1.0_dp
        ELSE IF (i == 3) THEN
          k = -k0
        END IF
        signk = k / ABS(k)
        sk = SQRT(k * k - ALPHA * ALPHA * projektil%z1 * projektil%z1) !Weaver, eqn. (30), top
        IF (k > 0.0_dp) THEN !Weaver, eqn. (30), bottom
          L = k
        ELSE
          L = -k - 1.0_dp
        END IF
        Cske = CMPLX(sk + 1.0_dp, eta, KIND=dp)
        Cexir = SQRT(CMPLX(k, -eta / gg, KIND=dp) / CMPLX(sk, -eta, KIND=dp)) !Weaver, eqn. (30), middle
        Cedr = Cexir * EXP(CMPLX(0.0_dp, -AIMAG(Clngamma(Cske)) + (PI / 2.0_dp) * (L - sk), KIND=dp)) !Weaver, eqn. (35)
        IF (FNSff == 'Y') THEN !finite nuclear size correction accounts for the finite size of atomic nuclei, TOTAL
          Cmske = CMPLX(-sk + 1.0_dp, eta, KIND=dp)
          Cexis = SQRT(CMPLX(k, -eta / gg, KIND=dp) / CMPLX(-sk, -eta, KIND=dp))
          Ceds = Cexis * EXP(CMPLX(0.0_dp, -AIMAG(Clngamma(Cmske)) + (PI / 2.0_dp) * (L + sk), KIND=dp)) !Weaver, eqn. (35)
          Caar = Cske
          Caas = Cmske !minus version of CSKE
          Cbbr = CMPLX(2.0_dp * sk + 1.0_dp, 0.0_dp, KIND=dp)
          Cbbs = CMPLX(-2.0_dp * sk + 1.0_dp, 0.0_dp, KIND=dp) !minus version of CBBR
          Czzr = CMPLX(0.0_dp, 2.0_dp * prh, KIND=dp)
          Clamr = Cexir * EXP(CMPLX(0.0_dp, -prh, KIND=dp)) * Chyperg(Caar, Cbbr, Czzr) !Weaver, eqn. (38)
          Clams = Cexis * EXP(CMPLX(0.0_dp, -prh, KIND=dp)) * Chyperg(Caas, Cbbs, Czzr)
          grgs = AIMAG(Clamr) / AIMAG(Clams)
          grgs = grgs * EXP(REAL(Clngamma(Cske)) - REAL(Clngamma(Cmske)) - REAL(Clngamma(Cbbr)) + REAL(Clngamma(Cbbs))             &
            + 2.0_dp * sk * LOG(2.0_dp * prh)) !Weaver, eqn. (39), underflow located with -ffpe-trap=underflow
          IF (COS(AIMAG(Clngamma(Cbbs))) < 1.0_dp) grgs = -1.0_dp * grgs
          IF (ABS(grgs) > 1.0D-09) THEN
            frgr = SQRT((gg - 1.0_dp) / (gg + 1.0_dp)) * REAL(Clamr) / AIMAG(Clamr) !Weaver, eqn. (37)
            fsgs = SQRT((gg - 1.0_dp) / (gg + 1.0_dp)) * REAL(Clams) / AIMAG(Clams)
            gz = -signk * (rho * gg + 1.5_dp * ALPHA * projektil%z1)
            z1 = -signk * projektil%z1
            b0 = 1.0_dp !Weaver, eqn. (41), 1
            a0 = ((2.0_dp * ABS(k) + 1.0_dp) * b0) / (rho - gz) !Weaver, eqn. (41), 2
            a1 = 0.5_dp * (gz + rho) * b0 !Weaver, eqn. (41), 3
            an = a1
            anm1 = a0
            bnm1 = b0
            asum = a0
            bsum = b0
            nn = 1.0_dp
            summation: DO
              IF (ABS(anm1 / asum) < 1.0D-06 .AND. ABS(bnm1 / bsum) < 1.0D-06) EXIT
              bn = ((rho - gz) * an + ALPHA * z1 * anm1 / 2.0_dp) / (2.0_dp * ABS(k) + 2.0_dp * nn + 1.0_dp) !Weaver, eqn. (41), 4
              anp1 = ((gz + rho) * bn - ALPHA * z1 * bnm1 / 2.0_dp) / (2.0_dp * (nn + 1.0_dp)) !Weaver, eqn. (41), 5
              asum = asum + an
              bsum = bsum + bn
              anm1 = an
              an = anp1
              bnm1 = bn
              nn = nn + 1.0_dp
            END DO summation
            IF (k > 0.0_dp) THEN
              figi = asum / bsum !Weaver, eqn. (40)
            ELSE
              figi = bsum / asum
            END IF
            H = ((frgr - figi) / (figi - fsgs)) * grgs !Weaver, eqns. (35) and (36)
          ELSE
            H = 0.0_dp !H provides the connection between the interior uniform spherically symmetric potential and the exterior Coulomb potential
          END IF
        ELSE
          H = 0.0_dp
          Ceds = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
        END IF !end FNS, back to LS
        dk(i) = Carg(Cedr + H * Ceds) !Weaver, eqn. (35)
      END DO maximum
      IF (n > 1) dk(3) = dmk
      sdm2 = SIN(dk(3) - dk(2))
      term1 = (k0 * (k0 + 1.0_dp) * sdm2 * sdm2) / (eta * eta * (2.0_dp * k0 + 1.0_dp)) !Weaver, eqn. (34), middle
      IF (n > 1) THEN
        sd2 = SIN(dk(1) - dkm1)
        term1 = term1 + (k0 * (k0 - 1.0_dp) * sd2 * sd2) / (eta * eta * (2.0_dp * k0 - 1.0_dp)) !Weaver, eqn. (34), top
      END IF
      sdd = SIN(dk(1) - dk(3))
      term2 = (k0 * sdd * sdd) / (eta * eta * (4.0_dp * k0 * k0 - 1.0_dp)) !Weaver, eqn. (34), bottom
      term3 = term1 - (1.0_dp / k0) !Weaver, eqn. (34), bottom
      total = total + term2 + term3
      n = n + 1
      dkm1 = dk(1)
      dmk = dk(2)
    END DO main
  ELSE
    total = -LOG(prh) - 0.2_dp !Weaver, eqn. (43), ultrarelativistic limit, asymptotic value of the LS correction
  END IF
  Lindhard = total + 0.5_dp * bf * bf !Weaver, eqn. (34)
END FUNCTION Lindhard

REAL(dp) FUNCTION Fbrems(x) !compute a mathematical function related to bremsstrahlung
  USE Shared, ONLY : dp, PI
  IMPLICIT NONE
  REAL(dp) x
  INTEGER n
  REAL(dp) s, t
  n = 1
  s = 0.0_dp
  t = 1.0_dp
  IF (NINT(x) == 1) THEN
    Fbrems = (PI * PI) / 12.0_dp !Weaver, eqn. (53)
  ELSE IF (x < 1.0_dp) THEN
    product: DO
      IF (ABS(t) < 1.0D-04 .OR. n > 10) EXIT
      t = t * (-x) !product accounts for (-x)^n
      s = s + t / FLOAT(n * n) !Weaver, eqn. (52)
      n = n + 1
    END DO product
    Fbrems = -s
  ELSE
    inverse: DO
      IF (ABS(t) < 1.0D-04 .OR. n > 10) EXIT
      t = t * (-1.0_dp / x)
      s = s + t / FLOAT(n * n)
      n = n + 1
    END DO inverse
    Fbrems = (PI * PI) / 12.0_dp + 0.5_dp * LOG(x) * LOG(x) + s !Weaver, eqn. (53)
  END IF
END FUNCTION Fbrems

REAL(dp) FUNCTION Carg(z) !arctangent intrinsic function with a complex argument
  USE Shared, ONLY : dp, PI
  IMPLICIT NONE
  COMPLEX(dp) z
  IF (ABS(AIMAG(z)) < 1.0D-03) THEN
    IF (REAL(z) > 0.0_dp) THEN
      Carg = 0.0_dp
    ELSE
      Carg = PI
    END IF
  ELSE IF (ABS(REAL(z)) < 1.0D-03) THEN
    IF (ABS(AIMAG(z)) < 1.0D-03) THEN
      Carg = 0.0_dp
    ELSE IF (AIMAG(z) > 0.0_dp) THEN
      Carg = PI / 2.0_dp
    ELSE
      Carg = -PI / 2.0_dp
    END IF
  ELSE
    Carg = ATAN2(AIMAG(z), REAL(z))
  END IF
END FUNCTION Carg

COMPLEX(dp) FUNCTION Clngamma(z) !complex gamma function (Appendix A)
  USE Shared, ONLY : dp, PI
  IMPLICIT NONE
  COMPLEX(dp) z
  INTEGER j
  REAL(dp) x, y, r, lterm1, lterm2, lterm3, aterm1, aterm2, aterm3, num, denom, cterm
  REAL(dp), DIMENSION(6) :: coeff
  COMPLEX(dp) result
  coeff = [76.18009172947146D+00, -86.50532032941677D+00, 24.01409824083091D+00, -1.231739572450155D+00, 0.1208650973866179D-02,   &
    -0.5395239384953D-05] !Weaver, eqn. (A.2)
  IF (REAL(z) > 0.0_dp) THEN
    x = REAL(z) - 1.0_dp
    y = AIMAG(z)
  ELSE
    x = -REAL(z)
    y = -AIMAG(z)
  END IF
  r = SQRT((x + 5.5_dp) * (x + 5.5_dp) + y * y) !rule of logarithms and powers
  lterm1 = (x + 0.5_dp) * LOG(r) !Weaver, eqn. (A.5), first line
  lterm2 = -y * ATAN2(y, x + 5.5_dp) - (x + 5.5_dp) + 0.5_dp * LOG(2.0_dp * PI) !Weaver, eqn. (A.5), second and third lines
  aterm1 = y * LOG(r) !Weaver, eqn. (A.6), first line
  aterm2 = (x + 0.5_dp) * ATAN2(y, x + 5.5_dp) - y !Weaver, eqn. (A.6), second line
  num = 0.0_dp !T
  denom = 1.000000000190015_dp !B, c_0
  parameters: DO j = 1, 6 !special choice of N = 6
    cterm = coeff(j) / ((x + FLOAT(j)) * (x + FLOAT(j)) + y * y) !Weaver, eqn. (A.3)
    num = num + cterm
    denom = denom + (x + FLOAT(j)) * cterm !Weaver, eqn. (A.4)
  END DO parameters
  num = -y * num !Weaver, eqn. (A.3)
  lterm3 = 0.5_dp * LOG(num * num + denom * denom) !Weaver, eqn. (A.5), third line
  aterm3 = ATAN2(num, denom) !Weaver, eqn. (A.6), third line
  result = CMPLX(lterm1 + lterm2 + lterm3, aterm1 + aterm2 + aterm3, KIND=dp) !rule of logarithms and addition
  IF (REAL(z) < 0.0_dp) THEN !reflection formula, Weaver, eqn. (A.7)
    result = CMPLX(LOG(PI), 0.0_dp, KIND=dp) - (result + LOG(SIN(PI * z)))
  END IF
  Clngamma = result !return the natural log
END FUNCTION Clngamma

COMPLEX(dp) FUNCTION Chyperg(a, b, z) !confluent hypergeometric function (Appendix B)
  USE Shared, ONLY : dp
  IMPLICIT NONE
  COMPLEX(dp) a, b, z
  INTEGER m
  REAL(dp) dm
  COMPLEX(dp) Cm, term, sum_term, previous_term
  sum_term = CMPLX(1.0_dp, 0.0_dp, KIND=dp) !Weaver, eqn. (B.3)
  term = CMPLX(1.0_dp, 0.0_dp, KIND=dp) !Weaver, eqn. (B.3)
  previous_term = term
  m = 1
  series: DO
    IF (ABS(term) < 1.0D-06 .AND. ABS(previous_term) < 1.0D-06) EXIT
    dm = FLOAT(m)
    Cm = CMPLX(dm - 1.0_dp, 0.0_dp, KIND=dp)
    term = previous_term * ((a + Cm) / (b + Cm) * (1.0_dp / dm) * z)
    sum_term = sum_term + term
    previous_term = term
    m = m + 1
  END DO series
  Chyperg = sum_term
END FUNCTION Chyperg

REAL(dp) FUNCTION Benton(e1f, projektil, ziel, absf)
  USE Shared, ONLY : dp, TAU, Projectile, Absorber !TAU is m_p/amu
  IMPLICIT NONE
  INTEGER absf
  REAL(dp) e1f
  TYPE(Projectile) projektil
  TYPE(Absorber) ziel
  INTEGER q, m, n
  REAL(dp) loglambda, term, b, g, x, bzz, logt, logi
  REAL(dp), DIMENSION(4) :: join
  REAL(dp), DIMENSION(2, 7), SAVE :: cjoin = RESHAPE([ &
       0.94D+00,  20.19D+00, -84.08D+00,  132.98D+00, -30.77D+00, -102.29D+00, 64.03D+00, &
      12.62D+00, -51.96D+00, 199.14D+00, -367.09D+00, 327.06D+00, -108.57D+00,  0.00D+00  & !eqns. (13) and (14), Dreamweaver_notes
  ], SHAPE(cjoin), order=[2, 1])
  REAL(dp), DIMENSION(3) :: prnglo
  REAL(dp), DIMENSION(3, 4, 4), SAVE :: amn = RESHAPE([ &
    -8.72500D+00,  1.88000D+00,  0.741900D+00,  0.752000D+00, &
     0.83090D+00,  0.11140D+00, -0.528800D+00, -0.555890D+00, &
    -0.13396D+00, -0.06481D+00,  0.126423D+00,  0.128431D+00, &
     0.01262D+00,  0.00540D+00, -0.009341D+00, -0.009306D+00, & !Table 1 of Dreamweaver_notes
    -7.6604D-01,  2.5398D+00, -2.4598D-01,  0.0000D+00, &
     7.3736D-02, -3.1200D-01,  1.1548D-01,  0.0000D+00, &
     4.0556D-02,  1.8664D-02, -9.9661D-03,  0.0000D+00, &
     0.0000D+00,  0.0000D+00,  0.0000D+00,  0.0000D+00, & !Table 2 of Dreamweaver_notes
    -8.0155D+00,  1.8371D+00,  4.5233D-02, -5.9898D-03, &
     3.6916D-01, -1.4520D-02, -9.5873D-04, -5.2315D-04, &
    -1.4307D-02, -3.0142D-02,  7.1303D-03, -3.3802D-04, &
     3.4718D-03,  2.3603D-03, -6.8538D-04,  3.9405D-05  & !Table 3 of Dreamweaver_notes
  ], SHAPE(amn), order=[3, 2, 1])
  REAL(dp), DIMENSION(4) :: cz
  REAL(dp), DIMENSION(4, 4), SAVE :: czmn = RESHAPE([ &
      -6.000D-05,  5.252D-02,  1.285D-01,  0.000D+00, &
      -1.850D-03,  7.355D-02,  7.171D-02, -2.723D-02, &
      -7.930D-02,  3.323D-01, -1.234D-01,  1.530D-02, &
       2.200D-01,  0.000D+00,  0.000D+00,  0.000D+00  & !eqns. (26)-(29), Dreamweaver_notes
  ], SHAPE(czmn), order=[2, 1])
  join = 0.0_dp !zero-out array
  verbinden: DO q = 1, 2
    join(q) = cjoin(q, 7) !initialize JOIN with CJOIN
    join_cjoin: DO m = 6, 1, -1
      join(q) = 1.0D-03 * ziel%I_pot(absf) * join(q) + cjoin(q, m) !units of 1000 eV
    END DO join_cjoin
  END DO verbinden
  prnglo = 0.0_dp
  logt = LOG(e1f * TAU)
  logi = LOG(ziel%I_pot(absf))
  prnglo = 0.0_dp !zero-out array
  proton_range: DO q = 1, 3 !get proton range, compute eqn. (9), Dreamweaver_notes
    loglambda = 0.0_dp
    four_amn: DO m = 4, 1, -1
      term = amn(q, m, 4)
      three_amn: DO n = 3, 1, -1
        term = term * logt + amn(q, m, n)
      END DO three_amn
      loglambda = loglambda + term
      loglambda = loglambda * logi
    END DO four_amn
    loglambda = loglambda / logi !cancel extra factor of LOGI
    loglambda = loglambda + LOG(ziel%a2(absf) / ziel%z2(absf))
    prnglo(q) = EXP(loglambda)
    IF (q == 2) prnglo(q) = prnglo(q) * 1.0D-03 !units of 1000 eV
  END DO proton_range
  g = 1.0_dp + (e1f / 931.49406121_dp)
  b = SQRT(1.0_dp - 1.0_dp / (g * g))
  x = 137.0_dp * b / projektil%z1
  cz = 0.0_dp !zero-out array
  range_extension: DO m = 1, 4 !range extension for emulsion
    cz(m) = 0.0_dp
    four_cz: DO n = 4, 1, -1
      cz(m) = cz(m) + czmn(m, n)
      cz(m) = cz(m) * x
    END DO four_cz
    cz(m) = cz(m) / x !cancel extra factor of X
  END DO range_extension
  q = 2 !medium
  IF (e1f < join(1)) q = 1 !low
  IF (e1f > join(2)) q = 3 !high
  IF (x <= 0.2_dp) THEN !eqn. (26), Dreamweaver_notes
    n = 1
  ELSE IF (x > 0.2_dp .AND. x <= 2.0_dp) THEN !eqn. (27), Dreamweaver_notes
    n = 2
  ELSE IF (x > 2.0_dp .AND. x <= 3.0_dp) THEN !eqn. (28), Dreamweaver_notes
    n = 3
  ELSE IF (x > 3.0_dp) THEN !eqn. (29), Dreamweaver_notes
    n = 4
  END IF
  bzz = (31.8_dp + 3.86_dp * EXP((5.0_dp / 8.0_dp) * logi)) * (ziel%a2(absf) / ziel%z2(absf)) * 1.0D-06                            &
    * EXP((8.0_dp / 3.0_dp) * LOG(projektil%z1)) !eqn. (25), Dreamweaver_notes
  Benton = ((projektil%a1 / TAU) / (projektil%z1 * projektil%z1)) * (prnglo(q) + bzz * cz(n))
END FUNCTION Benton