Initial commit from Polytechnique Montreal

1 files changed, 203 insertions, 0 deletions

diff --git a/Dragon/src/LIBSDC.f b/Dragon/src/LIBSDC.f
new file mode 100644
index 0000000..f231c23
--- /dev/null
+++ b/Dragon/src/LIBSDC.f
@@ -0,0 +1,203 @@
+*DECK LIBSDC
+      SUBROUTINE LIBSDC(NBMIX,NGROUP,NBISO,ISONRF,MIX,DEN,MASK,ENER,
+     1 KGAS,DENMAT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Apply Sternheimer density correction to the collision stopping power.
+*
+*Copyright:
+* Copyright (C) 2021 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert and A. Naceur
+*
+*Parameters: input
+* NBMIX   number of mixtures present in the calculation domain.
+* NGROUP  number of energy groups.
+* NBISO   number of isotopes present in the calculation domain.
+* ISONRF  character*12 reference names of isotopes.
+* MIX     mixture number of each isotope (can be zero).
+* DEN     density of each isotope.
+* MASK    mixture mask (=.true. if a mixture is to be made).
+* ENER    energy groups limits.
+* KGAS    state of each mixture (=0: solid/liquid; =1: gas).
+*
+*Parameters: input/output
+* DENMAT  Sterheimer density correction (delta).
+*
+*References:
+* [1]. L. J. Lorence Jr., J. E. Morel, and G. D. Valdez, "Physics guide
+*      to CEPXS: A multigroup coupled electron-photon cross section
+*      generating code," Technical report SAND89-1685, Sandia National
+*      Laboratories, Albuquerque, New Mexico 87185 and Livermore,
+*      California 94550.
+*
+* [2]. Sternheimer, R. M. (1956). Density effect for the ionization
+*      loss in various materials. Physical Review, 103(3), 511.
+*-----------------------------------------------------------------------
+*
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER :: NBMIX,NGROUP,NBISO,MIX(NBISO),KGAS(NBMIX)
+      CHARACTER(LEN=12) :: ISONRF(NBISO)
+      REAL :: DEN(NBISO),ENER(NGROUP+1) !,ESTOP(NBMIX,NGROUP+1)
+      REAL :: DENMAT(NBMIX,NGROUP+1)
+      LOGICAL :: MASK(NBMIX)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER :: I
+      DOUBLE PRECISION :: PITT(12)
+      CHARACTER(LEN=2) :: ELEMNT(100)
+      CHARACTER(LEN=131) :: HSMG
+      DOUBLE PRECISION, PARAMETER :: DM=3.0D0 ! Sternheimer exponent
+      DOUBLE PRECISION, PARAMETER :: CEMASS=0.510976D0 !Electron mass MeV
+      DOUBLE PRECISION, PARAMETER :: C2=0.249467D0 ! 3/8 Thomson xs
+      DOUBLE PRECISION :: XDRCST
+*----
+*  DATA STATEMENTS
+*----
+      DATA (ELEMNT(I),I=1,100) /
+     1 'h',  'he', 'li', 'be', 'b',  'c',  'n',
+     2 'o',  'f',  'ne', 'na', 'mg', 'al', 'si',
+     3 'p',  's',  'cl', 'ar', 'k',  'ca', 'sc',
+     4 'ti', 'v',  'cr', 'mn', 'fe', 'co', 'ni',
+     5 'cu', 'zn', 'ga', 'ge', 'as', 'se', 'br',
+     6 'kr', 'rb', 'sr', 'y',  'zr', 'nb', 'mo',
+     7 'tc', 'ru', 'rh', 'pd', 'ag', 'cd', 'in',
+     8 'sn', 'sb', 'te', 'i',  'xe', 'cs', 'ba',
+     9 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu',
+     1 'gd', 'tb', 'dy', 'ho', 'er', 'tm', 'yb',
+     2 'lu', 'hf', 'ta', 'w',  're', 'os', 'ir',
+     3 'pt', 'au', 'hg', 'tl', 'pb', 'bi', 'po',
+     4 'at', 'rn', 'fr', 'ra', 'ac', 'th', 'pa',
+     5 'u',  'np', 'pu', 'am', 'cm', 'bk', 'cf',
+     6 'es', 'fm' /
+      DATA (PITT(I),I=1,12) /
+     1  18.7D0, 42.0D0, 38.0D0, 60.0D0, 71.0D0, 78.0D0,
+     2  85.0D0, 89.0D0, 92.0D0, 131.0D0, 146.0D0, 156.0D0 /
+*----
+*  MAIN LOOP OVER MIXTURES
+*----
+      AVCON=1.0D-24*XDRCST('Avogadro','N/moles')
+      DO IBM=1,NBMIX
+        IF(MASK(IBM)) THEN
+*----
+*  CALCULATE THE MEAN IONIZATION ENERGY, EION, IN EV
+*----
+          EION=0.0D0
+          ZZA=0.0D0
+          DO ISOT=1,NBISO
+            IF((MIX(ISOT).NE.IBM).OR.(DEN(ISOT).EQ.0.0)) CYCLE
+            IZ=0
+            DO I=1,100
+              IF(ISONRF(ISOT)(:2).EQ.ELEMNT(I)) THEN
+                IZ=I
+                EXIT
+              ENDIF
+            ENDDO
+            IF(IZ.EQ.0) THEN
+              WRITE(HSMG,'(40HLIBSDC: UNABLE TO ASSIGN AN ATOMIC NUMBE,
+     1        5HR TO ,A,1H.)') ISONRF(ISOT)(:2)
+              CALL XABORT(HSMG)
+            ENDIF
+            WAZ=DEN(ISOT)*REAL(IZ)/AVCON
+            IF(IZ.GE.13) THEN
+*             for Z > 13, use definition of mean ionization given by
+*             Sternheimer
+              PIT=(9.76D0+58.8D0*(REAL(IZ)**(-1.19D0)))*REAL(IZ)
+            ELSE
+*             obtain ionization energy from the data statement
+              PIT=PITT(IZ)
+            ENDIF
+            EION=EION+WAZ*LOG(PIT)
+            ZZA=ZZA+WAZ
+          ENDDO
+          EION = EXP(EION/ZZA)  
+*----
+*  EVALUATE PLANCK'S CONSTANT TIMES THE PLASMA FREQUENCY IN EV
+*----
+          HNUP = 28.8D0*SQRT(ZZA)
+*----
+*  EVALUATE PARAMETERS IN THE STERNHEIMER FORMALISM
+*----
+          C = - (2.0D0*LOG(EION/HNUP) + 1.0D0)
+          IF(KGAS(IBM).EQ.0) THEN
+*           The material is a solid/liquid
+            IF(EION .GE. 100.D0) THEN
+              X1 = 3.0D0
+              IF(-C .GE. 5.215D0) THEN
+                X0 = -0.326D0 * C - 1.5D0
+              ELSE
+                X0 = 0.2D0
+              ENDIF
+            ELSE
+              X1 = 2.0D0
+              IF(-C .GE. 3.681D0) THEN
+                X0 = -0.326D0 * C - 1.0D0
+              ELSE
+                X0 = 0.2D0
+              ENDIF
+            ENDIF
+          ELSE
+*           The material is a gas
+            IF(-C .LT. 12.25D0) THEN
+              X1 = 4.0D0
+              X0 = 2.0D0
+              IF(-C .LT. 11.5D0) X0 = 1.9D0
+              IF(-C .LT. 11.0D0) X0 = 1.8D0
+              IF(-C .LT. 10.5D0) X0 = 1.7D0
+              IF(-C .LT. 10.0D0) X0 = 1.6D0
+            ELSE IF(-C .GE. 13.804D0) THEN
+              X1 = 5.0D0
+              X0 = -.326D0 * C - 2.5D0
+            ELSE
+              X1 = 5.0D0
+              X0 = 2.0D0
+            ENDIF
+          ENDIF
+*
+          IF(X1.LT.X0) THEN
+            WRITE(HSMG,*) 'LIBSDC: NEGATIVE REAL TO REAL POWER. HAVE ',
+     1      'YOU NEGLECTED THE "GAS" KEYWORD FOR A GASEOUS MIXTURE?'
+            CALL XABORT(HSMG)
+          ENDIF
+          B = (-C - 4.606D0*X0) / (X1 - X0)**DM
+          CONV = 2.0D0*AVCON*C2*CEMASS*ZZA
+*----
+*  CALCULATE THE DENSITY CORRECTION FACTOR
+*----
+          DO LLL=1,NGROUP+1
+            T = ENER(LLL)/CEMASS/1.0D6
+            T1 = T + 1.0D0
+            T2 = T + 2.0D0
+            BSQ = T * T2 / T1**2
+*----
+*  EVALUATE THE ELECTRON'S MOMENTUM / (MASS OF THE ELECTRON * C)
+*----
+            PMC = SQRT(2.0D0*T + T*T)
+*----
+*  EVALUTE THE ENERGY PARAMETER IN THE STERNHEIMER FORMALISM
+*----
+            X = LOG10(PMC)
+            IF(X.LE.X0) THEN
+              DS = 0.0D0
+            ELSE IF(X.LE.X1) THEN
+              DS = 4.606D0*X+C+B*((X1-X)**DM)
+            ELSE
+              DS = 4.606D0*X+C
+            ENDIF
+            IF(DS.LT.0.0) DS=0.0D0
+            DENMAT(IBM,LLL)= REAL(CONV*DS/BSQ)
+          ENDDO
+        ENDIF
+      ENDDO
+      RETURN
+      END