"""
The apec module contains routines crucial for the APEC code. This also
includes some interfaces to external C libraries (or will, eventually).
Version 0.1 - initial release
Adam Foster September 16th 2015
"""
import numpy, copy, pickle
import os, time
from . import util, atomdb, const, atomic
import scipy, ctypes
import astropy.io.fits as pyfits
#from joblib import Parallel, delayed
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_full_ionbal_old(Te, tau=False, init_pop='ionizing', Te_init=False, Zlist=False, teunit='K',\
extrap=True, cie=True, settings=False, datacache=False, allowmulti=True):
"""
Calculate the ionization balance for all the elements in Zlist.
One of init_pop or Te_init should be set. If neither is set, assume
all elements start from neutral.
Parameters
----------
Te : float
electron temperature in keV or K (default K)
tau : float
N_e * t for the non-equilibrium ioniziation (default False, i.e. off)
init_pop : dict of float arrays, indexed by Z
initial populations. E.g. init_pop[6]=[0.1,0.2,0.3,0.2,0.2,0.0,0.0]
Te_init : float
initial ionization balance temperature, same units as Te
Zlist : int array
array of nuclear charges to include in calculation (e.g. [8,26] for
oxygen and iron)
teunit : {'K' , 'keV'}
units of temperatures (default K)
extrap : bool
Extrappolate rates to values outside their given range. (default False)
cie : bool
If true, collisional ionization equilbrium calculation
(tau, init_pop, Te_init all ignored)
settings : dict
See description in atomdb.get_data
datacache : dict
Used for caching the data. See description in atomdb.get_data
allowmulti : bool (default True)
If set, include multiple stage ionization and recombination rates
Returns
-------
final_pop : dict of float arrays, indexed by Z
final populations. E.g. final_pop[6]=[0.1,0.2,0.3,0.2,0.2,0.0,0.0]
"""
# input checking
if teunit.lower() == 'kev':
kT = Te*1.0
elif teunit.lower() == 'ev':
kT = Te/1000.0
elif teunit.lower() == 'k':
kT = Te*const.KBOLTZ
else:
print("*** ERROR: unknown teunit %s, Must be keV or K. Exiting ***"%\
(teunits))
# input checking
if util.keyword_check(Te_init):
if teunit.lower() == 'kev':
kT_init = Te_init*1.0
elif teunit.lower() == 'ev':
kT_init = Te_init/1000.0
elif teunit.lower() == 'k':
kT_init = Te_init*const.KBOLTZ
else:
print("*** ERROR: unknown teunit %s, Must be keV or K. Exiting ***"%\
(teunits))
if not Zlist:
Zlist = list(range(1,29))
if (not util.keyword_check(Te_init)) & (not util.keyword_check(init_pop)) &\
(not util.keyword_check(cie)):
print("Warning: you have not specified an initial temperature or "+\
"ion population: assuming everything is neutral")
init_pop={}
for Z in Zlist:
init_pop[Z] = numpy.zeros(Z+1, dtype=float)
init_pop[Z][0] = 1.0
if (util.keyword_check(Te_init)!=False) & (util.keyword_check(init_pop)!=False):
print("Warning: you have specified both an initial temperature and "+\
"ion population: using ion population")
if cie:
# do some error checking
if util.keyword_check(Te_init):
print("Warning: you have specified both Te and Te_init for a CIE calculation. "+\
"Using Te for ionization balance calculation")
if (not init_pop) | (cie):
init_pop = {}
for Z in Zlist:
ionrate = numpy.zeros(Z, dtype=float)
recrate = numpy.zeros(Z, dtype=float)
multiionrate = []
if cie:
kT_init=kT
for z1 in range(1,Z+1):
tmp = \
atomdb.get_ionrec_rate(kT_init, False, Te_unit='keV', \
Z=Z, z1=z1, datacache=datacache, extrap=extrap,\
settings=settings, multiion=allowmulti)
if allowmulti:
ionrate[z1-1]=tmp[0]
recrate[z1-1]=tmp[1]
multiionrate.append(tmp[2])
else:
ionrate[z1-1], recrate[z1-1] =tmp
# now solve
init_pop[Z] = solve_ionbal(ionrate, recrate, multiionrate=multiionrate)
if cie:
# Exit here if the inital
return init_pop
# now solve the actual ionization balance we want.
pop = {}
for Z in Zlist:
# ionrate = numpy.zeros(Z, dtype=float)
# recrate = numpy.zeros(Z, dtype=float)
# dblionrate = numpy.zeros(Z, dtype=float)
# for z1 in range(1,Z+1):
# ionrate[z1-1], recrate[z1-1], dblionrate[z1-1]= \
# atomdb.get_ionrec_rate(kT, False, Te_unit='keV', \
# Z=Z, z1=z1, datacache=datacache, extrap=extrap, \
# doubleion=allowdouble)
# now solve
# pop[Z] = solve_ionbal(ionrate, recrate, doubleionrate = dblionrate, init_pop=init_pop[Z], tau=tau)
pop[Z] = calc_elem_ionbal(Z, kT, tau=tau, init_pop=init_pop,\
teunit='keV', extrap=extrap, settings=settings,\
datacache=datacache, multiion=allowmulti)
return pop
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_full_ionbal(Te, tau=False, init_pop='ionizing', Zlist=False, teunit='K',\
extrap=True, cie=True, settings=False, datacache=False, allowmulti=True):
"""
Calculate the ionization balance for all the elements in Zlist.
One of init_pop or Te_init should be set. If neither is set, assume
all elements start from neutral.
Parameters
----------
Te : float
electron temperature in keV or K (default K)
tau : float
N_e * t for the non-equilibrium ioniziation (default False, i.e. off)
init_pop : string or float
if 'ionizing' : all ionizing from neutral (so [1,0,0,0...])
if 'recombining': all recombining from ionized (so[...0,0,1])
if dict of arrays : the actual fractional populations (so init_pop[6] is the array for carbon)
if single float : the temperature (same units as Te)
Zlist : int array
array of nuclear charges to include in calculation (e.g. [8,26] for
oxygen and iron)
teunit : {'K' , 'keV'}
units of temperatures (default K)
extrap : bool
Extrappolate rates to values outside their given range. (default False)
cie : bool
If true, collisional ionization equilbrium calculation
(tau, init_pop, Te_init all ignored)
settings : dict
See description in atomdb.get_data
datacache : dict
Used for caching the data. See description in atomdb.get_data
allowmulti : bool (default True)
If set, include multiple stage ionization and recombination rates
Returns
-------
final_pop : dict of float arrays, indexed by Z
final populations. E.g. final_pop[6]=[0.1,0.2,0.3,0.2,0.2,0.0,0.0]
"""
# input checking
kT = util.convert_temp(Te, teunit, 'keV')
# input checking
if not cie:
# now we care about the initial population
if isinstance(init_pop, str):
if init_pop.lower() == 'ionizing':
init_pop_calc = {}
for Z in Zlist:
init_pop_calc[Z]=numpy.zeros(Z+1)
init_pop_calc[Z][0] =1.0
elif init_pop.lower() == 'recombining':
init_pop_calc = {}
for Z in Zlist:
init_pop_calc[Z]=numpy.zeros(Z+1)
init_pop_calc[Z][-1] =1.0
else:
raise util.OptionError("Error: init_pop is set as a string, must be 'ionizing' or 'recombining'. Currently %s."%\
(init_pop))
elif isinstance(init_pop, float):
# this is an initial temperature
kT_init = util.convert_temp(init_pop, teunit, 'keV')
# just provide temperature
init_pop_calc={}
for Z in Zlist:
init_pop_calc[Z] = kT_init
elif isinstance(init_pop, dict):
# user has provided all ion abundances
init_pop_calc = init_pop[Z]
if len(init_pop_calc) != Z+1:
print("Error: population supplied has %i ions, should have %i"%(len(init_pop_calc), Z+1))
else:
raise util.OptionError("Error: invalid type for init_pop")
if not Zlist:
Zlist = list(range(1,29))
if cie:
# do some error checking
if util.keyword_check(init_pop):
if init_pop.lower() != 'ionizing':
print("Warning: you have specified an initial population for a CIE calculation. "+\
"Ignoring initial population.")
pop = {}
for Z in Zlist:
pop[Z] = calc_elem_ionbal(Z, kT, tau=tau, init_pop=init_pop_calc[Z], teunit='keV',\
extrap=extrap, settings=settings, datacache=datacache, \
allowmulti=allowmulti)
return pop
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_elem_ionbal(Z, Te, tau=False, init_pop='ionizing', teunit='K',\
extrap=True, settings=False, datacache=False, \
allowmulti=True):
"""
Calculate the ionization balance for all the elements in Zlist.
If tau is set, then we perform an NEI calculation.
This code opens up the ionization and recombination rate files. It is
not the fast code based on the eigenvector files.
Parameters
----------
Z : int
nuclear charge to include in calculation (e.g. 8 for oxygen)
Te : float
electron temperature in keV or K (default K)
tau : float
N_e * t for the non-equilibrium ioniziation (default False, i.e. CIE)
init_pop : string or float
if 'ionizing' : all ionizing from neutral (so [1,0,0,0...])
if 'recombining': all recombining from ionized (so[...0,0,1])
if array: actual populations (e.g. [0, 0.1, 0.3, 0.5, 0.1, 0, 0])
if single float : the temperature (same units as Te)
teunit : {'K' , 'keV'}
units of temperatures (default K)
extrap : bool
Extrappolate rates to values outside their given range. (default False)
settings : dict
See description in atomdb.get_data
datacache : dict
Used for caching the data. See description in atomdb.get_data
allowmulti : bool (default True)
If set, include multiple stage ionization and recombination rates
Returns
-------
final_pop : array
final population. E.g. [0.1,0.2,0.3,0.2,0.2,0.0,0.0]
"""
kT = util.convert_temp(Te, teunit, 'keV')
if tau==False:
cie = True
init_pop_calc=False
else:
cie = False
if not cie:
# if it's not equilibrium, get the initial population
if isinstance(init_pop, str):
if init_pop.lower() == 'ionizing':
init_pop_calc = numpy.zeros(Z+1)
init_pop_calc[0] = 1.0
elif init_pop.lower() == 'recombining':
init_pop_calc = numpy.zeros(Z+1)
init_pop_calc[-1] = 1.0
else:
raise util.OptionError("Error: init_pop is set as a string, must be 'ionizing' or 'recombining'. Currently %s."%\
(init_pop))
elif isinstance(init_pop, float):
# this is an initial temperature
kT_init = util.convert_temp(init_pop, teunit, 'keV')
# rerun this routine in equilibrium mode to find the initial ion pop
init_pop_calc = calc_elem_ionbal(Z, kT_init, \
teunit='keV', \
datacache=datacache,\
settings = settings, extrap=extrap,\
allowmulti=allowmulti)
elif isinstance(init_pop, numpy.ndarray) or isinstance(init_pop, list):
init_pop_calc = init_pop
elif isinstance(init_pop, dict):
init_pop_calc = init_pop[Z]
else:
raise util.OptionError("Error: invalid type for init_pop")
if (cie):
A = numpy.zeros((Z+1, Z+1), dtype=float)
if cie:
kT_init=kT
for z1 in range(1,Z+1):
tmp = \
atomdb.get_ionrec_rate(kT_init, False, Te_unit='keV', \
Z=Z, z1=z1, datacache=datacache, extrap=extrap,\
settings=settings, multiion=allowmulti)
A[z1,z1-1] += tmp[0]
A[z1-1,z1] += tmp[1]
if allowmulti:
# in this case, multiple ionizations are returned as a dict
# indexed from 2 (for double ionization) in tmp[2].
for i in tmp[2].keys():
A[z1-1+i, z1-1]+=tmp[2][i]
init_pop_calc = solve_multi_ionbal(A)
if cie:
return(init_pop_calc)
# get the end point population
A = numpy.zeros((Z+1, Z+1), dtype=float)
for z1 in range(1,Z+1):
tmp = \
atomdb.get_ionrec_rate(kT, False, Te_unit='keV', \
Z=Z, z1=z1, datacache=datacache, extrap=extrap,\
settings=settings, multiion=allowmulti)
A[z1,z1-1] += tmp[0]
A[z1-1,z1] += tmp[1]
if allowmulti:
for i in tmp[2].keys():
A[z1-1+i, z1-1]+=tmp[2][i]
final_pop = solve_multi_ionbal(A, init_pop=init_pop_calc, tau=tau)
return final_pop
#------------------------------------------
#------------------------------------------
#------------------------------------------
[docs]
def solve_ionbal(ionrate, recrate, init_pop=False, tau=False, multiionrate=None, return_details=False):
"""
solve_ionbal: given a set of ionization and recombination rates, find
the equilibrium ionization balance. If init_pop and tau are set, do an
non-equilibrium calculation starting from init_pop and evolving for
n_e * t = tau (cm^-3 s)
Parameters
----------
ionrate : float array
the ionization rates, starting with neutral ionizing to +1
recrate : float array
the recombination rates, starting with singly ionized recombining to neutral
init_pop : float array
initial population of ions for non-equlibrium calculations. Will be renormalised to 1.
tau : float
N_e * t for the non-equilibrium ioniziation
multiionrate : list of dicts
The ionization rates for multiple stages, from each ion.
e.g. [ {2:1e-9,3:1e-10}, {2:2e-9}, {2:1.5e-9, 3:1.5e-10, 4:2.02e-11}]
would be double & triple ionization rates from neutral, double from 1+,
and double, triple and quadruple from the 2+ stage.
Returns
-------
final_pop : float array
final populations.
full_details : dict
If return_details is set, return parts of the calculation:
'eqpop' : equilibrium popn
'vl_out' : left eigenvector
'vr_out' : right eigenvector
'eig_out' : eigenvalues
Notes
-----
Note that init_pop & final_pop will have 1 more element than ionrate and recrate.
"""
#
# Version 0.1 Initial Release
# Adam Foster 16th September 2015
#
# first, calculate the equilibrium solution
# if (init_pop==False) & (tau==False): do_equilib=True
try:
if init_pop==False:
do_equilib=True
else:
do_equilib=False
except ValueError:
do_equilib=False
Z = len(ionrate)
b = numpy.zeros(Z+1, dtype=float)
a = numpy.zeros((Z+1,Z+1), dtype=float)
if multiionrate is None:
# set up dummy array of zeros
mcirate = numpy.zeros((1,len(ionrate)))
else:
# find the maximum index
maxjump = 0
for i in multiionrate:
if len(i.keys()) > 0:
maxjump = max([maxjump, max(i.keys())])
# create array of double, triple, quadruple etc ionization rates
if Z >=2:
mcirate = numpy.zeros((maxjump-1, Z), dtype=float)
for i in range(len(multiionrate)):
for j in multiionrate[i].keys():
mcirate[j-2,i] = multiionrate[i][j]
else:
mcirate = numpy.zeros((0,Z), dtype=float)
# Add the multi-ionization terms
for im in range(0, mcirate.shape[0]):
for iZ in range(0, Z-(1+im)):
a[iZ+im+2, iZ] += mcirate[im,iZ]
a[iZ,iZ] -= mcirate[im,iZ]
# Add the single ion & recomb terms
for iZ in range(0,Z):
a[iZ,iZ] -= ionrate[iZ]
a[iZ+1,iZ] += ionrate[iZ]
a[iZ,iZ+1] += (recrate[iZ])
a[iZ+1,iZ+1] -= (recrate[iZ])
# conservation of population
for iZ in range(0,Z+1):
a[0,iZ]=1.0
b[0]=1.0
eqpop=numpy.linalg.solve(a,b)
# remove any rounding error negatives
eqpop[eqpop<0] = 0.0
# set neutral popn = 1-sum(other pops)
eqpop[0] = max([1.0-sum(eqpop[1:]), 0])
if do_equilib == True:
return eqpop
# now the NEI part
#remake matrix a
Z=len(ionrate)+1
ndim=Z
AA = numpy.zeros((ndim-1,ndim-1), dtype=float)
# populate with stuff:
# Diagonal terms:
i = numpy.arange(ndim-1, dtype=int)
# add ionization rates
# These stop 1 early as you cannot ionize from final state
AA[i[:-1],i[:-1]] -= ionrate[1:]
# add multi ionization rates
# These stop, e.g. 2 early as you cannot double ionize from final 2 states
for im in range(0, mcirate.shape[0]):
AA[i[:-(im+2)],i[:-(im+2)]] -= mcirate[im, 1:-(im+1)]
# add recombination rates
# These cover the whole range as recrate[0] is recombination *from* the 2nd state
AA[i,i] -= recrate
# all other off diagonals
#[to,from] is the indexing
AA[i[1:],i[:-1]]+=ionrate[1:]
for im in range(0, mcirate.shape[0]):
AA[i[im+2:],i[:-(im+2)]] += mcirate[im, 1:-(im+1)]
AA[i[:-1],i[1:]]+=recrate[1:]
# specific to the 1 & 2 stages diagonal
# see Helfand and Hughes; removing the dF0/dt row leads to these extra terms.
# (double ionization is new but algebra isn't hard).
AA[0,:] -=ionrate[0]
if AA.shape[0] > 1:
for im in range(0, mcirate.shape[0]):
AA[im+1,:] -=mcirate[im,0]
w,v = numpy.linalg.eig(AA)
# now copy VR to AA
AA=v
# b_in is difference betwen initial and final populations
b_in=init_pop[1:]-eqpop[1:]
# solve for initial conditions
b_out=numpy.linalg.solve(AA,b_in)
# include exponential decay term
C = b_out*numpy.exp( w*tau )
# get change
G = numpy.dot(v,C)
# solve for population
Ion_pop=numpy.zeros(Z)
Ion_pop[1:]= G + eqpop[1:]
# make sure everything is > 0
Ion_pop[Ion_pop<0] = 0.0
# set neutral to be the residual (or zero)
Ion_pop[0]=max([1-sum(Ion_pop[1:]), 0.0])
if return_details:
vl = numpy.matrix(v)**-1
full_details={}
full_details['eqpop'] = eqpop
full_details['vr'] = v
full_details['vl'] = vl
full_details['eig'] = w
return full_details
# return the data
return Ion_pop
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def solve_multi_ionbal(ionrecrate, init_pop=False, tau=False, return_details=False):
"""
solve_ionbal: given a set of ionization and recombination rates, find
the equilibrium ionization balance. If init_pop and tau are set, do an
non-equilibrium calculation starting from init_pop and evolving for
n_e * t = tau (cm^-3 s)
Parameters
----------
ionrecrate : float array
A Z+1 x Z+1 element array, with ionization and recombination rates
indexed as [final, initial]
init_pop : float array
initial population of ions for non-equlibrium calculations. Will be renormalised to 1.
tau : float
N_e * t for the non-equilibrium ioniziation
Returns
-------
final_pop : float array
final populations.
full_details : dict
If return_details is set, return parts of the calculation:
'eqpop' : equilibrium popn
'vl_out' : left eigenvector
'vr_out' : right eigenvector
'eig_out' : eigenvalues
Notes
-----
Note that init_pop & final_pop will have 1 more element than ionrate and recrate.
"""
#
#
# first, calculate the equilibrium solution
try:
if init_pop==False:
do_equilib=True
else:
do_equilib=False
except ValueError:
do_equilib=False
# element is number of ions -1
Z = ionrecrate.shape[0]-1
# set up for equilibrium calculation
# lower case for equilibrium, upper case for non-equilibrium
b = numpy.zeros(Z+1, dtype=float)
a = numpy.zeros((Z+1,Z+1), dtype=float)
# add ionrec and diagonal elements into "a"
for iZ in range(0,Z+1):
a[iZ,iZ] -= sum(ionrecrate[:,iZ])
for iiZ in range(0,Z+1):
if iiZ==iZ: continue
a[iiZ,iZ] += ionrecrate[iiZ,iZ]
# conservation of population
for iZ in range(0,Z+1):
a[0,iZ]=1.0
b[0]=1.0
# solve
eqpop=numpy.linalg.solve(a,b)
# remove any rounding error negatives
eqpop[eqpop<0] = 0.0
# set neutral popn = 1-sum(other pops)
eqpop[0] = max([1.0-sum(eqpop[1:]), 0])
# if we are doing an equilibrium calculation, we are done.
if do_equilib == True:
return eqpop
# now the NEI part
# ndim=Z+1
# AA = numpy.zeros((ndim-1,ndim-1), dtype=float)
ndim = Z+1
# AAA array is smaller than ionrecrate, as we omit the neutral.
AAA = numpy.zeros((ndim-1, ndim-1), dtype=float)
# populate with stuff:
# DIAGONAL TERMS, omitting ionization from the neutral
for ifrom in range(1,ndim):
for ito in range(ndim):
AAA[ifrom-1,ifrom-1] -= ionrecrate[ito, ifrom]
# OFF DIAGONAL TERMS, again omitting the neutral
for ifrom in range(1,ndim):
for ito in range(1, ndim):
if ito == ifrom: continue
AAA[ito-1,ifrom-1] += ionrecrate[ito, ifrom]
# IONIZATION FROM NEUTRAL ION
# This gets subtracted from all ions.
for ito in range(1,ndim):
AAA[ito-1,:] -= ionrecrate[ito,0]
w,v = numpy.linalg.eig(AAA)
# now copy VR to AA
AA=v
# b_in is difference betwen initial and final populations
b_in=init_pop[1:]-eqpop[1:]
# solve for initial conditions
b_out=numpy.linalg.solve(AA,b_in)
# include exponential decay term
C = b_out*numpy.exp( w*tau )
# get change
G = numpy.dot(v,C)
# solve for population
Ion_pop=numpy.zeros(Z+1)
Ion_pop[1:]= G + eqpop[1:]
# make sure everything is > 0
Ion_pop[Ion_pop<0] = 0.0
# set neutral to be the residual (or zero)
Ion_pop[0]=max([1-sum(Ion_pop[1:]), 0.0])
if return_details:
vl = numpy.matrix(v)**-1
full_details={}
full_details['eqpop'] = eqpop
full_details['vr'] = v
full_details['vl'] = vl
full_details['eig'] = w
full_details['ionpop']=Ion_pop
return full_details
# return the data
return Ion_pop
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_brems_gaunt(E, T, z1, brems_type, datacache=False, \
settings=False):
"""
calculate the bremstrahulung free-free gaunt factor
Parameters
----------
E : float
Energy (in keV) to calculate gaunt factor
T : float
Temperature (in K) of plasma
z1 : int
Ion charge +1 of ion (e.g. 6 for C VI)
brems_type : int
Type of bremstrahlung requested:
1 = HUMMER = Non-relativistic: 1988ApJ...327..477H
2 = KELLOGG = Semi-Relativistic: 1975ApJ...199..299K
3 = RELATIVISTIC = Relativistic: 1998ApJ...507..530N
4 = BREMS_NONE = no bremstrahlung
settings : dict
See description in atomdb.get_data
datacache : dict
Used for caching the data. See description in atomdb.get_data
Returns
-------
gaunt_ff : float
The gaunt factor for the free-free process.
"""
Evec, Eisvec = util.make_vec(E)
gaunt_ff = numpy.zeros(len(Evec), dtype=float)
NUM_J=8
NUM_I=11
z0 = z1-1.0 # we need to use the actual ion charge
if brems_type==const.HUMMER:
# read in the hummer data
hdat = atomdb.get_data(False, False,'hbrems', settings = settings,\
datacache=datacache)
gaunt_D=hdat['BR_GAUNT'].data['COEFFICIENT']
gamma2 = z0**2 * const.RYDBERG/(const.KBOLTZ*T)
if ((gamma2 < 1e-3) | (gamma2 > 1e3)):
if (z0<10):
print("brems_hummer: Warning, gamma^2 = %e is out of range."%(gamma2))
gaunt_ff[:]=1.0
else:
u = Evec/(const.KBOLTZ*T)
j = numpy.where(u <1.e-4)[0]
if len(j) != 0:
print("brems_hummer: Warning, u is out of range: ", u[j])
gaunt_ff[j]=1.0
gaunt_ff[u>31.6227766]=1.0
# all out of range data points are set to 1.0 now. Do the good stuff.
j = numpy.where(gaunt_ff<1.0)[0]
if len(j) > 0:
x_u = (2*numpy.log10(u[j])+2.5)/5.5
x_g = numpy.log10(gamma2)/3.0
c_j = numpy.zeros(NUM_J, dtype=float)
for jj in range(NUM_J):
# We then sum the Chebyshev series, but only 0.5* the first value
tmpgaunt = gaunt_D[jj*NUM_I:(jj+1)*NUM_I]
tmpgaunt[0] *=0.5
c_j[jj]=numpy.polynomial.chebyshev.chebval(x_g, tmpgaunt)
# again, sum chebshev polynomial with first index * 0.5
c_j[0]*=0.5
for ij, jj in enumerate(j):
gaunt_ff[jj] = numpy.polynomial.chebyshev.chebval(x_u[ij], c_j)
if not Eisvec:
gaunt_ff = gaunt_ff[0]
return gaunt_ff
elif brems_type==const.KELLOGG:
gam2 = numpy.array([.7783,1.2217,2.6234, 4.3766,20.,70.])
gam3 = numpy.array([1.,1.7783,3.,5.6234, 10.,30. ])
a = numpy.array([1.001, 1.004, 1.017, 1.036, 1.056, 1.121, 1.001, 1.005,\
1.017, 1.046, 1.073, 1.115, .9991, 1.005, 1.03, 1.055,\
1.102, 1.176, .997, 1.005, 1.035, 1.069, 1.134, 1.186,\
.9962, 1.004, 1.042, 1.1, 1.193, 1.306,.9874, .9962,\
1.047, 1.156, 1.327, 1.485, .9681, .9755, 1.02, 1.208,\
1.525, 1.965, .3029, .1616, .04757, .013, .0049, -.0032,\
.4905, .2155, .08357, .02041, .00739, 2.9e-4, .654, \
.2833, .08057, .03257, .00759, -.00151, 1.029, .391, \
.1266, .05149, .01274, .00324, .9569, .4891, .1764, \
.05914, .01407, -2.4e-4, 1.236, .7579, .326, .1077, \
.028, .00548, 1.327, 1.017,.6017, .205, .0605, .00187, \
-1.323, -.254, -.01571, -.001, -1.84e-4, 8e-5, -4.762, \
-.3386, -.03571, -.001786, -3e-4, 1e-5, -6.349, -.4206, \
-.02571, -.003429, -2.34e-4, 5e-5, -13.231, -.59, \
-.04571, -.005714, -4.45e-4, -4e-5, -7.672, -.6852, \
-.0643, -.005857, -4.2e-4, 4e-5, -7.143, -.9947, \
-.12, -.01007, -8.51e-4, -4e-5, -3.175, -1.116, -.227, \
-.01821, -.001729, 2.3e-4])
kT = const.KBOLTZ*T
gam = z0 * z0 * const.RYDBERG / kT
gam1 = min(gam*1.e3,100.)
if (gam > .1):
gaunt_ff = kurucz(Evec/kT, gam)
elif (kT==0.0):
print("brems_kellog: Zero temperature!")
else:
u=Evec/kT
gaunt_ff[(u>50) | (u==0)]=0.0
i = numpy.where( (u<50)& (u !=0.0))[0]
g=1.0
u2=u[i]**2
u1=0.5*u[i]
t=u1/3.75
u12=u1*0.5
ak = numpy.zeros(len(i), dtype=float)
ii = numpy.where(u12<=1.0)[0]
t2=t[ii]**2
t4=t2**2
t8=t4**2
ai = t2 * 3.5156229 + 1. + t4 * 3.089942 + t2 * 1.2067492 * t4 +\
t8 * .2659732 + t8 * .0360768 * t2 + t8 * .0045813 * t4
u122 = u12[ii] * u12[ii]
ak[ii] = -numpy.log(u12[ii]) * ai - .57721566 + u122 *\
(u122 * \
(u122 * \
(u122 * \
(u122 * \
(u122 * 7.4e-6 + 1.075e-4) + \
.00262698) +\
.0348859) +\
.23069756) +\
.4227842)
ii = numpy.where(u12>1.0)[0]
type2=ii*1
u12im = -1. / u12[ii]
ak[ii] = u12im*\
(u12im*\
(u12im*\
(u12im*\
(u12im*\
(u12im*5.3208e-4 + .0025154) +\
.00587872) +\
.01062446) +\
.02189568) +\
.07832358) + \
1.25331414
ak[ii] /= numpy.exp(u1[ii]) * numpy.sqrt(u1[ii])
# see if born approx works
born = numpy.exp(u1) * .5513 * ak
if gam1<1.0:
gaunt_ff[i] = born
# print("FORCE BORN")
else:
# go to do polynomial expansion
u1 =u[i]
u1[u1<0.003]=0.003
u2=u1**2
n = numpy.zeros(len(i),dtype=int)
#m = numpy.zeros(len(i),dtype=int)
n[(u1<=0.03)]=1
n[(u1>0.03)&(u1<=0.3)]=2
n[(u1>0.3)&(u1<=1.0)]=3
n[(u1>1.0)&(u1<=5.0)]=4
n[(u1>5.0)&(u1<=15.0)]=5
n[(u1>15.0)]=6
if (gam1<=1.773):
m=1
elif ((gam1>1.773)&(gam1<=3.0)):
m=2
elif ((gam1>3.0)&(gam1<=5.6234)):
m=3
elif ((gam1>5.6234)&(gam1<=10.0)):
m=4
elif ((gam1>10.0)&(gam1<=30.0)):
m=5
else:
m=6
m1=m+1
g1= (a[n+(m+7)*6-49] +\
a[n+(m+14)*6-49]*u1+\
a[n+(m+21)*6-49]*u2)*born
g2 = (a[n+(m1+7)*6-49] +\
a[n+(m1+14)*6-49]*u1+\
a[n+(m1+21)*6-49]*u2)*born
p=(gam1-gam3[m-1])/gam2[m-1]
gaunt_ff[i]= (1.0-p)*g1 + p*g2
#aktype=numpy.ones(len(ak), dtype=int)
#aktype[type2]=2
# for i in range(len(E)):
# print "E[%i]=%e ak type %i, m=%i, n=%i, u=%e, g_ff=%e"%\
# (i,E[i], aktype[i],m,n[i], u[i],gaunt_ff[i])
#
# zzz=raw_input()
if not Eisvec:
gaunt_ff = gaunt_ff[0]
return gaunt_ff
elif brems_type==const.RELATIVISTIC:
#initialize
rdat = atomdb.get_data(False, False,'rbrems', settings = settings,\
datacache=datacache)
# gaunt_D=rdat['BR_GAUNT'].data['COEFFICIENT']
gaunt_U = rdat['GAUNT_FF'].data['LOG10_U']
gaunt_Z = rdat['GAUNT_FF'].data['Z']
gaunt_Ng = rdat['GAUNT_FF'].data['N_GAMMA2']
gaunt_g2 = rdat['GAUNT_FF'].data['LOG10_GAMMA2']
gaunt_gf = rdat['GAUNT_FF'].data['GAUNT_FF']
gamma2 = numpy.log10(z0*z0*const.RYDBERG/(const.KBOLTZ*T))
# extract the gaunt factors
if z0 in gaunt_Z:
# exact element is in file
Uvec, GauntFFvec = extract_gauntff(z0, gamma2, gaunt_U, gaunt_Z, gaunt_Ng, gaunt_g2, gaunt_gf)
# print "we are here"
# for iuvec in xrange(len(Uvec)):
# print "Uvec[%i]=%e, GauntFFvec[%i]=%e"%(iuvec, Uvec[iuvec],\
# iuvec, GauntFFvec[iuvec])
#zzz=raw_input()
else:
# find nearest elements in the file
zlo = gaunt_Z[numpy.where(gaunt_Z < z0)[0]]
if len(zlo) ==0:
zlo=0
else:
zlo = max(zlo)
zup = gaunt_Z[numpy.where(gaunt_Z > z0)[0]]
if len(zup)==0:
zup=100
else:
zup=min(zup)
# check for various cases:
if (zlo==0) & (zup<100):
Uvec, GauntFFvec = extract_gauntff(zup, gamma2, gaunt_U, gaunt_Z, gaunt_Ng, gaunt_g2, gaunt_gf)
if (zlo>0) & (zup==100):
Uvec, GauntFFvec = extract_gauntff(zlo, gamma2, gaunt_U, gaunt_Z, gaunt_Ng, gaunt_g2, gaunt_gf)
if (zlo==0) & (zup==100):
#no match found
print("brems relativistic: we should never be here")
if (zlo>0) & (zup<100):
# we are going to interpolate between the 2 of these
Uvecl, GauntFFvecl = extract_gauntff(zlo, gamma2, gaunt_U, gaunt_Z, gaunt_Ng, gaunt_g2, gaunt_gf)
Uvecu, GauntFFvecu = extract_gauntff(zup, gamma2, gaunt_U, gaunt_Z, gaunt_Ng, gaunt_g2, gaunt_gf)
if len(numpy.where(numpy.abs(Uvecl-Uvecu)>0.001)[0]) != 0 :
print("Error: brems_relativistic: U vector mismatch ", Uvecl, Uvecu)
Uvec = Uvecl
GauntFFvec = ((zup-z0)*GauntFFvecl + (z0-zlo)*GauntFFvecu)/(zup-zlo)
# print GauntFFvec
# ok, so now we have Uvec and GauntFFvec assigned
u = numpy.log10(Evec/(const.KBOLTZ*T))
#for (iBin=0;iBin < N; iBin++) {
# First, some special cases on u
# Long wavelength. Use approximation from Scheuer 1960, MNRAS, 120,231
# As noted in Hummer, 1988, ApJ, 325, 477
# low energy
i = numpy.where(u<Uvec[0])[0]
gaunt_ff[i]=-0.55133*(0.5*numpy.log(10.)*gamma2+numpy.log(10.)*u[i]+0.056745)
#high energy
i = numpy.where(u>Uvec[-1])[0]
gaunt_ff[i]=1.0
#other energy
i = numpy.where((u>=Uvec[0]) & (u <=Uvec[-1]))[0]
gaunt_ff[i] = numpy.interp(u[i],Uvec, GauntFFvec)
if not Eisvec:
gaunt_ff = gaunt_ff[0]
return gaunt_ff
else:
print("UNKNOWN BREMS TYPE: ", brems_type)
return -1
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def do_brems(Z, z1, T, abund, brems_type, eedges):
"""
Calculate the bremstrahlung emission in units of photon cm^3 s^-1 bin^-1
Parameters
----------
Z : int
nuclear charge for which result is required
z1 : int
ion charge +1
T : float
temperture (Kelvin)
abund : float
elemental abundance (should be between 1.0 and 0.0)
brems_type : int
Type of bremstrahlung requested:
1 = HUMMER = Non-relativistic: 1988ApJ...327..477H
2 = KELLOGG = Semi-Relativistic: 1975ApJ...199..299K
3 = RELATIVISTIC = Relativistic: 1998ApJ...507..530N
4 = BREMS_NONE = no bremstrahlung
eedges : array(float)
The energy bin edges for the spectrum (keV)
Returns
-------
array(float)
bremstrahlung emission in units of photon cm^3 s^-1 bin^-1
"""
# If neutral, there is no brems as no free electrons
if z1==1:
emission = numpy.zeros(len(eedges)-1, dtype=float)
return emission
kT = T*const.KBOLTZ # convert to keV
E = (eedges[1:]+eedges[:-1])/2.0
dE= (eedges[1:]-eedges[:-1])
EkT=E/kT
gaunt_ff = calc_brems_gaunt(E, T, z1, brems_type)
if z1 ==1:
# no brems
emission = numpy.zeros(len(dE))
else:
emission = const.BREMS_COEFF*abund*\
(z1-1)*(z1-1)*numpy.exp(-EkT)*(dE/numpy.sqrt(T))*gaunt_ff/(E*const.ERG_KEV)
return emission
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def kurucz(uin, gam):
"""
Correction factors to Kellogg bremstrahlung calculation by Bob Kurucz
Parameters
----------
uin : array(float)
energy grid, units of E/kT (both in keV)
gam : array(float)
Z**2/T, in units of Rydbergs
Returns
-------
array(float)
gaunt factors at high gam (> 0.1)
"""
ya = numpy.array([5.4,5.25, 5.,4.69,4.48,4.16,3.85,4.77,4.63,4.4,\
4.13,3.87,3.52, 3.27,4.15,4.02,3.8,3.57,3.27,\
2.98,2.7,3.54,3.41,3.22,2.97,2.7,2.45,2.2,2.94,\
2.81,2.65,2.44,2.21,2.01,1.81,2.41,2.32,2.19,2.02,\
1.84,1.67,1.5,1.95,1.9,1.8,1.68,1.52,1.41,1.3,\
1.55,1.56,1.51,1.42,1.33,1.25,1.17,1.17,1.3,1.32,\
1.3,1.2,1.15,1.11,.86,1.,1.15,1.18,1.15,1.11,1.08,\
.59,.76,.97,1.09,1.13,1.1,1.08,.38,.53,.76,.96,\
1.08,1.09,1.09])
gaunt = numpy.ones(len(uin), dtype=float)
rj = numpy.log10(gam) * 2. + 3.
j = numpy.array(numpy.floor(rj), dtype=int)
rj = numpy.array(j, dtype=float)
rk = numpy.log10(uin) * 2. + 9.
k = numpy.array(numpy.floor(rk), dtype=int)
k[k<1] = 1
rk = numpy.array(k, dtype=float)
gaunt[(j >= 7) | (j < 1) | (k>=12)]=0.0
i = numpy.where(gaunt>0)[0]
t = (numpy.log10(gam) - (rj - 3.) / 2.) / .5
u = (numpy.log10(uin[i]) - (rk[i] - 9.) / 2.) / .5
# so t is a scalar
# u is a vector
gaunt[i] = (1. - t) * (1. - u) * ya[j + k[i] * 7 - 8] + \
t * (1.-u)* ya[j + 1 + k[i]*7 - 8] + t*u*ya[j + 1 + (k[i] + 1) * 7 - 8] +\
(1. - t) * u * ya[j + (k[i] + 1) * 7 - 8]
return gaunt
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_ee_brems(E, T, N):
"""
calculate the electron-electron bremsstrahlung.
Parameters
----------
E : array (float)
energy grid (keV)
T : float
Electron temperature (keV)
N : float
electron density (cm^-3)
Returns
-------
array(float)
ee_brems in photons cm^s s^-1 keV-1 at each point E.
This should be multiplied by the bin width to get flux per bin.
References
----------
Need to check this!
"""
#
# T is the electron temperature (in keV)
# N is the electron density (in cm^-3)
# E is the list of energies (in keV)
#
# series of data constants
# Region I, k_BT<=1 keV
#print "T=%f"%(T)
aI1 = numpy.array([(3.15847E+0, -2.52430E+0, 4.04877E-1, 6.13466E-1, 6.28867E-1, 3.29441E-1),
(2.46819E-2, 1.03924E-1, 1.98935E-1, 2.18843E-1, 1.20482E-1, -4.82390E-2),
(-2.11118E-2, -8.53821E-2, -1.52444E-1, -1.45660E-1, -4.63705E-2, 8.16592E-2),
(1.24009E-2, 4.73623E-2, 7.51656E-2, 5.07201E-2, -2.25247E-2, -8.17151E-2),
(-5.41633E-3, -1.91406E-2, -2.58034E-2, -2.23048E-3, 5.07325E-2, 5.94414E-2),
(1.70070E-3, 5.39773E-3, 4.13361E-3, -1.14273E-2, -3.23280E-2, -2.19399E-2),
(-3.05111E-4, -7.26681E-4, 4.67015E-3, 1.24789E-2, -1.16976E-2, -1.13488E-2),
(-1.21721E-4, -7.47266E-4, -2.20675E-3, -2.74351E-3, -1.00402E-3, -2.38863E-3),
(1.77611E-4, 8.73517E-4, -2.67582E-3, -4.57871E-3, 2.96622E-2, 1.89850E-2),
(-2.05480E-5, -6.92284E-5, 2.95254E-5, -1.70374E-4, -5.43191E-4, 2.50978E-3),
(-3.58754E-5, -1.80305E-4, 1.40751E-3, 2.06757E-3, -1.23098E-2, -8.81767E-3)])
# column j=0-5; line i=0-10
aI2 = numpy.array([(-1.71486E-1, -3.68685E-1, -7.59200E-2, 1.60187E-1, 8.37729E-2),
(-1.20811E-1, -4.46133E-4, 8.88749E-2, 2.50320E-2, -1.28900E-2),
(9.87296E-2, -3.24743E-2, -8.82637E-2, -7.52221E-3, 1.99419E-2),
(-4.59297E-2, 5.05096E-2, 5.58818E-2, -9.11885E-3, -1.71348E-2),
(-2.11247E-2, -5.05387E-2, 9.20453E-3, 1.67321E-2, -3.47663E-3),
(1.76310E-2, 2.23352E-2, -4.59817E-3, -8.24286E-3, -3.90032E-4),
(6.31446E-2, 1.33830E-2, -8.54735E-2, -6.47349E-3, 3.72266E-2),
(-2.28987E-3, 7.79323E-3, 7.98332E-3, -3.80435E-3, -4.25035E-3),
(-8.84093E-2, -2.93629E-2, 1.02966E-1, 1.38957E-2, -4.22093E-2),
(4.45570E-3, -2.80083E-3, -5.68093E-3, 1.10618E-3, 2.33625E-3),
(3.46210E-2, 1.23727E-2, -4.04801E-2, -5.68689E-3, 1.66733E-2)])
# column j=6-10, line 0-10
# Region II (1 keV<=k_B<=300 keV)
abII = numpy.array([(-3.7369800E+1, -9.3647000E+0, 9.2170000E-1, -1.1628100E+1, -8.6991000E+0),
(3.8036590E+2, 9.5918600E+1, -1.3498800E+1, 1.2560660E+2, 6.3383000E+1),
(-1.4898014E+3, -3.9701720E+2, 7.6453900E+1, -5.3274890E+2, -1.2889390E+2),
(2.8614150E+3, 8.4293760E+2, -2.1783010E+2, 1.1423873E+3, -1.3503120E+2),
(-2.3263704E+3, -9.0730760E+2, 3.2097530E+2, -1.1568545E+3, 9.7758380E+2),
(-6.9161180E+2, 3.0688020E+2, -1.8806670E+2, 7.5010200E+1, -1.6499529E+3),
(2.8537893E+3, 2.9129830E+2, -8.2416100E+1, 9.9681140E+2, 1.2586812E+3),
(-2.0407952E+3, -2.9902530E+2, 1.6371910E+2, -8.8818950E+2, -4.0474610E+2),
(4.9259810E+2, 7.6346100E+1, -6.0024800E+1, 2.5013860E+2, 2.7335400E+1)])
# column a_0j-a2j, b0j,b1j; line j=0-8
cII = numpy.array([(-5.7752000E+0, 3.0558600E+1, -5.4327200E+1, 3.6262500E+1, -8.4082000E+0),
(4.6209700E+1, -2.4821770E+2, 4.5096760E+2, -3.1009720E+2, 7.4792500E+1),
(-1.6072800E+2, 8.7419640E+2, -1.6165987E+3, 1.1380531E+3, -2.8295400E+2),
(3.0500700E+2, -1.6769028E+3, 3.1481061E+3, -2.2608347E+3, 5.7639300E+2),
(-3.2954200E+2, 1.8288677E+3, -3.4783930E+3, 2.5419361E+3, -6.6193900E+2),
(1.9107700E+2, -1.0689366E+3, 2.0556693E+3, -1.5252058E+3, 4.0429300E+2),
(-4.6271800E+1, 2.6056560E+2, -5.0567890E+2, 3.8008520E+2, -1.0223300E+2)])
# column CII_2j-cII_6j; line j=0-6
# Region III (300 keV<= k_BT<=7 MeV)
abIII = numpy.array([(5.2163300E+1, 4.9713900E+1, 6.4751200E+1, -8.5862000E+0, 3.7643220E+2),
(-2.5703130E+2, -1.8977460E+2, -2.1389560E+2, 3.4134800E+1, -1.2233635E+3),
(4.4681610E+2, 2.7102980E+2, 1.7414320E+2, -1.1632870E+2, 6.2867870E+2),
(-2.9305850E+2, -2.6978070E+2, 1.3650880E+2, 2.9654510E+2, 2.2373946E+3),
(0.0000000E+0, 4.2048120E+2, -2.7148990E+2, -3.9342070E+2, -3.8288387E+3),
(7.7047400E+1, -5.7662470E+2, 8.9321000E+1, 2.3754970E+2, 2.1217933E+3),
(-2.3871800E+1, 4.3277900E+2, 5.8258400E+1, -3.0600000E+1, -5.5166700E+1),
(0.0000000E+0, -1.6053650E+2, -4.6080700E+1, -2.7617000E+1, -3.4943210E+2),
(1.9970000E-1, 2.3392500E+1, 8.7301000E+0, 8.8453000E+0, 9.2205900E+1)])
# column aII_0j-aIII_2j, bIII_0j, bIII_1j; line j=0-8
#print "E:",E
aI = numpy.hstack((aI1,aI2))
inum1 = numpy.arange(11)
jnum = numpy.resize(inum1,(11,11))
inum = numpy.transpose(jnum)
aII = numpy.zeros((9,3))
bII = numpy.zeros((9,2))
[aII, bII] = numpy.hsplit(abII,numpy.array([3]))
numII = numpy.arange(9)
aIIj = numpy.transpose(numpy.resize(numII,(3,9)))
aIIi = numpy.resize(numpy.arange(3),(9,3))
bIIj = numpy.transpose(numpy.resize(numpy.arange(9),(2,9)))
bIIi = numpy.resize(numpy.arange(2),(9,2))
cIIj = numpy.transpose(numpy.resize(numpy.arange(7),(5,7)))
cIIi = numpy.resize(numpy.arange(2,7),(7,5))
aIII = numpy.zeros((9,3))
bIII = numpy.zeros((9,2))
kTIII = 1.e3 #1 MeV, in unit of keV
[aIII, bIII] = numpy.hsplit(abIII,numpy.array([3]))
tao = T/const.ME_KEV
# find the length of the input
Earray, Eisvec = util.make_vec(E)
#if isinstance(E, (collections.Sequence, numpy.ndarray)):
# x = E/T
# (numx,) = x.shape
#else:
# Earray = numpy.array([E])
x = Earray/T
numx=len(Earray)
#print "x", x
#print "numx:", numx
GI = numpy.zeros((numx,))
AIIr = numpy.zeros((numx,))
BIIr = numpy.zeros((numx,))
GpwII = numpy.zeros((numx,))
GII = numpy.zeros((numx,))
FCCII = numpy.zeros((numx,))
Ei0 = numpy.zeros((numx,))
GpwIII = numpy.zeros((numx,))
if T<0.05:
ret = numpy.zeros(len(x), dtype=float)
elif 0.05<=T<70.:
# hmm
GI=numpy.zeros(len(x))
theta = (1/1.35) * ( numpy.log10(tao) + 2.65)
bigx = (1/2.5) * (numpy.log10(x) + 1.50)
for i in range(11):
for j in range(11):
GI += aI[i,j]*(theta**i)*(bigx**j)
GI *= numpy.sqrt(8/(3*numpy.pi))
ret = 1.455e-16*N**2*numpy.exp(-x)/(x*numpy.sqrt(tao))*GI
elif 70.<=T<300.:
taoII = tao
for k in range(numx):
def integrand(t):
return numpy.exp(-1.0*t)/t
[Ei0[k,],error] = scipy.integrate.quad(integrand,x[k,],\
numpy.inf,args=())
AIIr[k,] = numpy.sum(aII*taoII**(aIIj/8.)*x[k,]**(aIIi))
BIIr[k,] = numpy.sum(bII*taoII**(bIIj/8.)*x[k,]**(bIIi))
FCCII[k,] = 1.+numpy.sum(cII*taoII**(cIIj/6.)*x[k,]**(cIIi/8.))
GpwII[k,] = numpy.sum(aII*taoII**(aIIj/8.)*x[k,]**(aIIi))-\
numpy.exp(x[k,])*(-1.0)*Ei0[k,]*\
numpy.sum(bII*taoII**(bIIj/8.)*x[k,]**(bIIi))
GII[k,] = GpwII[k,]*FCCII[k,]
ret = 1.455e-16*N**2*numpy.exp(-x)/(x*numpy.sqrt(taoII))*GII
elif 300.<=T<7000.:
taoIII = tao
for k in range(numx):
GpwIII[k,] = numpy.sum(aIII*taoIII**(aIIj/8.)*x[k,]**(aIIi))-\
numpy.exp(x[k,])*(-1.0)*Ei0[k,]*\
numpy.sum(bIII*taoIII**(bIIj/8.)*x[k,]**(bIIi))
ret = 1.455e-16*N**2*numpy.exp(-x)/\
(x*numpy.sqrt(taoIII))*GpwIII
else:
taoIV = tao
GIV = 3./(4.*numpy.pi*numpy.sqrt(taoIV))*\
(28./3.+2.*x+x**2/2.+2.*(8./3.+4.*x/3.+x**2)*\
(numpy.log(2.*taoIV)-0.57721)-numpy.exp(x) \
*(-1.0*Ei0)*(8./3.-4.*x/3.+x**2))
ret = 1.455e-16*N**2*numpy.exp(-x)/(x*numpy.sqrt(taoIV))*GIV
if Eisvec==False:
# convert back from vector to scalar
ret = ret[0]
return ret/const.ME_KEV
[docs]
def parse_par_file(fname):
"""
Parse the apec.par input file for controlling APEC
Parameters
----------
fname : string
file name
Returns
-------
dict
The settings in "key:value" pairs.
"""
# Adam Foster 2016-07-27
# check the par file exists
f = open(fname,'r')
data = {}
rawdata = {}
for line in f:
# check for comment lines
if line[0] == '#': continue
lall = line.split('"')
l = []
for ill, ll in enumerate(lall):
if ill % 2 == 0:
# even, not a quote
ltmp = ll.split(',')
if ill > 0:
ltmp=ltmp[1:]
for iltmp in ltmp:
l.append(iltmp)
else:
l=l[:-1]
l.append(ll)
name = l[0]
dtype = l[1]
hidden = l[2]
value = l[3]
minval = l[4]
maxval= l[5]
descr = l[6]
rawdata[name]={}
rawdata[name]['dtype']=dtype
rawdata[name]['hidden']=hidden
rawdata[name]['value']=value
rawdata[name]['minval']=minval
rawdata[name]['maxval']=maxval
rawdata[name]['descr']=descr
# some checking
if dtype=='s':
value = value.strip('"')
minval = minval.strip('"')
#string
if '|' in minval:
# check the value is one of the options
options = minval.split('|')
if not value in options:
print("Error in paramater file %s. Item %s: %s is not in allowed parameter list %s"%\
(fname, name, value, minval))
print(options)
else:
data[name] = value
else:
data[name] = value
elif dtype =='b':
#boolean, yes or no.
if not value in ['yes','no']:
print("Error in paramater file %s. Item %s: should be boolean yes or no, supplied value is %s"%\
(fname, name, value))
else:
if value=='yes':
data[name]=True
else:
data[name]=False
elif dtype=='r':
# real
value = float(value)
if len(minval) > 0:
minval = float(minval)
if value < minval:
print("Error in paramater file %s. Item %s: should be > %e, supplied value is %e"%\
(fname, name, minval, value))
if len(maxval) > 0:
maxval = float(maxval)
if value > maxval:
print("Error in paramater file %s. Item %s: should be < %e, supplied value is %e"%\
(fname, name, maxval, value))
data[name] = value
elif dtype=='i':
# integer
value = int(value)
if len(minval) > 0:
minval = int(minval)
if value < minval:
print("Error in paramater file %s. Item %s: should be > %i, supplied value is %i"%\
(fname, name, minval, value))
if len(maxval) > 0:
maxval = int(maxval)
if value > maxval:
print("Error in paramater file %s. Item %s: should be < %i, supplied value is %i"%\
(fname, name, maxval, value))
data[name] = value
else:
print("Error in paramater file %s. Item %s: unknown data type %s"%\
(fname, name, dtype))
# now some massaging of the parameters
if 'BremsType' in list(data.keys()):
if data['BremsType']=='Hummer':
data['BremsType']=const.HUMMER
elif data['BremsType']=='Kellogg':
data['BremsType']=const.KELLOGG
elif data['BremsType']=='Relativistic':
data['BremsType']=const.RELATIVISTIC
else:
print("UNKNOWN BREMS TYPE: %s" %(data['BremsType']))
if 'NEIMinEpsilon' in list(data.keys()):
tmp = data['NEIMinEpsilon'].split(',')
tmp2 = []
for i in tmp:
tmp2.append(float(i))
data['NEIMinEpsilon']=tmp2
if 'NEIMinFrac' in list(data.keys()):
tmp = data['NEIMinFrac'].split(',')
tmp2 = []
for i in tmp:
tmp2.append(float(i))
data['NEIMinFrac']=tmp2
if 'IncAtoms' in list(data.keys()):
elsymblist = data['IncAtoms'].split(',')
Zlist=[]
for iel in elsymblist:
Zlist.append(atomic.elsymb_to_Z(iel))
data['Zlist']=Zlist
if not ('WriteIon' in list(data.keys())):
data['WriteIon'] = False
return data
[docs]
def make_vector(linear, minval, step, nstep):
"""
Create a vector from the given inputs
Parameters
----------
linear: boolean
Whether the array should be linear or log spaced
minval: float
initial value of the array. In dex if linear==False
step: float
step between points on the array. In dex if linear==False
nstep: int
number of steps
Returns
-------
array(float)
array of values spaced out use the above parameters
"""
arr = numpy.arange(nstep)*step+minval
if linear==False:
arr=10**arr
return arr
[docs]
def make_vector_nbins(linear, minval, maxval, nstep):
"""
Create a vector from the given inputs
Parameters
----------
linear: boolean
Whether the array should be linear or log spaced
minval: float
initial value of the array. In dex if linear==False
maxval: float
maximum value of the array. In dex if linear==False
nstep: int
number of steps
Returns
-------
array(float)
array of values spaced out use the above parameters
"""
arr = numpy.linspace(minval, maxval, nstep+1)
if linear==False:
arr=10**arr
return arr
[docs]
def run_apec(fname):
"""
Run the entire APEC code using the data in the parameter file fname
Parameters
----------
fname : string
file name
Returns
-------
None
"""
# get the input settings
settings = parse_par_file(fname)
# need to parse the IncAtoms parameter - this defines which atoms
# we will include
#
# we will transfer this to a "Zlist" parameter, in the form of
# nuclear charges
Zlist = settings['Zlist']
print("I will be running Z=", Zlist)
# run for each element, temperature, density
lhdulist = []
chdulist = []
seclHDUdat = numpy.zeros(settings['NumTemp']*settings['NumDens'], \
dtype=generate_datatypes('lineparams'))
seccHDUdat = numpy.zeros(settings['NumTemp']*settings['NumDens'], \
dtype=generate_datatypes('cocoparams'))
for iTe in range(settings['NumTemp']):
te = make_vector(settings['LinearTemp'], \
settings['TempStart'], \
settings['TempStep'], \
settings['NumTemp'])[iTe]
if settings['TempUnits']=='keV':
te /= const.KBOLTZ
for iDens in range(settings['NumDens']):
dens = make_vector(settings['LinearDens'], \
settings['DensStart'], \
settings['DensStep'], \
settings['NumDens'])[iDens]
# AT THIS POINT, GENERATE SHELL WHICH WILL GO IN THE HDU OF CHOICE
if settings['Ionization']=='CIE':
linedata = numpy.zeros(0,dtype=generate_datatypes('linelist_cie'))
cocodata = numpy.zeros(0,dtype=generate_datatypes('continuum', ncontinuum=0, npseudo=0))
if settings['Ionization']=='NEI':
linedata = numpy.zeros(0,dtype=generate_datatypes('linetype'))
cocodata = numpy.zeros(0,dtype=generate_datatypes('continuum', ncontinuum=0, npseudo=0))
for Z in Zlist:
print("Calling run_apec_element for Z=%i Te=%e dens=%e at %s"%(Z, te, dens, time.asctime()))
dat = run_apec_element(settings, te, dens, Z)
# append this data to the output
linedata = numpy.append(linedata, dat['lines'])
cocodata = continuum_append(cocodata, dat['cont'])
# now make an HDU for all of this
if settings['Ionization']=='CIE':
LHDUdat = create_lhdu_cie(linedata)
elif settings['Ionization']=='NEI':
LHDUdat = create_lhdu_nei(linedata)
# now update the headers
iseclHDUdat=iDens+iTe*settings['NumDens']
LHDUdat.header['EXTNAME']=("EMISSIVITY","name of this binary table extension")
LHDUdat.header['EXTVER']=(iseclHDUdat+1,"Index for this EMISSIVITY extension")
LHDUdat.header['HDUNAME'] = ("L%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
'Spectral emission data')
LHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
'Proposed OGIP standard')
LHDUdat.header['HDUCLAS1']=("LINE MODEL",\
'Line emission spectral model')
LHDUdat.header['HDUCLAS2']=("LINE",\
'Emission line data')
if settings['Ionization']=='CIE':
LHDUdat.header['HDUVERS1']=("1.1.0",\
'version of format')
elif settings['Ionization']=='NEI':
LHDUdat.header['HDUVERS1']=("2.0.0",\
'version of format')
LHDUdat.header['TEMPERATURE']=(te,\
'Electron temperature')
LHDUdat.header['DENSITY']=(dens,\
'Electron density')
LHDUdat.header['TIME']=(0,\
'IN EQUILIBRIUM')
if settings['Ionization']=='CIE':
tot_emiss = sum(linedata['epsilon']*const.HC_IN_ERG_A/linedata['lambda'])
else:
tot_emiss=0.0
LHDUdat.header['TOT_LINE']=(tot_emiss,\
'Total Line Emission (erg cm^3 s^-1)')
LHDUdat.header['N_LINES']=(len(linedata),\
'Number of emission lines')
seclHDUdat['kT'][iseclHDUdat]=te*const.KBOLTZ
seclHDUdat['EDensity'][iseclHDUdat]= dens
seclHDUdat['Time'][iseclHDUdat]= 0.0
seclHDUdat['Nelement'][iseclHDUdat]= len(Zlist)
seclHDUdat['Nline'][iseclHDUdat]= len(linedata)
lhdulist.append(LHDUdat)
# continuum data
# now make an HDU for all of this
CHDUdat = create_chdu_variable_cie(cocodata)
# now update the headers
iseccHDUdat=iDens+iTe*settings['NumDens']
CHDUdat.header['EXTNAME']=("EMISSIVITY","name of this binary table extension")
CHDUdat.header['EXTVER']=(iseccHDUdat+1,"Index for this EMISSIVITY extension")
CHDUdat.header['HDUNAME'] = ("C%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
'Spectral emission data')
CHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
'Proposed OGIP standard')
CHDUdat.header['HDUCLAS1']=("COMP CONT MODEL",\
'Compressed continua spectra')
CHDUdat.header['HDUCLAS2']=("COCO",\
'Compressed continuum data')
CHDUdat.header['HDUVERS1']=("2.0.0",\
'version of format')
CHDUdat.header['TEMPERATURE']=(te,\
'Electron temperature')
CHDUdat.header['DENSITY']=(dens,\
'Electron density')
CHDUdat.header['TIME']=("%.2e"%(0),\
'IN EQUILIBRIUM')
tot_emiss = calc_total_coco(cocodata, settings)
if settings['Ionization']=='CIE':
CHDUdat.header['TOT_COCO']=(tot_emiss,\
'Total Emission (erg cm^3 s^-1)')
else:
CHDUdat.header['TOT_COCO']=(0.0,\
'Total Emission (erg cm^3 s^-1)')
seccHDUdat['kT'][iseccHDUdat]=te*const.KBOLTZ
seccHDUdat['EDensity'][iseccHDUdat]= dens
seccHDUdat['Time'][iseccHDUdat]= 0.0
seccHDUdat['NElement'][iseccHDUdat]= len(cocodata)
seccHDUdat['NCont'][iseccHDUdat]= max(cocodata['N_Cont'])
seccHDUdat['NPseudo'][iseccHDUdat]= max(cocodata['N_Pseudo'])
chdulist.append(CHDUdat)
# make secHDUdat into a fits table
seclHDU = create_lparamhdu_cie(seclHDUdat)
seclHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
seclHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
seclHDU.header['HDUCLAS1']=('LINE MODEL','line emission spectra model')
seclHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
seclHDU.header['HDUVERS1']=('1.0.0','version of format')
seccHDU = create_cparamhdu_cie(seccHDUdat)
seccHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
seccHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
seccHDU.header['HDUCLAS1']=('COMP CONT MODEL','compressed continua spectra')
seccHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
seccHDU.header['HDUVERS1']=('1.0.0','version of format')
fileroot = settings['OutputFileStem']
if settings['Ionization']=='NEI':
fileroot+='_nei'
# create the Primary HDU
PrilHDU = pyfits.PrimaryHDU()
lhdulist.insert(0,PrilHDU)
lhdulist.insert(1,seclHDU)
tmplhdulist = pyfits.HDUList(lhdulist)
tmplhdulist.writeto('%s_line.fits'%(fileroot), overwrite=True, checksum=True)
PricHDU = pyfits.PrimaryHDU()
chdulist.insert(0,PricHDU)
chdulist.insert(1,seccHDU)
tmpchdulist = pyfits.HDUList(chdulist)
if settings['Ionization']=='CIE':
tmpchdulist.writeto('%s_coco.fits'%(fileroot), overwrite=True, checksum=True)
elif settings['Ionization']=='NEI':
tmpchdulist.writeto('%s_comp.fits'%(fileroot), overwrite=True, checksum=True)
[docs]
def calc_total_coco(cocodata, settings):
"""
Calculate the total emission in erg cm^3 s^-1
"""
emiss = 0.0
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
ecent = (ebins[1:]+ebins[:-1])/2
width = ebins[1:]-ebins[:-1]
for i in range(len(cocodata)):
E_Cont = cocodata[i]['E_Cont'][:cocodata[i]['N_Cont']]
Cont = cocodata[i]['Continuum'][:cocodata[i]['N_Cont']]
Cont_erg = Cont*E_Cont*const.ERG_KEV
e = numpy.interp(ecent, E_Cont, Cont_erg)
emiss += sum(e*width)
E_Cont = cocodata[i]['E_Pseudo'][:cocodata[i]['N_Pseudo']]
Cont = cocodata[i]['Pseudo'][:cocodata[i]['N_Pseudo']]
Cont_erg = Cont*E_Cont*const.ERG_KEV
e = numpy.interp(ecent, E_Cont, Cont_erg)
emiss += sum(e*width)
return emiss
[docs]
def continuum_append(a,b):
"""
Join two continuum arrays together, expanding arrays as necessary
Parameters
----------
a: numpy.array(dtype=continuum)
The first array
b: numpy.array(dtype=continuum)
The second array
Returns
c: numpy.array(dtype=continuum)
The two arrays combined, with continuum arrays resized as required.
"""
if len(a) == 0:
npseudomax = max(b['N_Pseudo'])
ncontmax = max(b['N_Cont'])
nlines = len(b)
else:
npseudomax = max([max(a['N_Pseudo']), max(b['N_Pseudo'])])
ncontmax = max([max(a['N_Cont']), max(b['N_Cont'])])
nlines = len(a) + len(b)
c = numpy.zeros(nlines, dtype=generate_datatypes('continuum',npseudo=npseudomax, ncontinuum=ncontmax))
ic = 0
if len(a) > 0:
for ia in range(len(a)):
c['Z'][ic] = a['Z'][ia]
c['rmJ'][ic] = a['rmJ'][ia]
c['N_Cont'][ic] = a['N_Cont'][ia]
c['E_Cont'][ic][:c['N_Cont'][ic]] = a['E_Cont'][ia][:a['N_Cont'][ia]]
c['Continuum'][ic][:c['N_Cont'][ic]] = a['Continuum'][ia][:a['N_Cont'][ia]]
c['N_Pseudo'][ic] = a['N_Pseudo'][ia]
c['E_Pseudo'][ic][:c['N_Pseudo'][ic]] = a['E_Pseudo'][ia][:a['N_Pseudo'][ia]]
c['Pseudo'][ic][:c['N_Pseudo'][ic]] = a['Pseudo'][ia][:a['N_Pseudo'][ia]]
ic +=1
for ib in range(len(b)):
c['Z'][ic] = b['Z'][ib]
c['rmJ'][ic] = b['rmJ'][ib]
c['N_Cont'][ic] = b['N_Cont'][ib]
c['E_Cont'][ic][:c['N_Cont'][ic]] = b['E_Cont'][ib][:b['N_Cont'][ib]]
c['Continuum'][ic][:c['N_Cont'][ic]] = b['Continuum'][ib][:b['N_Cont'][ib]]
c['N_Pseudo'][ic] = b['N_Pseudo'][ib]
c['E_Pseudo'][ic][:c['N_Pseudo'][ic]] = b['E_Pseudo'][ib][:b['N_Pseudo'][ib]]
c['Pseudo'][ic][:c['N_Pseudo'][ic]] = b['Pseudo'][ib][:b['N_Pseudo'][ib]]
ic +=1
return c
[docs]
def continuum_append_variable(a,b):
"""
Join two continuum arrays together, expanding arrays as necessary
Parameters
----------
a: numpy.array(dtype=continuum)
The first array
b: numpy.array(dtype=continuum)
The second array
Returns
c: numpy.array(dtype=continuum)
The two arrays combined, with continuum arrays resized as required.
"""
if len(a) == 0:
npseudomax = max(b['N_Pseudo'])
ncontmax = max(b['N_Cont'])
nlines = len(b)
else:
npseudomax = max([max(a['N_Pseudo']), max(b['N_Pseudo'])])
ncontmax = max([max(a['N_Cont']), max(b['N_Cont'])])
nlines = len(a) + len(b)
c = numpy.zeros(nlines, dtype=generate_datatypes('continuum',npseudo=npseudomax, ncontinuum=ncontmax))
ic = 0
if len(a) > 0:
for ia in range(len(a)):
c['Z'][ic] = a['Z'][ia]
c['rmJ'][ic] = a['rmJ'][ia]
c['N_Cont'][ic] = a['N_Cont'][ia]
c['E_Cont'][ic][:c['N_Cont'][ic]] = a['E_Cont'][ia][:a['N_Cont'][ia]]
c['Continuum'][ic][:c['N_Cont'][ic]] = a['Continuum'][ia][:a['N_Cont'][ia]]
c['N_Pseudo'][ic] = a['N_Pseudo'][ia]
c['E_Pseudo'][ic][:c['N_Pseudo'][ic]] = a['E_Pseudo'][ia][:a['N_Pseudo'][ia]]
c['Pseudo'][ic][:c['N_Pseudo'][ic]] = a['Pseudo'][ia][:a['N_Pseudo'][ia]]
ic +=1
for ib in range(len(b)):
c['Z'][ic] = b['Z'][ib]
c['rmJ'][ic] = b['rmJ'][ib]
c['N_Cont'][ic] = b['N_Cont'][ib]
c['E_Cont'][ic][:c['N_Cont'][ic]] = b['E_Cont'][ib][:b['N_Cont'][ib]]
c['Continuum'][ic][:c['N_Cont'][ic]] = b['Continuum'][ib][:b['N_Cont'][ib]]
c['N_Pseudo'][ic] = b['N_Pseudo'][ib]
c['E_Pseudo'][ic][:c['N_Pseudo'][ic]] = b['E_Pseudo'][ib][:b['N_Pseudo'][ib]]
c['Pseudo'][ic][:c['N_Pseudo'][ic]] = b['Pseudo'][ib][:b['N_Pseudo'][ib]]
ic +=1
return c
[docs]
def create_lhdu_cie(linedata):
# sort the data
tmp = numpy.zeros(len(linedata), dtype=numpy.dtype({'names':['negLambda','Element','Ion'],\
'formats':[float, float, float]}))
tmp['Element']= linedata['element']
tmp['Ion']= linedata['ion']
tmp['negLambda']= linedata['lambda']*-1
srt = numpy.argsort(tmp, order=['Element','Ion','negLambda'])
linedata = linedata[srt]
cols = []
cols.append(pyfits.Column(name='Lambda', format='1E', unit="A", array=linedata['lambda']))
cols.append(pyfits.Column(name='Lambda_Err', format='1E', unit="A", array=linedata['lambda_err']))
cols.append(pyfits.Column(name='Epsilon', format='1E', unit="photons cm^3 s^-1", array=linedata['epsilon']))
cols.append(pyfits.Column(name='Epsilon_Err', format='1E', unit="photons cm^3 s^-1", array=linedata['epsilon_err']))
cols.append(pyfits.Column(name='Element', format='1J', array=linedata['element']))
cols.append(pyfits.Column(name='Ion', format='1J', array=linedata['ion']))
cols.append(pyfits.Column(name='UpperLev', format='1J', array=linedata['upperlev']))
cols.append(pyfits.Column(name='LowerLev', format='1J', array=linedata['lowerlev']))
coldefs = pyfits.ColDefs(cols)
tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
print("Created a linelist HDU with %i lines"%(len(tbhdu.data)))
return tbhdu
[docs]
def create_lhdu_nei(linedata):
# sort the data, as required by XSPEC. Cannot sort in reverse order, so hackery ensues.
tmp = numpy.zeros(len(linedata), dtype=numpy.dtype({'names':['negLambda','Element','Ion'],\
'formats':[float, float, float]}))
tmp['Element']= linedata['element']
tmp['Ion']= linedata['ion']
tmp['negLambda']= linedata['lambda']*-1
srt = numpy.argsort(tmp, order=['Element','Ion','negLambda'])
linedata = linedata[srt]
cols = []
cols.append(pyfits.Column(name='Lambda', format='1E', unit="A", array=linedata['lambda']))
cols.append(pyfits.Column(name='Lambda_Err', format='1E', unit="A", array=linedata['lambda_err']))
cols.append(pyfits.Column(name='Epsilon', format='1E', unit="photons cm^3 s^-1", array=linedata['epsilon']))
cols.append(pyfits.Column(name='Epsilon_Err', format='1E', unit="photons cm^3 s^-1", array=linedata['epsilon_err']))
cols.append(pyfits.Column(name='Element', format='1J', array=linedata['element']))
cols.append(pyfits.Column(name='Elem_drv', format='1J', array=linedata['element']))
cols.append(pyfits.Column(name='Ion', format='1J', array=linedata['ion']))
cols.append(pyfits.Column(name='Ion_drv', format='1J', array=linedata['ion_drv']))
cols.append(pyfits.Column(name='UpperLev', format='1J', array=linedata['upperlev']))
cols.append(pyfits.Column(name='LowerLev', format='1J', array=linedata['lowerlev']))
coldefs = pyfits.ColDefs(cols)
tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
return tbhdu
[docs]
def create_chdu_cie(cocodata):
ncont = max(cocodata['N_Cont'])
npseudo = max(cocodata['N_Pseudo'])
cols = []
cols.append(pyfits.Column(name='Z', format='1J', array=cocodata['Z']))
cols.append(pyfits.Column(name='rmJ', format='1J', array=cocodata['rmJ']))
cols.append(pyfits.Column(name='N_Cont', format='1J', array=cocodata['N_Cont']))
cols.append(pyfits.Column(name='E_Cont', format='%iE'%(ncont), unit='keV', array=cocodata['E_Cont']))
cols.append(pyfits.Column(name='Continuum', format='%iE'%(ncont), unit='photons cm^3 s^-1 keV^-1', array=cocodata['Continuum']))
cols.append(pyfits.Column(name='Cont_Err', format='%iE'%(ncont), unit='photons cm^3 s^-1 keV^-1', array=cocodata['Cont_Err']))
cols.append(pyfits.Column(name='N_Pseudo', format='1J', array=cocodata['N_Pseudo']))
cols.append(pyfits.Column(name='E_Pseudo', format='%iE'%(npseudo), unit='keV', array=cocodata['E_Pseudo']))
cols.append(pyfits.Column(name='Pseudo', format='%iE'%(npseudo), unit='photons cm^3 s^-1 keV^-1', array=cocodata['Pseudo']))
cols.append(pyfits.Column(name='Pseudo_Err', format='%iE'%(npseudo), unit='photons cm^3 s^-1 keV^-1', array=cocodata['Pseudo_Err']))
coldefs = pyfits.ColDefs(cols)
tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
return tbhdu
[docs]
def create_chdu_variable_cie(cocodata):
# ncont = max(cocodata['N_Cont'])
# npseudo = max(cocodata['N_Pseudo'])
cols = []
cols.append(pyfits.Column(name='Z', format='1J', array=cocodata['Z']))
cols.append(pyfits.Column(name='rmJ', format='1J', array=cocodata['rmJ']))
cols.append(pyfits.Column(name='N_Cont', format='1J', array=cocodata['N_Cont']))
tmp = []
for i in range(len(cocodata['E_Cont'])):
tmp.append(numpy.array(cocodata['E_Cont'][i][:cocodata['N_Cont'][i]]))
cols.append(pyfits.Column(name='E_Cont', format='PE()', unit='keV', array=tmp))
tmp = []
for i in range(len(cocodata['Continuum'])):
tmp.append(numpy.array(cocodata['Continuum'][i][:cocodata['N_Cont'][i]]))
cols.append(pyfits.Column(name='Continuum', format='PE()', unit='keV', array=tmp))
tmp = []
for i in range(len(cocodata['Cont_Err'])):
tmp.append(numpy.array(cocodata['Cont_Err'][i][:cocodata['N_Cont'][i]]))
cols.append(pyfits.Column(name='Cont_Err', format='PE()', unit='keV', array=tmp))
cols.append(pyfits.Column(name='N_Pseudo', format='1J', array=cocodata['N_Pseudo']))
tmp = []
for i in range(len(cocodata['E_Pseudo'])):
tmp.append(numpy.array(cocodata['E_Pseudo'][i][:cocodata['N_Pseudo'][i]]))
cols.append(pyfits.Column(name='E_Pseudo', format='PE()', unit='keV', array=tmp))
tmp = []
for i in range(len(cocodata['Pseudo'])):
tmp.append(numpy.array(cocodata['Pseudo'][i][:cocodata['N_Pseudo'][i]]))
cols.append(pyfits.Column(name='Pseudo', format='PE()', unit='keV', array=tmp))
tmp = []
for i in range(len(cocodata['Pseudo_Err'])):
tmp.append(numpy.array(cocodata['Pseudo_Err'][i][:cocodata['N_Pseudo'][i]]))
cols.append(pyfits.Column(name='Pseudo_Err', format='PE()', unit='keV', array=tmp))
coldefs = pyfits.ColDefs(cols)
tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
return tbhdu
[docs]
def create_lparamhdu_cie(linedata):
cols = []
cols.append(pyfits.Column(name='kT', format='1E', unit="keV", array=linedata['kT']))
cols.append(pyfits.Column(name='EDensity', format='1E', unit="cm**-3", array=linedata['EDensity']))
cols.append(pyfits.Column(name='Time', format='1E', unit="s", array=linedata['Time']))
cols.append(pyfits.Column(name='Nelement', format='1J', array=linedata['Nelement']))
cols.append(pyfits.Column(name='Nline', format='1J', array=linedata['Nline']))
coldefs = pyfits.ColDefs(cols)
tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
return tbhdu
[docs]
def create_cparamhdu_cie(cocodata):
cols = []
cols.append(pyfits.Column(name='kT', format='1E', unit="keV", array=cocodata['kT']))
cols.append(pyfits.Column(name='EDensity', format='1E', unit="cm**-3", array=cocodata['EDensity']))
cols.append(pyfits.Column(name='Time', format='1E', unit="s", array=cocodata['Time']))
cols.append(pyfits.Column(name='NElement', format='1J', array=cocodata['NElement']))
cols.append(pyfits.Column(name='NCont', format='1J', array=cocodata['NCont']))
cols.append(pyfits.Column(name='NPseudo', format='1J', array=cocodata['NPseudo']))
coldefs = pyfits.ColDefs(cols)
tbhdu = pyfits.BinTableHDU.from_columns(coldefs)
return tbhdu
[docs]
def run_apec_element(settings, te, dens, Z):
"""
Run the APEC code using the settings provided for one element
Parameters
----------
settings: dictionary
The settings read from the apec.par file by parse_par_file
te: float
The electron temperature (K)
dens: float
The electron density (cm^-3)
Z: int
The nuclear charge of the element
Returns
-------
None
"""
if settings['Ionization']=='NEI':
z1list = list(range(1, Z+2))
ionfrac = numpy.ones(len(z1list), dtype=float)
elif settings['Ionization']=='CIE':
z1list = list(range(1, Z+2))
# time to calculate the ionization balance
if settings['UseIonBalanceTable']:
# read the ionization balance table
ionfrac = atomdb._get_precalc_ionfrac(os.path.expandvars(settings['IonBalanceTable']), Z, te)
else:
# calculate the ionization balance
ionftmp = calc_full_ionbal(te, 1e14, Te_init=te, Zlist=[Z], extrap=True)
ionfrac = ionftmp[Z]
else:
print("ERROR: settings['Ionization'] must be CIE or NEI, not %s"%(settings['Ionization']))
abundfile = atomdb.get_filemap_file('abund',\
Z,\
False,\
fmapfile=settings['FileMap'],\
atomdbroot=os.path.expandvars('$ATOMDB'),\
misc=True)
abundances = atomdb.get_abundance(abundfile = abundfile, abundset=settings['Abundances'])
abund=abundances[Z]
# create placeholders for all the data
linelist = numpy.zeros(0, dtype=generate_datatypes('linetype'))
contlist = {}
pseudolist = {}
# now go through each ion and assemble the data
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
ecent = (ebins[1:]+ebins[:-1])/2
z1_drv_list = numpy.arange(1,Z+2, dtype=int)
for z1_drv in range(1, Z+2):
tmplinelist, tmpcontinuum, tmppseudocont = run_apec_ion(settings, te, dens, Z, z1_drv, ionfrac, abund)
linelist = numpy.append(linelist, tmplinelist)
contlist[z1_drv] = tmpcontinuum
pseudolist[z1_drv] = tmppseudocont
# out = Parallel(n_jobs=4, max_nbytes=1e2)(delayed(run_apec_ion)(settings, te, dens,Z,z1_drv,ionfrac,abund) for z1_drv in range(1,Z+2))
#for ii,i in enumerate(out):
#linelist = numpy.append(linelist, i[0])
#contlist[ii+1] = i[1]
#pseudolist[ii+1] = i[2]
# now merge these together.
if settings['Ionization']=='CIE':
cieout = generate_cie_outputs(settings, Z, linelist, contlist, pseudolist)
return cieout
elif settings['Ionization']=='NEI':
ionftmp= calc_full_ionbal(te, 1e14, Te_init=te, Zlist=[Z], extrap=True)
ionfrac_nei = ionftmp[Z]
neiout = generate_nei_outputs(settings, Z, linelist, contlist, pseudolist, ionfrac_nei)
return neiout
[docs]
def generate_nei_outputs(settings, Z, linelist, contlist, pseudolist, ionfrac_nei):
"""
Convert a linelist and continuum values into a non-equilibrium AtomDB fits output
Parameters
----------
settings: dictionary
The settings read from the apec.par file by parse_par_file
Z: int
The nuclear charge of the element
linelist: numpy.array(dtype=linelisttype)
The list of lines, separated by ion
contlist: dict
Dictionary with the different continuum contributions from each ion.
Each is an array of ph cm^3 s^-1 bin^-1
pseudolist: dict
Dictionary with the different pseudocontinuum contributions from each ion.
Each is an array of ph cm^3 s^-1 bin^-1
Returns
-------
None
"""
# ok, so we are first up going to combine the lines from the different ions
linelist.sort(order=['element','ion','upperlev','lowerlev'])
igood = numpy.ones(len(linelist), dtype=bool)
linelist= linelist[igood]
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
pseudo = {}
cont = {}
# now do some weak line filtering
igood = numpy.ones(len(linelist), dtype=bool)
print("initially we have %i lines for Z =%i"%(len(linelist), Z))
for z1 in range(1, Z+2):
ionfrac = ionfrac_nei[z1-1]
mineps = settings['NEIMinEpsilon'][0]
pseudotmp = numpy.zeros(len(ebins)-1, dtype=float)
pseudotmp += pseudolist[z1]
for i in range(len(settings['NEIMinFrac'])):
if ionfrac < settings['NEIMinFrac'][i]:
mineps = settings['NEIMinEpsilon'][i+1]
print("z1 = %i. Ionfrac = %e. mineps = %e"%(z1,ionfrac, mineps))
weaklines = linelist[(linelist['element']==Z) &\
(linelist['ion_drv']==z1) &\
(linelist['epsilon']<mineps) &\
(linelist['lambda']>const.HC_IN_KEV_A /settings['GridMaximum']) &\
(linelist['lambda']<const.HC_IN_KEV_A /settings['GridMinimum'])]
print("identified %i weak lines"%(len(weaklines)))
for line in weaklines:
e = const.HC_IN_KEV_A /line['lambda']
ibin = numpy.where(ebins>e)[0][0] - 1
pseudotmp[ibin]+=line['epsilon']
igood[(linelist['element']==Z) &\
(linelist['ion_drv']==z1) &\
(linelist['epsilon']<mineps)] = False
print("Filtered by mineps %e: from %i to %i lines"%(mineps, len(igood), sum(igood)))
conttmp = contlist[z1]['rrc']+contlist[z1]['twophot']+contlist[z1]['brems']
econt, contin = compress_continuum(ebins, conttmp, const.TOLERANCE, minval=1e-38)
epseudo, pseudocont = compress_continuum(ebins, pseudotmp, const.TOLERANCE, minval=1e-38)
cont[z1] = {}
cont[z1]['E_Cont'] = econt
cont[z1]['Cont'] = contin
cont[z1]['E_Pseudo'] = epseudo
cont[z1]['Pseudo'] = pseudocont
igood[(linelist['lambda']<const.HC_IN_KEV_A /settings['GridMaximum']) |\
(linelist['lambda']>const.HC_IN_KEV_A /settings['GridMinimum'])] = False
print("Filtering lines on wavelength: keeping %i of %i lines"%\
(sum(igood), len(igood)))
linelist = linelist[igood]
ret={}
ret['lines']=linelist
maxnpseudo = 0
maxncont = 0
for i in list(cont.keys()):
if len(cont[i]['E_Cont'])>maxncont:
maxncont= len(cont[i]['E_Cont'])
if len(cont[i]['E_Pseudo'])>maxnpseudo:
maxnpseudo= len(cont[i]['E_Pseudo'])
ret['cont'] = numpy.zeros(len(list(cont.keys())), dtype=generate_datatypes('continuum', npseudo=maxnpseudo, ncontinuum=maxncont))
for iz1, z1 in enumerate(cont.keys()):
ret['cont']['Z'][iz1] = Z
ret['cont']['rmJ'][iz1] = z1
ret['cont']['N_Cont'][iz1] = len(cont[z1]['E_Cont'])
ret['cont']['E_Cont'][iz1][:ret['cont']['N_Cont'][iz1]] = cont[z1]['E_Cont']
ret['cont']['Continuum'][iz1][:ret['cont']['N_Cont'][iz1]] = cont[z1]['Cont']
ret['cont']['N_Pseudo'][iz1] = len(cont[z1]['E_Pseudo'])
ret['cont']['E_Pseudo'][iz1][:ret['cont']['N_Pseudo'][iz1]] = cont[z1]['E_Pseudo']
ret['cont']['Pseudo'][iz1][:ret['cont']['N_Pseudo'][iz1]] = cont[z1]['Pseudo']
print("returning ret['lines'] with length %i"%(len(ret['lines'])))
return ret
[docs]
def generate_cie_outputs(settings, Z, linelist, contlist, pseudolist):
"""
Convert a linelist and continuum values into an equilibrium AtomDB fits output
Parameters
----------
settings: dictionary
The settings read from the apec.par file by parse_par_file
Z: int
The nuclear charge of the element
linelist: numpy.array(dtype=linelisttype)
The list of lines, separated by ion
contlist: dict
Dictionary with the different continuum contributions from each ion.
Each is an array of ph cm^3 s^-1 bin^-1
pseudolist: dict
Dictionary with the different pseudocontinuum contributions from each ion.
Each is an array of ph cm^3 s^-1 bin^-1
Returns
-------
None
"""
# ok, so we are first up going to combine the lines from the different ions
linelist.sort(order=['element','ion','upperlev','lowerlev'])
igood = numpy.ones(len(linelist), dtype=bool)
# pickle.dump(linelist, open('whole_linedebugZ%i.'%(Z), 'wb'))
for i in range(1,len(linelist)):
if ((linelist['element'][i]==linelist['element'][i-1]) &\
(linelist['ion'][i]==linelist['ion'][i-1]) &\
(linelist['upperlev'][i]==linelist['upperlev'][i-1]) &\
(linelist['lowerlev'][i]==linelist['lowerlev'][i-1])):
linelist['epsilon'][i] += linelist['epsilon'][i-1]
igood[i-1] = False
linelist= linelist[igood]
linelist_equ = numpy.zeros(len(linelist), dtype=generate_datatypes('linelist_cie'))
for i in linelist_equ.dtype.names:
linelist_equ[i] = linelist[i]
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
pseudo = numpy.zeros(len(ebins)-1, dtype=float)
# now do some weak line filtering
MinEpsilon = settings['MinEpsilon']
weaklines = linelist_equ[(linelist_equ['epsilon']< MinEpsilon) &\
(linelist_equ['lambda']>const.HC_IN_KEV_A /settings['GridMaximum']) &\
(linelist_equ['lambda']<const.HC_IN_KEV_A /settings['GridMinimum'])]
for line in weaklines:
e = const.HC_IN_KEV_A /line['lambda']
ibin = numpy.where(ebins>e)[0][0] - 1
pseudo[ibin]+=line['epsilon']
linelist_equ = linelist_equ[linelist_equ['epsilon'] > MinEpsilon]
# now do the same for the continuum:
cont_rrc = numpy.zeros(len(ebins)-1, dtype=float)
cont_2ph = numpy.zeros(len(ebins)-1, dtype=float)
cont_bre = numpy.zeros(len(ebins)-1, dtype=float)
for z1 in list(contlist.keys()):
cont_rrc += contlist[z1]['rrc']
cont_2ph += contlist[z1]['twophot']
cont_bre += contlist[z1]['brems']
for z1 in list(pseudolist.keys()):
pseudo += pseudolist[z1]
# compress things
# ret = {}
# ret['lines']= linelist_equ
# ret['cont_rrc'] = cont_rrc
# ret['cont_2ph'] = cont_2ph
# ret['cont_bre'] = cont_bre
# ret['pseudo'] = pseudo
# ret['ebins'] = ebins
# ret['Z'] = Z
# create the real outputs
econt, cont = compress_continuum(ebins, cont_rrc+cont_2ph+cont_bre, const.TOLERANCE, minval=1e-38)
epseudo, pseudo = compress_continuum(ebins, pseudo, const.TOLERANCE, minval=1e-38)
# make adjustments for "zero"
pseudo[pseudo < 0] = 0
igood = numpy.ones(len(pseudo), dtype=bool)
for i in range(1, len(pseudo)-1):
if ((pseudo[i]<=0) &\
(pseudo[i-1]<=0) &\
(pseudo[i+1]<=0)):
igood[i] = False
pseudo = pseudo[igood]
epseudo = epseudo[igood]
ret={}
ret['lines']=linelist_equ
maxnpseudo = len(pseudo)
maxncont = len(cont)
ret['cont'] = numpy.zeros(1, dtype=generate_datatypes('continuum', npseudo=maxnpseudo, ncontinuum=maxncont))
ret['cont']['Z'][0] = Z
ret['cont']['rmJ'][0] = 0
ret['cont']['N_Cont'][0] = maxncont
ret['cont']['E_Cont'][0][:maxncont] = econt
ret['cont']['Continuum'][0][:maxncont] = cont
ret['cont']['N_Pseudo'][0] = maxnpseudo
ret['cont']['E_Pseudo'][0][:maxnpseudo] = epseudo
ret['cont']['Pseudo'][0][:maxnpseudo] = pseudo
return ret
[docs]
def gather_rates(Z, z1, te, dens, datacache=False, settings=False,\
do_la=True, do_ai=True, do_ec=True, do_pc=True,\
do_ir=True):
"""
fetch the rates for all the levels of Z, z1
Parameters
----------
Z: int
The nuclear charge of the element
z1 : int
ion charge +1
te : float
temperture (Kelvin)
dens: float
electron density (cm^-3)
settings : dict
See description in atomdb.get_data
datacache : dict
Used for caching the data. See description in atomdb.get_data
Returns
-------
up: numpy.array(float)
Initial level of each transition
lo: numpy.array(float)
Final level of each transition
rate: numpy.array(float)
Rate for each transition (in s-1)
"""
Te_arr, dummy = util.make_vec(te)
print("Starting Gather Rates Z=%i, z1=%iat %s"%(Z, z1,time.asctime()))
lvdat = atomdb.get_data(Z, z1, 'LV', datacache=datacache, \
settings = settings)
nlev = len(lvdat[1].data)
diagterms = numpy.zeros(nlev)
# get the LA data:
if 'AAUT_TOT' in lvdat[1].data.names:
has_sum_lv = True
else:
has_sum_lv = False
print("has_sum_lv is ", has_sum_lv)
laup = numpy.zeros(0, dtype=int)
lalo = numpy.zeros(0, dtype=int)
larate = numpy.zeros(0, dtype=float)
if do_la:
print("Starting Gather Rates do_la at %s"%(time.asctime()))
t1=time.time()
if has_sum_lv:
diagterms+=lvdat[1].data['ARAD_TOT']
ladat = atomdb.get_data(Z, z1, 'LA', datacache=datacache, \
settings = settings)
if ladat != False:
laup = ladat[1].data['Upper_Lev'] - 1
lalo = ladat[1].data['Lower_Lev'] - 1
larate = ladat[1].data['Einstein_A']
if not(has_sum_lv):
# diagterms[laup] +=larate
for i in range(len(laup)):
diagterms[laup[i]] +=larate[i]
t2 = time.time()
print("Finished Gather Rates do_la at %s: took %f seconds"%(time.asctime(),t2-t1))
# create dummy results
# get the AI data:
aiup = numpy.zeros(0, dtype=int)
ailo = numpy.zeros(0, dtype=int)
airate = numpy.zeros(0, dtype=float)
if do_ai:
if has_sum_lv:
diagterms+= lvdat[1].data['AAUT_TOT']
print("Starting Gather Rates do_ai at %s"%(time.asctime()))
t1=time.time()
aidat = atomdb.get_data(Z, z1, 'AI', datacache=datacache, \
settings = settings)
if aidat != False:
aiup = aidat[1].data['Level_Init'][:] - 1
ailo = numpy.zeros(len(aidat[1].data), dtype=int)
airate = aidat[1].data['Auto_Rate']
if not(has_sum_lv):
for i in range(len(aiup)):
diagterms[aiup[i]] +=airate[i]
t2=time.time()
print("Finished Gather Rates do_ai at %s: took %f seconds"%(time.asctime(),t2-t1))
# get the EC data:
ecup = numpy.zeros(0, dtype=int)
eclo = numpy.zeros(0, dtype=int)
ecrate = numpy.zeros(0, dtype=float)
if do_ec:
print("Starting Gather Rates do_ec at %s"%(time.asctime()))
t1 = time.time()
ecdat = atomdb.get_data(Z, z1, 'EC', datacache=datacache, \
settings = settings)
if ecdat != False:
lvdat = atomdb.get_data(Z, z1, 'LV', datacache=datacache, \
settings = settings)
ecup = numpy.zeros(len(ecdat[1].data), dtype=int)
eclo = numpy.zeros(len(ecdat[1].data), dtype=int)
ecrate = numpy.zeros(len(ecdat[1].data), dtype=float)
decup = numpy.zeros(len(ecdat[1].data), dtype=int)
declo = numpy.zeros(len(ecdat[1].data), dtype=int)
decrate = numpy.zeros(len(ecdat[1].data), dtype=float)
idex = 0
# need to loop and calculate each result
# Te_arr = numpy.array([te])
deglarr = lvdat[1].data['LEV_DEG'][ecdat[1].data['lower_lev']-1]
deguarr = lvdat[1].data['LEV_DEG'][ecdat[1].data['upper_lev']-1]
deltaearr = lvdat[1].data['ENERGY'][ecdat[1].data['upper_lev']-1]-\
lvdat[1].data['ENERGY'][ecdat[1].data['lower_lev']-1]
ladat = atomdb.get_data(Z, z1, 'LA', datacache=datacache, \
settings = settings)
for i in range(len(ecdat[1].data)):
# check if we need to swap things
if deltaearr[i] >=0:
degl = deglarr[i]
degu = deguarr[i]
ilo = ecdat[1].data['Lower_Lev'][i] - 1
iup = ecdat[1].data['Upper_Lev'][i] - 1
else:
degl = deguarr[i]
degu = deglarr[i]
ilo = ecdat[1].data['Upper_Lev'][i] - 1
iup = ecdat[1].data['Lower_Lev'][i] - 1
deltaearr[i] *= -1
exc,dex, tmp = atomdb._calc_maxwell_rates(ecdat[1].data['coeff_type'][i],\
ecdat[1].data['min_temp'][i],\
ecdat[1].data['max_temp'][i],\
ecdat[1].data['temperature'][i],\
ecdat[1].data['effcollstrpar'][i],\
deltaearr[i]/1e3, Te_arr, Z, \
degl, degu,\
force_extrap=True,\
levdat=lvdat,ladat=ladat, \
lolev=ilo+1, uplev=iup+1)
ecup[i] = ilo
eclo[i] = iup
ecrate[i] = exc*dens
if dex > 0:
decup[idex] = iup
declo[idex] = ilo
decrate[idex] = dex*dens
idex += 1
# now merge the de-excitation data
decup = decup[:idex]
declo = declo[:idex]
decrate = decrate[:idex]
ecup = numpy.append(ecup, decup)
eclo = numpy.append(eclo, declo)
ecrate = numpy.append(ecrate, decrate)
# create dummy results
# if not(has_sum_lv):
for i in range(len(ecup)):
diagterms[ecup[i]] +=ecrate[i]
t2 = time.time()
print("Finished Gather Rates do_ec at %s: took %f seconds"%(time.asctime(),t2-t1))
# get the PC data:
pcup = numpy.zeros(0, dtype=int)
pclo = numpy.zeros(0, dtype=int)
pcrate = numpy.zeros(0, dtype=float)
if do_pc:
print("Starting Gather Rates do_pc at %s"%(time.asctime()))
t1 = time.time()
pcdat = atomdb.get_data(Z, z1, 'PC', datacache=datacache, \
settings = settings)
if pcdat != False:
pcup = numpy.zeros(len(pcdat[1].data), dtype=int)
pclo = numpy.zeros(len(pcdat[1].data), dtype=int)
pcrate = numpy.zeros(len(pcdat[1].data), dtype=float)
dpcup = numpy.zeros(len(pcdat[1].data), dtype=int)
dpclo = numpy.zeros(len(pcdat[1].data), dtype=int)
dpcrate = numpy.zeros(len(pcdat[1].data), dtype=float)
idex = 0
# need to loop and calculate each result
# Te_arr = numpy.array([te])
deglarr = lvdat[1].data['LEV_DEG'][pcdat[1].data['lower_lev']-1]
deguarr = lvdat[1].data['LEV_DEG'][pcdat[1].data['upper_lev']-1]
deltaearr = lvdat[1].data['ENERGY'][pcdat[1].data['upper_lev']-1]-\
lvdat[1].data['ENERGY'][pcdat[1].data['lower_lev']-1]
for i in range(len(pcdat[1].data)):
exc,dex, tmp = atomdb._calc_maxwell_rates(pcdat[1].data['coeff_type'][i],\
pcdat[1].data['min_temp'][i],\
pcdat[1].data['max_temp'][i],\
pcdat[1].data['temperature'][i],\
pcdat[1].data['effcollstrpar'][i],\
deltaearr[i]/1e3, Te_arr, Z, \
deglarr[i], deguarr[i],\
force_extrap=True)
pcup[i] = pcdat[1].data['Lower_Lev'][i] - 1
pclo[i] = pcdat[1].data['Upper_Lev'][i] - 1
pcrate[i] = exc*dens
if dex > 0:
dpcup[idex] = pcdat[1].data['Upper_Lev'][i] - 1
dpclo[idex] = pcdat[1].data['Lower_Lev'][i] - 1
dpcrate[idex] = dex*dens
idex += 1
# now merge the de-excitation data
dpcup = dpcup[:idex]
dpclo = dpclo[:idex]
dpcrate = dpcrate[:idex]
pcup = numpy.append(pcup, dpcup)
pclo = numpy.append(pclo, dpclo)
pcrate = numpy.append(pcrate, dpcrate)
# if not(has_sum_lv):
for i in range(len(pcup)):
diagterms[pcup[i]] +=pcrate[i]
t2=time.time()
print("Finished Gather Rates do_pc at %s: took %f seconds"%(time.asctime(),t2-t1))
# get the IR data for colln ionization:
irup = numpy.zeros(0, dtype=int)
irlo = numpy.zeros(0, dtype=int)
irrate = numpy.zeros(0, dtype=float)
if do_ir:
print("Starting Gather Rates do_ir at %s"%(time.asctime()))
t1 = time.time()
irdat = atomdb.get_data(Z, z1, 'IR', datacache=datacache, \
settings = settings)
if irdat != False:
irup = numpy.zeros(len(irdat[1].data), dtype=int)
irlo = numpy.zeros(len(irdat[1].data), dtype=int)
irrate = numpy.zeros(len(irdat[1].data), dtype=float)
iir = 0
# need to loop and calculate each result
Te_arr = numpy.array(te)
irtmp = irdat[1].data[(irdat[1].data['TR_TYPE']=='XI') | \
(irdat[1].data['TR_TYPE']=='CI')]
deglarr = lvdat[1].data['LEV_DEG'][irtmp['level_init']-1]
lvdatp1 = atomdb.get_data(Z,z1+1,'LV', settings=settings, datacache=datacache)
if not (lvdatp1==False) :
deglarr = lvdatp1[1].data['LEV_DEG'][irtmp['level_final']-1]
else:
deglarr = numpy.ones(len(irtmp['level_final']), dtype=int)
ionpot = float(irdat[1].header['ionpot'])
for i in range(len(irdat[1].data)):
if not irdat[1].data['TR_TYPE'][i] in ['XI','CI']:
continue
else:
rate=atomdb.get_maxwell_rate(Te_arr, irdat, i, lvdat, \
lvdatap1=lvdatp1, ionpot=ionpot, \
exconly=True)
irup[iir] = irdat[1].data['level_init'][i] - 1
irlo[iir] = 0
irrate[iir] = rate*dens
iir += 1
# now merge the de-excitation data
irup = irup[:iir]
irlo = irlo[:iir]
irrate = irrate[:iir]
# for i in range(len(irup)):
# print "IR %i %i = %e"%(irup[i], irlo[i], irrate[i])
# if not(has_sum_lv):
for i in range(len(irup)):
diagterms[irup[i]] +=irrate[i]
t2=time.time()
print("Finished Gather Rates do_ir at %s: took %f seconds"%(time.asctime(),t2-t1))
print("Gather Rates: Starting combining rates into arrays at %s"%(time.asctime()))
tmp={}
tmp['up']={}
tmp['lo']={}
tmp['rate']={}
tmp['up']['ec'] = ecup
tmp['lo']['ec'] = eclo
tmp['rate']['ec'] = ecrate
tmp['up']['pc'] = pcup
tmp['lo']['pc'] = pclo
tmp['rate']['pc'] = pcrate
tmp['up']['ir'] = irup
tmp['lo']['ir'] = irlo
tmp['rate']['ir'] = irrate
tmp['up']['ai'] = aiup
tmp['lo']['ai'] = ailo
tmp['rate']['ai'] = airate
tmp['up']['la'] = laup
tmp['lo']['la'] = lalo
tmp['rate']['la'] = larate
tmp['up']['diag'] = numpy.arange(nlev, dtype=int)
tmp['lo']['diag'] = numpy.arange(nlev, dtype=int)
tmp['rate']['diag'] = diagterms*-1
# pickle.dump(tmp, open('tmp_gather_%i_%i.pkl'%(Z,z1),'wb'))
t1=time.time()
up_out = numpy.append(laup, numpy.append(aiup, numpy.append(ecup, numpy.append(pcup, irup))))
lo_out = numpy.append(lalo, numpy.append(ailo, numpy.append(eclo, numpy.append(pclo, irlo))))
rate_out = numpy.append(larate, numpy.append(airate, numpy.append(ecrate, numpy.append(pcrate, irrate))))
up_out = numpy.append(up_out, numpy.arange(nlev, dtype=int))
lo_out = numpy.append(lo_out, numpy.arange(nlev, dtype=int))
rate_out = numpy.append(rate_out, diagterms*-1)
t2=time.time()
print("Gather Rates: Finished combining rates into arrays at %s: took %f seconds"%(time.asctime(), t2-t1))
print("Finished Gather Rates at %s"%(time.asctime()))
return up_out, lo_out, rate_out
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def generate_datatypes(dtype, npseudo=0, ncontinuum=0):
"""
returns the various data types needed by apec
Parameters
----------
dtype : string
One of "linetype", "cielinetype", "continuum"
npseudo : int (default=0)
Number of pseudocontinuum points for "continuum" type
ncontinuum : int (default=0)
Number of continuum points for "continuum" type
Returns
-------
numpy.dtype
The data dtype in question
"""
if dtype == 'linetype':
ret = numpy.dtype({'names':['lambda',\
'lambda_err',\
'epsilon',\
'epsilon_err',\
'element',\
'ion', \
'elem_drv',\
'ion_drv', \
'upperlev',\
'lowerlev'],\
'formats':[float,\
float,\
float,\
float,\
int,\
int,\
int,\
int,\
int,\
int]})
elif dtype =='linelist_cie':
ret = numpy.dtype({'names':['lambda',\
'lambda_err',\
'epsilon',\
'epsilon_err',\
'element',\
'ion', \
'upperlev',\
'lowerlev'],\
'formats':[float,\
float,\
float,\
float,\
int,\
int,\
int,\
int]})
elif dtype =='linelist_cie_spectrum':
ret = numpy.dtype({'names':['Lambda',\
'Lambda_Err',\
'Epsilon',\
'Epsilon_Err',\
'Element',\
'Ion', \
'UpperLev',\
'LowerLev'],\
'formats':[numpy.float32,\
numpy.float32,\
numpy.float32,\
numpy.float32,\
numpy.int32,\
numpy.int32,\
numpy.int32,\
numpy.int32]})
elif dtype =='linelist_nei_spectrum':
ret = numpy.dtype({'names':['Lambda',\
'Lambda_Err',\
'Epsilon',\
'Epsilon_Err',\
'Element',\
'Elem_drv',\
'Ion', \
'Ion_drv', \
'UpperLev',\
'LowerLev'],\
'formats':[numpy.float32,\
numpy.float32,\
numpy.float32,\
numpy.float32,\
numpy.int32,\
numpy.int32,\
numpy.int32,\
numpy.int32,\
numpy.int32,\
numpy.int32]})
elif dtype == 'linetype_cap':
ret = numpy.dtype({'names':['Lambda',\
'Lambda_Err',\
'Epsilon',\
'Epsilon_Err',\
'Element',\
'Ion', \
'Elem_Drv',\
'Ion_Drv', \
'UpperLev',\
'LowerLev'],\
'formats':[float,\
float,\
float,\
float,\
int,\
int,\
int,\
int,\
int,\
int]})
elif dtype == 'linelist_cie_cap':
ret = numpy.dtype({'names':['Lambda',\
'Lambda_Err',\
'Epsilon',\
'Epsilon_Err',\
'Element',\
'Ion', \
'UpperLev',\
'LowerLev'],\
'formats':[float,\
float,\
float,\
float,\
int,\
int,\
int,\
int]})
elif dtype == 'continuum':
if ncontinuum==0:
ncontinuum+=1
if npseudo==0:
npseudo+=1
ret = numpy.dtype({'names':['Z','rmJ','N_Cont','E_Cont','Continuum','Cont_Err','N_Pseudo','E_Pseudo','Pseudo','Pseudo_Err'],\
'formats':[int, int, \
int, (float, ncontinuum), (float, ncontinuum),(float, ncontinuum),\
int, (float, npseudo), (float, npseudo),(float, npseudo)]})
elif dtype == 'lineparams':
ret = numpy.dtype({'names':['kT','EDensity','Time','Nelement','Nline'],\
'formats':[float, float, float, int, int]})
elif dtype == 'cocoparams':
ret = numpy.dtype({'names':['kT','EDensity','Time','NElement','NCont', 'NPseudo'],\
'formats':[float, float, float, int, int, int]})
elif dtype == 'ecdat':
# Electron collisional data
ret = numpy.dtype({'names':['lower_lev','upper_lev','coeff_type',\
'min_temp','max_temp', 'temperature', \
'effcollstrpar','inf_limit','reference'],\
'formats':[int, int, int, \
float, float, (float,20), \
(float,20), float, '|S20']})
else:
print("Unknown dtype %s in generate_datatypes"%(dtype))
return ret
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def solve_level_pop(init,final,rates,settings):
"""
Solve the level population
Parameters
----------
init : array(int)
The initial level for each transition
final : array(int)
The initial level for each transition
rates : array(float)
The rate for each transition
settings: dictionary
The settings read from the apec.par file by parse_par_file
Returns
-------
array(float)
The level population
"""
import scipy.sparse as sparse
from scipy.sparse.linalg import spsolve
# stuff = {}
# stuff['init'] = init
# stuff['final'] = final
# stuff['rates'] = rates
# pickle.dump(stuff, open('tmpA.pkl','wb'))
# pickle.dump(stuff, open('tmpB.pkl','wb'))
# zzz=raw_input()
# creat the arrays
nlev = max([max(init), max(final)])+1
#-- print "Starting Solve Level Pop at %s"%(time.asctime())
tstart = time.time()
k=numpy.where((final==0) & (init==0))[0]
if len(k)==1:
rates[k[0]] = 0.0
k=numpy.isfinite(rates)
rates=rates[k]
init=init[k]
final=final[k]
# print "first test"
# for i in range(len(final)):
# if final[i] == init[i]: print init[i],final[i], rates[i]
# print "end first test"
nlev_old = max([max(init), max(final)])+1
irate = numpy.where(rates>1e-40)[0]
# print irate
nlev = max([max(init[irate]), max(final[irate])])
# print "nlev =", nlev
irate = numpy.where((init<=nlev) & (final<=nlev))[0]
# print "keeping %i of %i rates"%(len(irate), len(init))
# print len(irate), irate[-1]
init=init[irate]
final=final[irate]
rates=rates[irate]
nlev = max([max(init), max(final)])+1
print("nlev = %i, nlev_old =%i"%(nlev, nlev_old))
# print "second test"
# for i in range(len(final)):
# if final[i] == init[i]: print init[i],rates[i]
# print "end second test"
if nlev <= const.NLEV_NOSPARSE:
print("Using regular solver.")
print("Starting generation of matrixA at %s"%(time.asctime()))
t1=time.time()
# convert to a regular solver
matrixA = numpy.zeros([nlev,nlev], dtype=float)
matrixB = numpy.zeros([nlev], dtype=float)
t0 = time.time()
t1=time.time()
for i in range(len(init)):
matrixA[final[i], init[i]] += rates[i]
t2 = time.time()
print("time differences: %f vs %f seconds"%(t1-t0, t2-t1))
# popn conservation
matrixB[0] = 1.0
matrixA[0,:] = 1.0
print("Starting check of diagonal terms at %s"%(time.asctime()))
# bug-u-fix
for i in range(1, len(matrixB)):
if matrixA[i,i] >= 0:
matrixA[i,i]=-1e10
print("ATieing level %i to ground with rate 1e10"%(i))
print("Finished check of diagonal terms at %s"%(time.asctime()))
# a = {}
# a['A'] = matrixA
# a['B'] = matrixB
# fname='dump_matrixa_1.pkl'
# ifile = 1
# while os.path.isfile(fname):
# ifile +=1
# fname = 'dump_matrixa_%i.pkl'%(ifile)
# pickle.dump(a, open(fname,'wb'))
# print "Wrote %s"%(fname)
# print "a['A'].shape", a['A'].shape
t2=time.time()
print("Finished generation of matrixA at %s: took %f seconds"%(time.asctime(), t2-t1))
print("Starting calling solver at %s"%(time.asctime()))
try:
popn = numpy.linalg.solve(matrixA, matrixB)
except numpy.linalg.linalg.LinAlgError:
raise
t3 = time.time()
print("Finished calling solver at %s: took %f seconds"%(time.asctime(), t3-t2))
else:
print("using sparse solver")
matrixA={}
matrixB = numpy.zeros(nlev, dtype=float)
# matrixA['init'] = numpy.append(init, init)
# matrixA['final'] = numpy.append(final, init)
# matrixA['rate'] = numpy.append(rates, -1.0*rates)
matrixA['init'] = init
matrixA['final'] = final
matrixA['rate'] = rates
maxlev = max(matrixA['final'][matrixA['rate']>1e-40])
k = numpy.where(matrixA['rate']<1e-40)[0]
# filter for the ground state levels
i = matrixA['final']>0
matrixA['final'] = matrixA['final'][i]
matrixA['init'] = matrixA['init'][i]
matrixA['rate'] = matrixA['rate'][i]
# add in continuity eqn
i = (matrixA['final']<=maxlev) & (matrixA['init']<=maxlev)
matrixA['final'] = matrixA['final'][i]
matrixA['init'] = matrixA['init'][i]
matrixA['rate'] = matrixA['rate'][i]
matrixA['final'] = numpy.append(matrixA['final'], \
numpy.zeros(maxlev, dtype=int))
matrixA['init'] = numpy.append(matrixA['init'], \
numpy.arange(maxlev, dtype=int))
matrixA['rate'] = numpy.append(matrixA['rate'], \
numpy.ones(maxlev, dtype=float))
# A = sparse.coo_matrix((matrixA['rate'],\
# (matrixA['final'],matrixA['init'])), \
# shape=(maxlev+1,maxlev+1)).tocsr()
A = numpy.zeros([maxlev+1,maxlev+1])
for i in range(len(matrixA['final'])):
A[matrixA['final'][i], matrixA['init'][i]]+=matrixA['rate'][i]
# print "A.shape", A.shape, " maxlev", maxlev
# for i in range(len(matrixA['init'])):
# A[matrixA['final'][i], matrixA['init'][i]] += matrixA['rate'][i]
# if matrixA['final'][i] == matrixA['init'][i]: print matrixA['init'][i],matrixA['rate'][i]
matrixB[0] = 1.0
A[0,:] = 1.0
hasdat = numpy.zeros(len(matrixB), dtype=bool)
for i in range(1,maxlev+1):
if A[i,i]>=0.0:
A[i,i] = -1e10
print("BTieing level %i to ground with rate 1e10"%(i))
matrixB[0] = 1.0
# tmp={}
#tmp['A'] = matrixA
#tmp['B'] = matrixB
# pickle.dump(tmp, open('dump_tmp.pkl','wb'))
popn = numpy.zeros(nlev_old)
matrixA=0
popn[:maxlev+1] = numpy.linalg.solve(A, matrixB[:maxlev+1])
# popn[:maxlev+1] = spsolve(A, matrixB[:maxlev+1])
# sparse solver sometimes returns small negative pop. set to 0.
popn[popn<0] = 0.0
# check for low population levels which are not
# correctly handled by the sparse solver
# anything under 1e-28 is suspect.
# tocheck = numpy.where((popn>= 0) &\
# (popn<= 1e-28))[0]
# if len(tocheck) > 0:
# for i in tocheck:
# itot_in = numpy.where((matrixA['init']==0) &\
# (matrixA['final']==i))[0]
#
# itot_out = numpy.where((matrixA['init']==i) &\
# (matrixA['final']==i))[0]
#
# tot_in = -sum(matrixA['rate'][itot_in])
# tot_out = sum(matrixA['rate'][itot_out])
# p =tot_in/tot_out
#
# if (popn[i] < 1e-30):
# if tot_in < 1e-21:
# popn[i] = p
# elif (popn[i] < 1e-28):
# if (tot_in < 1e-30):
# popn[i] = p
tfinish=time.time()
print("Finished Solve Level Pop at %s: took %i seconds"%(time.asctime(), tfinish-tstart))
return popn
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def do_lines(Z, z1, lev_pop, N_e, datacache=False, settings=False, z1_drv_in=-1):
"""
Convert level populations into line lists
Parameters
----------
Z: int
The nuclear charge of the element
z1 : int
Ion charge +1 of ion (e.g. 6 for C VI)
lev_pop : array(float)
The level population for the ion. Should already have elemental abundance
and ion fraction multiplied in.
N_e : float
Electron Density (cm^-3)
datacache : dict
Used for caching the data. See description in atomdb.get_data
settings : dict
See description in atomdb.get_data
z1_drv_in : int
the driving ion for this calculation, if not z1 (defaults to z1)
Returns
-------
linelist: numpy.dtype(linetype)
The list of lines and their emissivities. see generate_datatypes
twophot: array(float)
The two-photon continuum on the grid specified by the settings
If settings['TwoPhoton'] is False, then returns a grid of zeros.
"""
print("starting do_lines at %s"%(time.asctime()))
tstart=time.time()
ladat = atomdb.get_data(Z,z1,'LA', datacache=datacache, settings=settings)
lvdat = atomdb.get_data(Z,z1,'LV', datacache=datacache, settings=settings)
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
twoph = numpy.zeros(settings['NumGrid'], dtype=float)
linelist = numpy.zeros(len(ladat[1].data), \
dtype= generate_datatypes('linetype'))
goodlines = numpy.ones(len(linelist), dtype=bool)
if z1_drv_in < 0:
z1_drv = z1
else:
z1_drv = z1_drv_in
# now do this
linelist['epsilon'] = ladat[1].data.field('einstein_a') * \
lev_pop[ladat[1].data.field('upper_lev')-1]/N_e
linelist['lambda'] = ladat[1].data.field('wavelen')
igood = numpy.isfinite(ladat[1].data['wave_obs'])
if sum(igood) > 0 :
igood = numpy.where(igood==True)[0]
igood = igood[numpy.where(ladat[1].data['wave_obs'][igood]>0)[0]]
if len(igood) >0:
linelist['lambda'][igood]=ladat[1].data['wave_obs'][igood]
# for iline, line in enumerate(ladat[1].data):
# if numpy.isfinite(line['wave_obs']):
# if line['wave_obs'] > 0:
# linelist[iline]['lambda'] = line['wave_obs']
# else:
# linelist[iline]['lambda'] = line['wavelen']
# else:
# linelist[iline]['lambda'] = line['wavelen']
linelist['lambda_err'] = numpy.nan
linelist['epsilon_err'] = numpy.nan
linelist['element'] = Z
linelist['elem_drv'] = Z
linelist['ion'] = z1
linelist['ion_drv']= z1_drv
linelist['upperlev'] = ladat[1].data['upper_lev']
linelist['lowerlev'] = ladat[1].data['lower_lev']
# I have a linelist. Yay.
t1=time.time()
print("finished making linelist at %s: took %f seconds"%(time.asctime(), t1-tstart))
# now check for 2 photon transitions
print("starting check for 2 photon transitions at %s"%(time.asctime()))
if (Z-z1==1):
# He-like:
nup = lvdat[1].data['n_quan'][linelist['upperlev']-1]
lup = lvdat[1].data['l_quan'][linelist['upperlev']-1]
degup = lvdat[1].data['lev_deg'][linelist['upperlev']-1]
deggd = lvdat[1].data['lev_deg'][0]
ila = numpy.where((nup==2) &\
(lup==0) &\
(degup==deggd) &\
(linelist['lowerlev']==1))[0]
if len(ila)>0:
# do the 2 photon fun
twoph = numpy.zeros(len(ebins)-1, dtype=float)
if settings['TwoPhoton']:
# if len(ila) > 1:
for iila in ila:
tmp2ph={}
tmp2ph['lambda'] = linelist['lambda'][iila]
tmp2ph['lev_pop'] = lev_pop[ladat[1].data.field('upper_lev')[iila]-1]/N_e
tmp2ph['einstein_a'] = ladat[1].data.field('einstein_a')[iila]
twoph += atomdb.calc_two_phot(tmp2ph['lambda'],\
tmp2ph['einstein_a'],\
tmp2ph['lev_pop'], ebins)
goodlines[ila]=False
elif (Z-z1==0):
# H-like
nup = lvdat[1].data['n_quan'][linelist['upperlev']-1]
lup = lvdat[1].data['l_quan'][linelist['upperlev']-1]
degup = lvdat[1].data['lev_deg'][linelist['upperlev']-1]
deggd = lvdat[1].data['lev_deg'][0]
ila = numpy.where((nup==2) &\
(lup==0) &\
(degup==deggd)&\
(linelist['lowerlev']==1))[0]
if len(ila)>0:
twoph = numpy.zeros(len(ebins)-1, dtype=float)
if settings['TwoPhoton']:
for iila in ila:
tmp2ph={}
tmp2ph['lambda'] = linelist['lambda'][iila]
tmp2ph['lev_pop'] = lev_pop[ladat[1].data.field('upper_lev')[iila]-1]/N_e
tmp2ph['einstein_a'] = ladat[1].data.field('einstein_a')[iila]
twoph += atomdb.calc_two_phot(tmp2ph['lambda'],\
tmp2ph['einstein_a'],\
tmp2ph['lev_pop'], ebins)
goodlines[ila]=False
elif (Z-z1==3):
nup = lvdat[1].data['n_quan'][linelist['upperlev']-1]
lup = lvdat[1].data['l_quan'][linelist['upperlev']-1]
degup = lvdat[1].data['lev_deg'][linelist['upperlev']-1]
deggd = lvdat[1].data['lev_deg'][0]
ila = numpy.where((linelist['upperlev']==2) &\
(linelist['lowerlev']==1))[0]
if len(ila)>0:
twoph = numpy.zeros(len(ebins)-1, dtype=float)
if settings['TwoPhoton']:
for iila in ila:
tmp2ph={}
tmp2ph['lambda'] = linelist['lambda'][iila]
tmp2ph['lev_pop'] = lev_pop[ladat[1].data.field('upper_lev')[iila]-1]/N_e
tmp2ph['einstein_a'] = ladat[1].data.field('einstein_a')[iila]
twoph += atomdb.calc_two_phot(tmp2ph['lambda'],\
tmp2ph['einstein_a'],\
tmp2ph['lev_pop'], ebins)
goodlines[ila]=False
t2=time.time()
print("finished checking two photon transitions at %s: took %f seconds"%(time.asctime(), t2-t1))
linelist = linelist[goodlines]
tfinish=time.time()
print("finished do_lines at %s, took %f seconds"%(time.asctime(), tfinish-tstart))
return linelist, twoph
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_satellite(Z, z1, T, datacache=False, settings=False):
"""
Calcaulate DR satellite lines
Parameters
----------
Z: int
The nuclear charge of the element Z: int
z1 : int
Recombined Ion charge +1 of ion (e.g. 5 for C VI -> C V)
te: float
The electron temperature (K)
settings: dictionary
The settings read from the apec.par file by parse_par_file
Returns
-------
array(linelist)
List of DR lines
array(levlistin)
Rates into each lower level, driven by DR
"""
# recombining ion charge
z1_drv=z1+1
drdat = atomdb.get_data(Z,z1_drv,'DR', settings=settings, datacache=datacache)
lvdat = atomdb.get_data(Z,z1_drv,'LV', settings=settings, datacache=datacache)
lvdatrec = atomdb.get_data(Z,z1,'LV', settings=settings, datacache=datacache)
# pseudocont = numpy.zeros(settings['NumGrid']-1, dtype=float)
# print "start DR"
if drdat==False:
linelist = numpy.zeros(0,dtype=generate_datatypes('linetype'))
lev_rates_in = 0.0
else:
if not(lvdatrec):
# no level data for recombined ion.
lomax = max(drdat[1].data['lower_lev'])
lev_rates_in=numpy.zeros(lomax+1, dtype=float)
else:
lev_rates_in = numpy.zeros(len(lvdatrec[1].data), dtype=float)
linelist = numpy.zeros(len(drdat[1].data),dtype=generate_datatypes('linetype'))
kT = const.KBOLTZ*T
for iline in range(len(linelist)):
epsilson = 0.0
ll = drdat[1].data['lower_lev'][iline]
lu = drdat[1].data['upper_lev'][iline]
lam = drdat[1].data['wavelen'][iline]
if (numpy.isfinite(drdat[1].data['wave_obs'][iline]) &\
(drdat[1].data['wave_obs'][iline]>0.0)):
lam = drdat[1].data['wave_obs'][iline]
lamerr = drdat[1].data['wave_err'][iline]
e_excite = drdat[1].data['E_excite'][iline]
if drdat[1].data['dr_type'][iline]==const.ROMANIK:
q_exc = drdat[1].data['SatelInt'][iline]
epsilon = q_exc*numpy.exp(-e_excite/kT)/(T**1.5)
elif drdat[1].data['dr_type'][iline]==const.SAFRANOVA:
try:
gl = lvdat[1].data['lev_deg'][0]
except TypeError:
if lvdat==False:
gl = 1
else: raise
q_exc = drdat[1].data['SatelInt'][iline]
epsilon = const.SAF_COEFF* (const.RYDBERG/kT)**1.5 * \
(q_exc/gl) * numpy.exp(-e_excite/kT)
else:
print("Error in calc_satellite: unknown DR type %i"%\
(drdat[1].data['type'][iline]))
epsilon = numpy.nan
# if there is a translation table from DR SATELLITE to AtomDB levels, use it.
if len(drdat) > 2:
itmp = drdat[2].data[drdat[2].data['DRLEVID']==ll]['APEDID']
if len(itmp) ==1:
if itmp[0] != 0:
ll = itmp[0]
# if ll >= len(lev_rates_in):
# print("warning: DR satellite line recombining into non existant level %i of ion Z=%i, z1=%i"%\
# (ll, Z, z1))
# else:
# lev_rates_in[ll-1] += epsilon
linelist[iline]['lambda'] = lam
linelist[iline]['lambda_err'] = lamerr
linelist[iline]['epsilon'] = epsilon
linelist[iline]['epsilon_err'] = numpy.nan
linelist[iline]['element'] = Z
linelist[iline]['ion'] = z1
linelist[iline]['elem_drv'] = Z
linelist[iline]['ion_drv'] = z1_drv
linelist[iline]['upperlev'] = lu
linelist[iline]['lowerlev'] = ll
linelist = linelist[numpy.where(linelist['epsilon'] > 0)[0]]
return linelist, lev_rates_in
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_recomb_popn(levpop, Z, z1, z1_drv,T, dens, drlevrates, rrlevrates,\
settings=False, datacache=False, dronly=False,\
rronly=False):
"""
Calculate the level population of a recombined ion
Parameters
----------
levpop: array(float)
Level populations, already taking into account elemental abundance
and ion fraction of z1_drv
Z: int
z1: int
z1_drv: int
T: electron temperature (K)
dens: electron density (cm^-3)
drlevrates: array(float)
Rates into each level from DR calculations
rrlevrates: array(float)
Rates into each level from RR calculations
Returns
-------
array(float)
Level population
"""
import scipy.sparse as sparse
from scipy.sparse.linalg import spsolve
# levpop at this point should alread have the corrected abundance in
# there
lvdat = atomdb.get_data(Z,z1,'LV', settings=settings, datacache=datacache)
if not lvdat:
nlev = 1
levpop = numpy.zeros(1, dtype=float)
return levpop
nlev = len(lvdat[1].data)
Tarr, dummy = util.make_vec(T)
if nlev > const.NLEV_NOSPARSE:
print("using sparse solver for recomb")
# sort the levels
aidat = atomdb.get_data(Z,z1,'AI', settings=settings, datacache=datacache)
if aidat:
ailev = numpy.array(util.unique(aidat[1].data['level_init']))-1
nailev = len(ailev)
isbound = numpy.ones(nlev, dtype=bool)
isbound[ailev]=False
else:
nailev = 0
isbound = numpy.ones(nlev, dtype=bool)
#sortdat = sort_levels(isbound, nlev, nailev, Z, z1, settings)
recombrate = numpy.zeros(nlev, dtype=float)
# get the recomb data
irdat = atomdb.get_data(Z, z1, 'IR', settings=settings, datacache=datacache)
for iir, ir in enumerate(irdat[1].data):
# check we have the right data types
if ir['TR_TYPE'] in ['RR','DR','XR']:
recrate = atomdb.get_maxwell_rate(Tarr, irdat, iir, lvdat)*dens
if not (numpy.isfinite(recrate)):
print("iir=%i, recrate is not finite!"%(iir))
else:
recombrate[ir['level_final']-1] += recrate*levpop[ir['level_init']-1]
maxlev = numpy.where(recombrate > 0)[0]
if len(maxlev) == 0: # so recombrate has all the influxes
return numpy.zeros(nlev)
maxlev=maxlev[-1]
matrixB = recombrate
# print "Sum recomb rates in level>1:", sum(matrixB[1:])
matrixA = {}
matrixA['init'], matrixA['final'], matrixA['rate']=\
gather_rates(Z, z1, T, dens, datacache=datacache, settings=settings,\
do_la=True, do_ai=False, do_ec=False, do_pc=False,\
do_ir=False)
# matrixA['init'] = numpy.array(ladat[1].data['upper_lev'], dtype=int)-1
# matrixA['final'] = numpy.array(ladat[1].data['lower_lev'], dtype=int)-1
# matrixA['rate'] = numpy.array(ladat[1].data['einstein_a'], dtype=float)
### FIXME HERE
matrixB *= -1
nlev = len(matrixB)
# print "recombination rates"
# for i in range(len(matrixB)):
# print i, matrixB[i]
# no recombination to speak of. Return zeros
# append extras onto matrixA:
# matrixA['final'] = numpy.append(matrixA['final'],matrixA['init'])
# matrixA['init'] = numpy.append(matrixA['init'],matrixA['init'])
# matrixA['rate'] = numpy.append(matrixA['rate'],-1*matrixA['rate'])
# remove all matrix A from and to level 0
i = (matrixA['final']>0) & (matrixA['init']>0)
matrixA['final'] = matrixA['final'][i]
matrixA['init'] = matrixA['init'][i]
matrixA['rate'] = matrixA['rate'][i]
# remove all matrix A from and to high levels
i = (matrixA['final']<maxlev+1) & (matrixA['init']<maxlev+1)
matrixA['final'] = matrixA['final'][i]
matrixA['init'] = matrixA['init'][i]
matrixA['rate'] = matrixA['rate'][i]
# subtract 1 from the levels
matrixA['init']-= 1
matrixA['final']-= 1
# solve
# A = sparse.coo_matrix((matrixA['rate'],\
# (matrixA['final'],matrixA['init'])), \
# shape=(nlev-1,nlev-1)).tocsr()
A = numpy.zeros([maxlev,maxlev])
for i in range(len(matrixA['final'])):
A[matrixA['final'][i], matrixA['init'][i]]+=matrixA['rate'][i]
# A = numpy.linalg.solve((matrixA['rate'],\
# (matrixA['final'],matrixA['init'])), \
# shape=(nlev-1,nlev-1)).tocsr()
levpop_this = numpy.zeros(len(matrixB))
# for i in numpy.where(matrixB != 0)[0]:
# print i, matrixB[i]
if sum(matrixB[1:] < 0):
levpop_this[1:maxlev+1] = numpy.linalg.solve(A, matrixB[1:maxlev+1])
# levpop_this = calc_cascade(recombrate, Z, z1, isbound, sortdat, T, settings, noauto=True)
else:
print("using regular solver for recomb")
rrrecombrate = numpy.zeros(nlev, dtype=float)
drrecombrate = numpy.zeros(nlev, dtype=float)
irdat = atomdb.get_data(Z, z1, 'IR', settings=settings, datacache=datacache)
havedrrate=False
haverrrate=False
for iir, ir in enumerate(irdat[1].data):
# check we have the right data types
if ir['TR_TYPE'] in ['RR','XR']:
recrate = atomdb.get_maxwell_rate(Tarr, irdat, iir, lvdat)
rrrecombrate[ir['level_final']-1] += recrate*levpop[ir['level_init']-1]*dens
if ((ir['TR_TYPE'] in ['RR','XR']) & (ir['level_final']>1)):
haverrrate=True
if ir['TR_TYPE'] in ['DR','XD']:
recrate = atomdb.get_maxwell_rate(Tarr, irdat, iir, lvdat)
drrecombrate[ir['level_final']-1] += recrate*levpop[ir['level_init']-1]*dens
if ((ir['TR_TYPE'] in ['DR','XD']) & (ir['level_final']>1)):
havedrrate=True
if havedrrate:
tmpdrlevrates=0.0
else:
tmpdrlevrates=drlevrates
if haverrrate:
tmprrlevrates=0.0
else:
tmprrlevrates=rrlevrates
if isinstance(rrlevrates, numpy.ndarray):
sumrrlevrates = sum(rrlevrates)
else:
sumrrlevrates = 0.0
if isinstance(drlevrates, numpy.ndarray):
sumdrlevrates = sum(drlevrates)
else:
sumdrlevrates = 0.0
print("DR: sum from satellite lines: %e, sum from IR file: %e" %\
(sumdrlevrates, sum(drrecombrate)))
print("RR: sum from PI xsections: %e, sum from IR file: %e" %\
(sumrrlevrates, sum(rrrecombrate)))
matrixB = rrrecombrate+drrecombrate+tmpdrlevrates+tmprrlevrates
if dronly:
matrixB = drrecombrate+tmpdrlevrates
if rronly:
matrixB = rrrecombrate+tmprrlevrates
matrixA = numpy.zeros([nlev,nlev],dtype=float)
ladat = atomdb.get_data(Z, z1, 'LA', settings=settings, datacache=datacache)
matrixA_in = {}
matrixA_in['init'], matrixA_in['final'], matrixA_in['rate']=\
gather_rates(Z, z1, T, dens, datacache=datacache, settings=settings,\
do_la=True, do_ai=True, do_ec=False, do_pc=False,\
do_ir=False)
datacache={}
for i in range(len(matrixA_in['init'])):
matrixA[matrixA_in['final'][i], matrixA_in['init'][i]]+=matrixA_in['rate'][i]
# solve unless matrixB ==0
if sum(matrixB[1:])>0:
matrixB = -1*matrixB
levpop_this = calc_cascade_population(matrixA, matrixB)
else:
levpop_this = numpy.zeros(nlev)
#print("level population for recombination into Z=%i, z1=%i, z1_drv=%i, levpop=%i"%\
#(Z, z1, z1_drv,levpop))
for i in range(len(levpop_this)):
print(i, levpop_this[i])
return levpop_this
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_cascade_population(matrixA, matrixB):
# replace the first line in the matrix with the population
# conservation equation (sum of level pops = 1)
matrixB[0] = 1.0
matrixA[0,:] = 1.0
mb =matrixB[1:]
ma =matrixA[1:,1:]
# solve
# print('ma', ma)
# print('mb', mb)
# print('ma.shape', ma.shape)
# print('mb.shape', mb.shape)
try:
# print('TRY1')
popn = numpy.linalg.solve(ma,mb)
except numpy.linalg.linalg.LinAlgError:
print('TRY2')
if ma[0,0]==0.0:
print('hacking ma[0,0]=1.0')
ma[0,0]=-1.0
try:
popn = numpy.linalg.solve(ma,mb)
except numpy.linalg.linalg.LinAlgError:
print("failed again")
# look for levels with no way to ground. Put in a -1rate
for i in range(ma.shape[0]):
if ma[i,i]>= 0.0:
ma[i,i]=-1.0
try:
popn = numpy.linalg.solve(ma,mb)
except:
print('triple fail')
print(ma)
print(mb)
raise
else:
print("tie levels to ground")
# look for levels with no way to ground. Put in a -1rate
for i in range(ma.shape[0]):
if ma[i,i]>= 0.0:
ma[i,i]=-1.0
print("tieing level %i to ground"%(i+2))
try:
popn = numpy.linalg.solve(ma,mb)
except:
print('triple fail')
print(ma)
print(mb)
raise
#check
soln = numpy.allclose(numpy.dot(ma, popn), mb)
if soln==False:
print("ERROR Solving population matrix!")
popn=numpy.append(numpy.array([0.0]), popn)
return popn
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def calc_ioniz_popn(levpop, Z, z1, z1_drv,T, Ne, settings=False, \
datacache=False, do_xi=False):
"""
Calculate the level population due to ionization into the ion
Parameters
----------
levpop: array(float)
The level population of the parent ion. Should already have abundance
and ion fraction built in.
Z: int
z1: int
z1_drv: int
T: float
Ne: float
settings: dict
datacache: dict
do_xi: bool
Include collisional ionization
Returns
-------
levpop_out: array(float)
The level populations of the Z,z1 ion
"""
# levpop at this point should alread have the corrected abundance in
# there
# This calculates the population of levels of z1, given a previous
# ion's (z1-1) population of levpop.
import scipy.sparse as sparse
from scipy.sparse.linalg import spsolve
print("Starting calc_ioniz_popn at %s"%(time.asctime()))
lvdat = atomdb.get_data(Z,z1,'LV', settings=settings, datacache=datacache)
# if we have no lv data, ignore.
if not util.keyword_check(lvdat):
nlev = 1
return numpy.array([0.0])
nlev = len(lvdat[1].data)
# get populating rate from previous ion
aidat = atomdb.get_data(Z, z1-1, 'AI', settings=settings, datacache=datacache)
ionizrateai=numpy.zeros(nlev, dtype=float)
ionizrateir=numpy.zeros(nlev, dtype=float)
print("Starting calc_ioniz_popn aidat loop at %s"%(time.asctime()))
if aidat:
tmp_pop = levpop[aidat[1].data['level_init']-1]
for iai in range(len(aidat[1].data)):
ionizrateai[aidat[1].data['level_final'][iai]-1] += \
tmp_pop[iai]*aidat[1].data['auto_rate'][iai]
#aidat.close()
print("Finished calc_ioniz_popn aidat loop at %s"%(time.asctime()))
print("Starting calc_ioniz_popn xidat loop at %s"%(time.asctime()))
if do_xi:
irdat = atomdb.get_data(Z, z1-1, 'IR', settings=settings, datacache=datacache)
ionpot = float(irdat[1].header['ionpot'])
if z1 >1:
lvdatm1 = atomdb.get_data(Z, z1-1, 'LV', settings=settings, datacache=datacache)
# go through each excitation, have fun
for iir, ir in enumerate(irdat[1].data):
if ir['TR_TYPE'] in ['XI']:
Te = numpy.array([T])
ionrate=atomdb.get_maxwell_rate(Te, irdat, iir, lvdatm1, \
lvdatap1=lvdat, ionpot=ionpot)
ionizrateir[ir['level_final']-1] += levpop[ir['level_init']-1]*\
ionrate
#get_maxwell_rate(Te, colldata=False, index=-1, lvdata=False, Te_unit='K', \
# lvdatap1=False, ionpot = False, \
# force_extrap=False, silent=True,\
# finallev=False, initlev=False,\
# Z=-1, z1=-1, dtype=False, exconly=False,\
# datacache=False, settings=False):
# print
print("Finished calc_ioniz_popn xidat loop at %s"%(time.asctime()))
ionizrate=ionizrateir+ionizrateai
matrixB = ionizrate
# save some time if there is nothing to ionize.
if sum(matrixB[1:]) ==0:
levpop_this = numpy.zeros(len(matrixB))
return levpop_this
maxlev = numpy.where(matrixB > 1e-40)[0]
if len(maxlev)==0:
print("No significant ionization found")
popn = numpy.zeros(len(matrixB))
return popn
maxlev=maxlev[-1]
print("maxlev=", maxlev)
matrixA_in={}
matrixA_in['init'], matrixA_in['final'], matrixA_in['rate'] = \
gather_rates(Z, z1, T, Ne, datacache=datacache, settings=settings,\
do_la=True, do_ai=True, do_ec=False, do_pc=False,\
do_ir=False)
i = (matrixA_in['init']<=maxlev) & (matrixA_in['final']<=maxlev)
matrixA_in['init']=matrixA_in['init'][i]
matrixA_in['final']=matrixA_in['final'][i]
matrixA_in['rate']=matrixA_in['rate'][i]
# fix the rates
for i in range(len(matrixA_in['init'])):
if matrixA_in['init'][i]==matrixA_in['final'][i]:
if matrixA_in['rate'][i] >=0.0:
matrixA_in['rate'][i] -=1e10
print("CTieing level %i to ground with rate 1e10"%(i))
if (maxlev <= const.NLEV_NOSPARSE):
# convert to a regular solver
print("regular solver")
matrixA = numpy.zeros([maxlev+1,maxlev+1], dtype=float)
for i in range(len(matrixA_in['init'])):
matrixA[matrixA_in['final'][i], matrixA_in['init'][i]] += matrixA_in['rate'][i]
# matrixA[matrixA_in['init'][i], matrixA_in['init'][i]] -= matrixA_in['rate'][i]
# popn conservation
# matrixB[0] = 1.0
# matrixA[0,:] = 1.0
# bug-u-fix
for i in range(1, maxlev):
if matrixA[i,i] >= 0:
matrixA[i,i]=-1e10
print("FIXING matrixA[%i,%i] = -1.0"%(i,i))
popn = numpy.zeros(nlev)
matrixB*=-1
try:
popn[1:maxlev] = numpy.linalg.solve(matrixA[1:maxlev,1:maxlev], matrixB[1:maxlev])
except numpy.linalg.linalg.LinAlgError:
"EEK ERROR!"
raise
#if sum(matrixB[1:])<1e-40:
#levpop_this = numpy.zeros(len(matrixB))
#else:
#matrixB = -1*matrixB
#levpop_this = numpy.zeros(len(matrixB))
# check for zeros
#for iLev in range(1,len(matrixB)):
#if not(matrixA[iLev,iLev]< 0.0):
# matrixA[iLev,iLev]=-1e10
#levpop_this[1:] = numpy.linalg.solve(matrixA[1:,1:], matrixB[1:])
# calc_cascade_population(matrixA, matrixB)
else:
# add into sparse solver
print("Using sparse solver")
matrixA={}
matrixB *= -1
nlev = len(matrixB)
if sum(matrixB)>=0:
return numpy.zeros(len(matrixB))
# add in the reverse processes
# matrixA['init'] = numpy.append(matrixA['init'], matrixA['init'])
# matrixA['final'] = numpy.append(matrixA['final'], matrixA['init'])
# matrixA['rate'] = numpy.append(matrixA['rate'], -1*matrixA['rate'])
# remove ground level
i = (matrixA_in['init']>0) & (matrixA_in['final']>0)
matrixA['init'] = matrixA_in['init'][i]
matrixA['final'] = matrixA_in['final'][i]
matrixA['rate'] = matrixA_in['rate'][i]
i = (matrixA['init']<=maxlev+1) & (matrixA['final']<=maxlev+1)
matrixA['init'] = matrixA['init'][i]
matrixA['final'] = matrixA['final'][i]
matrixA['rate'] = matrixA['rate'][i]
# subtract 1 from the levels
matrixA['init']-=1
matrixA['final']-=1
# A = sparse.coo_matrix((matrixA['rate'],\
# (matrixA['final'],matrixA['init'])), \
# shape=(maxlev,maxlev)).tocsr()
A = numpy.zeros([maxlev, maxlev])
for i in range(len(matrixA['final'])):
A[matrixA['final'][i], matrixA['init'][i]]+=matrixA['rate'][i]
popn = numpy.zeros(len(matrixB))
# popn[1:maxlev+1] = spsolve(A, matrixB[1:maxlev+1])
popn[1:maxlev+1] = numpy.linalg.solve(A, matrixB[1:maxlev+1])
popn_bak = popn*1.0
popn[popn<0] = 0.0
# check for low population levels which are not
# correctly handled by the sparse solver
# anything under 1e-28 is suspect.
# tocheck = numpy.where((popn>= 1e-40) &\
# (popn<= 1e-28))[0]
#
# if len(tocheck) > 0:
# for i in tocheck:
# # lowest lying levels are generally OK.
# if i < 100: continue
# tot_in = matrixB[i-1]
#
# itot_out = numpy.where((matrixA['init']==i) &\
# (matrixA['final']==i))[0]
#
# #tot_in = -sum(matrixA['rate'][itot_in])
# tot_out = sum(matrixA['rate'][itot_out])
# p =tot_in/tot_out
#
# if (popn[i] < 1e-30):
# if tot_in < 1e-21:
# popn[i] = p
# elif (popn[i] < 1e-28):
# if (tot_in < 1e-30):
# popn[i] = p
#
#
#
print("level population for Z=%i, z1=%i, z1_drv=%i"%(Z,z1,z1_drv))
for i in range(len(popn)):
print("%6i %e"%(i, popn[i]))
return popn
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def run_apec_ion(settings, te, dens, Z, z1, ionfrac, abund):
"""
Run the APEC code using the settings provided for an individual ion.
Parameters
----------
settings: dictionary
The settings read from the apec.par file by parse_par_file
te: float
The electron temperature (K)
dens: float
The electron density (cm^-3)
Z: int
The nuclear charge of the element
z1: int
The ion charge +1 of the ion
ionfrac: float
The fractional abundance of this ion (between 0 and 1)
abund: float
The elemental abundance of the element (normalized to H)
Returns
-------
linelist : numpy array
List of line details and emissivities
continuum : array
Continuum emission in photons bin-1 s-1.
This is a 3-item dict, with "rrc", "twophot", "brems" entries
for each continuum source
pseudocont : array
Pseudo Continuum emission in photons bin-1 s-1
"""
# get the data.
z1_drv=z1*1
# First, get the elemental abundance
# make some compatilibty copies of entries in the settings file
settings['filemap'] = settings['FileMap']
settings['atomdbroot'] = os.path.expandvars('$ATOMDB')
# get the output energy bins for the continuum
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
## FIXME CUTOFF FOR MIN IONPOP
linelist = numpy.zeros(0, dtype=generate_datatypes('linetype'))
pseudo = numpy.zeros(settings['NumGrid'], dtype=float)
continuum = {}
continuum['brems'] = numpy.zeros(settings['NumGrid'], dtype=float)
continuum['twophot'] = numpy.zeros(settings['NumGrid'], dtype=float)
continuum['rrc'] = numpy.zeros(settings['NumGrid'], dtype=float)
## FIXME CUTOFF FOR MIN IONPOP
if ionfrac[z1_drv-1] < const.MIN_IONPOP:
print("OMITTING Z=%i, z1=%i as ionfrac %e is below threshold of %e"%\
(Z, z1_drv, ionfrac[z1_drv-1], const.MIN_IONPOP))
return linelist, continuum, pseudo
print("NOT OMITTING Z=%i, z1=%i as ionfrac %e is above threshold of %e"%\
(Z, z1_drv, ionfrac[z1_drv-1], const.MIN_IONPOP))
# set up the datacache
datacache = {}
# Find the number of levels
lvdat = atomdb.get_data(Z,z1,'LV', datacache=datacache, settings=settings)
linelist_exc = numpy.zeros(0, dtype= generate_datatypes('linetype'))
linelist_dr = numpy.zeros(0, dtype= generate_datatypes('linetype'))
linelist_rec = numpy.zeros(0, dtype= generate_datatypes('linetype'))
if lvdat!= False:
# find the number of levels
nlev = len(lvdat[1].data)
# check if we need to do any of the line-related calculations
if (settings['EmissionLines'] or settings['TwoPhoton']):
# gather all the level to level rates
print("calling gather_rates from run_apec_ion")
up, lo, rates = gather_rates(Z, z1, te, dens, datacache=datacache, settings=settings)
print("finished calling gather_rates from run_apec_ion")
# purge the datacache here, as this often requires heavy memory use
# datacache={}
# solve everything
print("calling solve_level_pop from run_apec_ion")
lev_pop = solve_level_pop(up,lo,rates, settings)
print("finished calling solve_level_pop from run_apec_ion")
# just in case, add zeros to lengthen the lev_pop appropriately
if len(lev_pop) < nlev:
lev_pop = numpy.append(lev_pop, numpy.zeros(nlev-len(lev_pop), dtype=float))
# fix any sub-zero level populations
lev_pop[lev_pop<0] = 0.0
# now we have the level populations, make a line list for each ion
# scale lev_pop by the ion and element abundance.
#print("lev_pop Z=%i, z1=%i,z1_drv=%i, abund*ionfrac=%e, sum(pop)=%e:"%(Z,z1, z1, abund*ionfrac[z1-1], sum(lev_pop)*abund*ionfrac[z1-1]))
lev_pop *= abund*ionfrac[z1-1]
# for i in range(len(lev_pop)):
# print(i, lev_pop[i])
print("calling do_lines from run_apec_ion")
linelist_exc, continuum['twophot'] = do_lines(Z, z1, lev_pop, dens, datacache=datacache, settings=settings, z1_drv_in=z1_drv)
print("finished calling do_lines from run_apec_ion")
print("Excitation Z=%i z1=%i z1_drv=%i created %i lines"%(Z, z1, z1_drv, len(linelist_exc)))
else:
# skipping the exact level calculation, fill it with zeros, ground state with 1.
lev_pop = numpy.zeros(nlev, dtype=float)
lev_pop[0] = 1.0*abund*ionfrac[z1-1]
else:
lev_pop=numpy.ones(1, dtype=float)*abund*ionfrac[z1-1]
print("lev_pop Z=%i, z1=%i,z1_drv=%i, abund*ionfrac=%e, sum(pop)=%e: (ZEROS)"%(Z,z1, z1, abund*ionfrac[z1-1], sum(lev_pop)*abund*ionfrac[z1-1]))
for i in range(len(lev_pop)):
print(i, lev_pop[i])
# calculate some continuum processes: brems
if settings['Bremsstrahlung'] ==True:
print("Calling do_brems from run_apec_ion at %s"%(time.asctime()))
brems = do_brems(Z, z1, te, 1.0, settings['BremsType'], ebins)
# scale for ion and element abundance.
continuum['brems']=brems*abund*ionfrac[z1-1]
print("Finished Calling do_brems from run_apec_ion at %s"%(time.asctime()))
else:
continuum['brems']=numpy.zeros(len(ebins)-1, dtype=float)
# now look at the neighbouring ions
#datacache={}
# z1_drv=z1*1
if z1_drv>1:
z1=z1_drv-1
print("Start processing Recombination from run_apec_ion at %s, Z=%i, z1=%i, z1_drv=%i"%(time.asctime(),Z,z1,z1_drv))
# do the DR satellite lines
if settings['DRSatellite']:
print("Start calc_satellte run_apec_ion at %s"%(time.asctime()))
linelist_dr, drlevrates = calc_satellite(Z, z1, te, datacache=datacache, settings=settings)
linelist_dr['epsilon']*=ionfrac[z1_drv-1]*abund
drlevrates *=ionfrac[z1_drv-1]*abund
print("Finished calc_satellte run_apec_ion at %s"%(time.asctime()))
else:
linelist_dr = numpy.zeros(0, dtype= generate_datatypes('linetype'))
drlevrates = 0.0
# Radiative Recombination
if settings['RRC']:
print("Start calc_rad_rec_cont at %s"%(time.asctime()))
rrc, rrlevrates = atomdb.calc_rad_rec_cont(Z, z1, z1_drv, te, ebins, settings=settings, datacache=datacache)
continuum['rrc'] = rrc*ionfrac[z1_drv-1]*abund
rrlevrates*=ionfrac[z1_drv-1]*abund
# print "rrlevrates from Z=%i, z1_drv=%i to z1=%i"%(Z, z1_drv, z1)
# for ilev in range(len(rrlevrates)):
# print ilev, rrlevrates[ilev]
print("Finished calc_rad_rec_cont at %s"%(time.asctime()))
else:
continuum['rrc'] = numpy.zeros(len(ebins)-1, dtype=float)
rrlevrates=0.0
# for i in range(len(rrlevrates)):
# print i+1, rrlevrates[i]
# zzz=raw_input()
# get the level specific recombination rates
# print len(drlevrates)
# if there is recombination to process:
tmpdrlevrates,xxx = util.make_vec(drlevrates)
tmprrlevrates,xxx = util.make_vec(rrlevrates)
if sum(tmpdrlevrates) + sum(tmprrlevrates)>0:
print("Start calc_recomb_popn at %s"%(time.asctime()))
levpop_recomb=calc_recomb_popn(lev_pop, Z, z1,\
z1_drv, te, dens, drlevrates,\
rrlevrates,\
datacache=datacache, settings=settings)
print("Finish calc_recomb_popn at %s"%(time.asctime()))
#zzz=raw_input()
print("Start do_lines Z=%i, z1=%i, z1_drv=%i at %s"%(Z,z1,z1_drv,time.asctime()))
linelist_rec, tmptwophot = \
do_lines(Z, z1, levpop_recomb , dens, datacache=datacache, settings=settings, z1_drv_in=z1_drv)
continuum['twophot']+= tmptwophot
print("linelist_rec Z=%i, z1=%i,z1_drv=%i, nlines=%i:"%(Z,z1, z1_drv, len(linelist_rec)))
print("Finish do_lines Z=%i, z1=%i, z1_drv=%i at %s"%(Z,z1,z1_drv,time.asctime()))
# now do the ionizing cases
linelist_ion = numpy.zeros(0,dtype= generate_datatypes('linetype'))
if z1_drv < Z:
#datacache={}
z1=z1_drv+1
lev_pop_parent = lev_pop*1.0
print("Sum lev_pop_parent[1:] = %e"%(sum(lev_pop_parent[1:])))
while (sum(lev_pop_parent[1:]) > 1e-40) &\
(z1 <= Z):
print("Processing ionization into Z=%i, z1=%i, z1_drv=%i at %s"%(Z,z1,z1_drv,time.asctime()))
# do we need to include collisional ionzation here?
# only in first ion
if z1== z1_drv+1:
do_xi = True
else:
do_xi = False
print("Calling calc_ioniz_popn at %s"%(time.asctime()))
lev_pop=calc_ioniz_popn(lev_pop_parent, Z, z1, z1_drv, te, dens, \
settings=settings, datacache=datacache, \
do_xi=do_xi)
print("Finished calc_ioniz_popn at %s"%(time.asctime()))
lev_pop[lev_pop<const.MIN_LEVPOP] = 0.0
if sum(lev_pop[1:]) > 0:
print("Start do_lines Z=%i, z1=%i, z1_drv=%i at %s"%(Z,z1,z1_drv,time.asctime()))
linelist_ion_tmp, tmptwophot = \
do_lines(Z, z1, lev_pop, dens, datacache=datacache, settings=settings, z1_drv_in=z1_drv)
print("linelist_ion Z=%i, z1=%i,z1_drv=%i, nlines=%i:"%(Z,z1, z1_drv, len(linelist_ion_tmp)))
print("Finished do_lines Z=%i, z1=%i, z1_drv=%i at %s"%(Z,z1,z1_drv,time.asctime()))
linelist_ion = numpy.append(linelist_ion, linelist_ion_tmp)
continuum['twophot']+=tmptwophot
lev_pop_parent = lev_pop
z1+=1
# generate return data
print("Start merging linelist at %s"%(time.asctime()))
linelist = numpy.append(linelist_exc, numpy.append(linelist_dr, numpy.append(linelist_ion, linelist_rec)))
print("Finished merging linelist at %s"%(time.asctime()))
# filter line list
print("Start filtering linelist at %s"%(time.asctime()))
MinEpsilon = settings['MinEpsilon']
if settings['Ionization']=='CIE':
MinEpsilon*=0.001
# istrong = linelist>=MinEpsilon
pseudocont = numpy.zeros(len(ebins)-1, dtype=float)
if len(linelist) > 0:
weaklines = linelist[(linelist['epsilon']< MinEpsilon) &\
(linelist['lambda']>const.HC_IN_KEV_A /settings['GridMaximum']) &\
(linelist['lambda']<const.HC_IN_KEV_A /settings['GridMinimum'])]
# print "filtered out % i weak lines"%(len(weaklines))
for line in weaklines:
e = const.HC_IN_KEV_A /line['lambda']
ibin = numpy.where(ebins>e)[0][0] - 1
pseudocont[ibin]+=line['epsilon']
linelist = linelist[linelist['epsilon'] > MinEpsilon]
print("Finish filtering linelist at %s"%(time.asctime()))
# print "kept % i strong lines"%(len(weaklines))
# ret = {}
# ret['lines'] = linelist
# ret['Z'] = Z
# ret['z1'] = z1
# ret['pseudocont'] = pseudocont
# ret['continuum'] = continuum
# ret['ionfrac'] = ionfrac
# fname = 'dump_%i_%i.pkl'%(Z,z1_drv)
# pickle.dump(ret, open(fname, 'wb'))
# print "wrote %s"%(fname)
if settings['WriteIon']==True:
ret = {}
ret['lines'] = linelist
ret['continuum'] = continuum
ret['pseudocont'] = pseudocont
ret['ionfrac'] = ionfrac
ret['te'] = te
ret['dens'] = dens
ret['settings'] = settings
ret['abund'] = abund
fname = settings['WriteIonFname']+'.pkl'
pickle.dump(ret, open(fname, 'wb'))
return linelist, continuum, pseudocont
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def compress_continuum(xin, yin, tolerance, minval = 0.0):
"""
Compress the continuum into linear interpolatable grids
Parameters
----------
xin : array(float)
The bin edges (keV)
yin : array(float)
The continuum in photons (or ergs) cm^3 s^-1 bin^-1.
Should be 1 element shorter then xin
tolerance : float
The tolerance of the final result (if 0.01, the result will always
be within 1% of the original value)
Returns
-------
xout : array (float)
The energy points of the compressed energy grid (keV)
yout : array (float)
The continuum, in photons(or ergs) cm^3 s^-1 keV^-1
"""
try:
from pyatomdb import liblinapprox
except OSError:
print("Unable to open liblinapprox. Hopefully you are on Readthedocs")
return
npts = len(yin)
xin_tmp = ctypes.c_double * npts
yin_tmp = ctypes.c_double * npts
ein_tmp = ctypes.c_double * npts
xout_tmp = ctypes.c_double * npts
xout_tmp = ctypes.c_double * npts
yout_tmp = ctypes.c_double * npts
wout_tmp = ctypes.c_double * npts
x = xin_tmp()
y = yin_tmp()
e = ein_tmp()
xout= xout_tmp()
yout= yout_tmp()
wout= wout_tmp()
m = ctypes.c_int(npts)
for i in range(npts):
x[i] = (xin[i]+xin[i+1])/2.0
y[i] = yin[i]/(xin[i+1]-xin[i]) # convert to per keV
e[i] = tolerance*y[i]
k = ctypes.c_int(0)
ip = ctypes.c_int(-6)
liblinapprox.linear_approx(ctypes.byref(x),\
ctypes.byref(y),\
ctypes.byref(e),\
ctypes.byref(m),\
ctypes.byref(xout),\
ctypes.byref(yout),\
ctypes.byref(wout),\
ctypes.byref(k),ctypes.byref(ip))
nout = k.value+1
xout = numpy.array(xout[:nout], dtype=float)
yout = numpy.array(yout[:nout], dtype=float)
# now trim if required
i = numpy.where(yout < minval)[0]
if len(i) > 0:
# we have something to do
yout[i] = 0.0
isgood = numpy.ones(len(yout), dtype=bool)
for i in range(1, len(yout)-1):
if ( (yout[i]<=minval) &\
(yout[i-1]<=minval) &\
(yout[i+1]<=minval)):
isgood[i] = False
xout = xout[isgood]
yout = yout[isgood]
return xout, yout
[docs]
def wrap_ion_directly(fname, ind, Z, z1):
import time
# read in the file
settings = parse_par_file(fname)
te = make_vector(settings['LinearTemp'], \
settings['TempStart'], \
settings['TempStep'], \
settings['NumTemp'])
if settings['TempUnits']=='keV':
te /= const.KBOLTZ
dens = make_vector(settings['LinearDens'], \
settings['DensStart'], \
settings['DensStep'], \
settings['NumDens'])
ite = ind //len(dens)
idens = ind%len(dens)
print(ite, idens)
Te = te[ite]
Dens = dens[idens]
if settings['Ionization']=='NEI':
z1list = list(range(1, Z+2))
ionfrac = numpy.ones(len(z1list), dtype=float)
elif settings['Ionization']=='CIE':
z1list = list(range(1, Z+2))
# time to calculate the ionization balance
if settings['UseIonBalanceTable']:
# read the ionization balance table
ionfrac = atomdb._get_precalc_ionfrac(os.path.expandvars(settings['IonBalanceTable']), Z, te)
else:
# calculate the ionization balance
ionftmp = calc_full_ionbal(Te, 1e14, Te_init=Te, Zlist=[Z], extrap=True)
ionfrac = ionftmp[Z]
else:
print("ERROR: settings['Ionization'] must be CIE or NEI, not %s"%(settings['Ionization']))
abundfile = atomdb.get_filemap_file('abund',\
Z,\
False,\
fmapfile=settings['FileMap'],\
atomdbroot=os.path.expandvars('$ATOMDB'),\
misc=True)
abundances = atomdb.get_abundance(abundfile=abundfile, abundset=settings['Abundances'])
#abundances = atomdb.get_abundance()
abund=abundances[Z]
# update the output filename
settings['WriteIonFname'] ="Z_%i_z1_%i_iT_%iiN_%i.pkl"%(Z,z1,ite,idens)
settings['HDUIndex'] = ind
settings['Te_used'] = Te
settings['Ne_used'] = Dens
x=run_apec_ion(settings, Te, Dens, Z, z1, ionfrac, abund)
ret = {}
ret['settings'] = settings
ret['data'] = x
fname = settings['OutputFileStem']+'_'+settings['WriteIonFname']
pickle.dump(ret, open(fname,'wb'))
print("wrote file %s"%(fname))
print("Finished cleanly at %s"%(time.asctime()))
[docs]
def wrap_run_apec(fname, readpickle=False, writepickle=False, override_iTe=False):
"""
After running the APEC code ion by ion, use this to combine into
FITS files.
Parameters
----------
fname : string
file name
readpickle : bool
Load apec results by element from pickle files, instead of regenerating
Returns
-------
None
"""
# get the input settings
settings = parse_par_file(fname)
# need to parse the IncAtoms parameter - this defines which atoms
# we will include
#
# we will transfer this to a "Zlist" parameter, in the form of
# nuclear charges
# if len(Zlist)==0:
Zlist = settings['Zlist']
print("I will be running Z=", Zlist)
# run for each element, temperature, density
lhdulist = []
chdulist = []
seclHDUdat = numpy.zeros(settings['NumTemp']*settings['NumDens'], \
dtype=generate_datatypes('lineparams'))
seccHDUdat = numpy.zeros(settings['NumTemp']*settings['NumDens'], \
dtype=generate_datatypes('cocoparams'))
listiTe = range(settings['NumTemp'])
if override_iTe != False:
listiTe = [override_iTe]
# for iTe in range(settings['NumTemp']):
for iTe in listiTe:
te = make_vector(settings['LinearTemp'], \
settings['TempStart'], \
settings['TempStep'], \
settings['NumTemp'])[iTe]
if settings['TempUnits']=='keV':
te /= const.KBOLTZ
for iDens in range(settings['NumDens']):
dens = make_vector(settings['LinearDens'], \
settings['DensStart'], \
settings['DensStep'], \
settings['NumDens'])[iDens]
# AT THIS POINT, GENERATE SHELL WHICH WILL GO IN THE HDU OF CHOICE
if settings['Ionization']=='CIE':
linedata = numpy.zeros(0,dtype=generate_datatypes('linelist_cie'))
cocodata = numpy.zeros(0,dtype=generate_datatypes('continuum', ncontinuum=0, npseudo=0))
if settings['Ionization']=='NEI':
linedata = numpy.zeros(0,dtype=generate_datatypes('linetype'))
cocodata = numpy.zeros(0,dtype=generate_datatypes('continuum', ncontinuum=0, npseudo=0))
for Z in Zlist:
print("Calling run_apec_element for Z=%i Te=%e dens=%e at %s"%(Z, te, dens, time.asctime()))
dat = wrap_run_apec_element(settings, te, dens, Z,iTe,iDens, readpickle=readpickle, writepickle=writepickle)
# append this data to the output
#pickle.dump(dat, open('dump_%i.pkl'%(Z),'wb'))
linedata = numpy.append(linedata, dat['lines'])
cocodata = continuum_append(cocodata, dat['cont'])
print("Z=%i, nlines=%i"%(Z, len(dat['lines'])))
# now make an HDU for all of this
if settings['Ionization']=='CIE':
LHDUdat = create_lhdu_cie(linedata)
elif settings['Ionization']=='NEI':
# pickle.dump(linedata,open('argh.pkl','wb'))
LHDUdat = create_lhdu_nei(linedata)
# now update the headers
iseclHDUdat=iDens+iTe*settings['NumDens']
LHDUdat.header['EXTNAME']=("EMISSIVITY","name of this binary table extension")
LHDUdat.header['EXTVER']=(iseclHDUdat+1,"Index for this EMISSIVITY extension")
LHDUdat.header['HDUNAME'] = ("L%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
'Spectral emission data')
LHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
'Proposed OGIP standard')
LHDUdat.header['HDUCLAS1']=("LINE MODEL",\
'Line emission spectral model')
LHDUdat.header['HDUCLAS2']=("LINE",\
'Emission line data')
if settings['Ionization']=='CIE':
LHDUdat.header['HDUVERS1']=("1.0.0",\
'version of format')
elif settings['Ionization']=='NEI':
LHDUdat.header['HDUVERS1']=("2.0.0",\
'version of format')
LHDUdat.header['TEMPERATURE']=(te,\
'Electron temperature')
LHDUdat.header['DENSITY']=(dens,\
'Electron density')
LHDUdat.header['TIME']=(0,\
'IN EQUILIBRIUM')
if settings['Ionization']=='CIE':
tot_emiss = sum(linedata['epsilon']*const.HC_IN_ERG_A/linedata['lambda'])
else:
tot_emiss=0.0
LHDUdat.header['TOT_LINE']=(tot_emiss,\
'Total Line Emission (erg cm^3 s^-1)')
LHDUdat.header['N_LINES']=(len(linedata),\
'Number of emission lines')
seclHDUdat['kT'][iseclHDUdat]=te*const.KBOLTZ
seclHDUdat['EDensity'][iseclHDUdat]= dens
seclHDUdat['Time'][iseclHDUdat]= 0.0
seclHDUdat['Nelement'][iseclHDUdat]= len(Zlist)
seclHDUdat['Nline'][iseclHDUdat]= len(linedata)
lhdulist.append(LHDUdat)
# continuum data
# now make an HDU for all of this
CHDUdat = create_chdu_variable_cie(cocodata)
# now update the headers
iseccHDUdat=iDens+iTe*settings['NumDens']
CHDUdat.header['EXTNAME']=("EMISSIVITY","name of this binary table extension")
CHDUdat.header['EXTVER']=(iseccHDUdat+1,"Index for this EMISSIVITY extension")
CHDUdat.header['HDUNAME'] = ("C%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
'Spectral emission data')
CHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
'Proposed OGIP standard')
CHDUdat.header['HDUCLAS1']=("COMP CONT MODEL",\
'Compressed continua spectra')
CHDUdat.header['HDUCLAS2']=("COCO",\
'Compressed continuum data')
CHDUdat.header['HDUVERS1']=("1.0.0",\
'version of format')
CHDUdat.header['TEMPERATURE']=(te,\
'Electron temperature')
CHDUdat.header['DENSITY']=(dens,\
'Electron density')
CHDUdat.header['TIME']=("%.2e"%(0),\
'IN EQUILIBRIUM')
tot_emiss = calc_total_coco(cocodata, settings)
if settings['Ionization']=='CIE':
CHDUdat.header['TOT_COCO']=(tot_emiss,\
'Total Emission (erg cm^3 s^-1)')
else:
CHDUdat.header['TOT_COCO']=(0.0,\
'Total Emission (erg cm^3 s^-1)')
seccHDUdat['kT'][iseccHDUdat]=te*const.KBOLTZ
seccHDUdat['EDensity'][iseccHDUdat]= dens
seccHDUdat['Time'][iseccHDUdat]= 0.0
seccHDUdat['NElement'][iseccHDUdat]= len(cocodata)
seccHDUdat['NCont'][iseccHDUdat]= max(cocodata['N_Cont'])
seccHDUdat['NPseudo'][iseccHDUdat]= max(cocodata['N_Pseudo'])
chdulist.append(CHDUdat)
# make secHDUdat into a fits table
seclHDU = create_lparamhdu_cie(seclHDUdat)
seclHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
seclHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
seclHDU.header['HDUCLAS1']=('LINE MODEL','line emission spectra model')
seclHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
seclHDU.header['HDUVERS1']=('1.0.0','version of format')
seccHDU = create_cparamhdu_cie(seccHDUdat)
seccHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
seccHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
seccHDU.header['HDUCLAS1']=('COMP CONT MODEL','compressed continua spectra')
seccHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
seccHDU.header['HDUVERS1']=('1.0.0','version of format')
fileroot = settings['OutputFileStem']
if settings['Ionization']=='NEI':
fileroot+='_nei'
# create the Primary HDU
PrilHDU = pyfits.PrimaryHDU()
lhdulist.insert(0,PrilHDU)
lhdulist.insert(1,seclHDU)
tmplhdulist = pyfits.HDUList(lhdulist)
PricHDU = pyfits.PrimaryHDU()
chdulist.insert(0,PricHDU)
chdulist.insert(1,seccHDU)
tmpchdulist = pyfits.HDUList(chdulist)
# now add header blurb
generate_apec_headerblurb(settings, tmplhdulist, tmpchdulist)
# write out results
tmplhdulist.writeto('%s_line.fits'%(fileroot), overwrite=True, checksum=True)
if settings['Ionization']=='CIE':
tmpchdulist.writeto('%s_coco.fits'%(fileroot), overwrite=True, checksum=True)
elif settings['Ionization']=='NEI':
tmpchdulist.writeto('%s_comp.fits'%(fileroot), overwrite=True, checksum=True)
[docs]
def wrap_update_one_ion(fname, linefile, cocofile, fnameout=None, noWrite=False):
"""
After running the APEC code ion by ion, use this to combine into
FITS files.
Parameters
----------
fname : string
file name of pickle file
linefile : string
file name of line file to update
cocofile : string
file name of cont file to update
Returns
-------
None
"""
# get the input settings
dat = pickle.load(open(fname, 'rb'))
#settings = parse_par_file(fname)
if type(linefile)==str:
ldat = pyfits.open(linefile)
else:
ldat = linefile
if type(cocofile)==str:
cdat = pyfits.open(cocofile)
else:
cdat = cocofile
# now, get the ion
lines = dat['data'][0]
if len(lines) > 0:
Z = lines[0]['element']
z1_drv = lines[0]['ion_drv']
else:
xxx1 = fname.split('_')
xxx = xxx1.index('Z')
Z=int(xxx1[xxx+1])
z1_drv = int(xxx1[xxx+3])
hduindex= dat['settings']['HDUIndex'] # the HDU in the ldat and cdat is this +2
# lines are the easy part:
ldat_dat_new = ldat[hduindex+2].data
cdat_dat_new = cdat[hduindex+2].data
# remove the lines for this ion
i = (ldat_dat_new['Element'] ==Z) & (ldat_dat_new['Ion_drv']==z1_drv)
ldat_dat_new =ldat_dat_new[~i]
# get the ion fraction
ionftmp= calc_full_ionbal(dat['settings']['Te_used'], 1e14, Te_init=dat['settings']['Te_used'], Zlist=[Z], extrap=True)
ionfrac_nei = ionftmp[Z]
ionfrac = ionfrac_nei[z1_drv-1]
print("Ionfrac is ", ionfrac)
# make a pseudocontinuum
ebins = make_vector_nbins(dat['settings']['LinearGrid'], \
dat['settings']['GridMinimum'], \
dat['settings']['GridMaximum'], \
dat['settings']['NumGrid'])
pseudo = {}
cont = {}
# now do some weak line filtering
mineps = dat['settings']['NEIMinEpsilon'][0]
pseudotmp = dat['data'][2]
for i in range(len(dat['settings']['NEIMinFrac'])):
if ionfrac < dat['settings']['NEIMinFrac'][i]:
mineps = dat['settings']['NEIMinEpsilon'][i+1]
linelist = dat['data'][0]
if len(dat['data'][0])> 0:
igood = numpy.ones(len(linelist), dtype=bool)
weaklines = linelist[(linelist['element']==Z) &\
(linelist['ion_drv']==z1_drv) &\
(linelist['epsilon']<mineps) &\
(linelist['lambda']>const.HC_IN_KEV_A /dat['settings']['GridMaximum']) &\
(linelist['lambda']<const.HC_IN_KEV_A /dat['settings']['GridMinimum'])]
print("identified %i weak lines"%(len(weaklines)))
for line in weaklines:
e = const.HC_IN_KEV_A /line['lambda']
ibin = numpy.where(ebins>e)[0][0] - 1
pseudotmp[ibin]+=line['epsilon']
igood[(linelist['element']==Z) &\
(linelist['ion_drv']==z1_drv) &\
(linelist['epsilon']<mineps)] = False
print("Filtered by mineps %e: from %i to %i lines"%(mineps, len(igood), sum(igood)))
linelist = linelist[igood]
conttmp =dat['data'][1]['rrc']+dat['data'][1]['twophot']+dat['data'][1]['brems']
econt, contin = compress_continuum(ebins, conttmp, const.TOLERANCE, minval=1e-38)
epseudo, pseudocont = compress_continuum(ebins, pseudotmp, const.TOLERANCE, minval=1e-38)
# now let's update!
ldat_new = numpy.zeros(len(linelist)+len(ldat_dat_new), dtype= generate_datatypes("linetype"))
for i in range(len(ldat_dat_new)):
for key in ldat_new.dtype.names:
ldat_new[key][i] = ldat_dat_new[key][i]
for i in range(len(linelist)):
for key in ldat_new.dtype.names:
ldat_new[key][i+len(ldat_dat_new)] = linelist[key][i]
LHDUdat = create_lhdu_nei(ldat_new)
# trim out the continuum data
i = (cdat_dat_new['Z'] ==Z) & (cdat_dat_new['rmJ']==z1_drv)
cdat_dat_new =cdat_dat_new[~i]
# gather the continuum data
cdat_new = numpy.zeros(1, dtype = generate_datatypes('continuum',
npseudo = len(pseudocont),
ncontinuum= len(contin)))
cdat_new['Z'][0] = Z
cdat_new['rmJ'][0] = z1_drv
cdat_new['N_Cont'][0] = len(econt)
cdat_new['E_Cont'][0][:len(econt)] = econt
cdat_new['Continuum'][0][:len(econt)] = contin
cdat_new['N_Pseudo'][0] = len(epseudo)
cdat_new['E_Pseudo'][0][:len(epseudo)] = epseudo
cdat_new['Pseudo'][0][:len(epseudo)] = pseudocont
cdat_new = continuum_append_variable(cdat_dat_new,cdat_new)
cdat_new = numpy.sort(cdat_new, order=['Z','rmJ'])
CHDUdat = create_chdu_variable_cie(cdat_new)
ldat[hduindex+2].data=LHDUdat.data
ldat[1].data['Nline'][hduindex]=len(ldat[hduindex+2].data)
print(CHDUdat.data.dtype)
print(cdat[hduindex+2].data.dtype)
cdat[hduindex+2].data=CHDUdat.data
cdat[1].data['NCont'][hduindex]=len(cdat[hduindex+2].data['E_Cont'][0])
cdat[1].data['NPseudo'][hduindex]=len(cdat[hduindex+2].data['E_Pseudo'][0])
if fnameout is None:
fnameout = linefile[:-8]
if noWrite==False:
ldat.writeto(fnameout+'line.fits', overwrite=True)
cdat.writeto(fnameout+'comp.fits', overwrite=True)
# LHDUdat.header['EXTVER']=(iseclHDUdat+1,"Index for this EMISSIVITY extension")
# LHDUdat.header['HDUNAME'] = ("L%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
# 'Spectral emission data')
# LHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
# 'Proposed OGIP standard')
# LHDUdat.header['HDUCLAS1']=("LINE MODEL",\
# 'Line emission spectral model')
# LHDUdat.header['HDUCLAS2']=("LINE",\
# 'Emission line data')
# if settings['Ionization']=='CIE':
# LHDUdat.header['HDUVERS1']=("1.0.0",\
# 'version of format')
# elif settings['Ionization']=='NEI':
# LHDUdat.header['HDUVERS1']=("2.0.0",\
# 'version of format')
# LHDUdat.header['TEMPERATURE']=(te,\
# 'Electron temperature')
# LHDUdat.header['DENSITY']=(dens,\
# 'Electron density')
# LHDUdat.header['TIME']=(0,\
# 'IN EQUILIBRIUM')
# if settings['Ionization']=='CIE':
# tot_emiss = sum(linedata['epsilon']*const.HC_IN_ERG_A/linedata['lambda'])
# else:
# tot_emiss=0.0
# LHDUdat.header['TOT_LINE']=(tot_emiss,\
# 'Total Line Emission (erg cm^3 s^-1)')
# LHDUdat.header['N_LINES']=(len(linedata),\
# 'Number of emission lines')
# seclHDUdat['kT'][iseclHDUdat]=te*const.KBOLTZ
# seclHDUdat['EDensity'][iseclHDUdat]= dens
# seclHDUdat['Time'][iseclHDUdat]= 0.0
# seclHDUdat['Nelement'][iseclHDUdat]= len(Zlist)
# seclHDUdat['Nline'][iseclHDUdat]= len(linedata)
# lhdulist.append(LHDUdat)
# # continuum data
# # now make an HDU for all of this
# CHDUdat = create_chdu_cie(cocodata)
# # now update the headers
# iseccHDUdat=iDens+iTe*settings['NumDens']
# CHDUdat.header['EXTNAME']=("EMISSIVITY","name of this binary table extension")
# CHDUdat.header['EXTVER']=(iseccHDUdat+1,"Index for this EMISSIVITY extension")
# CHDUdat.header['HDUNAME'] = ("C%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
# 'Spectral emission data')
# CHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
# 'Proposed OGIP standard')
# CHDUdat.header['HDUCLAS1']=("COMP CONT MODEL",\
# 'Compressed continua spectra')
# CHDUdat.header['HDUCLAS2']=("COCO",\
# 'Compressed continuum data')
# CHDUdat.header['HDUVERS1']=("1.0.0",\
# 'version of format')
# CHDUdat.header['TEMPERATURE']=(te,\
# 'Electron temperature')
# CHDUdat.header['DENSITY']=(dens,\
# 'Electron density')
# CHDUdat.header['TIME']=("%.2e"%(0),\
# 'IN EQUILIBRIUM')
# tot_emiss = calc_total_coco(cocodata, settings)
# if settings['Ionization']=='CIE':
# CHDUdat.header['TOT_COCO']=(tot_emiss,\
# 'Total Emission (erg cm^3 s^-1)')
# else:
# CHDUdat.header['TOT_COCO']=(0.0,\
# 'Total Emission (erg cm^3 s^-1)')
# seccHDUdat['kT'][iseccHDUdat]=te*const.KBOLTZ
# seccHDUdat['EDensity'][iseccHDUdat]= dens
# seccHDUdat['Time'][iseccHDUdat]= 0.0
# seccHDUdat['NElement'][iseccHDUdat]= len(cocodata)
# seccHDUdat['NCont'][iseccHDUdat]= max(cocodata['N_Cont'])
# seccHDUdat['NPseudo'][iseccHDUdat]= max(cocodata['N_Pseudo'])
# chdulist.append(CHDUdat)
# # make secHDUdat into a fits table
# seclHDU = create_lparamhdu_cie(seclHDUdat)
# seclHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
# seclHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
# seclHDU.header['HDUCLAS1']=('LINE MODEL','line emission spectra model')
# seclHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
# seclHDU.header['HDUVERS1']=('1.0.0','version of format')
# seccHDU = create_cparamhdu_cie(seccHDUdat)
# seccHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
# seccHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
# seccHDU.header['HDUCLAS1']=('COMP CONT MODEL','compressed continua spectra')
# seccHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
# seccHDU.header['HDUVERS1']=('1.0.0','version of format')
# fileroot = settings['OutputFileStem']
# if settings['Ionization']=='NEI':
# fileroot+='_nei'
# # create the Primary HDU
# PrilHDU = pyfits.PrimaryHDU()
# lhdulist.insert(0,PrilHDU)
# lhdulist.insert(1,seclHDU)
# tmplhdulist = pyfits.HDUList(lhdulist)
# PricHDU = pyfits.PrimaryHDU()
# chdulist.insert(0,PricHDU)
# chdulist.insert(1,seccHDU)
# tmpchdulist = pyfits.HDUList(chdulist)
# # now add header blurb
# generate_apec_headerblurb(settings, tmplhdulist, tmpchdulist)
# write out results
#tmplhdulist.writeto('%s_line.fits'%(fileroot), overwrite=True, checksum=True)
#if settings['Ionization']=='CIE':
# tmpchdulist.writeto('%s_coco.fits'%(fileroot), overwrite=True, checksum=True)
#elif settings['Ionization']=='NEI':
# tmpchdulist.writeto('%s_comp.fits'%(fileroot), overwrite=True, checksum=True)
[docs]
def wrap_update_multi_ion(fnamelist, linefile, cocofile, fnameout=None, noWrite=False):
"""
After running the APEC code ion by ion, use this to combine into
FITS files.
Parameters
----------
fname : string
file name of pickle file
linefile : string
file name of line file to update
cocofile : string
file name of cont file to update
Returns
-------
None
"""
# get the input settings
datacache={}
if type(linefile)==str:
ldat = pyfits.open(linefile)
else:
ldat = linefile
if type(cocofile)==str:
cdat = pyfits.open(cocofile)
else:
cdat = cocofile
for ifname, fname in enumerate(fnamelist):
print(fname)
dat = pickle.load(open(fname, 'rb'))
print('file read in')
#settings = parse_par_file(fname)
if ifname == 0:
hduindex= dat['settings']['HDUIndex'] # the HDU in the ldat and cdat is this +2
ldat_tmp = ldat[hduindex+2].data
ldat_dat_new = numpy.zeros(len(ldat_tmp), dtype= generate_datatypes("linetype"))
for i in range(len(ldat_dat_new)):
for key in ldat_dat_new.dtype.names:
ldat_dat_new[key][i] = ldat_tmp[key][i]
cdat_dat_new = cdat[hduindex+2].data
# now, get the ion
lines = dat['data'][0]
if len(lines) > 0:
Z = lines[0]['element']
z1_drv = lines[0]['ion_drv']
else:
xxx1 = fname.split('_')
xxx = xxx1.index('Z')
Z=int(xxx1[xxx+1])
z1_drv = int(xxx1[xxx+3])
# lines are the easy part:
print("Removing Lines")
# remove the lines for this ion
i = (ldat_dat_new['element'] ==Z) & (ldat_dat_new['ion_drv']==z1_drv)
ldat_dat_new =ldat_dat_new[~i]
print("Calc ion frac")
# get the ion fraction
ionftmp= calc_full_ionbal(dat['settings']['Te_used'], 1e14, cie=True, Zlist=[Z], extrap=True, datacache=datacache)
ionfrac_nei = ionftmp[Z]
ionfrac = ionfrac_nei[z1_drv-1]
print("Ionfrac is ", ionfrac)
# make a pseudocontinuum
ebins = make_vector_nbins(dat['settings']['LinearGrid'], \
dat['settings']['GridMinimum'], \
dat['settings']['GridMaximum'], \
dat['settings']['NumGrid'])
pseudo = {}
cont = {}
# now do some weak line filtering
mineps = dat['settings']['NEIMinEpsilon'][0]
pseudotmp = dat['data'][2]
for i in range(len(dat['settings']['NEIMinFrac'])):
if ionfrac < dat['settings']['NEIMinFrac'][i]:
mineps = dat['settings']['NEIMinEpsilon'][i+1]
linelist = dat['data'][0]
if len(dat['data'][0])> 0:
igood = numpy.ones(len(linelist), dtype=bool)
weaklines = linelist[(linelist['element']==Z) &\
(linelist['ion_drv']==z1_drv) &\
(linelist['epsilon']<mineps) &\
(linelist['lambda']>const.HC_IN_KEV_A /dat['settings']['GridMaximum']) &\
(linelist['lambda']<const.HC_IN_KEV_A /dat['settings']['GridMinimum'])]
print("identified %i weak lines"%(len(weaklines)))
for line in weaklines:
e = const.HC_IN_KEV_A /line['lambda']
ibin = numpy.where(ebins>e)[0][0] - 1
pseudotmp[ibin]+=line['epsilon']
igood[(linelist['element']==Z) &\
(linelist['ion_drv']==z1_drv) &\
(linelist['epsilon']<mineps)] = False
print("Filtered by mineps %e: from %i to %i lines"%(mineps, len(igood), sum(igood)))
linelist = linelist[igood]
conttmp =dat['data'][1]['rrc']+dat['data'][1]['twophot']+dat['data'][1]['brems']
econt, contin = compress_continuum(ebins, conttmp, const.TOLERANCE, minval=1e-38)
epseudo, pseudocont = compress_continuum(ebins, pseudotmp, const.TOLERANCE, minval=1e-38)
# now let's update!
# print(type(linelist))
# print(type(numpy.array(ldat_dat_new)))
# print("Mergingin lines")
# ldat_new = numpy.zeros(len(linelist)+len(ldat_dat_new), dtype= generate_datatypes("linetype"))
# for i in range(len(ldat_dat_new)):
# for key in ldat_new.dtype.names:
# ldat_new[key][i] = ldat_dat_new[key][i]
# for i in range(len(linelist)):
# for key in ldat_new.dtype.names:
# ldat_new[key][i+len(ldat_dat_new)] = linelist[key][i]
ldat_dat_new =numpy.append(ldat_dat_new, linelist)
print("Lines merged")
# for i in range(len(ldat_dat_new)):
# for key in ldat_new.dtype.names:
# ldat_new[key][i] = ldat_dat_new[key][i]
#
# for i in range(len(linelist)):
# for key in ldat_new.dtype.names:
# ldat_new[key][i+len(ldat_dat_new)] = linelist[key][i]
print("now doing the comp")
i = (cdat_dat_new['Z'] ==Z) & (cdat_dat_new['rmJ']==z1_drv)
print('cdat_dat_new[0] before:', cdat_dat_new[0])
cdat_dat_new =cdat_dat_new[~i]
# gather the continuum data
cdat_new = numpy.zeros(1, dtype = generate_datatypes('continuum',
npseudo = len(pseudocont),
ncontinuum= len(contin)))
print("Generated new CDAT!")
cdat_new['Z'][0] = Z
cdat_new['rmJ'][0] = z1_drv
cdat_new['N_Cont'][0] = len(econt)
cdat_new['E_Cont'][0][:len(econt)] = econt
cdat_new['Continuum'][0][:len(econt)] = contin
cdat_new['N_Pseudo'][0] = len(epseudo)
cdat_new['E_Pseudo'][0][:len(epseudo)] = epseudo
cdat_new['Pseudo'][0][:len(epseudo)] = pseudocont
print('cdat_new', cdat_new)
print('cdat_dat_new[0] at call:', cdat_dat_new[0])
cdat_dat_new = continuum_append_variable(cdat_dat_new,cdat_new)
print("Continuum Appended")
cdat_dat_new = numpy.sort(cdat_dat_new, order=['Z','rmJ'])
# dat_dat_new = numpy.sort(cdat_new, order=['Z','rmJ'])
LHDUdat = create_lhdu_nei(ldat_dat_new)
# trim out the continuum data
CHDUdat = create_chdu_variable_cie(cdat_dat_new)
ldat[hduindex+2].data=LHDUdat.data
ldat[1].data['Nline'][hduindex]=len(ldat[hduindex+2].data)
cdat[hduindex+2].data=CHDUdat.data
cdat[1].data['NCont'][hduindex]=max(cdat[hduindex+2].data['N_Cont'])
cdat[1].data['NPseudo'][hduindex]=max(cdat[hduindex+2].data['N_Pseudo'])
if fnameout is None:
fnameout = linefile[:-8]
if noWrite==False:
ldat.writeto(fnameout+'line.fits', overwrite=True)
cdat.writeto(fnameout+'comp.fits', overwrite=True)
# LHDUdat.header['EXTVER']=(iseclHDUdat+1,"Index for this EMISSIVITY extension")
# LHDUdat.header['HDUNAME'] = ("L%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
# 'Spectral emission data')
# LHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
# 'Proposed OGIP standard')
# LHDUdat.header['HDUCLAS1']=("LINE MODEL",\
# 'Line emission spectral model')
# LHDUdat.header['HDUCLAS2']=("LINE",\
# 'Emission line data')
# if settings['Ionization']=='CIE':
# LHDUdat.header['HDUVERS1']=("1.0.0",\
# 'version of format')
# elif settings['Ionization']=='NEI':
# LHDUdat.header['HDUVERS1']=("2.0.0",\
# 'version of format')
# LHDUdat.header['TEMPERATURE']=(te,\
# 'Electron temperature')
# LHDUdat.header['DENSITY']=(dens,\
# 'Electron density')
# LHDUdat.header['TIME']=(0,\
# 'IN EQUILIBRIUM')
# if settings['Ionization']=='CIE':
# tot_emiss = sum(linedata['epsilon']*const.HC_IN_ERG_A/linedata['lambda'])
# else:
# tot_emiss=0.0
# LHDUdat.header['TOT_LINE']=(tot_emiss,\
# 'Total Line Emission (erg cm^3 s^-1)')
# LHDUdat.header['N_LINES']=(len(linedata),\
# 'Number of emission lines')
# seclHDUdat['kT'][iseclHDUdat]=te*const.KBOLTZ
# seclHDUdat['EDensity'][iseclHDUdat]= dens
# seclHDUdat['Time'][iseclHDUdat]= 0.0
# seclHDUdat['Nelement'][iseclHDUdat]= len(Zlist)
# seclHDUdat['Nline'][iseclHDUdat]= len(linedata)
# lhdulist.append(LHDUdat)
# # continuum data
# # now make an HDU for all of this
# CHDUdat = create_chdu_cie(cocodata)
# # now update the headers
# iseccHDUdat=iDens+iTe*settings['NumDens']
# CHDUdat.header['EXTNAME']=("EMISSIVITY","name of this binary table extension")
# CHDUdat.header['EXTVER']=(iseccHDUdat+1,"Index for this EMISSIVITY extension")
# CHDUdat.header['HDUNAME'] = ("C%.2f_%.2f"%(numpy.log10(te), numpy.log10(dens)),\
# 'Spectral emission data')
# CHDUdat.header['HDUCLASS'] = ("Proposed OGIP",\
# 'Proposed OGIP standard')
# CHDUdat.header['HDUCLAS1']=("COMP CONT MODEL",\
# 'Compressed continua spectra')
# CHDUdat.header['HDUCLAS2']=("COCO",\
# 'Compressed continuum data')
# CHDUdat.header['HDUVERS1']=("1.0.0",\
# 'version of format')
# CHDUdat.header['TEMPERATURE']=(te,\
# 'Electron temperature')
# CHDUdat.header['DENSITY']=(dens,\
# 'Electron density')
# CHDUdat.header['TIME']=("%.2e"%(0),\
# 'IN EQUILIBRIUM')
# tot_emiss = calc_total_coco(cocodata, settings)
# if settings['Ionization']=='CIE':
# CHDUdat.header['TOT_COCO']=(tot_emiss,\
# 'Total Emission (erg cm^3 s^-1)')
# else:
# CHDUdat.header['TOT_COCO']=(0.0,\
# 'Total Emission (erg cm^3 s^-1)')
# seccHDUdat['kT'][iseccHDUdat]=te*const.KBOLTZ
# seccHDUdat['EDensity'][iseccHDUdat]= dens
# seccHDUdat['Time'][iseccHDUdat]= 0.0
# seccHDUdat['NElement'][iseccHDUdat]= len(cocodata)
# seccHDUdat['NCont'][iseccHDUdat]= max(cocodata['N_Cont'])
# seccHDUdat['NPseudo'][iseccHDUdat]= max(cocodata['N_Pseudo'])
# chdulist.append(CHDUdat)
# # make secHDUdat into a fits table
# seclHDU = create_lparamhdu_cie(seclHDUdat)
# seclHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
# seclHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
# seclHDU.header['HDUCLAS1']=('LINE MODEL','line emission spectra model')
# seclHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
# seclHDU.header['HDUVERS1']=('1.0.0','version of format')
# seccHDU = create_cparamhdu_cie(seccHDUdat)
# seccHDU.header['EXTNAME']=('PARAMETERS','name of this binary table extension')
# seccHDU.header['HDUCLASS']=('Proposed OGIP','Proposed OGIP standard')
# seccHDU.header['HDUCLAS1']=('COMP CONT MODEL','compressed continua spectra')
# seccHDU.header['HDUCLAS2']=('Parameters','extension containing parameter info')
# seccHDU.header['HDUVERS1']=('1.0.0','version of format')
# fileroot = settings['OutputFileStem']
# if settings['Ionization']=='NEI':
# fileroot+='_nei'
# # create the Primary HDU
# PrilHDU = pyfits.PrimaryHDU()
# lhdulist.insert(0,PrilHDU)
# lhdulist.insert(1,seclHDU)
# tmplhdulist = pyfits.HDUList(lhdulist)
# PricHDU = pyfits.PrimaryHDU()
# chdulist.insert(0,PricHDU)
# chdulist.insert(1,seccHDU)
# tmpchdulist = pyfits.HDUList(chdulist)
# # now add header blurb
# generate_apec_headerblurb(settings, tmplhdulist, tmpchdulist)
# write out results
#tmplhdulist.writeto('%s_line.fits'%(fileroot), overwrite=True, checksum=True)
#if settings['Ionization']=='CIE':
# tmpchdulist.writeto('%s_coco.fits'%(fileroot), overwrite=True, checksum=True)
#elif settings['Ionization']=='NEI':
# tmpchdulist.writeto('%s_comp.fits'%(fileroot), overwrite=True, checksum=True)
[docs]
def wrap_run_apec_element(settings, te, dens, Z, ite, idens, writepickle=False, readpickle=False):
"""
Combine wrap_run_apec_ion results for an element
Parameters
----------
settings: dictionary
The settings read from the apec.par file by parse_par_file
te: float
The electron temperature (K)
dens: float
The electron density (cm^-3)
Z: int
The nuclear charge of the element
ite: int
The temperature index
idens: int
The density index
writepickle: bool
Dump data into a pickle file. Useful for rapidly combining data
after runs.
readpickle: bool
Read data from a pickle file. Useful for rapidly combining data
after runs. Usually the result of a previous call using
writepickle=True
Returns
-------
None
"""
if readpickle:
setpicklefname = "%s_Z_%i_elem_iT_%iiN_%i.pkl"%(settings['OutputFileStem'],Z,ite,idens)
print("Trying to read in the whole element pickle file %s"%(setpicklefname))
if os.path.exists(setpicklefname):
ret = pickle.load(open(setpicklefname,'rb'))
print("read %s"%(setpicklefname))
return ret
else:
print("Not found. Going to assemble it.")
linelist = numpy.zeros(0, dtype=generate_datatypes('linetype'))
contlist = {}
pseudolist = {}
# now go through each ion and assemble the data
ebins = make_vector_nbins(settings['LinearGrid'], \
settings['GridMinimum'], \
settings['GridMaximum'], \
settings['NumGrid'])
ecent = (ebins[1:]+ebins[:-1])/2
for z1_drv in range(1,Z+2):
setpicklefname = "%s_Z_%i_z1_%i_iT_%iiN_%i.pkl"%(settings['OutputFileStem'],Z,z1_drv,ite,idens)
print("loading %s"%(setpicklefname))
if not os.path.exists(setpicklefname):
print("Warning: no such file: %s"%(setpicklefname))
contlist[z1_drv]={}
contlist[z1_drv]['rrc']=numpy.zeros(settings['NumGrid'], dtype=float)
contlist[z1_drv]['twophot']=numpy.zeros(settings['NumGrid'], dtype=float)
contlist[z1_drv]['brems']=numpy.zeros(settings['NumGrid'], dtype=float)
pseudolist[z1_drv] = numpy.zeros(settings['NumGrid'], dtype=float)
else:
dat= pickle.load(open(setpicklefname,'rb'))
tmplinelist, tmpcontinuum, tmppseudocont = dat['data']
# check for NAN
nlines = len(tmplinelist)
print(nlines)
ngoodlines = sum(numpy.isfinite(tmplinelist['epsilon']))
if nlines != ngoodlines:
print("Bad lines found in %s"%(setpicklefname))
linelist = numpy.append(linelist, tmplinelist)
for key in list(tmpcontinuum.keys()):
tmpncont = len(tmpcontinuum[key])
if tmpncont != sum(numpy.isfinite(tmpcontinuum[key])):
print("Bad continuum found in %s %s"%(key, setpicklefname), end=' ')
# if key=='rrc':
# tmpcontinuum['rrc'][numpy.isnan(tmpcontinuum['rrc'])]=0.0
print("")
contlist[z1_drv] = tmpcontinuum
tmpncont = len(tmppseudocont)
if tmpncont != sum(numpy.isfinite(tmppseudocont)):
print("Bad pseudocont found in %s"%( setpicklefname))
pseudolist[z1_drv] = tmppseudocont
# now merge these together.
if settings['Ionization']=='CIE':
cieout = generate_cie_outputs(settings, Z, linelist, contlist, pseudolist)
if writepickle:
setpicklefname = "%s_Z_%i_elem_iT_%iiN_%i.pkl"%(settings['OutputFileStem'],Z,ite,idens)
pickle.dump(cieout, open(setpicklefname,'wb'))
print("wrote %s"%(setpicklefname))
return cieout
elif settings['Ionization']=='NEI':
ionftmp= calc_full_ionbal(te, extrap=True, cie=True, settings=settings, Zlist=[Z])
ionfrac_nei = ionftmp[Z]
neiout = generate_nei_outputs(settings, Z, linelist, contlist, pseudolist, ionfrac_nei)
if writepickle:
setpicklefname = "%s_Z_%i_elem_iT_%iiN_%i.pkl"%(settings['OutputFileStem'],Z,ite,idens)
pickle.dump(neiout, open(setpicklefname,'wb'))
print("wrote %s"%(setpicklefname))
return neiout
#-------------------------------------------------------------------------------
[docs]
def run_wrap_run_apec(fname, Z, iTe, iDens):
"""
After running the APEC code ion by ion, use this to combine into
FITS files.
Parameters
----------
fname : string
file name of par file
Z: int
The atomic numbers
iTe: int
The temperature index
iDens: int
The density index
Returns
-------
None
"""
# get the input settings
settings = parse_par_file(fname)
# need to parse the IncAtoms parameter - this defines which atoms
# we will include
#
# we will transfer this to a "Zlist" parameter, in the form of
# nuclear charges
print("I will be running Z=", Z)
# run for each element, temperature, density
lhdulist = []
chdulist = []
seclHDUdat = numpy.zeros(settings['NumTemp']*settings['NumDens'], \
dtype=generate_datatypes('lineparams'))
seccHDUdat = numpy.zeros(settings['NumTemp']*settings['NumDens'], \
dtype=generate_datatypes('cocoparams'))
te = make_vector(settings['LinearTemp'], \
settings['TempStart'], \
settings['TempStep'], \
settings['NumTemp'])[iTe]
if settings['TempUnits']=='keV':
te /= const.KBOLTZ
dens = make_vector(settings['LinearDens'], \
settings['DensStart'], \
settings['DensStep'], \
settings['NumDens'])[iDens]
# AT THIS POINT, GENERATE SHELL WHICH WILL GO IN THE HDU OF CHOICE
if settings['Ionization']=='CIE':
linedata = numpy.zeros(0,dtype=generate_datatypes('linelist_cie'))
cocodata = numpy.zeros(0,dtype=generate_datatypes('continuum', ncontinuum=0, npseudo=0))
if settings['Ionization']=='NEI':
linedata = numpy.zeros(0,dtype=generate_datatypes('linetype'))
cocodata = numpy.zeros(0,dtype=generate_datatypes('continuum', ncontinuum=0, npseudo=0))
print("Calling run_apec_element for Z=%i Te=%e dens=%e at %s"%(Z, te, dens, time.asctime()))
dat = wrap_run_apec_element(settings, te, dens, Z,iTe,iDens, writepickle=True)
print("Done safely")
#-------------------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
[docs]
def return_ionbal(Z, Te, init_pop=False, tau=False,\
teunit='K', \
filename=False, datacache=False, fast=True,
settings= False, debug=False, extrap=True,\
allowmulti=True):
"""
Solve the ionization balance for a element Z.
Parameters
----------
Z : int
atomic number of element
Te : float or array
electron temperature(s), default in K
init_pop : float array
initial population of ions for non-equlibrium calculations. Will be renormalised to 1.
tau : float or array
N_e * t for the non-equilibrium ioniziation, in cm^3 s^-1.
Te_init : float
initial ionization balance temperature, same units as Te
teunit : {'K' , 'keV'}
units of temperatures (default K)
filename : string
Can optionally point directly to the file in question, i.e. to look at older data
look at $HEADAS/../spectral/modelData/eigenELSYMB_v3.0.fits.
If not set, download from AtomDB FTP site.
datacache : dict
Used for caching the data. See description in atomdb.get_data
fast : bool
If true, use precalculated eigenvector files to obtain CIE and NEI results
Returns
-------
final_pop : float array
final populations.
"""
if fast:
ionbal = _solve_ionbal_eigen(Z, Te, init_pop=init_pop, tau=tau, \
teunit=teunit, \
filename=filename, datacache=datacache, debug=debug)
return ionbal
else:
# if tau!=False:
# cie=False
# else:
# cie=True
ionbal = calc_elem_ionbal(Z, Te, tau=tau, init_pop=init_pop, teunit=teunit,\
extrap=extrap, settings=settings, datacache=datacache,\
allowmulti=allowmulti)
return ionbal
def _solve_ionbal_eigen(Z, Te, init_pop=False, tau=False, \
teunit='K', \
filename=False, datacache=False, debug=False):
"""
Solve the ionization balance for element Z using the eigenvector
approach and files as distributed in XSPEC.
Parameters
----------
Z : int
atomic number of element
Te : float or array
electron temperature(s), default in K
init_pop : float array
initial population of ions for non-equlibrium calculations. Will be renormalised to 1.
tau : float or array
N_e * t for the non-equilibrium ioniziation, in cm^3 s^-1.
teunit : {'K' , 'keV'}
units of temperatures (default K)
filename : string
Can optionally point directly to the file in question, i.e. to look at older data
look at $HEADAS/../spectral/modelData/eigenELSYMB_v3.0.fits.
If not set, download from AtomDB FTP site.
datacache : dict
Used for caching the data. See description in atomdb.get_data
Returns
-------
final_pop : float array
final populations.
"""
#
# Version 0.1 Initial Release
# Adam Foster 16th September 2015
#
kT = util.convert_temp(Te, teunit, 'keV')
if type(tau)==bool:
if tau==False:
cie = True
init_pop_calc=False
else:
raise ValueError("Error: tau should be False, a float, or an array of floats. Received "+repr(tau))
else:
cie = False
if not cie:
# if it's not equilibrium, get the initial population
if isinstance(init_pop, str):
if init_pop.lower() == 'ionizing':
init_pop_calc = numpy.zeros(Z+1)
init_pop_calc[0] = 1.0
elif init_pop.lower() == 'recombining':
init_pop_calc = numpy.zeros(Z+1)
init_pop_calc[-1] = 1.0
else:
raise util.OptionError("Error: init_pop is set as a string, must be 'ionizing' or 'recombining'. Currently %s."%\
(init_pop))
elif isinstance(init_pop, float):
# this is an initial temperature
kT_init = util.convert_temp(init_pop, teunit, 'keV')
init_pop_calc = return_ionbal(Z, kT_init, \
teunit='keV', \
filename=filename,\
datacache=datacache,fast=True)
elif isinstance(init_pop, numpy.ndarray) or isinstance(init_pop, list):
init_pop_calc = init_pop
elif isinstance(init_pop, dict):
init_pop_calc = init_pop[Z]
else:
raise util.OptionError("Error: invalid type for init_pop")
# open the eigenvector data file
if util.keyword_check(filename):
# we have a filename specified!
if type(filename)==str:
fname = os.path.expandvars(filename)
if not os.path.isfile(fname):
print("Specified file %s does not exist. Exiting"%(fname))
return
d = pyfits.open(fname)
else:
d = filename
else:
d = atomdb.get_data(Z, False, 'eigen', datacache=datacache)
telist = numpy.logspace(4,9,1251)
kTlist=telist*const.KBOLTZ
try:
dd = d['EIGEN']
except KeyError:
elsymb = atomic.Ztoelsymb(Z).lower()
dd=d[elsymb]
# if we are looking for equilibrium, return the nearest data
if cie:
ikTlist = numpy.argsort(numpy.abs(kTlist-kT))
ite = [min(ikTlist[:2]), max(ikTlist[:2])]
Tdiff = kTlist[ite[1]] - kTlist[ite[0]]
if Tdiff > 0.0:
factorlow = (kTlist[ite[1]]-kT)/Tdiff
factorhigh = (kT-kTlist[ite[0]])/Tdiff
equilib = factorlow * dd.data['FEQB'][ite[0]]+\
factorhigh * dd.data['FEQB'][ite[1]]
else:
equilib = dd.data['FEQB'][ite[0]]
#renormalize
equilib /= sum(equilib)
return equilib
# now do the non-equilibrium data
# renormalize
# make Te into a vector
kT_vec, kT_isvec = util.make_vec(kT)
tau_vec, tau_isvec = util.make_vec(tau)
frac_out = numpy.zeros([len(kT_vec),len(tau_vec),Z+1], dtype=float)
for ikT, kT in enumerate(kT_vec):
kTindex = numpy.argmin(numpy.abs(kTlist-kT))
lefteigenvec = numpy.zeros([Z,Z], dtype=float)
righteigenvec = numpy.zeros([Z,Z], dtype=float)
if Z==1:
for i in range(Z):
for j in range(Z):
lefteigenvec[i,j] = dd.data['VL'][kTindex]
righteigenvec[i,j] = dd.data['VR'][kTindex]
else:
for i in range(Z):
for j in range(Z):
lefteigenvec[i,j] = dd.data['VL'][kTindex][i*Z+j]
righteigenvec[i,j] = dd.data['VR'][kTindex][i*Z+j]
work = numpy.array(init_pop_calc[1:] - dd.data['FEQB'][kTindex][1:], dtype=float)
fspectmp = numpy.matrix(lefteigenvec) * numpy.matrix(work).transpose()
delt = 1.0
worktmp = numpy.zeros(Z)
for itau, ttau in enumerate(tau_vec):
if Z >1:
for i in range(Z):
worktmp[i] = fspectmp[i]*numpy.exp(dd.data['EIG'][kTindex,i]*delt*ttau)
else:
worktmp[0] = fspectmp[0]*numpy.exp(dd.data['EIG'][kTindex]*delt*ttau)
frac = numpy.zeros(Z+1)
for i in range(Z):
for j in range(Z):
frac[i+1] += worktmp[j]*righteigenvec[j][i]
frac[i+1] += dd.data['FEQB'][kTindex][i+1]
if debug:
frac_out[ikT, itau,:] = frac
frac[frac<0.0] = 0.0
if sum(frac)> 1.0:
frac = frac/sum(frac)
frac[0] = 1-sum(frac[1:])
frac[frac<0.0] = 0.0
if not(debug):
frac_out[ikT,itau,:]=frac
if not tau_isvec:
frac_out = frac_out.sum(1)
if not kT_isvec:
frac_out = frac_out.sum(0)
return frac_out