#!/usr/bin/env python
import os
import numpy
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
import matplotlib.gridspec as gridspec
from astropy import units as u
from astropy.table import QTable, Table
from astropy.visualization import quantity_support
quantity_support()
from . import pyspex_f2py
from . import model
[docs]class PlotModel:
"""Class containing the model spectra for plotting purposes.
:ivar nsector: Number of sectors available.
:vartype nsector: int
:ivar sector: List of plot objects for each sector.
:vartype sector: list of objects
"""
def __init__(self):
self.nsector = 0 # Number of sectors available.
self.sector = [] # List of plot objects for each sector.
[docs] def update(self, wave=False):
"""Update the model spectrum information for all sectors."""
pyspex_f2py.api_model_sector.api_sector_update()
self.nsector = int(pyspex_f2py.api_model_sector.sect_nsector)
self.sector = []
for i in numpy.arange(self.nsector):
pls = PlotModelSector()
pls.update(i + 1, wave)
self.sector.append(pls)
[docs] def plot(self, xlog=False, ylog=False, wave=False, title='SPEX Model Spectrum', show=True):
"""Plot the model spectrum for all sectors.
:param xlog: (Optional) Set the X-axis to be logarithmic.
:type xlog: bool
:param ylog: (Optional) Set the Y-axis to be logarithmic.
:type ylog: bool
:param wave: (Optional) Plot in wavelength (Angstrom)
:type wave: bool
:param title: (Optional) Set the title of the plot.
:type title: str
:param show: (Optional) Show the plot (True) or return the plot object (False).
:type show: bool
:return: Matplotlib plt object.
:rtype: matplotlib.pyplot
"""
self.update(wave)
fig, ax = plt.subplots()
for i in numpy.arange(self.nsector):
plt.step(self.sector[i].tabmodel['X_ctr'], self.sector[i].tabmodel['Model'], where='pre', c='r')
if xlog:
plt.xscale('log')
formatter = FuncFormatter(lambda y, _: '{:.16g}'.format(y))
ax.xaxis.set_major_formatter(formatter)
if ylog:
plt.yscale('log')
if wave:
plt.xlabel('Wavelength ({0})'.format(self.sector[0].tabmodel['X_ctr'].unit.to_string('latex_inline')))
else:
plt.xlabel('Energy ({0})'.format(self.sector[0].tabmodel['X_ctr'].unit.to_string('latex_inline')))
plt.ylabel('Flux ({0})'.format(self.sector[0].tabmodel['Model'].unit.to_string('latex_inline')))
plt.title(title)
if show:
plt.show()
else:
return plt
[docs] def adum_xspec_qdp(self, isect, qdpfile, overwrite=False):
"""Write the model to a QDP file in Xspec output format. This is useful to use SPEX models as input
for various simulation programs, like SIMX, SIXTE and HEASIM.
:param isect: Sector number to write (Usually 1).
:type isect: int
:param qdpfile: File name of the QDP file to write.
:type qdpfile: str
:param overwrite: Overwrite existing files (default=False).
:type overwrite: bool
"""
self.update(wave=False)
if isect > self.nsector or isect <= 0:
return 1
# Get the bin centers
xc = self.sector[isect-1].xc
# Calculate the 'error' in energy (e.g. half the bin width)
dx = self.sector[isect-1].xb - self.sector[isect-1].xc
# Calculate the flux in Photons / cm**2 / s / keV
m = self.sector[isect-1].m.to(u.ph / (u.cm ** 2 * u.s * u.keV))
xspectable = Table([xc, dx, m], names=('X_ctr', 'X_err', 'Model'),
meta={'name': 'Model spectrum (Xspec format)'})
xspectable.write(qdpfile, format='ascii.qdp', overwrite=overwrite)
return 0
[docs]class PlotModelSector:
"""Class containing the spectral model for a particular sector for plotting purposes.
:ivar xc: Central bin energy (keV).
:vartype xc: astropy.units.quantity.Quantity
:ivar xb: Upper boundary of bin energy (keV).
:vartype xb: astropy.units.quantity.Quantity
:ivar m: Model spectrum (Photons/m**2/s/keV).
:vartype m: astropy.units.quantity.Quantity
:ivar n: Number of bins.
:vartype n: int
:ivar tabmodel: Astropy QTable containing the model spectrum
:vartype tabmodel: astropy.table.table.QTable
"""
def __init__(self):
self.tabmodel = QTable()
self.xc = numpy.array([0.], dtype=object) * u.keV # Central bin energy
self.xb = numpy.array([0.], dtype=object) * u.keV # Upper boundary of bin energy
self.m = numpy.array([0.], dtype=object) * u.ph / (u.m ** 2 * u.s * u.keV) # Model spectrum
self.n = 0 # Number of bins
[docs] def update(self, isect, wave=False):
"""Update the values of the model spectrum for this sector.
:param isect: Sector number.
:type isect: int
"""
pyspex_f2py.api_plot_model.plot_model_update(float(isect))
self.n = pyspex_f2py.api_plot_model.model_nbins
self.xc = numpy.zeros(self.n, dtype=object) * u.keV
self.xb = numpy.zeros(self.n, dtype=object) * u.keV
self.m = numpy.zeros(self.n, dtype=object) * u.ph / (u.m ** 2 * u.s * u.keV)
for i in numpy.arange(self.n):
self.xc[i] = float(pyspex_f2py.api_plot_model.model_xc[i]) * u.keV
self.xb[i] = float(pyspex_f2py.api_plot_model.model_xb[i]) * u.keV
self.m[i] = float(pyspex_f2py.api_plot_model.model_m[i]) * u.ph / (u.m ** 2 * u.s * u.keV)
if wave:
self._convert_to_ang()
self.tabmodel = QTable([self.xc, self.xb, self.m], names=('X_ctr', 'X_upp', 'Model'),
meta={'name': 'Model spectrum'})
pyspex_f2py.api_plot_model.plot_model_final()
def _convert_to_ang(self):
"""Convert the spectrum in keV to Angstrom."""
# Calculate the bin width in keV
delta_e = 2. * (self.xb - self.xc)
# Convert the individual arrays
self.xc = self.xc.to(u.Angstrom, equivalencies=u.spectral())
self.xb = self.xb - delta_e
self.xb = self.xb.to(u.Angstrom, equivalencies=u.spectral())
delta_w = 2. * (self.xb - self.xc)
self.m = self.m * delta_e / delta_w
# Reverse the order of the arrays
self.xc = numpy.flip(self.xc)
self.xb = numpy.flip(self.xb)
self.m = numpy.flip(self.m)
[docs]class PlotArea:
"""Class for plotting effective area.
:ivar ninst: Number of instruments.
:vartype ninst: int
:ivar inst: List of instrument objects.
:vartype inst: list of objects
"""
def __init__(self):
self.ninst = 0 # Number of instruments
self.inst = [] # List of instrument objects
[docs] def update(self, wave=False):
"""Update the effective area information for the plot."""
self.ninst = int(pyspex_f2py.api_instrument.api_instrument_n())
for i in numpy.arange(self.ninst):
ins = PlotAreaInst()
ins.update(i + 1, wave)
self.inst.append(ins)
[docs] def plot(self, xlog=False, ylog=False, wave=False, title='Effective Area', show=True):
"""Plot the effective area.
:param xlog: (Optional) Set the X-axis to be logarithmic.
:type xlog: bool
:param ylog: (Optional) Set the Y-axis to be logarithmic.
:type ylog: bool
:param wave: (Optional) Plot the effective area as a function of wavelength.
:type wave: bool
:param title: (Optional) Set the title of the plot.
:type title: str
:param show: (Optional) Show the plot (True) or return the plot object (False).
:type show: bool
:return: Matplotlib plt object.
:rtype: matplotlib.pyplot
"""
self.update(wave)
fig, ax = plt.subplots()
# The emin and emax retrieve the minimum and maximum energy to plot on the x-axis.
emin = 1.E+8
emax = 0.
for i in numpy.arange(self.ninst):
nreg = self.inst[i].nreg
for j in numpy.arange(nreg):
pyspex_f2py.api_region.api_reg_update(float(i + 1), float(j + 1))
rmin = float(pyspex_f2py.api_region.reg_edmin) * u.keV
emin = min(rmin.value, emin) * u.keV
rmax = float(pyspex_f2py.api_region.reg_edmax) * u.keV
emax = max(rmax.value, emax) * u.keV
plt.step(self.inst[i].reg[j].tabarea['X_ctr'], self.inst[i].reg[j].tabarea['Area'], where='pre', c='r')
if xlog:
plt.xscale('log')
formatter = FuncFormatter(lambda y, _: '{:.16g}'.format(y))
ax.xaxis.set_major_formatter(formatter)
if ylog:
plt.yscale('log')
if wave:
wmax = emin.to(u.Angstrom, equivalencies=u.spectral())
wmin = emax.to(u.Angstrom, equivalencies=u.spectral())
plt.xlim(wmin.value, wmax.value)
plt.xlabel('Wavelength ({0})'.format(self.inst[0].reg[0].tabarea['X_ctr'].unit.to_string('latex_inline')))
else:
plt.xlim(emin.value, emax.value)
plt.xlabel('Energy ({0})'.format(self.inst[0].reg[0].tabarea['X_ctr'].unit.to_string('latex_inline')))
plt.ylabel('Effective area ({0})'.format(self.inst[0].reg[0].tabarea['Area'].unit.to_string('latex_inline')))
plt.title(title)
if show:
plt.show()
else:
return plt
[docs]class PlotAreaInst:
"""Instrument definition and region list for area plots.
:ivar nreg: Number of regions.
:vartype nreg: int
:ivar reg: List of region area plots.
:vartype reg: list
"""
def __init__(self):
self.nreg = 0 # Number of regions
self.reg = [] # List of region area plots
[docs] def update(self, inst, wave):
"""Update the area information in the regions of this instrument.
:param inst: Instrument number to update.
:type inst: int
"""
pyspex_f2py.api_instrument.api_instrument_update(float(inst))
self.nreg = int(pyspex_f2py.api_instrument.inst_nregion)
for i in numpy.arange(self.nreg):
r = PlotAreaReg()
r.update(inst, i + 1, wave)
self.reg.append(r)
[docs]class PlotAreaReg:
"""Effective area information for each region.
:ivar xc: Energy bin center value (keV)
:vartype xc: astropy.units.quantity.Quantity
:ivar xb: Energy bin upper boundary value (keV)
:vartype xb: astropy.units.quantity.Quantity
:ivar area: Effective area (:math:`m^2`)
:vartype area: astropy.units.quantity.Quantity
:ivar tabarea: Astropy QTable containing the effective area plot
:vartype tabarea: astropy.table.QTable
"""
def __init__(self):
self.xc = numpy.array([0.], dtype=object) * u.keV # Energy bin center value (keV)
self.xb = numpy.array([0.], dtype=object) * u.keV # Energy bin upper boundary value (keV)
self.area = numpy.array([0.], dtype=object) * u.m ** 2 # Effective area (m**2)
self.tabarea = QTable()
[docs] def update(self, inst, reg, wave=False):
"""Update effective area information for this region (reg) in instrument (inst).
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
"""
pyspex_f2py.api_plot_area.plot_area_update(float(inst), float(reg))
self.n = pyspex_f2py.api_plot_area.area_nbins
self.xc = numpy.zeros(self.n, dtype=object) * u.keV
self.xb = numpy.zeros(self.n, dtype=object) * u.keV
self.area = numpy.zeros(self.n, dtype=object) * u.m ** 2
for i in numpy.arange(self.n):
self.xc[i] = float(pyspex_f2py.api_plot_area.area_xc[i]) * u.keV
self.xb[i] = float(pyspex_f2py.api_plot_area.area_xb[i]) * u.keV
self.area[i] = float(pyspex_f2py.api_plot_area.area_a[i]) * u.m ** 2
if wave:
self._convert_to_ang()
self.tabarea = QTable([self.xc, self.xb, self.area], names=('X_ctr', 'X_upp', 'Area'),
meta={'name': 'Effective Area'})
pyspex_f2py.api_plot_area.plot_area_final()
def _convert_to_ang(self):
"""Convert the spectrum in keV to Angstrom."""
# Calculate the bin width in keV
delta_e = 2. * (self.xb - self.xc)
# Convert the individual arrays
self.xc = self.xc.to(u.Angstrom, equivalencies=u.spectral())
self.xb = self.xb - delta_e
self.xb = self.xb.to(u.Angstrom, equivalencies=u.spectral())
# Reverse the order of the arrays
self.xc = numpy.flip(self.xc)
self.xb = numpy.flip(self.xb)
self.area = numpy.flip(self.area)
[docs]class PlotData:
"""Class containing the observed spectrum and the convolved model spectrum for plotting.
:ivar ninst: Number of instruments
:vartype ninst: int
:ivar inst: List of instruments with data plot information
:vartype inst: list
:ivar colors: List of colors.
:vartype colors: tuple
"""
def __init__(self):
self.ninst = 0 # Number of instruments
self.inst = [] # List of instruments with data plot information
self.colors = ('g', 'b', 'c', 'm', 'y')
[docs] def update(self, chi='dchi', wave=False):
"""Update data plot information for this instrument.
:param chi: (Optional) Type of residual. Either 'dchi' (data-model)/error or 'rel' (data-model)/model.
:type chi: str
"""
self.inst = []
self.ninst = int(pyspex_f2py.api_instrument.api_instrument_n())
for i in numpy.arange(self.ninst):
ins = PlotDataInst()
ins.update(i + 1, tchi=chi, wave=wave)
self.inst.append(ins)
[docs] def plot_data(self, ax, xlog=False, ylog=False):
"""Generate plot of the data, model and background spectrum. This is a subfunction of plot() and chiplot().
:param ax: Axis object from Matplotlib.
:type ax: matplotlib.axes.Axes
:param xlog: (Optional) Set the X-axis to be logarithmic.
:type xlog: bool
:param ylog: (Optional) Set the Y-axis to be logarithmic.
:type ylog: bool
"""
for i in numpy.arange(self.ninst):
nreg = self.inst[i].nreg
for j in numpy.arange(nreg):
ax.step(self.inst[i].reg[j].tabdata['X_upp'], self.inst[i].reg[j].tabdata['Bkg'], where='pre',
linestyle='--', c='b')
xerror = self.inst[i].reg[j].tabdata['X_ctr'] - self.inst[i].reg[j].tabdata['X_low']
ax.errorbar(self.inst[i].reg[j].tabdata['X_ctr'], self.inst[i].reg[j].tabdata['Data'], xerr=xerror,
yerr=self.inst[i].reg[j].tabdata['Error'], marker='', ls='', c='k')
ax.step(self.inst[i].reg[j].tabdata['X_upp'], self.inst[i].reg[j].tabdata['Model'], where='pre',
c='r')
if xlog:
ax.set_xscale('log')
if ylog:
ax.set_yscale('log')
[docs] def plot_chi(self, ax):
"""Generate plot data for a plot of the residuals.
:param ax: Axis object from Matplotlib.
:type ax: matplotlib.axes.Axes
"""
for i in numpy.arange(self.ninst):
nreg = self.inst[i].nreg
for j in numpy.arange(nreg):
xerror = self.inst[i].reg[j].tabdata['X_ctr'] - self.inst[i].reg[j].tabdata['X_low']
ax.errorbar(self.inst[i].reg[j].tabdata['X_ctr'], self.inst[i].reg[j].tabdata['Chi'], xerr=xerror,
yerr=self.inst[i].reg[j].tabdata['Chi_err'], marker='', ls='', c='k')
[docs] def plot(self, xlog=False, ylog=False, wave=False, title='SPEX', show=True):
"""Plot the observed data, model and background spectra for all instruments.
:param xlog: (Optional) Set the X-axis to be logarithmic.
:type xlog: bool
:param ylog: (Optional) Set the Y-axis to be logarithmic.
:type ylog: bool
:param wave: (Optional) Plot the data on a wavelength scale (Angstrom).
:type wave: bool
:param title: (Optional) Set the title of the plot.
:type title: str
:param show: (Optional) Show the plot (True) or return the plot object (False).
:type show: bool
:return: Matplotlib plt object.
:rtype: matplotlib.pyplot
"""
self.update(wave=wave)
fig, ax = plt.subplots()
self.plot_data(ax, xlog=xlog, ylog=ylog)
if wave:
plt.xlabel('Wavelength ({0})'.format(self.inst[0].reg[0].tabdata['X_ctr'].unit.to_string('latex_inline')))
else:
plt.xlabel('Energy ({0})'.format(self.inst[0].reg[0].tabdata['X_ctr'].unit.to_string('latex_inline')))
plt.ylabel('Counts ({0})'.format(self.inst[0].reg[0].tabdata['Data'].unit.to_string('latex_inline')))
if xlog:
formatter = FuncFormatter(lambda y, _: '{:.16g}'.format(y))
ax.xaxis.set_major_formatter(formatter)
plt.title(title)
if show:
plt.show()
else:
return plt
[docs] def chiplot(self, xlog=False, ylog=False, chi='dchi', wave=False, title='SPEX', show=True):
"""Plot the observed data, model and background spectra for all instruments, and include a subwindow
with the residuals. Chi options are 'dchi': (observed - model)/error and 'rel': (observed - model)/model.
:param xlog: (Optional) Set the X-axis to be logarithmic.
:type xlog: bool
:param ylog: (Optional) Set the Y-axis to be logarithmic.
:type ylog: bool
:param chi: (Optional) Type of residual. Either 'dchi' (data-model)/error or 'rel' (data-model)/model.
:type chi: str
:param wave: (Optional) Plot residuals as a function of wavelength (Angstrom)
:type wave: bool
:param title: (Optional) Set the title of the plot.
:type title: str
:param show: (Optional) Show the plot (True) or return the plot object (False).
:type show: bool
:return: Matplotlib plt object.
:rtype: matplotlib.pyplot
"""
self.update(chi=chi, wave=wave)
fig = plt.figure()
gs = gridspec.GridSpec(ncols=1, nrows=4)
gs.update(hspace=0.05)
ax1 = fig.add_subplot(gs[:3])
ax2 = fig.add_subplot(gs[3], sharex=ax1)
self.plot_data(ax1, xlog=xlog, ylog=ylog)
ax1.set_ylabel('Counts ({0})'.format(self.inst[0].reg[0].tabdata['Data'].unit.to_string('latex_inline')))
ax1.set_title(title)
plt.setp(ax1.get_xticklabels(), visible=False)
self.plot_chi(ax2)
if chi == 'dchi':
ax2.set_ylabel('(Observed - Model) / Error')
elif chi == 'rel':
ax2.set_ylabel('(Observed - Model) / Model')
if xlog:
formatter = FuncFormatter(lambda y, _: '{:.16g}'.format(y))
ax2.xaxis.set_major_formatter(formatter)
if wave:
ax2.set_xlabel(
'Wavelength ({0})'.format(self.inst[0].reg[0].tabdata['X_ctr'].unit.to_string('latex_inline')))
else:
ax2.set_xlabel('Energy ({0})'.format(self.inst[0].reg[0].tabdata['X_ctr'].unit.to_string('latex_inline')))
if show:
plt.show()
else:
return plt
[docs] def plot_comp(self, xlog=False, ylog=False, wave=False, title='SPEX', show=True):
"""Plot the total spectrum and the individual additive model components.
:param xlog: (Optional) Set the X-axis to be logarithmic.
:type xlog: bool
:param ylog: (Optional) Set the Y-axis to be logarithmic.
:type ylog: bool
:param wave: (Optional) Plot as function of wavelenght (Angstrom).
:type wave: bool
:param title: (Optional) Set the title of the plot.
:type title: str
:param show: (Optional) Show the plot (True) or return the plot object (False).
:type show: bool
:return: Matplotlib plt object.
:rtype: matplotlib.pyplot
"""
self.update(wave=wave)
fig, ax = plt.subplots()
self.plot_data(ax, xlog=xlog, ylog=ylog)
# Add the individual model components
mod = model.Model()
mod.update()
# Find additive components and store their normalisations.
norms = numpy.array([], dtype=float)
s = numpy.array([], dtype=numpy.uint8)
c = numpy.array([], dtype=numpy.uint8)
n = 0
for i in numpy.arange(mod.nsector):
for j in numpy.arange(mod.sect[i].ncomp):
if mod.sect[i].comp[j].add:
par = mod.par_get(i + 1, j + 1, 'norm')
for k in numpy.arange(mod.sect[i].comp[j].npar):
# print("{:d} {:d} {:s}".format(i+1, j+1, par[k].name))
if mod.sect[i].comp[j].par[k].name == 'norm':
norms = numpy.append(norms, mod.sect[i].comp[j].par[k].value)
s = numpy.append(s, i + 1)
c = numpy.append(c, j + 1)
n = n + 1
print("Number of additive components: {:d}".format(n))
# Clear the model components
for i in numpy.arange(self.ninst):
nreg = self.inst[i].nreg
for j in numpy.arange(nreg):
self.inst[i].reg[j].comp = []
self.inst[i].reg[j].ncomp = 0
# Calculate model for each component
for k in numpy.arange(n):
print("Calculate component {:d} {:d}...".format(s[k], c[k]))
self._comp_norm_select(mod, k, norms, s, c)
ier = pyspex_f2py.api_command('calculate')
ier = pyspex_f2py.api_command('plot')
self.update(wave=wave)
for i in numpy.arange(self.ninst):
nreg = self.inst[i].nreg
for j in numpy.arange(nreg):
if wave:
com = numpy.zeros(self.inst[i].reg[j].n, dtype=float) * u.ct / (u.s * u.Angstrom)
else:
com = numpy.zeros(self.inst[i].reg[j].n, dtype=float) * u.ct / (u.s * u.keV)
for l in numpy.arange(self.inst[i].reg[j].n):
com[l] = self.inst[i].reg[j].pmodel[l]
if self.inst[i].reg[j].n != 0:
self.inst[i].reg[j].comp.append(com)
ax.step(self.inst[i].reg[j].tabdata['X_upp'],
self.inst[i].reg[j].comp[self.inst[i].reg[j].ncomp], where='pre', c=self.colors[k])
self.inst[i].reg[j].ncomp = self.inst[i].reg[j].ncomp + 1
# Reset model to original values
self._comp_norm_select(mod, -1, norms, s, c)
if xlog:
ax.set_xscale('log')
if ylog:
ax.set_yscale('log')
if wave:
plt.xlabel('Wavelength ({0})'.format(self.inst[0].reg[0].tabdata['X_ctr'].unit.to_string('latex_inline')))
else:
plt.xlabel('Energy ({0})'.format(self.inst[0].reg[0].tabdata['X_ctr'].unit.to_string('latex_inline')))
plt.ylabel('Counts ({0})'.format(self.inst[0].reg[0].tabdata['Data'].unit.to_string('latex_inline')))
if xlog:
formatter = FuncFormatter(lambda y, _: '{:.16g}'.format(y))
ax.xaxis.set_major_formatter(formatter)
plt.title(title)
if show:
plt.show()
else:
return plt
def _comp_norm_select(self, mod, inorm, norms, s, c):
"""Select normalisation.
:param mod: PYSPEX model object containing the SPEX model.
:type mod: model.Model
:param inorm: Component number to calculate (inorm=-1 restores the model to the original values).
:type inorm: int
:param norms: Array with the normalisations of the additive components.
:type norms: numpy.ndarray
:param s: Array with sector numbers of the additive components.
:type s: numpy.ndarray
:param c: Array with the component numbers of the additive components.
:type c: numpy.ndarray
"""
if inorm == -1:
for i in numpy.arange(norms.size):
mod.par_set(s[i], c[i], 'norm', norms[i])
else:
for i in numpy.arange(norms.size):
if i == inorm:
mod.par_set(s[i], c[i], 'norm', norms[i])
else:
mod.par_set(s[i], c[i], 'norm', 0.)
[docs]class PlotDataInst:
"""Class containing the data plot information for this instrument.
:ivar nreg: Number of regions.
:vartype nreg: int
:ivar reg: List of regions with data plot information
:vartype reg: list
"""
def __init__(self):
self.nreg = 0 # Number of regions
self.reg = [] # List of regions with data plot information
[docs] def update(self, inst, tchi='dchi', wave=False):
"""Update the data plot information for this instrument.
:param inst: Instrument number.
:type inst: int
:param tchi: (Optional) Type of residual. Either 'dchi' (data-model)/error or 'rel' (data-model)/model.
:type tchi: str
"""
self.reg = []
pyspex_f2py.api_instrument.api_instrument_update(float(inst))
self.nreg = int(pyspex_f2py.api_instrument.inst_nregion)
for i in numpy.arange(self.nreg):
r = PlotDataReg()
status = r.update(inst, i + 1, tchi=tchi, wave=wave)
if status == 0:
self.reg.append(r)
[docs]class PlotDataReg:
"""Class containing the data plot information for a region.
:ivar n: Number of bins
:vartype n: int
:ivar elow: Bin lower boundary (keV)
:vartype elow: astropy.units.quantity.Quantity
:ivar ectr: Bin center value (keV)
:vartype ectr: astropy.units.quantity.Quantity
:ivar eupp: Bin upper boundary (keV)
:vartype eupp: astropy.units.quantity.Quantity
:ivar ewid: Bin width (keV)
:vartype ewid: astropy.units.quantity.Quantity
:ivar pdata: Data value (Counts/s/keV)
:vartype pdata: astropy.units.quantity.Quantity
:ivar pmodel: Model value (Counts/s/keV)
:vartype pmodel: astropy.units.quantity.Quantity
:ivar pback: Background value (Counts/s/keV)
:vartype pback: astropy.units.quantity.Quantity
:ivar perror: Error value
:vartype perror: astropy.units.quantity.Quantity
:ivar tints: Exposure time (s)
:vartype tints: astropy.units.quantity.Quantity
:ivar tarea: Inverse of the effective area value (1/m**2)
:vartype tarea: astropy.units.quantity.Quantity
:ivar comp: List of model components for a component plot
:vartype comp: list
:ivar ncomp: Number of model components
:vartype ncomp: int
:ivar tabdata: Astropy QTable containing the plot data
:vartype tabdata: astropy.table.QTable
:ivar chi: Residual value
:vartype chi: numpy.ndarray
:ivar chierr: Residual error
:vartype chierr: numpy.ndarray
"""
def __init__(self):
self.n = 0 # Number of bins
self.elow = numpy.array([0.], dtype=object) * u.keV # Bin lower boundary
self.ectr = numpy.array([0.], dtype=object) * u.keV # Bin center value
self.eupp = numpy.array([0.], dtype=object) * u.keV # Bin upper boundary
self.ewid = numpy.array([0.], dtype=object) * u.keV # Bin width
self.pdata = numpy.array([0.], dtype=object) * u.ct / (u.s * u.keV) # Data value
self.pmodel = numpy.array([0.], dtype=object) * u.ct / (u.s * u.keV) # Model value
self.pback = numpy.array([0.], dtype=object) * u.ct / (u.s * u.keV) # Background value
self.perror = numpy.array([0.], dtype=object) * u.ct / (u.s * u.keV) # Error value
self.tints = numpy.array([0.], dtype=object) * u.s # Exposure time
self.tarea = numpy.array([0.], dtype=object) / u.m ** 2 # Inverse Effective Area value
self.tabdata = QTable()
# Model components
self.comp = [] # List of model components for a component plot
self.ncomp = 0 # Number of model components
# Residuals
self.chi = numpy.array([0.], dtype=object) # Residual value
self.chierr = numpy.array([0.], dtype=object) # Residual error
[docs] def update(self, inst, reg, tchi='dchi', wave=False):
"""Update the data plot information for this region.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param tchi: (Optional) Type of residual. Either 'dchi' (data-model)/error or 'rel' (data-model)/model.
:type tchi: str
"""
pyspex_f2py.api_plot_data.plot_data_update(float(inst), float(reg))
self.n = int(pyspex_f2py.api_plot_data.data_nbins)
if self.n == 0:
return 1
self.elow = numpy.zeros(self.n, dtype=object) * u.keV
self.ectr = numpy.zeros(self.n, dtype=object) * u.keV
self.eupp = numpy.zeros(self.n, dtype=object) * u.keV
self.ewid = numpy.zeros(self.n, dtype=object) * u.keV
self.pdata = numpy.zeros(self.n, dtype=object) * u.ct / (u.s * u.keV)
self.pmodel = numpy.zeros(self.n, dtype=object) * u.ct / (u.s * u.keV)
self.pback = numpy.zeros(self.n, dtype=object) * u.ct / (u.s * u.keV)
self.perror = numpy.zeros(self.n, dtype=object) * u.ct / (u.s * u.keV)
self.tints = numpy.zeros(self.n, dtype=object) * u.s
self.tarea = numpy.zeros(self.n, dtype=object) / u.m ** 2
for i in numpy.arange(self.n):
self.elow[i] = float(pyspex_f2py.api_plot_data.data_elow[i]) * u.keV
self.ectr[i] = float(pyspex_f2py.api_plot_data.data_ectr[i]) * u.keV
self.eupp[i] = float(pyspex_f2py.api_plot_data.data_eupp[i]) * u.keV
self.ewid[i] = float(pyspex_f2py.api_plot_data.data_ewid[i]) * u.keV
self.pdata[i] = float(pyspex_f2py.api_plot_data.data_data[i]) * u.ct / (u.s * u.keV)
self.pmodel[i] = float(pyspex_f2py.api_plot_data.data_model[i]) * u.ct / (u.s * u.keV)
self.pback[i] = float(pyspex_f2py.api_plot_data.data_back[i]) * u.ct / (u.s * u.keV)
self.perror[i] = float(pyspex_f2py.api_plot_data.data_error[i]) * u.ct / (u.s * u.keV)
self.tints[i] = float(pyspex_f2py.api_plot_data.data_tints[i]) * u.s
self.tarea[i] = float(pyspex_f2py.api_plot_data.data_tarea[i]) / u.m ** 2
pyspex_f2py.api_plot_data.plot_data_final()
if wave:
self._convert_to_ang()
self.chitype(tchi, wave=wave)
self.tabdata = QTable([self.elow, self.ectr, self.eupp, self.pdata, self.perror, self.pmodel, self.pback,
self.tints, self.tarea, self.chi, self.chierr],
names=('X_low', 'X_ctr', 'X_upp', 'Data', 'Error', 'Model', 'Bkg', 'Exposure',
'1/Area', 'Chi', 'Chi_err'),
meta={'name': 'Plot data'})
return 0
[docs] def chitype(self, tchi, wave=False):
"""Calculate the residuals for the chi plot window.
:param tchi: (Optional) Type of residual. Either 'dchi' (data-model)/error or 'rel' (data-model)/model.
:type tchi: str
"""
self.chi = numpy.zeros(self.n, dtype=object) * u.dimensionless_unscaled
self.chierr = numpy.zeros(self.n, dtype=object) * u.dimensionless_unscaled
if tchi == 'dchi':
for i in numpy.arange(self.n):
if self.perror[i] != 0:
self.chi[i] = (self.pdata[i] - self.pmodel[i]) / self.perror[i]
else:
if wave:
self.chi[i] = (self.pdata[i] - self.pmodel[i]) / (u.ct / (u.s * u.Angstrom))
else:
self.chi[i] = (self.pdata[i] - self.pmodel[i]) / (u.ct / (u.s * u.keV))
self.chierr[i] = 1. * u.dimensionless_unscaled
elif tchi == 'rel':
for i in numpy.arange(self.n):
self.chi[i] = (self.pdata[i] - self.pmodel[i]) / self.pmodel[i]
self.chierr[i] = self.perror[i] / self.pmodel[i]
def _convert_to_ang(self):
"""Convert the data plots to Angstrom"""
# Convert the individual arrays
wlow = numpy.copy(self.eupp)
wupp = numpy.copy(self.elow)
delta_e = numpy.copy(self.ewid)
self.ectr = self.ectr.to(u.Angstrom, equivalencies=u.spectral())
self.elow = wlow.to(u.Angstrom, equivalencies=u.spectral())
self.eupp = wupp.to(u.Angstrom, equivalencies=u.spectral())
self.ewid = self.eupp - self.elow
delta_w = self.ewid
self.pdata = self.pdata * delta_e / delta_w
self.pmodel = self.pmodel * delta_e / delta_w
self.pback = self.pback * delta_e / delta_w
self.perror = self.perror * delta_e / delta_w
# Reverse the order of the arrays
self.ectr = numpy.flip(self.ectr)
self.elow = numpy.flip(self.elow)
self.eupp = numpy.flip(self.eupp)
self.ewid = numpy.flip(self.ewid)
self.pdata = numpy.flip(self.pdata)
self.pmodel = numpy.flip(self.pmodel)
self.pback = numpy.flip(self.pback)
self.perror = numpy.flip(self.perror)
self.tints = numpy.flip(self.tints)
self.tarea = numpy.flip(self.tarea)
