Source code for pyspex.plot

#!/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)