PyAtomDB Example Scripts¶
Module Structure¶
These are examples of using the PyAtomDB module in your projects. They can all be found in the examples subdirectory of your PyAtomDB installation.
The examples here are separated by the module in which they are implemented. Generally speaking, the roles of these modules can be split into the following areas:
- spectrum
Obtaining emissivities from the AtomDB emissivity files (e.g. the
apecmodel). This module is used to create spectra as required.- apec
Related to the APEC code - creating the AtomDB emissivity files from the underlying atomic database
- atomdb
Accessing the atomic database - returning rates and coefficients, fetching files
- atomic
Basic atomic functions - getting atomic masses, converting to element symbols etc.
- util
File utilities (not generally interesting to end users)
- const
List of physical constants and code-related constants used throughout PyAtomDB (again, not generally interesting to end users)
Spectra¶
The spectrum.py module contains routines for taking the spectrum generated from the apec model and extracting line emissivities or continuum emission, applying responses, changing abundances, etc. In these examples, we will use the Chandra ACIS-S MEG +1 order grating responses as examples where one is required, but you can use any.
Similarly, we use pylab for plotting. You can of course use whatever system you
like for plotting, this is just what we did. It is included in matplotlib and
therefore hopefully present on most systems.
CIESession Class¶
The heart of the spectral analysis is the spectrum.py class. This reads in the results of an apec run (by default, $ATOMDB/apec_line.fits and $ATOMDB/apec_coco.fits) and allows the user to obtain spectra at a range of temperatures accounting for instrument responses, thermal and velocity broadening, abundance changes and other issues.
Making a Spectrum¶
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
# create a set of energy bins (in keV) for the response. Note these are
# the n edges of the n-1 bins.
ebins = numpy.linspace(0.3,1.0,10000)
# set the response (raw keyword tells pyatomdb it is not a real response file)
sess.set_response(ebins, raw=True)
kT = 0.4 # temperature in keV
spec = sess.return_spectrum(kT)
# we now have a spectrum and the energy bins,
# prepend a zero to the spectrum so that energy bins and spectrum have
# the same length. If you use a different plotting system you may
# need to add this to the end.
spec = numpy.append(0, spec)
# Returned spectrum has units of photons cm^5 s^-1 bin^-1
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(211)
ax.plot(sess.ebins_out, spec, drawstyle='steps', label='dummy response')
# label the figure
ax.set_xlabel('Energy (keV)')
ax.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
# zoom in a bit
#ax.set_xlim([0.8,0.85])
# make plot appear
# now repeat the process with a real response
# set response
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# get spectrum, prepend zero
spec = sess.return_spectrum(kT)
spec = numpy.append(0, spec)
ax2 = fig.add_subplot(212)
ax2.plot(sess.ebins_out, spec, drawstyle='steps', label='Chandra MEG')
# add legends
ax.legend(loc=0)
ax2.legend(loc=0)
# zoom in on small sections of the spectrum
ax.set_xlim([0.7, 0.9])
ax2.set_xlim([0.7, 0.9])
# set axes
ax2.set_xlabel('Energy (keV)')
ax2.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
# adjust plot spacing so labels are visible
pylab.matplotlib.pyplot.subplots_adjust(hspace = 0.34)
# draw graphs
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_1_1.pdf')
fig.savefig('spectrum_session_examples_1_1.svg')
A kT=0.4keV simple spectrum created with and without an instrument response¶
Handling Large Response Matrices¶
Some missions, for example XRISM, have a large number of bins in their response files (XRISM has 50,000). Creating a 50,000x50,000 response matrix would require of order 16GB of memory, most of which would be zeros. Setting sparse=True when setting the response will use sparse matrices to save memory. The code below shows how to do this and also how to compare any innacuracies caused (testing so far has shown nothing larger than numerical rounding errors, 1 part in 10 15)
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
kT = 0.4 # temperature in keV
# Returned spectrum has units of photons cm^5 s^-1 bin^-1
# create figure
gs = pylab.matplotlib.gridspec.GridSpec(2, 1, height_ratios=[2, 1])
fig = pylab.figure()
fig.show()
ax1 = fig.add_subplot(gs[0]) # for data
ax2 = fig.add_subplot(gs[1], sharex=ax1) # for ratio
# label the figure
ax2.set_xlabel('Energy (keV)')
ax1.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
ax2.set_ylabel('Ratio')
# set response - using a regular matrix
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# get spectrum, prepend zero
spec = sess.return_spectrum(kT)
spec = numpy.append(spec[0], spec)
ax1.plot(sess.ebins_out, spec, drawstyle='steps', label='Regular')
# change the response to a sparse one
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf',
sparse=True)
spec_sparse = sess.return_spectrum(kT)
spec_sparse = numpy.append(spec_sparse[0], spec_sparse)
ax1.plot(sess.ebins_out, spec_sparse, '--', drawstyle='steps', label='Sparse')
# plot the ratio
ax2.plot(sess.ebins_out, spec/spec_sparse, drawstyle='steps',\
label='Regular/Sparse')
# add legend
ax1.legend(loc=0)
ax2.legend(loc=0)
# zoom in on small sections of the spectrum
ax1.set_xlim([0.5, 2.0])
# adjust plot spacing so labels are visible
pylab.matplotlib.pyplot.subplots_adjust(hspace = 0)
# draw graphs
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_1b_1.pdf')
fig.savefig('spectrum_session_examples_1b_1.svg')
A kT=0.4keV simple spectrum created with regular and sparse instrument responses. The same code can be used to compare the accuracy of the two methods for different response matrices, but will require large amounts of memory if using large responses.¶
Line Broadening¶
By default, line broadening is off. The command session.set_broadening allows you to turn on thermal broadening and, if desired, add additional turbulent velocity broadening too.
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
# create a set of energy bins (in keV) for the response. Note these are
# the n edges of the n-1 bins.
ebins = numpy.linspace(1.0,1.2,501)
# set the response (raw keyword tells pyatomdb it is not a real response file)
sess.set_response(ebins, raw=True)
kT = 4.0 # temperature in keV
vel = 400.0 # velocity to broaden with when using velocity broadening (km/s)
spec = sess.return_spectrum(kT)
# we now have a spectrum and the energy bins,
# prepend a zero to the spectrum so that energy bins and spectrum have
# the same length. If you use a different plotting system you may
# need to add this to the end.
spec = numpy.append(0, spec)
# Returned spectrum has units of photons cm^5 s^-1 bin^-1
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(211)
ax.plot(sess.ebins_out, spec, drawstyle='steps', label='dummy response')
# label the figure
ax.set_xlabel('Energy (keV)')
ax.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
sess.set_broadening(True)
spec = sess.return_spectrum(kT)
spec = numpy.append(0, spec)
ax.plot(sess.ebins_out, spec, drawstyle='steps', label='thermal')
sess.set_broadening(True, velocity_broadening=vel)
spec = sess.return_spectrum(kT)
spec = numpy.append(0, spec)
ax.plot(sess.ebins_out, spec, drawstyle='steps', label='%.0f km/s'%(vel))
# make plot appear
# now repeat the process with a real response
# reset broadening to off.
sess.set_broadening(False)
# set response
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# get spectrum, prepend zero
spec = sess.return_spectrum(kT)
spec = numpy.append(0, spec)
ax2 = fig.add_subplot(212)
ax2.plot(sess.ebins_out, spec, drawstyle='steps', label='Chandra MEG')
# add legends
ax.legend(loc=0)
ax2.legend(loc=0)
# zoom in on small sections of the spectrum
ax.set_xlim([1.150, 1.180])
ax2.set_xlim([1.150, 1.180])
# set axes
ax2.set_xlabel('Energy (keV)')
ax2.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
sess.set_broadening(True)
spec = sess.return_spectrum(kT)
spec = numpy.append(0, spec)
ax2.plot(sess.ebins_out, spec, drawstyle='steps', label='thermal')
sess.set_broadening(True, velocity_broadening=vel)
spec = sess.return_spectrum(kT)
spec = numpy.append(0, spec)
ax2.plot(sess.ebins_out, spec, drawstyle='steps', label='%.0f km/s'%(vel))
# adjust plot spacing so labels are visible
pylab.matplotlib.pyplot.subplots_adjust(hspace = 0.34)
# draw graphs
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_2_1.pdf')
fig.savefig('spectrum_session_examples_2_1.svg')
A kT=3.0keV spectrum unbroadened, thermally broadened and then additionally velocity broadend.¶
Changing Abundances¶
There are several ways to change the abundances. By default, all are set to 1.0 times the solar value of Anders and Grevesse.
The session.set_abund command implements this.
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
# set response
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# for plotting
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)
# make a spectrum at 1keV
s1 = sess.return_spectrum(1.0)
# change the abundance of Iron
sess.set_abund(26,2.0)
s2 = sess.return_spectrum(1.0)
# change the abundance of Magnesium and Silicon
sess.set_abund([12,14],0.5)
s3 = sess.return_spectrum(1.0)
# change the abundance of Neon and Iron separately
sess.set_abund([10,26],[2.0,0.5])
s4 = sess.return_spectrum(1.0)
# plot all this
offset = max(s1)
ax.plot(sess.ebins_out, numpy.append(0,s1), drawstyle='steps', label='Original')
ax.plot(sess.ebins_out, numpy.append(0,s2)+offset, drawstyle='steps', label='Fex2')
ax.plot(sess.ebins_out, numpy.append(0,s3)+offset*2, drawstyle='steps', label='Mg, Six0.5')
ax.plot(sess.ebins_out, numpy.append(0,s4)+offset*3, drawstyle='steps', label='Nex2, Fex0.5')
ax.legend(ncol=4, loc=0)
ax.set_xlim([0.8,1.5])
ax.set_xlabel('Energy (keV)')
ax.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_3_1.pdf')
fig.savefig('spectrum_session_examples_3_1.svg')
A kT=1.0keV spectrum with assorted different abundances applied.¶
You can also change the entire abundance set in use using session.set_abundset.
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
# set response
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# for plotting
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)
# get a list of the available abundance sets
sess.set_abundset()
# this prints them to screen.
for iabset, abset in enumerate(['Allen', 'AG89', 'GA88', 'Feldman', 'GA10', 'Lodd09', 'AE82']):
# set abundset
sess.set_abundset(abset)
# get spectrum at 1keV
s1 = sess.return_spectrum(1.0)
if iabset == 0:
baseoffset = max(s1)
offset = iabset*baseoffset
#plot
ax.plot(sess.ebins_out, numpy.append(0,s1)+offset, drawstyle='steps', label=abset)
ax.set_xlim([0.8,2.1])
ax.set_xlabel('Energy (keV)')
ax.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
ax.legend(loc=0, ncol=4)
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_4_1.pdf')
fig.savefig('spectrum_session_examples_4_1.svg')
A kT=1.0keV spectrum with assorted different abundance sets applied.¶
Return line list¶
To obtain a list of lines and their emissivities in a spectral range, use session.return_linelist. This
returns the data as a numpy array of the lines in the region, applying all current abundance information
to the linelist. It also interpolates in temperature between the two nearest points. Note
that since AtomDB cuts off lines with emissivities below \(10^{-20}\) \(ph ~ cm^3 s^{-1}\), lines
will diappear below this emissivity. Lines beyond the last temperature point at which they are tabulated will
not be included in the returned linelist.
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
# set response
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# get a linelist between 1.85 and 1.88 Angstroms for a kT=4keV plasma
llist = sess.return_linelist(4.0,[1.85,1.88])
# select only the strongest 10 lines
llist.sort(order='Epsilon')
llist = llist[-10:]
# show the column names
print(llist.dtype.names)
# expected output:
# ('Lambda', 'Lambda_Err', 'Epsilon', 'Epsilon_Err', 'Element', 'Ion', 'UpperLev', 'LowerLev')
# print the data
print(llist)
# expected output:
#[(1.8635 , nan, 8.8202775e-19, nan, 26, 24, 47, 1)
# (1.8532 , nan, 1.0113579e-18, nan, 26, 24, 10089, 6)
# (1.8622 , nan, 1.0681659e-18, nan, 26, 24, 10242, 3)
# (1.863 , nan, 2.5785468e-18, nan, 26, 24, 10247, 2)
# (1.8611 , nan, 2.8751725e-18, nan, 26, 24, 48, 1)
# (1.8659 , nan, 3.8414809e-18, nan, 26, 24, 10248, 3)
# (1.8554125, nan, 7.2097873e-18, nan, 26, 25, 6, 1)
# (1.8595169, nan, 7.7184372e-18, nan, 26, 25, 5, 1)
# (1.8681941, nan, 1.2001933e-17, nan, 26, 25, 2, 1)
# (1.8503995, nan, 3.5732330e-17, nan, 26, 25, 7, 1)]
# we can then do the same thing, but including the effective area of the instrument.
# So the Epsilon columnw will be multiplied by the effective area at the line's energy
llist = sess.return_linelist(4.0,[1.85,1.88], apply_aeff=True)
llist.sort(order='Epsilon')
llist = llist[-10:]
print('After apply effective area:')
print(llist)
#[(1.8635 , nan, 1.6238816e-18, nan, 26, 24, 47, 1)
# (1.8532 , nan, 1.7191134e-18, nan, 26, 24, 10089, 6)
# (1.8622 , nan, 1.9665764e-18, nan, 26, 24, 10242, 3)
# (1.863 , nan, 4.7473049e-18, nan, 26, 24, 10247, 2)
# (1.8611 , nan, 5.2934158e-18, nan, 26, 24, 48, 1)
# (1.8659 , nan, 7.3193854e-18, nan, 26, 24, 10248, 3)
# (1.8554125, nan, 1.2779108e-17, nan, 26, 25, 6, 1)
# (1.8595169, nan, 1.3680673e-17, nan, 26, 25, 5, 1)
# (1.8681941, nan, 2.2867946e-17, nan, 26, 25, 2, 1)
# (1.8503995, nan, 6.0738072e-17, nan, 26, 25, 7, 1)]
Return line emissivity¶
To calculate the emissvitiy of a specific line, you can use session.return_line_emissivity along
with the ion, element, upper and lower levels of the transition. You can supply a single or a range
of temperatures, and a dictionary will be returned containing much of the information along with
emissivity (epsilon) you requested.
import pyatomdb
import numpy
import pylab
# declare the Collisional Ionization Equilibrium session
sess = pyatomdb.spectrum.CIESession()
# set response
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
# get the emissivity of the resonance line of O VII at 5 different temperatures
kTlist = numpy.array([0.1,0.3,0.5,0.7,0.9])
Z=8
z1=7
up=7
lo = 1
ldata = sess.return_line_emissivity(kTlist, Z, z1, up, lo)
fig = pylab.figure()
fig.show()
ax= fig.add_subplot(111)
ax.semilogy(ldata['Te'], ldata['epsilon'])
ax.set_xlabel("Temperature (keV)")
ax.set_ylabel("Epsilon ph cm$^3$ s$^{-1}$")
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_6_1.pdf')
fig.savefig('spectrum_session_examples_6_1.svg')
Emissivity of the resonance line of O VII.¶
NEISession Class¶
Derived from the CIESession class, this handles non-equilibrium spectra. As such, all of the calls to it are exactly the same. Adding response, setting abundances obtaining spectra, etc all work the same way. Therefore I will only outline what is different.
The main difference is that a non equilibrium plasma requires 2 more pieces of information to define it. The current electron temperature is as with the CIE case. The other two are:
- tau
The ionization timescale (\(n_e t\)) in cm-3 s1 since the plasma left its previous equilibrium
- init_pop
The initial population the plasma was in at tau=0.
Tau is a single number in all cases. init_pop can be be defined in a range of ways:
- float
If it is a single number, this is assumed to be an electron temperature. The ionization fraction will be calculated at this \(T_e\).
- dict
If a dictionary is provided, it should be an explicit ionization fraction for each element. So dict[6]= [0.0, 0.1, 0.2, 0.3, 0.3, 0.1, 0.0] would be the ionization fraction for carbon, dict[8]= [0.0, 0.0, 0.1, 0.1, 0.1, 0.15,0.3,0.2,0.05] would be for oxygen etc.
- ‘ionizing’
If the string ionizing is provided, set all elements to be entirely neutral. This is the default.
- ‘recombining’
If the string recombining is provided, set all elements to be fully ionized.
Making a Spectrum¶
As an example, this will plot a simple recombining spectrum within initial temperature of 1.0keV.
import pyatomdb
import numpy
import pylab
# declare the Non-Equilibrium Ionization session. I am only using
# oxygen and neon for simplicity
sess = pyatomdb.spectrum.NEISession(elements=[8,10])
# set the response (raw keyword tells pyatomdb it is not a real response file)
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
kT = 0.4 # temperature in keV
tau = 1e10 # n_e * t in cm^-3 s
spec = sess.return_spectrum(kT, tau, init_pop=1.0)
# Returned spectrum has units of photons cm^5 s^-1 bin^-1
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)
offset = max(spec)
ax.plot(sess.ebins_out, numpy.append(0,spec), drawstyle='steps', label='mild recombining')
spec = sess.return_spectrum(kT, tau, init_pop='ionizing')
ax.plot(sess.ebins_out, offset+numpy.append(0,spec), drawstyle='steps', label='ionizing')
spec = sess.return_spectrum(kT, tau, init_pop='recombining')
ax.plot(sess.ebins_out, 2*offset+numpy.append(0,spec), drawstyle='steps', label='recombining')
# manually specified initial population
init_pop_dict = {}
init_pop_dict[8] = numpy.zeros(9) # 9 stages of oxygen (include neutral an fully stripped)
init_pop_dict[10] = numpy.zeros(11)
#put everything in the H-like stage
init_pop_dict[8][8] = 1.0
init_pop_dict[10][10] = 1.0
spec = sess.return_spectrum(kT, tau, init_pop=init_pop_dict)
ax.plot(sess.ebins_out, 3*offset+numpy.append(0,spec), drawstyle='steps', label='all H-like')
# label the figure
ax.set_xlabel('Energy (keV)')
ax.set_ylabel('Intensity (ph cm$^5$ s$^{-1}$ bin$^{-1}$)')
# zoom in a bit
ax.set_xlim([0.6,1.2])
ax.legend(loc=0)
# make plot appear
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('spectrum_session_examples_7_1.pdf')
fig.savefig('spectrum_session_examples_7_1.svg')
Recombining spectrum specified by an initial temperature¶
Making a Linelist¶
This is again exactly the same as the CIESession case, except with extra parameters tau and init_pop.
The only additional option is the ability to separate the linelist by driving ion. This is the ion which gives
rise to the emission. For example, if the emission is largely from ionization into an excited state and then
subsequent cascade, then that part of the emissivity has a different driving ion than that driven by
excitation of the ion’s ground state.
Separating out these features can be turned on and off using the by_ion_drv keyword,
import pyatomdb
import numpy
import pylab
# declare the Non-Equilibrium Ionization session. I am only using
# oxygen and neon for simplicity
sess = pyatomdb.spectrum.NEISession()
# set the response (raw keyword tells pyatomdb it is not a real response file)
sess.set_response('aciss_meg1_cy22.grmf', arf = 'aciss_meg1_cy22.garf')
kT = 4.0 # temperature in keV
tau = 1e11 # n_e * t in cm^-3 s
llist = sess.return_linelist(kT, tau, [1.7,1.9],init_pop=1.0)
# filter lines to strongest 20
llist.sort(order='Epsilon')
llist= llist[-20:]
print("strongest 20 lines")
print(llist.dtype.names)
print(llist)
llist = sess.return_linelist(kT, tau, [1.7,1.9],by_ion_drv=True,init_pop=1.0)
# filter lines to strongest 20
llist.sort(order='Epsilon')
llist= llist[-20:]
print("Strongest 20 lines, separated by driving ion")
print(llist.dtype.names)
print(llist)
# Expected output
# strongest 20 lines
# ('Lambda', 'Lambda_Err', 'Epsilon', 'Epsilon_Err', 'Element', 'Ion', 'UpperLev', 'LowerLev')
# [(1.8827 , nan, 4.3892664e-19, nan, 26, 22, 20009, 1)
# (1.8622 , nan, 4.5014698e-19, nan, 26, 24, 10242, 3)
# (1.8824 , nan, 4.8421002e-19, nan, 26, 22, 483, 1)
# (1.8825 , nan, 5.8088691e-19, nan, 26, 22, 482, 1)
# (1.8851 , nan, 5.9669826e-19, nan, 26, 22, 20008, 2)
# (1.8747 , nan, 6.6686722e-19, nan, 26, 24, 44, 1)
# (1.8737 , nan, 6.9597891e-19, nan, 26, 23, 20126, 3)
# (1.863 , nan, 1.0866564e-18, nan, 26, 24, 10247, 2)
# (1.8756001, nan, 1.0874650e-18, nan, 26, 23, 20122, 4)
# (1.8706 , nan, 1.3853237e-18, nan, 26, 24, 46, 1)
# (1.8738 , nan, 1.4878691e-18, nan, 26, 24, 45, 1)
# (1.8659 , nan, 1.6188786e-18, nan, 26, 24, 10248, 3)
# (1.857 , nan, 1.9076071e-18, nan, 26, 24, 50, 1)
# (1.8554125, nan, 2.3919818e-18, nan, 26, 25, 6, 1)
# (1.8595169, nan, 2.8080335e-18, nan, 26, 25, 5, 1)
# (1.8635 , nan, 3.6176967e-18, nan, 26, 24, 47, 1)
# (1.8704 , nan, 7.2666667e-18, nan, 26, 23, 309, 1)
# (1.8681941, nan, 8.8379641e-18, nan, 26, 25, 2, 1)
# (1.8611 , nan, 1.1892281e-17, nan, 26, 24, 48, 1)
# (1.8503995, nan, 1.4463351e-17, nan, 26, 25, 7, 1)]
# Strongest 20 lines, separated by driving ion
# ('Lambda', 'Lambda_Err', 'Epsilon', 'Epsilon_Err', 'Element', 'Elem_drv', 'Ion', 'Ion_drv', 'UpperLev', 'LowerLev')
# [(1.8622 , nan, 4.5014698e-19, nan, 26, 26, 24, 25, 10242, 3)
# (1.8824 , nan, 4.8421002e-19, nan, 26, 26, 22, 22, 483, 1)
# (1.8825 , nan, 5.6233666e-19, nan, 26, 26, 22, 22, 482, 1)
# (1.8851 , nan, 5.9669826e-19, nan, 26, 26, 22, 23, 20008, 2)
# (1.8747 , nan, 6.5364458e-19, nan, 26, 26, 24, 24, 44, 1)
# (1.8737 , nan, 6.9597891e-19, nan, 26, 26, 23, 24, 20126, 3)
# (1.863 , nan, 1.0866564e-18, nan, 26, 26, 24, 25, 10247, 2)
# (1.8756001, nan, 1.0874650e-18, nan, 26, 26, 23, 24, 20122, 4)
# (1.8706 , nan, 1.3512313e-18, nan, 26, 26, 24, 24, 46, 1)
# (1.8738 , nan, 1.4619381e-18, nan, 26, 26, 24, 24, 45, 1)
# (1.8659 , nan, 1.6188786e-18, nan, 26, 26, 24, 25, 10248, 3)
# (1.857 , nan, 1.8978882e-18, nan, 26, 26, 24, 24, 50, 1)
# (1.8554125, nan, 2.3303760e-18, nan, 26, 26, 25, 25, 6, 1)
# (1.8595169, nan, 2.7651213e-18, nan, 26, 26, 25, 25, 5, 1)
# (1.8681941, nan, 3.3107048e-18, nan, 26, 26, 25, 25, 2, 1)
# (1.8635 , nan, 3.6059715e-18, nan, 26, 26, 24, 24, 47, 1)
# (1.8681941, nan, 5.3853784e-18, nan, 26, 26, 25, 24, 2, 1)
# (1.8704 , nan, 7.2068723e-18, nan, 26, 26, 23, 23, 309, 1)
# (1.8611 , nan, 1.1865456e-17, nan, 26, 26, 24, 24, 48, 1)
# (1.8503995, nan, 1.4405425e-17, nan, 26, 26, 25, 25, 7, 1)]
#
Getting Atomic Data¶
The atomdb module is designed to read data from the database and either
return raw opened files (e.g. all the energy levels of an ion) or
individual useful data (e.g. a rate coefficient).
AtomDB stores data sorted by ion in a series of files in the $ATOMDB/APED. There are several files for
each ion, covering different data types:
‘IR’ - ionization and recombination
‘LV’ - energy levels
‘LA’ - radiative transition data (lambda and A-values)
‘EC’ - electron collision data
‘PC’ - proton collision data
‘DR’ - dielectronic recombination satellite line data
‘PI’ - XSTAR photoionization data
‘AI’ - autoionization data
‘ALL’ - reads all of the above. Does not return anything. Used for bulk downloading.
Or, for non-ion-specific data (abundances and bremstrahlung coeffts):
‘ABUND’ - abundance tables
‘HBREMS’ - Hummer bremstrahlung coefficients
‘RBREMS’ - relativistic bremstrahlung coefficitients
‘IONBAL’ - ionization balance tables
‘EIGEN’ - eigenvalue files (for non-equilibrium ionization)
The raw atomic data files can be downloaded using atomdb.get_data.
Note
Datacache
Often, when using AtomDB data you will want to use the same file over and over again. However,
the whole database is too large to be loaded into memory as a matter of course. To get around this,
the datacache keyword is used in several routines. This is a dictionary which is populated
with any data files that you have to open. The atomdb function which ultimately opens the raw
data files, get_data, will always check for an already existing data file here before
going to the disk to open the file anew. If there is no data there, then the file is opened
and its data is stored in the datacache so that next time the routine is called it will
already be there.
If you need to save memory, you can empty the cache by declaring the datacache as a
new dictionary, i.e. datacache={}.
import pyatomdb, time
# Example of using get data. I am going to get a range of different atomic data
#
# This routine returns the raw opened FITS data. No other processing is done.
Z = 8
z1 = 7
# get the energy levels of O VII
a = pyatomdb.atomdb.get_data(Z,z1,'LV')
# print out the first few lines
# note that for all the atomic data files, the useful data is in HDU 1.
print(a[1].data.names)
print(a[1].data[:5])
# get radiative transition information
b = pyatomdb.atomdb.get_data(Z,z1,'LA')
# find the transitions from level 7
bb = b[1].data[ (b[1].data['UPPER_LEV']==7)]
print()
print(bb.names)
print(bb)
# and so on.
# For non ion-specific data, Z and z1 are ignored.
abund = pyatomdb.atomdb.get_data(0,0,'ABUND')
print()
print(abund[1].data.names)
print(abund[1].data[:2])
# How the datacache works
# The datacache variable stores opened files in memory, preventing
# having to go and re-open them each time. This reduces disk access.
# make sure file is downloaded for a fair test:
a = pyatomdb.atomdb.get_data(Z,z1,'EC')
# case 1: no datacache
t1 = time.time()
for i in range(10):
a = pyatomdb.atomdb.get_data(Z,z1,'EC')
t2 = time.time()
# case 2: using datacache
datacache = {}
t3 = time.time()
for i in range(10):
a = pyatomdb.atomdb.get_data(Z,z1,'EC', datacache=datacache)
t4 = time.time()
print("Time without / with datacache: %f / %f seconds"%(t2-t1, t4-t3))
Getting Rate Coefficients¶
Ionization, recombination, and excitation rate coefficients can all be
obtained using get_maxwell_rate.
import pyatomdb, numpy, pylab
# We will get some maxwellian rates for a O 6+ 7->1
Z = 8
z1 = 7
up = 7
lo = 1
# excitation
Te = numpy.logspace(6, 7, 15)
# >>> # (1) Get excitation rates for row 12 of an Fe XVII file
# >>> colldata = pyatomdb.atomdb.get_data(26,17,'EC')
# >>> exc, dex = get_maxwell_rate(Te, colldata=colldata, index=12)
# >>> # (2) Get excitation rates for row 12 of an Fe XVII file
# >>> exc, dex = get_maxwell_rate(Te, Z=26,z1=17, index=12)
# >>> (3) Get excitation rates for transitions from level 1 to 15 of FE XVII
# >>> exc, dex = get_maxwell_rate(Te, Z=26, z1=17, dtype='EC', finallev=15, initlev=1)
datacache = {}
# get data by specifying ion, upper and lower levels
exc, dexc = pyatomdb.atomdb.get_maxwell_rate(Te, Z=Z, z1=z1, initlev = lo, finallev=up, dtype='EC', datacache=datacache)
fig= pylab.figure()
fig.show()
ax = fig.add_subplot(111)
ax.loglog(Te, exc, label='excitation')
ax.loglog(Te, dexc, label = 'deexcitation')
# preload data and find a specific transition
ecdat = pyatomdb.atomdb.get_data(8,7,'EC', datacache=datacache)
i = numpy.where( (ecdat[1].data['Upper_Lev']==up) &\
(ecdat[1].data['Lower_Lev']==lo))[0][0]
exc, dex = pyatomdb.atomdb.get_maxwell_rate(Te,colldata=ecdat, index=i, datacache=datacache)
ax.loglog(Te, exc, 'o', label='excitation')
ax.loglog(Te, dexc, 'o', label = 'deexcitation')
ax.legend(loc=0)
ax.set_xlabel("Temperature (K)")
ax.set_ylabel("Rate Coefficient (cm$^3$ s$^{-1}$)")
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('atomdb_examples_2_1.pdf')
fig.savefig('atomdb_examples_2_1.svg')
# you can also obtain ionization or recombination rates (see get_maxwell_rates writeup for options)
# in most cases you need to specify initlev and finallev as 1 to get the ion to ion rates.
ion= pyatomdb.atomdb.get_maxwell_rate(Te,Z=Z, z1=z1, initlev=1, finallev=1, dtype='CI', datacache=datacache)
print(ion)
#
# However it is recommended that you use get_ionrec_rate instead - it gives you everythign in one go
# combining the different types
ion, rec = pyatomdb.atomdb.get_ionrec_rate(Te, Z=Z, z1=z1, datacache=datacache)
# as separate entities
CI, EA, RR, DR=pyatomdb.atomdb.get_ionrec_rate(Te, Z=Z, z1=z1, datacache=datacache, separate=True)
print("Collisional ionization: ",CI)
print("Excitation Autoionization: ",EA)
print("Radiative Recombination: ",RR)
print("Dielectroni Recombination: ",DR)
ax.cla()
ax.loglog(Te, ion, label='Ionization')
ax.loglog(Te, rec, label='Recombination')
ax.loglog(Te, CI, '--',label='CI')
ax.loglog(Te, EA, '--',label='EA')
ax.loglog(Te, RR, '--',label='RR')
ax.loglog(Te, DR, '--',label='DR')
ax.legend(loc=0)
ax.set_xlabel("Temperature (K)")
ax.set_ylabel("Rate Coefficient (cm$^3$ s$^{-1}$)")
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('atomdb_examples_2_2.pdf')
fig.savefig('atomdb_examples_2_2.svg')
#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, ladat=False):
Excitation rates for O6+.¶
Ionization and recombination rates for O6+ to O7+. Note that for the components breakdown, there is no excitation autionization contribution in the files.¶
Calculating Plasma Conditions¶
The apec module is used for calculating plasma related phenomena. For example,
ionization fractions, non-equilibrium ionization, level populations, emission spectra.
Currently, it is mostly written in such a way that it is difficult to use it outside
of running the entire APEC code. A rewrite will happen soon. However, for now,
here are a couple of useful routines which should be useable.
Getting Charge State Distribution¶
import pyatomdb, numpy, pylab
# Calculate an equilibrium ionization balance for an element
kT = 0.5
Z = 8
# going to declare a datacache to speed up the calculations as we
# will be repeating things
datacache = {}
ionfrac_fast = pyatomdb.apec.return_ionbal(Z, kT, teunit='keV', datacache=datacache)
# The "fast" keyword makes the code open up a precalculated set of
# ionization balance data, a grid os 1251 values from 10^4K <= T <= 10^9K.
# By default it is true. Here, we turn it off for comparison
ionfrac_slow = pyatomdb.apec.return_ionbal(Z, kT, teunit='keV', datacache=datacache,
fast=False)
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)
ax.semilogy(range(0,Z+1), ionfrac_fast, 'o', label='fast')
ax.semilogy(range(0,Z+1), ionfrac_slow, '^', label='slow')
ax.set_ylim([1e-6,1])
ax.set_xlabel('Ion Charge')
ax.set_ylabel('Fraction')
ax.legend(loc=0)
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('apec_examples_1_1.pdf')
fig.savefig('apec_examples_1_1.svg')
# now to do the same for several different ways of defining non-equilbrium:
tau = 1e11
# 1: from neutral
init_pop='ionizing' # this is the default
nei1 = pyatomdb.apec.return_ionbal(Z, kT, tau = tau, init_pop = init_pop,
teunit='keV', datacache=datacache)
# 2: from fully stripped
init_pop='recombining'
nei2 = pyatomdb.apec.return_ionbal(Z, kT, tau = tau, init_pop = init_pop,
teunit='keV', datacache=datacache)
# 3: from a given Temperature
init_pop=0.1
nei3 = pyatomdb.apec.return_ionbal(Z, kT, tau = tau, init_pop = init_pop,
teunit='keV', datacache=datacache)
# 4: from a given array
init_pop=numpy.array([0.1, #O0+
0.1, #O1+
0.1, #O2+
0.1, #etc
0.3,
0.1,
0.1,
0.1,
0.0])
nei4 = pyatomdb.apec.return_ionbal(Z, kT, tau = tau, init_pop = init_pop,
teunit='keV', datacache=datacache)
# 5: from a dict of arrays
# (this is useful if you have lots of elements to do)
init_pop = {}
init_pop[8] = numpy.array([0.5, #O0+
0.0, #O1+
0.0, #O2+
0.0, #etc
0.0,
0.0,
0.0,
0.0,
0.5])
nei5 = pyatomdb.apec.return_ionbal(Z, kT, tau = tau, init_pop = init_pop,
teunit='keV', datacache=datacache)
ax.cla()
ax.semilogy(range(Z+1), nei1, 'o', label='ionizing')
ax.semilogy(range(Z+1), nei2, '<', label='recombining')
ax.semilogy(range(Z+1), nei3, 's', label='0.1keV')
ax.semilogy(range(Z+1), nei4, '^', label='array')
ax.semilogy(range(Z+1), nei5, 'x', label='dict')
ax.set_ylim([1e-6,1])
ax.set_xlabel('Ion Charge')
ax.set_ylabel('Fraction')
ax.legend(loc=0)
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('apec_examples_1_2.pdf')
fig.savefig('apec_examples_1_2.svg')
Ionization balance for oxygen¶
Charge state distribution calculated many different ways for non equilibrium ionization.¶
Including APEC models in PyXSPEC¶
from apec_xspec import *
# declare a new model
# inital import creates pyapec, pyvapec, pyvvapec, analagous to apec, vapec, vvapec
m1 = xspec.Model('pyapec')
m1.show()
m1.pyapec.kT=5.0
# let's plot a spectrum
xspec.Plot.device='/xs'
xspec.Plot('model')
# you can fiddle with some aspects of the model directly
cie.set_eebrems(False) # turn off electron-electron bremsstrahlung
xspec.Plot('model')
# it is also possible to go in and change broadening, etc. The interface can be
# adjust quite easily by looking at apec_xspec.py to add extra parameters
There are other wrapper files which work in the same way for a model with
a separate temperature set for line broadening, and others will be added in
the wrapper directory.
Individual Use Cases¶
These are typically scripts created to answer user’s questions which might be interesting. As a result, sometimes the exact response files etc may not be available to you. Please swap in what you need to.
Get PI Cross Sections¶
Extract the PI cross section data: photoionization_data.py
import pyatomdb, numpy,os, pylab
try:
import astropy.io.fits as pyfits
except:
import pyfits
# This is a sample routine that reads in the photoionization data
# It also demonstrates using get_data, which should download the data you
# need automagically from the AtomDB site.
#
# It also shows how to get the raw XSTAR PI cross sections.
# going to get PI cross section from iron 16+ to 17+ (Fe XVII-XVIII)
Z = 26
z1 = 17
# get the AtomDB level data
lvdata = pyatomdb.atomdb.get_data(Z, z1, 'LV')
# get the XSTAR PI data from AtomDB
pidata = pyatomdb.atomdb.get_data(Z, z1, 'PI')
# set up the figure
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)
# to calculate the cross section (in cm^2) at a given single energy E (in keV)
# does not currently work with vector input, so have to call in a loop if you
# want multiple energies [I will fix this]
E = 10.
# get the ground level (the 0th entry in LV file) data
lvd = lvdata[1].data[0]
# This is the syntax for calculating the PI cross section of a given line
# This will work for non XSTAR data too.
sigma = pyatomdb.atomdb.sigma_photoion(E, Z, z1,lvd['phot_type'], lvd['phot_par'], \
xstardata=pidata, xstarfinallev=1)
# To get the raw XSTAR cross sections (units: energy = keV, cross sections = Mb)
# for level 1 -> 1 (ground to ground)
pixsec = pyatomdb.atomdb.__sort_pi_data(pidata, 1,1)
ax.loglog(pixsec['energy'], pixsec['pi_param']*1e-18, label='raw xstar data')
# label the plot
ax.set_title('Plotting raw XSTAR PI cross sections. Fe XVII gnd to Fe XVIII gnd')
ax.set_xlabel("Energy (keV)")
ax.set_ylabel("PI cross section (cm$^{2}$)")
pylab.draw()
zzz=input('press enter to continue')
Make Cooling Curve¶
Make a cooling curve, total emissivity in keV cm3 s-1, for each element in a specfied spectral range (e.g. 2 to 10 keV).
import pyatomdb, numpy, os, pylab
"""
This code is an example of generating a cooling curve: the total power
radiated in keV cm3 s-1 by each element at each temperature. It will
generate a text file with the emission per element at each temperature
from 1e4 to 1e9K.
This is similar to the atomdb.lorentz_power function, but with a few
normalizations removed to run a little quicker.
Note that the Anders and Grevesse (1989) abundances are built in to
this. These can be looked up using atomdb.get_abundance(abundset='AG89'),
or the 'angr' column of the table at
https://heasarc.nasa.gov/xanadu/xspec/xspec11/manual/node33.html#SECTION00631000000000000000
Adjustable parameters (energy range, element choice) are in the block
marked #### ADJUST THINGS HERE
Note that any warnings of the nature:
"kT = 0.000862 is below minimum range of 0.000862. Returning lowest kT spectrum available"
should be ignored. This is returning the temperature at our lowest tabulated value but is
a rounding error making it think it is outside our range.
Usage: python3 calc_power.py
"""
def calc_power(Zlist, cie, Tlist):
"""
Zlist : [int]
List of element nuclear charges
cie : CIESession
The CIEsession object with all the relevant data predefined.
Tlist : array(float)
The temperatures at which to calculate the power (in K)
"""
res = {}
res['power'] = {}
res['temperature'] = []
cie.set_abund(numpy.arange(1,31), 0.0)
kTlist = Tlist*pyatomdb.const.KBOLTZ
en = (cie.ebins_out[1:]+cie.ebins_out[:-1])/2
for i, kT in enumerate(kTlist):
print("Doing temperature iteration %i of %i"%(i, len(kTlist)))
T = Tlist[i]
res['temperature'].append(T)
res['power'][i] = {}
for Z in Zlist:
if Z==0:
# This is the electron-electron bremstrahlung component alone
#set all abundances to 1 (I need a full census of electrons in the plasma for e-e brems)
cie.set_abund(Zlist[1:], 1.0)
# turn on e-e bremsstrahlung
cie.set_eebrems(True)
spec = cie.return_spectrum(kT, dolines=False, docont=False, dopseudo=False)
else:
# This is everything else, element by element.
# turn off all the elements
cie.set_abund(Zlist[1:], 0.0)
# turn back on this element
cie.set_abund(Z, 1.0)
# turn off e-e bremsstrahlung (avoid double counting)
cie.set_eebrems(False)
spec = cie.return_spectrum(kT)
# if Z = 1, do the eebrems (only want to calculate this once)
# if Z == 1:
# cie.set_eebrems(True)
# else:
# cie.set_eebrems(False)
# get spectrum in ph cm3 s-1
#spec = cie.return_spectrum(kT)
# convert to keV cm3 s-1, sum
res['power'][i][Z] = sum(spec*en)
return res
if __name__=='__main__':
############################
#### ADJUST THINGS HERE ####
############################
# Elements to include
#Zlist = range(31) <- all the elements
Zlist = [0,1,2,6,7,8,10,12,13,14,16,18,20,26,28] #<- just a few
# Note that for this purpose, Z=0 is the electron-electron bremsstrahlung
# continuum. This is not a general AtomDB convention, just what I've done here
# to make this work.
# specify photon energy range you want to integrate over (min = 0.001keV, max=100keV)
Elo = 0.001 #keV
Ehi = 100.0 #
# temperatures at which to calculate curve (K)
Tlist = numpy.logspace(4,9,51)
# specify output file name (default output.txt)
outfile = 'output.txt'
################################
#### END ADJUST THINGS HERE ####
################################
# set up the spectrum
cie = pyatomdb.spectrum.CIESession()
ebins = numpy.linspace(Elo, Ehi, 10001)
cie.set_response(ebins, raw=True)
cie.set_eebrems(True)
# crunch the numbers
k = calc_power(Zlist, cie, Tlist)
# output generation
o = open(outfile, 'w')
# header row
s = '# Temperature log10(K)'
for i in range(len(Zlist)):
s += ' %12i'%(Zlist[i])
o.write(s+'\n')
# for each temperature
k['totpower'] = numpy.zeros(len(k['temperature']))
for i in range(len(k['temperature'])):
s = '%22e'%(numpy.log10(k['temperature'][i]))
for Z in Zlist:
s+=' %12e'%(k['power'][i][Z])
k['totpower'][i]+=k['power'][i][Z]
o.write(s+'\n')
# notes
o.write("# Total Emissivity in keV cm^3 s^-1 for each element with AG89 abundances, between %e and %e keV\n"%(Elo, Ehi))
o.write("# To get cooling power, multiply by Ne NH")
o.write("# Z=0 component is electron-electron bremsstrahlung component, only significant at high T")
o.close
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)
ax.loglog(k['temperature'], k['totpower']*pyatomdb.const.ERG_KEV)
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Radiated Power (erg cm$^3$ s$^{-1}$)')
ax.set_xlim(min(Tlist), max(Tlist))
# draw graphs
pylab.draw()
zzz=input("Press enter to continue")
# save image files
fig.savefig('calc_power_examples_1_1.pdf')
fig.savefig('calc_power_examples_1_1.svg')
The total emissivity of the plasma between 0.001 and 100 keV.¶