Update population files

The warmabs model uses pre-calculated ion population and auxiliary files for efficient computation of synthetic spectra. Default popoulation files are calculated for an incident spectral energy distribution (SED) of a powerlaw with index 2. For any other SEDs, the population files should be re-calculated to ensure accurate results. This example illustrates the calculation of population files from and exisiting xspec SED model.

1. Tabulate the SED in a text file

In this exmaple we use the pyxspec to tabulate the comptt model. For a given source, the incident SED will have to be obtained by fitting an appropriate unabsorbed continuum model.

The first line of the text file must be the number of energies listed in the table. The remaining lines are the energy channel (in eV) and the flux in units of \(\mathsf{photons}\,\mathsf{cm}^{-2}\,\mathsf{s}^{-1}\,\mathsf{erg}^{-1}\) or \(\mathsf{erg}\,\mathsf{cm}^{-2}\,\mathsf{s}^{-1}\,\mathsf{erg}^{-1}\). Here we use the latter (default units). Note that XSTAR will appropriately renormalize the luminosity.

>>> import xspec
>>> xspec.Model("comptt")

========================================================================
Model compTT<1> Source No.: 1   Active/Off
Model Model Component  Parameter  Unit     Value
par  comp
1    1   compTT     Redshift            0.0          frozen
2    1   compTT     T0         keV      0.100000     +/-  0.0
3    1   compTT     kT         keV      50.0000      +/-  0.0
4    1   compTT     taup                1.00000      +/-  0.0
5    1   compTT     approx              1.00000      frozen
6    1   compTT     norm                1.00000      +/-  0.0
________________________________________________________________________

>>> xspec.AllData.dummyrsp(.01,100.,1000,"log")
>>> xspec.plot.device = "/null"
>>> xspec.Plot("model")
>>> yVals= xspec.Plot.model(1) # Photon bin density
>>> xVals = xspec.Plot.x() # Bin center in keV

Inspect the photon spectrum.

>>> import matplotlib.pyplot as plt
>>> plt.plot(xVals, yVals)
>>> plt.xscale("log")
>>> plt.yscale("log")
>>> plt.xlabel("Energy (keV)")
>>> plt.ylabel("phs/keV")
>>> plt.show()
../_images/ised-comptt.png

Write spectrum to file.

>>> of = open("comptt.txt", "w")
>>> of.write("%d\n" % len(yVals)) # First line contains number of energy bins
>>> for ii in range(0, len(yVals)):
>>>    # energy in eV, flux in erg/cm^2/s/erg
>>>    of.write("%.4E\t%.4E\n" % (xVals[ii]*1e3, yVals[ii]*xVals[ii]))
>>> of.close()

unix> head comptt.txt
1000
1.0046E+01   6.9607E-03
1.0139E+01   7.1498E-03
1.0233E+01   7.3440E-03
1.0328E+01   7.5434E-03
1.0423E+01   7.7481E-03
1.0520E+01   7.9584E-03
1.0617E+01   8.1742E-03
1.0715E+01   8.3959E-03
1.0814E+01   8.6235E-03

2. Run xstar with this incident SED

We start in a fresh directory to prevent previous population files from being overwritten.

unix> mkdir popfiles
unix> cd popfiles
unix> mv <path/to/sed>/comptt.txt .

Now run xstar.

unix> xstar \
cfrac=1.0 \
temperature=1.0E+04 \
lcpres=0 \
pressure=3.0E-2 \
density=1.0E+4 \
spectrum='file' \
spectrum_file='comptt.txt' \
spectun=0 \
trad=-1. \
rlrad38=10.E-12 \
column=10.E+16 \
rlogxi=5.0 \
nsteps=20 \
niter=99 \
lwrite=1 \
lprint=1 \
lstep=0 \
habund=1.0E+00 \
heabund=1.0E+00 \
liabund=0.0E+00 \
beabund=0.0E+00 \
babund=0.0E+00 \
cabund=1.0E+00 \
nabund=1.0E+00 \
oabund=1.0E+00 \
fabund=1.0E+00 \
neabund=1.0E+00 \
naabund=1.0E+00 \
mgabund=1.0E+00 \
alabund=1.0E+00 \
siabund=1.0E+00 \
pabund=1.0E+00 \
sabund=1.0E+00 \
clabund=1.0E+00 \
arabund=1.0E+00 \
kabund=1.0E+00 \
caabund=1.0E+00 \
scabund=1.0E+00 \
tiabund=1.0E+00 \
vabund=1.0E+00 \
crabund=1.0E+00 \
mnabund=1.0E+00 \
feabund=1.0E+00 \
coabund=1.0E+00 \
niabund=1.0E+00 \
cuabund=0.0E+00 \
znabund=0.0E+00 \
emult=5.0E-01 \
taumax=1.5E+01 \
xeemin=4.0E-02 \
critf=1.0E-06 \
vturbi=1.0E+00 \
npass=1.0E+00 \
modelname='custom_pops' \
loopcontrol=0

xstar version 2.59b
Loading Atomic Database...
initializng database...

pass number=           1          -1
log(r) delr/r log(N) log(xi) x_e   log(n) log(t) h-c(%) h-c(%) log(tau)
                                                              fwd    rev
9.00 -36.00 -10.00   5.00   1.21   4.00   7.79  -0.30   0.00 -10.00 -10.00  8
9.00 -36.00 -10.00   5.00   1.21   4.00   7.79  -0.30  -0.00 -10.00 -10.00  1
9.02  -1.32  11.70   4.96   1.21   4.00   7.79  -0.66  -0.00 -10.00 -10.00  1
9.04  -1.03  12.01   4.92   1.21   4.00   7.78  -0.02  -0.00 -10.00 -10.00  3
9.06  -0.87  12.20   4.87   1.21   4.00   7.78  -0.46  -0.00 -10.00 -10.00  1
9.08  -0.75  12.33   4.83   1.21   4.00   7.78  -0.94  -0.00 -10.00 -10.00  1
[...]

Upon completion, the following files will be created:

unix> ls

comptt.txt
xo01_detal2.fits
xo01_detal4.fits
xout_cont1.fits
xout_rrc1.fits
xout_step.log
xo01_detail.fits
xo01_detal3.fits
xout_abund1.fits
xout_lines1.fits
xout_spect1.fits
xstar.par

Note that some of these files can be several GB in size.