Python |
def fnagex2(pt, pcf, kncf):
"""
This function calculates the electron-impact excitation cross sections
as a function of the electron energy in eV.
pt: electron energy in eV
pcf: parameter data array
pcf[0]: type of the fit, either 1 or 2
pcf[1]: excitation energy
pcf[2]: lower limit of fitting of reduced energy x
pcf[3]: upper limit of fitting of reduced energy x
pcf[4]: statistical weight of initial state (2s+1)*(2l+1)
pcf[5:14]: fitting parameters, depending upon the type of fit
"""
itype = pcf[0]
pia02 = 0.87973554e-16
rydberg = 13.6056981
wi = pcf[4]
vif = pcf[1]
A = pcf[5]
B = pcf[6]
C = pcf[7]
D = pcf[8]
E = pcf[9]
X = pt / vif
if itype == 2:
F = pcf[10]
ioffset = 1
else:
ioffset = 0
if kncf >= 12:
P = pcf[11 + ioffset]
Q = pcf[12 + ioffset]
X1 = pcf[13 + ioffset]
itype = 3
else:
X1 = 0.0
if itype == 3:
cs = P * X + Q
elif itype == 1:
cs = A + B / X + C / (X * X) + D / (X * X * X) + E * np.log(X)
elif itype == 2:
cs = A / (X * X) + B * np.exp(-F * X) + C * np.exp(-2 * F * X) + \
D * np.exp(-3 * F * X) + E * np.exp(-4 * F * X)
else:
raise ValueError("Invalid integer for fit type in fnagex (1 or 2)")
pfit = (cs * pia02 * rydberg) / (wi * pt)
if pfit <= 0.0:
raise ValueError("Error: Reaction rate is negative. Check data and temperature range.")
return pfit |
Fortran |
c
c##############################################################################
c
c
c###################################################################
c
subroutine fnagex2(pt, pcf, kncf, pfit, kermsg)
c
c this is an iaea subroutine to calculate the electron impact
c excitation cross sections as a function of the electron
c energy in eV.
c
c Cross sections are calculated from the collision strengths:
c xs (pi*a0**2) = cs / wi / Ee (Ryd) where
c Ee energy of the impact electron in Rydberg
c wi statical weight of the initial state
c
c pt = electron energy in ev
c pcf is the coefficient data array, where
c pcf(1) = itype, index for type of fit, either 1 or 2
c pcf(2) = excitation energy , v (also referred to as delta e)
c pcf(3) = lower limit of fitting of reduced energy x
c pcf(4) = upper limit of fitting of reduced energy x
c pcf(5) = statistical weight of initial state (2s+1)*(2l+1)
c pcf(6) = parameter A
c pcf(7) = parameter B
c pcf(8) = parameter C
c pcf(9) = parameter D
c pcf(10) = parameter E
c type 1: cs = A + B/X + C/(X*X) + D/(X*X*X) + E*log(X)
c
c if ftype = 2, the pcf array element 10 can be followed by
c pcf(11) = parameter F
c type 2: cs = A/(X*X) + B*exp(-F*X) + C*exp(-2*F*X) + D*exp(-3*F*X) + E*exp(-4*F*X)
c
c This can be followed by three more parameters for
c the region which contains resonances represented by a linear term
c resonnance region: cs = PX + Q where 1 < X < X1
c type 1:
c pcf(11) = parameter P
c pcf(12) = parameter Q
c pcf(13) = parameter X1, the upper limit of the range over which
c the collision strength is represented by a linear
c approximation.
c type 2:
c pcf(12) = parameter P
c pcf(13) = parameter Q
c pcf(14) = parameter X1, the upper limit of the range over which
c the collision strength is represented by a linear
c approximation.
c
c kermsg = blank if no errors
c pfit = cross section in cm2
c
c written by D. Humbert , iaea atomic and molecular data unit
c fit function from ADNDT,33,149,(1985) y. itikawa et al, nagoya
c 22 August 2006
c
c------------------------------------------------------------------------
c
double precision pt, pcf, pfit
dimension pcf(14)
character*(*) kermsg
data s/0.5/
c
kermsg = ' '
itype = pcf(1)
pia02 = .87973554e-16
rydberg = 13.6056981
wi = pcf(5)
vif = pcf(2)
A = pcf(6)
B = pcf(7)
C = pcf(8)
D = pcf(9)
E = pcf(10)
c
X = pt/vif
c
if (itype .eq. 2) then
F = pcf(11)
ioffset = 1
else
ioffset = 0
endif
if (kncf .ge. 12) then
P = pcf(11 + ioffset)
Q = pcf(12 + ioffset)
X1 = pcf(13 + ioffset)
itype = 3
else
X1 = 0.
endif
c
c--- rate coefficient with resonances: itype = 3
if (itype .eq. 3) then
cs = P*X + Q
c--- rate coefficient without resonances
else if (itype .eq. 1) then
cs = A + B/X + C/(X*X) + D/(X*X*X) + E*log(X)
else if (itype .eq. 2) then
cs = A/(X*X) + B*exp(-F*X) + C*exp(-2*F*X) +
1 D*exp(-3*F*X) + E*exp(-4*F*X)
else
kermsg = ' invalid integer for fit type in fnagex (1 or 2)'
return
endif
c
c collision strengh to cross section
pfit = (cs * pia02 * rydberg) / (wi * pt)
c
if (pfit .le. 0.0d0) then
kermsg =
1 ' error reaction rate is negative check data and temp. range '
return
endif
c
end |