diff options
Diffstat (limited to 'Dragon/src/LIBSDC.f')
| -rw-r--r-- | Dragon/src/LIBSDC.f | 203 |
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 |