Spectrum Calculation with ATOMDB
Data in the ATOMDB refers can be used to calculate the spectrum emitted from a optically-thin collisionally ionized plasma using a number of tools. Some sample scripts (using various languages) are shown here.
Spectrum Calculation with ATOMDB in Sherpa
Sherpa can use the ATOMDB in a number of ways, and can generate spectra from it in two ways. The first method is to use the xsapec model, as shown here used with an linear grid from 10 eV to 10 keV to calculate the spectrum from a 2 keV plasma:
sherpa> dataspace (0.01:10:0.01) histogram
sherpa> source = xsapec[a1]
a1.kT parameter value [1] 2.0
a1.Abundanc parameter value [1] 1.0
a1.Redshift parameter value [0] 0.0
a1.norm parameter value [1] 1.0
sherpa> set plot lin log
sherpa> lplot source
which gives the following result: (click for larger version)
In addition to this method, it is also possible to use an S-lang routine to calculate a spectrum. This code uses the importable version of ISIS and is available here. Currently, it runs only in wavelength mode, and so the dataspace must be given in Angstroms, not keV. It can be used in the following manner:
sherpa> import("isis")
Imported ISIS module version X.X
sherpa> plasma(aped)
Read 129 energy level files
Read 129 wavelength files
Tables have 163598 lines between 1.198 and 2.381e+07 Angstrom
Initializing: emissivity data
Scanning: line emissivity tables [50 hdus]
hdu: 50/50
Tables have 169354 lines between 1.198 and 2.381e+07 Angstrom
Loading: line emissivity tables [50 hdus]
hdu: 50/50
sherpa> evalfile("apec.sl")
1
sherpa> dataspace (1:20:0.01) histogram
sherpa> source = apec[a1]
a1.kT parameter value [1] 2.0
a1.abund parameter value [1] 1.0
a1.norm parameter value [1] 1.0
sherpa> set plot lin log
sherpa> lp source
which gives the following result: (click for larger version)
Spectrum Calculation with ATOMDB in XSPEC
XSPEC version 11.0 includes a number of optically-thin collisional plasma models, including apec (which uses the ATOMDB data), raymond (an implementation of the Raymond-Smith code), and mekal (the Mewe-Kaastra-Liedahl code). To use the apec model in XSPEC with an linear grid from 10 eV to 10 keV to calculate a 2 keV spectrum, use the following commands:
XSPEC> dummyrsp 0.01 10 1000
XSPEC> model apec
Input parameter value, delta, min, bot, top, and max values for ...
Current: 1 0.01 0.008 0.008 64 64
apec:kT> 2.0
Current: 1 -0.001 0 0 5 5
apec:Abundanc> 1.0
Current: 0 -0.001 0 0 2 2
apec:Redshift> 1.0
Current: 1 0.01 0 0 1E+24 1E+24
apec:norm> 1.0
---------------------------------------------------------------------------
---------------------------------------------------------------------------
Model: apec[1]
Model Fit Model Component Parameter Unit Value
par par comp
1 1 1 apec kT keV 2.000 +/- 0.
2 2 1 apec Abundanc 1.000 frozen
3 3 1 apec Redshift 0. frozen
4 4 1 apec norm 1.000 +/- 0.
---------------------------------------------------------------------------
---------------------------------------------------------------------------
XSPEC> plot model
At this point, we get the following result (click for larger view):
Spectrum Calculation with ATOMDB in IDL
A number of IDL scripts have been written to read and use data from the ATOMDB. The demonstration that follows assumes that these files have been downloaded and installed; directions can be found here.
IDL> read_linelist,'APEC_v1.3.2_line.fits',line,t1,n1,l1
IDL> read_coco,'APEC_v1.3.2_coco.fits',coco,t3,n3,l3
IDL> Ebin = 0.01 + findgen(1000)*(10.0 - 0.01)/999.
IDL> spectrum = calc_spectrum(line,coco,Ebin,2.0/8.617e-8)
IDL> plot_io, Ebin, spectrum
At this point, we get the following result (click for larger view):
Spectrum Calculation with ATOMDB in ISIS
The ISIS package has been written by the CXC MIT team to analyze grating data using data in the ATOMDB. This sample script shows how to use ISIS with the ATOMDB to calculate the same 2 keV plasma model as shown in the previous examples:
isis> plasma(aped);
isis> state = default_plasma_state();
isis> state.temperature = 2.0/8.617e-8;
isis> define_model([state]);
isis> binlo = [1:20:0.01];
isis> binhi = binlo + 0.01;
isis> spectrum = model_spectrum(binlo, binhi);
isis> plot(binlo, spectrum);
At this point, we get the following result (click for larger view):