MODULE Shared USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION IMPLICIT NONE INTEGER, PARAMETER :: MAXE = 200, Target_dat = 10, dEdx_dat = 11 !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 / 931.49406121, 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)' CHARACTER(*), PARAMETER :: dEdx_format = '(6f22.15)' !.15 is dp 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 dEdx USE Shared, ONLY : dp, Target_dat, dEdx_dat, my_format, dEdx_format, Projectile, Absorber IMPLICIT NONE INTEGER :: k, abs, count, line, MAXAB REAL(dp) :: second, second_read, third, depth, depth_read, length, length_read, rel, dE, dx, resultat, en, Benton, stoppingPower REAL(dp) :: column1, column2, column3, column4, column5, column6 TYPE(Projectile) :: projektil TYPE(Absorber) :: ziel CHARACTER(LEN=1) :: blank EXTERNAL stoppingPower, Fbrems, Lindhard, Carg, Chyperg, Clngamma, Benton OPEN(Target_dat, file='Target.dat', access='sequential', status='old', action='read') OPEN(dEdx_dat, file='dEdx.dat', access='sequential', status='unknown', action='write', asynchronous='yes') 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, asynchronous='yes') 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) CALL Stuff(MAXAB, projektil, ziel, abs) PRINT *, " " !leave these for readability WRITE(6, '(a16, a9)') "Target material:", ziel%tname(abs) PRINT *, " " WRITE(6, '(a37)', advance='no') "Initial kinetic energy, E_i [MeV/n]: " READ *, second_read second = second_read PRINT *, " " WRITE(6, '(a32)', advance='no') "Depth of the target, x [g/cm2]: " READ *, depth_read depth = depth_read PRINT *, " " WRITE(6, '(a30)', advance='no') "Unit path length, dx [g/cm2]: " READ *, length_read length = length_read PRINT *, " " column1 = 0.0_dp !depth in g/cm2 go: DO IF (column1 > depth) EXIT third = 0.0_dp !REL is activated by setting this to a positive number resultat = stoppingPower(second, third, projektil, ziel, abs) !get LET column2 = column1 / ziel%rho(abs) !depth in cm column3 = (resultat * FLOAT(NINT(projektil%a1)) * ziel%rho(abs)) / 10.0_dp !LET in keV/um column4 = resultat * FLOAT(NINT(projektil%a1)) !LET in MeV-cm2/g column5 = second !E in MeV/n column6 = second * FLOAT(NINT(projektil%a1)) !E in MeV WRITE(dEdx_dat, dEdx_format) column1, column2, column3, column4, column5, column6 !write en = second !pass column1 = column1 + length !depth in g/cm2 dx = column1 / 20.0_dp !slices KE_f: DO count = 1, 20 !get final kinetic energy rel = 0.0_dp !REL is activated by setting this to a positive number dE = stoppingPower(en, rel, projektil, ziel, abs) !get LET IF (en < 1.0_dp .OR. dE < 0.0_dp) EXIT !less than 1 MeV/n en = en - (dx * dE) !kinetic energy decreases END DO KE_f IF (en < 1.0_dp .OR. dE < 0.0_dp) THEN !less than 1 MeV/n en = 0.0_dp STOP END IF second = en !E in MeV/n END DO go 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) CLOSE(dEdx_dat) END PROGRAM dEdx 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, '(a24)') "Beams: 0. User Specified" WRITE(6, '(a15)') " 1. 56-Fe" WRITE(6, '(a15)') " 2. 28-Si" WRITE(6, '(a15)') " 3. 16-O " WRITE(6, '(a15)') " 4. 12-C " WRITE(6, '(a15)') " 5. 4-He" WRITE(6, '(a15)') " 6. 1-H " PRINT *, " " WRITE(6, '(a21)', advance='no') "Enter the beam type: " READ *, beam_read beam = beam_read IF (beam == 0) THEN PRINT *, " " WRITE(6, '(a4)', advance='no') "Z = " READ *, projektil%z1 WRITE(6, '(a4)', 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, '(a9, "1.", a9)') "Targets: ", ziel%tname(1) available: DO i = 2, MAXABs PRINT '(i10, ".", a9)', i, ziel%tname(i) END DO available PRINT *, " " WRITE(6, '(a23)', advance='no') "Enter the target type: " READ *, abss END SUBROUTINE Stuff REAL(dp) FUNCTION stoppingPower(secondf, thirdf, projektil, ziel, absf) 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 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 CHARACTER(LEN=1) FD, FE, FSH, FLE, FKIN, FRAD, FLS, FBA, FBR EXTERNAL Fbrems, Lindhard g = 1.0_dp + (secondf / 931.49406121_dp) !gamma FD = 'Y' !density effect correction IF (FD == 'Y') THEN X0 = ziel%X0(absf) X1 = ziel%X1(absf) cbar = 2.0_dp * LOG(ziel%I_pot(absf) / ziel%pla(absf)) + 1.0_dp !Seltzer, eq. (14) b = SQRT(1.0_dp - 1.0_dp / (g * g)) !beta 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) 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)) FE = 'Y' !electron capture IF (FE == 'Y') THEN 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 !everything else: Weaver, eq. (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 FSH = 'Y' !shell corrections FLE = 'N' !Leung relativistic shell correction (only place "bindf") IF (FSH == 'Y') THEN 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)) !eq. (64) f2 = f2 - (cadj / ziel%z2(absf)) !Weaver, eq. (62) IF (FLE == 'Y') THEN 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, eq.(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 FRAD = 'N' !radiative correction IF (FRAD == 'Y') THEN 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, eq.(46) END IF !end FRAD; not able to account for "8.0314" FKIN = 'N' !kinematic correction IF (FKIN == 'Y') THEN 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, eq. (48) END IF !end FKIN FBR = 'N' !projectile bremsstrahlung correction IF (FBR == 'Y') THEN 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)) - (2.0_dp & * 2.414_dp * ALPHA * ALPHA * ziel%z2(absf) * ziel%z2(absf)) !Weaver, eq. (50); Weaver, eq. (54) accounts for the inadequacy of the first Born approximation for large Z_2 IF (Blim < Bbr) THEN Bbr = Blim END IF Sbr = 9.7822719D-08 * ((z1 * z1 * z1 * z1 * ziel%z2(absf) * ziel%z2(absf)) / (a1 * a1 * ziel%a2(absf))) * g * Bbr !Weaver, eq. (49) END IF !end FBR FLS = 'Y' !LS replaces BMA IF (FLS == 'Y') THEN f3 = Lindhard(projektil, b) END IF !end FLS FBA = 'Y' !Barkas correction IF (FBA == 'Y') THEN v = (b * g) / (ALPHA * SQRT(ziel%z2(absf))) !Weaver, eq. (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, eq. (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) stoppingPower = REL !Restricted Energy Loss ELSE stoppingPower = S END IF END FUNCTION stoppingPower REAL(dp) FUNCTION Lindhard(projektil, bf) !Lindhard-Sorenson correction replaces Bloch-Mott-Ahlen USE Shared, ONLY : dp, PI, ALPHA, COMPTON, Projectile IMPLICIT NONE TYPE(Projectile) :: projektil REAL(dp) :: bf 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 CHARACTER(LEN=1) FNS EXTERNAL Chyperg, Clngamma dk = 0.0_dp !zero-out this array (good housekeeping) a3 = EXP(LOG(projektil%a1) / 3.0_dp) eta = projektil%z1 * ALPHA / bf gg = 1.0_dp / SQRT(1.0_dp - bf * bf) !gamma rho = a3 * COMPTON prh = bf * gg * rho total = 0.0_dp !the whole point of this function is to get this sum term1 = 0.0_dp !terms in Weaver, eq. (34) term2 = 0.0_dp term3 = 1.0_dp !why? FNS = 'Y' !finite nuclear size correction accounts for the finite size of atomic nuclei IF (gg < 10.0_dp / rho .OR. FNS == '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, eq. (30), top IF (k > 0.0_dp) THEN !Weaver, eq. (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, eq. (30), middle Cedr = Cexir * EXP(CMPLX(0.0_dp, -AIMAG(Clngamma(Cske)) + (PI / 2.0_dp) * (L - sk), KIND=dp)) !Weaver, eq. (35), I.; solve (29) for xi_k, then put into left side of middle (30) then a^(n + m) = a^n * a^m IF (FNS == 'Y') THEN 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, eq. (35), II.; same way 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, eq. (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, eq. (39); underflow located with -ffpe-trap=underflow IF (COS(AIMAG(Clngamma(Cbbs))) < 1.0_dp) grgs = grgs * (-1.0_dp) IF (ABS(grgs) > 1.0D-09) THEN frgr = SQRT((gg - 1.0_dp) / (gg + 1.0_dp)) * REAL(Clamr) / AIMAG(Clamr) !Weaver, eq. (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, eq. (41), 1 a0 = ((2.0_dp * ABS(k) + 1.0_dp) * b0) / (rho - gz) !eq. (41), 2 a1 = 0.5_dp * (gz + rho) * b0 !eq. (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) !eq. (41), 4 anp1 = ((gz + rho) * bn - ALPHA * z1 * bnm1 / 2.0_dp) / (2.0_dp * (nn + 1.0_dp)) !eq. (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, eq. (40) ELSE figi = bsum / asum END IF H = ((frgr - figi) / (figi - fsgs)) * grgs !Weaver, eq. (35), III.; Weaver, eq. (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, eq. (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, eq. (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, eq. (34), top END IF sdd = SIN(dk(1) - dk(3)) term2 = (k0 * sdd * sdd) / (eta * eta * (4.0_dp * k0 * k0 - 1.0_dp)) !doesn't match eq. (34), bottom term3 = term1 - (1.0_dp / k0) !doesn't match eq. (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, eq. (43), ultrarelativistic limit; asymptotic value of the LS correction END IF Lindhard = total + 0.5_dp * bf * bf !Weaver, eq. (34) END FUNCTION Lindhard REAL(dp) FUNCTION Fbrems(x) 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 !why not PI^2/6.0? 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, eq. (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, eq. (53); why not PI^2/6.0? 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 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, previous_term, term, sum_term sum_term = CMPLX(1.0_dp, 0.0_dp, KIND=dp) !Weaver, eq. (B.3) term = CMPLX(1.0_dp, 0.0_dp, KIND=dp) !Weaver, eq. (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) !each term in the series can be obtained from the previous one by multiplying by (a+n-1)z/(b+n-1)n sum_term = sum_term + term previous_term = term m = m + 1 END DO series Chyperg = sum_term END FUNCTION Chyperg 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, cterm, aterm1, aterm2, aterm3, lterm1, lterm2, lterm3, num, denom REAL(dp), DIMENSION(6) :: coeff COMPLEX(dp) :: result coeff = [76.180092D+00, -86.505320D+00, 24.014098D+00, -1.2317396D+00, 0.1208651D-02, -0.5395239D-05] !Weaver, eq. (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) !y/2.0 in front? aterm1 = y * LOG(r) !Weaver, eq. (A.6), first two "arg" lines aterm2 = (x + 0.5_dp) * ATAN2(y, x + 5.5_dp) - y !0.5* in front? lterm1 = (x + 0.5_dp) * LOG(r) !Weaver, eq. (A.5), first two "ln" lines lterm2 = -y * ATAN2(y, x + 5.5_dp) - (x + 5.5_dp) + 0.5_dp * LOG(2.0_dp * PI) num = 0.0_dp !"T" denom = 1.0_dp !"B", c_o A3: DO j = 1, 6 cterm = coeff(j) / ((x + FLOAT(j)) * (x + FLOAT(j)) + y * y) !Weaver, eq. (A.3) num = num + cterm denom = denom + (x + FLOAT(j)) * cterm END DO A3 num = (-y) * num !Weaver, eq. (A.3) aterm3 = ATAN2(num, denom) !Weaver, eq. (A.6), last line lterm3 = 0.5_dp * LOG(num * num + denom * denom) !Weaver, eq. (A.5), last line result = CMPLX(lterm1 + lterm2 + lterm3, aterm1 + aterm2 + aterm3, KIND=dp) IF (REAL(z) < 0.0_dp) THEN !reflection formula; Weaver, eq. (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