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