CHEB| Fit Function | \[\begin{align*} \text{pfit} (pet) &= \exp\left(\frac{1}{2} \text{pcf}(1) + \sum_{i=1}^{8} \text{pcf}(i+1) \cdot T_i(X)\right) \\ \text{where} \quad X &= \frac{\Big(\ln(pet) - \ln(\text{pcf}(10))\Big) - \Big(\ln(\text{pcf}(11)) - \ln(pet)\Big)}{{\ln(\text{pcf}(11)) - \ln(\text{pcf}(10))}}\\ T_0(X) &=1; \; T_1(X) = X; T_{i+1} = 2XT_i(X) - T_{i-1}(X) \\ \end{align*}\] |
| Comments | Python code requires NumPy imported as `np`. |
| Arguments |
|
||||||||||||||||||||||||
| Return values |
|
||||||||||||||||||||||||
| Code | c
c######################################################################
c
subroutine alcheb(pet, pcf, kncf, pfit, kermsg)
c
c this is an ornl:cfadc subroutine to calculate cross sections in
c (cm[2]) versus energy in (ev/amu) or rate coefficients in
c (cm[3]/s) versus maxwellian temperature in (ev) from chebyshev
c polynomial fitting coefficients
c
c these fits are valid only between the limits emin and emax,
c which are coefficients pcf(10) and pcf(11) in the entry data field
c
c pet = collision energy in ev/amu or maxwellian temperature in ev
c
c kermsg = blank if no errors
c
c pfit = cross section in cm[2] or rate coefficient in cm[3]/s
c
c written by h. hunter, cfadc oak ridge national laboratory
c (modified to aladdin calling structure 4/21/88 r.a. hulse)
c
c------------------------------------------------------------------------
c
double precision pet, pcf, pfit
double precision emin, emax, cheb, eminl, emaxl, enl, xnorm
double precision twox, prev, prev2
dimension pcf(11)
character*(*) kermsg
emin = pcf(10)
emax = pcf(11)
if(pet .ge. emin .and. pet .le. emax) then
kermsg = ' '
else
kermsg = 'outside range of fit in alcheb'
return
endif
c
c calculate polynomial using recursion relation
c
k = 9
cheb = pcf(k)
eminl = dlog(emin)
emaxl = dlog(emax)
enl= dlog(pet)
k = k-1
xnorm = (enl-eminl-(emaxl-enl)) / (emaxl-eminl)
twox = 2.0d0 * xnorm
prev2 = 0.0d+00
10 prev = cheb
if(k .ne. 1) then
cheb = pcf(k) + twox*prev - prev2
prev2 = prev
k = k-1
go to 10
endif
cheb = 0.5d0*pcf(1) + xnorm*prev - prev2
pfit = dexp(cheb)
100 return
c
end |
| Arguments |
|
||||||||||||
| Return values |
|
||||||||||||
| Code | def cheb(pet, pcf):
"""
This function calculates cross sections in cm^2 versus energy in eV/amu
or rate coefficients in cm^3/s versus electron temperature in eV
using Chebyshev polynomial fitting coefficients.
pe: collision energy in eV/amu or electron temperature in eV
pcf: parameter data array
pcf[0:9]: parameters for fit to the cross section
pcf[9]: emin
pcf[10]: emax, the fit is valid between the limits emin and emax
"""
emin = pcf[9]
emax = pcf[10]
if not (emin <= pet <= emax):
raise ValueError('Energy outside range of validity of fit in alcheb')
k = 9
cheb = pcf[k]
eminl = np.log(emin)
emaxl = np.log(emax)
enl = np.log(pet)
k -= 1
xnorm = (enl - eminl - (emaxl - enl)) / (emaxl - eminl)
twox = 2.0 * xnorm
prev2 = 0.0
while k != 0:
prev = cheb
cheb = pcf[k] + twox * prev - prev2
prev2 = prev
k -= 1
cheb = 0.5 * pcf[0] + xnorm * prev - prev2
pfit = np.exp(cheb)
return pfit |