5. Plotting the Fitted Spectra Outside SPEX

After restoring the data file, in order to plot the spectra, we first need to compute the photon flux, error bars, background flux, model flux and so on for a specific interval.

IDL>interval = 0
 

To compute the photon flux (=(obsi-backi)/convi) :

IDL> phflux = (obsi-backi)/convi
 

To compute the error of the photon flux (=eobsi/convi):

IDL>errflux = sqrt(eobsi^2 + ebacki^2) /convi
 

To compute the background flux (=backi/convi):

IDL> backflux=backi/convi
IDL> help, phflux, errflux, backflux
PHFLUX FLOAT = Array[82, 37]
ERRFLUX FLOAT = Array[82, 37]
BACKFLUX FLOAT = Array[82, 37]
 

Since the “edges†is the edges of energy channel, we need to compute the mean value of each channel:

IDL> print, edges[*,0:2]
3.00000 4.00000
4.00000 5.00000

IDL> edge_products, edges, mean = eph
IDL> print, eph[0:2]
3.50000 4.50000 5.50000
 

So far, we have all the spectral data for plots. We also need the f_vth_bpow model flux so that we can overplot on the data like SPEX does. By using xdoc introduced in "Spectroscopy - First Steps", we know f_vth_pow model uses f_vth.pro and f_bpow.pro to compute the thermal and nonthermal photon flux. The total photon flux is the sum of the two fluxes. Let us compute the model flux for the first interval using the saved fitting parameters from apar_arr.

First, we get the fitting parameters for the first spectral fit:

IDL> fitpara=apar_arr[*, interval]
IDL> print, fitpara
0.756227 3.20582 9.46661 3.31592 75.3729 3.83593
 

Secondly, we can check the UT time of that interval to see whether they are correct by comparing the time history window:

IDL> uttime = anytim(xselect[*, interval] + utbase, /ecs)
IDL> print, uttime
2002/07/23 00:27:24.000 2002/07/23 00:27:43.999
 

Thirdly, compute the thermal photon flux:

IDL> thermflux = f_vth(edges, fitpara[0:1])
IDL> help, thermflux
THERMFLUX FLOAT = Array[82]
 

Finally, we compute the nonthermal photon flux. But we should do something beforehand. There is always a low energy cut-off for the power-law photon distribution. From Spectroscopy - First Steps, we know the default is 1 keV cut-off, and spectral index below 1 keV is -1.5. The nonthermal f_ bpow model uses these two parameters, but in the saved fitting parameters, they are not saved. So we have to add them at the end of the variable array "fitpara":

IDL> a_cutoff = [1, 1.5]
IDL> fitpara = [fitpara, a_cutoff]
IDL> print, fitpara
0.756227 3.20582 9.46661 3.31592 75.3729 3.83593 1.00000 1.50000

IDL> nonthermflux = f_bpow(eph, fitpara[2:7])
IDL> help, nonthermflux
NONTHERMFLUX FLOAT = Array[82]
 

Now we have all the data we want to plot the spectrum of the first interval. The spectrum can be plotted as follows:

IDL> loadct, 4
IDL> window, 10, xsize = 600, ysize = 500
IDL> plot, eph, thermflux, /xlog, /ylog, xtitle = 'Photon Energy (keV)', ytitle = 'Flux (photons s!u-1!n cm!u-2!n keV!u-1!n)', xrange = [10.0, 300.], yrange=[1.0e-3, 1.0e6], xstyle = 1, ystyle = 1, /nodata
IDL> oplot, eph, thermflux, lineStyle = 1, color = 47, thick = 4 ;plot thermal flux
IDL> oplot, eph, nonthermFlux, lineStyle = 2, color = 110, thick = 4 ;plot nonthermal flux
IDL> oplot, eph, thermflux + nonthermflux, linestyle = 0, color =255 , thick = 4 ;plot total flux
IDL> oplot, eph, phflux[*, interval], psym = 1, color = 144, symsize = 1, thick = 2 ;plot spectrum data
IDL> oplot_err, eph, phflux[*, interval], yerr = errFlux[*,interval], psym = 1, symsize = 1, color = 144 ;plot spectrum data error bar
IDL> oplot, eph, backflux[*, interval], psym = 5, symsize = 0.5, color = 208 ;plot the background flux
 

Here is an IDL procedure that does all of the above steps. You can use it directly by cutting and pasting it into an IDL/SSW session.

The following is the spectrum we just plotted:

This page last updated: June 27, 2011