FNAGEX| Fit Function | \[p_{\text{fit}} (pt) = \text{rnor} + \text{rres}\] |
| Comments | Python code requires NumPy imported as `np`. In the above equation, rres and rnor represent rate coefficients for the resonance and non-resonance regions, respectively. These rate coefficients can be calculated using a power-law or exponential type of fit based upon the input fitting parameters. For more details, please refer to the source code. |
| Arguments |
|
||||||||||||||||||||||||
| Return values |
|
||||||||||||||||||||||||
| Code | c
c###################################################################
c
function eiexp(x)
double precision x, x2, x3, x4, x5
double precision a0, a1, a2, a3, a4, a5, b1, b2, b3, b4
c
c eiexp=e1(x)*exp(x)
c handbook of mathematical functions, page 231
c by m. abramowitz and i. a. stegun
c (routine copied from ippj-am-27, y. itikawa et al 1983)
x2=x*x
x3=x2*x
x4=x3*x
x5=x4*x
c
if (x.le.1.0) then
a0=-0.57721566
a1=0.99999193
a2=-0.24991055
a3=0.05519968
a4=-0.00976004
a5=0.00107857
eiexp=(a0 +a1*x +a2*x2 +a3*x3 +a4*x4 +a5*x5 - log(x))*exp(x)
c
else
a1=8.5733287401
a2=18.0590169730
a3=8.6347608925
a4=.2677737343
b1=9.5733223454
b2=25.6329561486
b3=21.0996530827
b4=3.9584969228
eiexp=(x4 +a1*x3+ a2*x2+ a3*x+ a4) /
1 (x4 +b1*x3 +b2*x2 +b3*x+ b4)/x
endif
c
return
end function eiexp
c###################################################################
c
subroutine fnagex(pt, pcf, kncf, pfit, kermsg)
c
c this is an iaea subroutine to calculate the electron impact
c excitation rate coefficients as a function of the electron
c temperature.
c
c pt = electron temperature in ev
c
c pcf is the coefficient data array, where
c
c pcf(1) = itype, index for type of fit, either 1 or 2
c
c pcf(2) = excitation energy , v (also referred to as delta e)
c
c pcf(3) = lower limit of fitting of reduced energy x
c
c pcf(4) = upper limit of fitting of reduced energy x
c
c pcf(5) = statistical weight of initial state (2s+1)*(2l+1)
c
c pcf(6) = parameter a
c
c pcf(7) = parameter b
c
c pcf(8) = parameter c
c
c pcf(9) = parameter d
c
c pcf(10) = parameter e
c
c if ftype = 1 this can be followed by three more parameters for
c the region which contains resonances represented by a linear term
c
c pcf(11) = parameter p
c
c pcf(12) = parameter q
c
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
c if ftype = 2, the pcf array element 10 can be followed by
c
c pcf(11) = parameter f
c
c and possibly three more parameters follow for the region which
c contains resonances represented by a linear term
c
c pcf(12) = parameter p
c
c pcf(13) = parameter q
c
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
c pfit = rate coefficient in cm[3]/s
c
c written by j. j. smith , iaea atomic and molecular data unit
c (taken from report ippj-am-27, y. itikawa et al, nagoya,
c institute of plasma physics, nsgoya univ., (1983))
c
c------------------------------------------------------------------------
c
double precision pt, pcf, pfit, eiexp
dimension pcf(14)
character*(*) kermsg
data s/0.5/
c
c DH, Dec, 16, 2003 ---- kermsg = ' ' to avoid loop on error message ----
kermsg = ' '
itype = pcf(1)
const = 8.010e-8 / (dsqrt(pt)*pcf(5))
y = pcf(2) / pt
c
if (kncf .lt. 12) then
c
c--- rate coefficient without resonances
c
if (itype .eq. 1) then
c
c--- rate coefficient represented by power-log type of fit
c
ra = pcf(6)/y + pcf(8) + s * pcf(9) * (1.0 - y)
rb = pcf(7) - pcf(8)*y + s * pcf(9) * y * y + pcf(10)/y
rnor = const * exp(-y) * y * (ra + eiexp(y) * rb)
rres = 0.0
c
else if (itype .eq. 2) then
c
c--- rate coefficient represented by exponential type of fit
c
ra = pcf(6) * (1.0 - eiexp(y)*y)
f = pcf(11)
rb = pcf(7) * exp(-f) / (f+y)
1 + pcf(8) * exp (-2.*f) / (2.*f+y)
2 + pcf(9) * exp (-3.*f) / (3.*f+y)
3 + pcf(10) * exp (-4.*f) / (4.*f+y)
rnor = const * exp(-y) * y * (ra + rb)
rres = 0.0
c
else
kermsg =
1 ' invalid integer for fit type in fnagex - must be 1 or 2'
return
endif
c
else
c
c--- rate coefficient with resonances
c
if (itype .eq. 1) then
c
c--- rate coefficient represented by power-log type of fit
c
x1 = pcf(13)
y1 = y * x1
ra = pcf(6)/y + pcf(8)/x1
1 + s * pcf(9) * (1./(x1 * x1) - y/x1)
2 + pcf(10) * log(x1) / y
rb = pcf(7) - pcf(8)*y + s * pcf(9) * y * y + pcf(10)/y
rnor = const * y * exp(-y1) * (ra + eiexp(y1) * rb)
rres = const * exp (-y) * (pcf(11) * (1. + 1./y)
1 * (1. - exp ((1. - x1) * y) * (x1 + 1./y)
2 / (1. + 1./y)) + pcf(12) * (1. - exp ((1. - x1) * y )))
c
else if (itype .eq. 2) then
c
c--- rate coefficient represented by exponential type of fit
c
x1 = pcf(14)
y1 = y * x1
ra = pcf(6) * (1. / x1 - eiexp(y1) * y)
f = pcf(11)
rb = pcf(7) * exp(-f*x1) / (f+y)
1 + pcf(8) * exp (-2.*f*x1) / (2.*f+y)
2 + pcf(9) * exp (-3.*f*x1) / (3.*f+y)
3 + pcf(10) * exp (-4.*f*x1) / (4.*f+y)
rnor = const * y * exp(-y1) * (ra + rb)
rres = const * exp (-y) * (pcf(12) * (1. + 1./y)
1 * (1. - exp ((1. - x1) * y) * (x1 + 1./y)
2 / (1. + 1./y)) + pcf(13) * (1. - exp ((1. - x1) * y )))
else
kermsg =
1 ' invalid integer for fit type in fnagex - must be 1 or 2'
return
endif
c
endif
c
pfit = rnor + rres
c
if (pfit .le. 0.0d0) then
kermsg =
1 ' error reaction rate is negative check data and temp. range '
return
endif
c
end |
| Arguments |
|
||||||||||||
| Return values |
|
||||||||||||
| Code | def eiexp(x):
x2 = x * x
x3 = x2 * x
x4 = x3 * x
x5 = x4 * x
if x <= 1.0:
a0 = -0.57721566
a1 = 0.99999193
a2 = -0.24991055
a3 = 0.05519968
a4 = -0.00976004
a5 = 0.00107857
eiexp = (a0 + a1 * x + a2 * x2 + a3 * x3 + a4 * x4 + a5 * x5 - np.log(x)) * np.exp(x)
else:
a1 = 8.5733287401
a2 = 18.0590169730
a3 = 8.6347608925
a4 = 0.2677737343
b1 = 9.5733223454
b2 = 25.6329561486
b3 = 21.0996530827
b4 = 3.9584969228
eiexp = (x4 + a1 * x3 + a2 * x2 + a3 * x + a4) / (x4 + b1 * x3 + b2 * x2 + b3 * x + b4) / x
return eiexp
def fnagex(pt, pcf):
"""
This function calculates the electron-impact excitation rate coefficients
as a function of the electron temperature.
pt: electron temperature in eV
pcf: parameter data array
pcf[0]: type of the fit, to represent resonace and non-resonace region
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
"""
const = 8.010e-8 / (np.sqrt(pt) * pcf[5])
y = pcf[2] / pt
if len(pcf) < 12:
if pcf[1] == 1:
ra = pcf[6] / y + pcf[8] + 0.5 * pcf[9] * (1.0 - y)
rb = pcf[7] - pcf[8] * y + 0.5 * pcf[9] * y * y + pcf[10] / y
rnor = const * np.exp(-y) * y * (ra + eiexp(y) * rb)
rres = 0.0
elif pcf[1] == 2:
ra = pcf[6] * (1.0 - eiexp(y) * y)
f = pcf[11]
rb = pcf[7] * np.exp(-f) / (f + y) + pcf[8] * np.exp(-2.0 * f) / (2.0 * f + y) \
+ pcf[9] * np.exp(-3.0 * f) / (3.0 * f + y) + pcf[10] * np.exp(-4.0 * f) / (4.0 * f + y)
rnor = const * np.exp(-y) * y * (ra + rb)
rres = 0.0
else:
raise ValueError('Invalid integer for fit type in fnagex - must be 1 or 2')
else:
if pcf[1] == 1:
x1 = pcf[13]
y1 = y * x1
ra = pcf[6] / y + pcf[8] / x1 + s * pcf[9] * (1.0 / (x1 * x1) - y / x1) + pcf[10] * np.log(x1) / y
rb = pcf[7] - pcf[8] * y + pcf[9] * y * y + pcf[10] / y
rnor = const * y * np.exp(-y1) * (ra + eiexp(y1) * rb)
rres = const * np.exp(-y) * (pcf[11] * (1.0 + 1.0 / y) *
(1.0 - np.exp((1.0 - x1) * y) * (x1 + 1.0 / y) / (1.0 + 1.0 / y))
+ pcf[12] * (1.0 - np.exp((1.0 - x1) * y)))
elif pcf[1] == 2:
x1 = pcf[14]
y1 = y * x1
ra = pcf[6] * (1.0 / x1 - eiexp(y1) * y)
f = pcf[11]
rb = pcf[6] * np.exp(-f * x1) / (f + y) \
+ pcf[7] * np.exp(-2.0 * f * x1) / (2.0 * f + y) \
+ pcf[8] * np.exp(-3.0 * f * x1) / (3.0 * f + y) \
+ pcf[9] * np.exp(-4.0 * f * x1) / (4.0 * f + y)
rnor = const * y * np.exp(-y1) * (ra + rb)
rres = const * np.exp(-y) * (pcf[12] * (1.0 + 1.0 / y) \
* (1.0 - np.exp((1.0 - x1) * y) * (x1 + 1.0 / y) / (1.0 + 1.0 / y)) \
+ pcf[13] * (1.0 - np.exp((1.0 - x1) * y)))
else:
raise ValueError('Invalid integer for fit type in fnagex - must be 1 or 2')
pfit = rnor + rres
if pfit <= 0.0:
raise ValueError('Error: Reaction rate is negative. Check data and temperature range.')
return pfit |