#!/usr/bin/env python
from . import pyspex_f2py
from . import data
from . import model
from . import optimize
from . import log
from . import ascdump
from . import plot
import copy
import numpy
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
from astropy.visualization import quantity_support
quantity_support()
pyspex_f2py.spexapi.api_version()
__version__ = pyspex_f2py.spexapi.version.astype('|S1').tobytes().decode('utf-8')
class _Singleton(type):
"""This class is introduced to make sure that a PYSPEX session is only
initialized once. F2PY is not capable of creating truly separate
sessions, unfortunately. This class provides a common Python Singleton
implementation as a metaclass for Session."""
instance = None
def __call__(cls, *args, **kwargs):
if not cls.instance:
cls.instance = super(_Singleton, cls).__call__(*args, **kwargs)
else:
print("Warning: You can only start one session of PYSPEX within Python.")
print("This function will return the instance of the original session.")
return cls.instance
class Singleton(_Singleton(str("SingletonMeta"), (object,), {})): pass
[docs]class Session(Singleton):
"""This class contains all the classes and commands that are available in a pySPEX session."""
[docs] def __init__(self):
"""Function to initialize the PYSPEX session. It starts an instance of SPEX in the background.
Furthermore, it initializes a set of objects containing information about the loaded data,
model values, optimization and logs.
:ivar version: The SPEX version number
:vartype version: str
:ivar dataset: Class containing the data information
:vartype dataset: data.Data
:ivar mod_abundance: Model class for the abundance setting
:vartype mod_abundance: model.Abundance
:ivar mod_distance: Model class containing distance tools
:vartype mod_distance: model.Distance
:ivar mod_egrid: Model class for the definition of energy grids
:vartype mod_egrid: model.Egrid
:ivar mod_flux: Model class containing flux and luminosity information
:vartype mod_flux: model.Fluxes
:ivar mod_ibal: Model class for the ionisation balance setting
:vartype mod_ibal: model.Ibal
:ivar mod_ions: Model class for setting the used ions
:vartype mod_ions: model.Ions
:ivar mod_var: Model class for plasma model settings
:vartype mod_var: model.Var
:ivar mod_dem: Model class for DEM modeling
:vartype mod_dem: model.Dem
:ivar mod_spectrum: Model class containing commands to extract model results from SPEX
:vartype mod_spectrum: model.Spectrum
:ivar mod: Model class containing sectors, components and parameters.
:vartype mod: model.Model
:ivar opt_fit: Initialize the Fit class for spectral fitting
:vartype opt_fit: optimize.Fit
:ivar asc: Class for the Ascdump output
:vartype asc: ascdump.Ascdump
:ivar logs: Class for Log saving and execution
:vartype logs: log.Log
"""
pyspex_f2py.spexapi.api_init() # Initialize Fortran to Python interface
self.version = __version__ #: The SPEX version number
# Initialize class for data related commands
self.dataset = data.Data() #: Class containing the data information
# Initialize classes from the Model group
self.mod_abundance = model.Abundance()
self.mod_distance = model.Distance()
self.mod_egrid = model.Egrid() #: Model class for the definition of energy grids
self.mod_flux = model.Fluxes() #: Model class containing flux and luminosity information
self.mod_ibal = model.Ibal() #: Model class for the ionisation balance setting
self.mod_ions = model.Ions() #: Model class for setting the used ions
self.mod_var = model.Var() #: Model class for plasma model settings
self.mod_dem = model.Dem() #: Model class for DEM modeling
self.mod_spectrum = model.Spectrum() #: Model class containing commands to extract model results from SPEX.
self.mod = model.Model() #: Model class containing sectors, components and parameters.
# Initialize classes from the Optimize group
self.opt_fit = optimize.Fit() #: Initialize the Fit class for spectral fitting
# Initialize ascdump output
self.asc = ascdump.Ascdump()
# Initialize class for log saving and executing
self.logs = log.Log() #: Class for Log saving and execution
[docs] def command(self, comm):
"""Send a command to SPEX. The command string uses the same syntax as the SPEX commands on the command line.
:param comm: SPEX command.
:type comm: str
"""
ier = pyspex_f2py.api_command(comm)
return ier
# ======================================================
# Abundance commands
# ======================================================
[docs] def abun(self, abunset):
"""Set the abundance table in SPEX. The abunset parameter is one of 'reset', 'ag', 'allen', 'asplund',
'ra', 'grevesse', 'gs', 'lodders', and 'solar'.
:param abunset: String representing the abundance set used in SPEX.
:type abunset: str
:return: Object containing the abundance class.
:rtype: pyspex.model.Abundance
"""
ab = model.Abundance()
ab.set(abunset)
return ab
[docs] def abun_show(self):
"""Show the current ionisation balance.
:return: Reference to the used ionisation balance.
:rtype: str
"""
ab = self.mod_abundance.get()
print(ab)
return self.mod_abundance.ref
# ======================================================
# Aerror command
# ======================================================
[docs] def aerror(self, isect, icomp, name, shell=0):
"""Command to calculate the uncertainty in a model parameter due to the uncertainties
in the atomic data.
:param isect: Sector number of the component.
:type isect: int
:param icomp: Component number.
:type icomp: int
:param name: Parameter name.
:type name: str
:param shell: Shell number for calculation (L-shell = 1, K-shell = 2)
:type shell: int
:return: The atomic uncertainty for the selected parameter.
:rtype: float
"""
aerror = model.Aerror()
aerr = aerror.get(isect, icomp, name, shell=shell)
return aerr
# ======================================================
# Ascdump command
# ======================================================
[docs] def ascdump(self, isect, icomp, atype):
"""Generic ascdump method to obtain the various numeric outputs of the spectral models.
(see https://spex-xray.github.io/spex-help/reference/commands/ascdump.html for
more information).
:param isect: Sector number of the component.
:type isect: int
:param icomp: Component number.
:type icomp: int
:param atype: Ascdump type
:type atype: str
"""
asc = ascdump.Ascdump()
asc.get(isect, icomp, atype)
return asc
# ======================================================
# Bin command
# ======================================================
[docs] def bin(self, inst, reg, elow, ehigh, factor, unit=None):
"""Bin the spectrum using a fixed factor.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param factor: Binning factor
:type factor: int
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.bin(inst, reg, elow, ehigh, factor, unit=unit)
# ======================================================
# Calculate command
# ======================================================
[docs] def calc(self):
"""Calculate the current model."""
self.command('calculate')
self.mod.update()
# ======================================================
# Comp command
# ======================================================
[docs] def com(self, name, isect=1):
"""Add model component.
:param name: Name of the model component, for example 'reds', 'hot', 'cie', etc.
:type name: str
:param isect: Sector number to add component to (default is sector 1).
:type isect: int
"""
self.mod.comp_new(name, isect=isect)
[docs] def com_rel(self, isect, icomp, rel):
"""Relate model components.
:param isect: Sector number of the component to relate.
:type isect: int
:param icomp: Component number to relate.
:type icomp: int
:param rel: Array containing the multiplicative component numbers counted from the source to the observer.
:type rel: numpy.ndarray
"""
self.mod.comp_set_rel(isect, icomp, rel)
[docs] def com_del(self, isect, icomp):
"""Delete component from model.
:param isect: Sector number of the component to delete.
:type isect: int
:param icomp: Component number to delete.
:type icomp: int
"""
self.mod.comp_delete(isect, icomp)
# ======================================================
# Data commands
# ======================================================
[docs] def data(self, resfile, spofile):
"""Load instrument data into SPEX using a .res and a .spo file.
:param resfile: SPEX response file name (including .res extension).
:type resfile: str
:param spofile: SPEX spectrum file name (including .spo extension).
:type spofile: str
"""
ier = self.dataset.load(resfile, spofile)
return ier
[docs] def data_del(self, ins):
"""Delete an instrument from the dataset.
:param ins: Instrument number to delete.
:type ins: int
"""
self.dataset.delete(ins)
[docs] def data_save(self, ins, spofile, overwrite=False):
"""Save the (simulated) spectrum to a .spo file.
:param ins: Instrument number to save.
:type ins: int
:param spofile: Output filename for SPEX output file (including .spo extension).
:type spofile: str
:param overwrite: (Optional) Overwrite an existing .spo file with the same name.
:type overwrite: bool
"""
self.dataset.save(ins, spofile, overwrite=overwrite)
def data_show(self):
"""Show a summary of the loaded data."""
self.command('data show')
# ======================================================
# Dem commands
# ======================================================
[docs] def dem_lib(self):
"""Create the basic DEM library of isothermal spectra."""
self.command('dem lib')
[docs] def dem_reg_auto(self, s=1):
"""Do DEM analysis using the regularization method using an automatic search of the optimum regularization
parameter. It determines the regularisation parameter R in such a way that
:math:`\chi^2(R) = \chi^2(0) [1 + s \sqrt{2/(n-n_{\mathrm T}}]` where the scaling factor s=1, n is the number of
spectral bins in the data set and :math:`n_{\mathrm T}` is the number of temperature components in the DEM library.
:param s: Scaling factor in the automatic search for R (optional).
:type s: int
:return chisq: Best fit :math:`\chi^2(R)`.
:rtype chisq: float
:return dempen: DEM penalty (number of zero emission components)
:rtype dempen: float
"""
self.command('dem reg auto {0}'.format(s))
self.mod_dem.update()
return self.mod_dem.chisq, self.mod_dem.dempen
[docs] def dem_reg(self, r):
"""Do DEM analysis using the regularization method, using a
fixed regularization parameter R = r.
:param r: Regularization parameter R
:type r: int
:return chisq: Best fit :math:`\chi^2`.
:rtype chisq: float
:return dempen: DEM penalty (number of zero emission components)
:rtype dempen: float
"""
self.command('dem reg {0}'.format(r))
self.mod_dem.update()
return self.mod_dem.chisq, self.mod_dem.dempen
[docs] def dem_chireg(self, r1, r2, n):
"""Do a grid search over the regularization parameter R, with ``n`` steps and R distributed logarithmically
between ``r1`` and ``r2``. Useful to scan the :math:`\chi^2(R)` curve whenever it is complicated and to see how
much `penalty` (negative DEM values) there are for each value of R.
:param r1: Lower boundary for R.
:type r1: float
:param r2: Upper boundary for R.
:type r2: float
:param n: Number of steps.
:type n: int
:return chi_table: Table with the Chi^2 and DEM penalty values for each regularization parameter.
:rtype chi_table: astropy.table.QTable
"""
self.command('dem chireg {0}:{1} {2}'.format(r1, r2, n))
self.mod_dem.update()
return self.mod_dem.chi_table
[docs] def dem_clean(self):
"""Do DEM analysis using the clean method.
:return chisq: Best fit :math:`\chi^2`.
:rtype chisq: float
:return dempen: DEM penalty (number of zero emission components)
:rtype dempen: float
"""
self.command('dem clean')
self.mod_dem.update()
return self.mod_dem.chisq, self.mod_dem.dempen
[docs] def dem_poly(self, degree):
"""Do DEM analysis using the polynomial method, where ``degree`` is the degree of the polynomial.
:param degree: The degree of the polynomial.
:type degree: int
:return chisq: Best fit :math:`\chi^2`.
:rtype chisq: float
:return dempen: DEM penalty (number of zero emission components)
:rtype dempen: float
"""
self.command('dem poly {0}'.format(degree))
self.mod_dem.update()
return self.mod_dem.chisq, self.mod_dem.dempen
[docs] def dem_mult(self, ncomp):
"""Do DEM analysis using the multi-temperature method, where ``ncomp`` is the number of broad components.
:param ncomp: Number of broad temperature components.
:type ncomp: int
:return mult_table: Table containing the best temperature components and their emission measure.
:rtype mult_table: astropy.table.QTable
:return chisq: Best fit :math:`\chi^2`.
:rtype chisq: float
"""
self.command('dem mult {0}'.format(ncomp))
self.mod_dem.update()
return self.mod_dem.chisq, self.mod_dem.mult_table
[docs] def dem_gene(self, pop, gen):
"""Do DEM analysis using the genetic algorithm, using a population size given by ``pop`` (maximum value 1024)
and ``gen`` is the number of generations (no limit, in practice after about 100 generations not much change in
the solution. Experiment with these numbers for your practical case).
:param pop: Population size
:type pop: int
:param gen: Number of generations
:type gen: int
:return mult_table: Table containing the best temperature components and their emission measure.
:rtype mult_table: astropy.table.QTable
:return chisq: Best fit :math:`\chi^2`.
:rtype chisq: float
"""
self.command('dem gene {0} {1}'.format(pop, gen))
self.mod_dem.update()
return self.mod_dem.chisq, self.mod_dem.dempen
[docs] def dem_read(self, file):
""" Read a DEM distribution from a file named #a which automatically gets the extension ``.dem``. It is an ascii
file with at least two columns, the first column is the temperature in keV and the second column the
differential emission measure, in units of :math:`10^{64} m^{-3} keV^{-1}`. The maximum number of data points
in this file is 8192. Temperature should be in increasing order. The data will be interpolated to match
the temperature grid defined in the dem model (which is set by the user).
:param file: Filename without ``.dem`` extension.
:type file: str
"""
self.command('dem read {0}'.format(file))
self.mod_dem.update()
[docs] def dem_write(self, file):
"""Save the DEM to a file ``file`` with extension “.dem”. The same format as above is used for the file. A third
column has the corresponding error bars on the DEM as determined by the DEM method used (not always relevant
or well defined, exept for the regularization method).
:param file: Filename without ``.dem`` extension.
:type file: str
"""
self.command('dem write {0}'.format(file))
[docs] def dem_smooth(self, width):
"""Smoothes a DEM previously determined by any DEM method using a block filter/ Here ``width`` is the full width
of the filter expressed in :math:`^{10}\log T`. Note that this smoothing will in principle worsen the
:math:`\chi^2` of the solution, but it is sometimes useful to “wash out” some residual noise in the DEM
distribution, preserving total emission measure.
:param width: The full width of the filter expressed in :math:`^{10}\log T`.
:type width: float
"""
self.command('dem smooth {0}'.format(width))
self.mod_dem.update()
[docs] def dem_get(self):
"""Get the best-fit DEM distribution.
:return table: Table with the DEM distribution
:rtype table: astropy.table.QTable
"""
self.mod_dem.update()
return self.mod_dem.table
[docs] def dem_plot(self, xlog=False, ylog=False):
"""Plot the DEM distribution.
:param xlog: Make the X-axis logarithmic.
:type xlog: bool
:param ylog: Make the Y-axis logarithmic.
:type ylog: bool
"""
fig, ax = plt.subplots()
ax.errorbar(self.mod_dem.table['kT'], self.mod_dem.table['DY'], yerr=self.mod_dem.table['DY_Err'],
marker='o', ls='', c='k')
if xlog:
plt.xscale('log')
formatter = FuncFormatter(lambda y, _: '{:.16g}'.format(y))
ax.xaxis.set_major_formatter(formatter)
if ylog:
plt.yscale('log')
plt.xlabel('Temperature ({0})'.format(self.mod_dem.table['kT'].unit.to_string('latex_inline')))
plt.ylabel('Differential Emission Measure ({0})'.format(self.mod_dem.table['DY'].unit.to_string('latex_inline')))
plt.show()
# ======================================================
# Distance commands
# ======================================================
[docs] def dist(self, isect, dist, unit):
"""Set the distance in SPEX.
:param isect: Sector number.
:type isect: int
:param dist: Distance value.
:type dist: float
:param unit: The unit of the distance value, for example: 'm', 'au', 'ly', 'pc', 'kpc', 'mpc', 'z', 'cz'
:type unit: str
"""
self.mod_distance.set(isect, dist, unit)
return self.mod_distance
[docs] def dist_cosmo(self, h0, omega_m, omega_l, omega_r):
"""Set the cosmology for the distance calculation.
:param h0: Hubble constant (km/s/Mpc).
:type h0: float
:param omega_m: Omega matter.
:type omega_m: float
:param omega_l: Omega lambda.
:type omega_l: float
:param omega_r: Omega R.
:type omega_r: float
"""
self.mod_distance.set_cosmo(h0, omega_m, omega_l, omega_r)
[docs] def dist_get(self, isect):
"""Get the current distances in SPEX.
:param isect: Sector number.
:type isect: int
"""
self.mod_distance.get(isect)
return self.mod_distance
# ======================================================
# Egrid commands
# ======================================================
[docs] def egrid(self, elow, ehigh, nbins, unit, log):
"""Egrid command when the number of desired bins is known.
:param elow: Lowest energy/wavelength for energy grid.
:type elow: float
:param ehigh: Highest energy/wavelength for energy grid.
:type ehigh: float
:param nbins: Number of bins for energy grid.
:type nbins: int
:param unit: Unit of the energy/wavelength range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
:param log: Make the energy grid logarithmic (True or False)
:type log: bool
"""
self.mod_egrid.set(elow, ehigh, nbins, unit, log)
[docs] def egrid_step(self, elow, ehigh, step, unit, log):
"""Egrid command when the stepsize for the grid is known.
:param elow: Lowest energy/wavelength for energy grid.
:type elow: float
:param ehigh: Highest energy/wavelength for energy grid.
:type ehigh: float
:param step: Step size for the energy grid.
:type step: float
:param unit: Unit of the energy/wavelength range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
:param log: Make the energy grid logarithmic (True or False)
:type log: bool
"""
self.mod_egrid.set_step(elow, ehigh, step, unit, log)
[docs] def egrid_save(self, savefile):
"""Save energy grid to file.
:param savefile: Filename to save the energy grid to (including .egr extension)
:type savefile: str
"""
self.mod_egrid.save(savefile)
[docs] def egrid_read(self, readfile):
"""Read energy grid from file.
:param readfile: Filename to read the energy grid from (including .egr extension)
:type readfile: str
"""
self.mod_egrid.read(readfile)
[docs] def egrid_get(self):
"""Get the energy grid and return the grid as an object.
:return: An Egrid object from pyspex.
:rtype: pyspex.model.Egrid
"""
self.mod_egrid.get()
return self.mod_egrid
[docs] def egrid_set(self, ebounds):
"""Provide a grid to SPEX by providing a numpy array with the bin boundaries.
Please note that the length of this array is the number of bins + 1!
:param ebounds: Array containing the energy boundaries of the new energy grid.
:type ebounds: numpy.ndarray
"""
self.mod_egrid.grid(ebounds)
# ======================================================
# Elim command
# ======================================================
[docs] def elim(self, elow, ehigh, unit):
"""Set the energy limits for flux output.
:param elow: Lowest energy/wavelength for energy interval.
:type elow: float
:param ehigh: Highest energy/wavelength for energy interval.
:type ehigh: float
:param unit: Unit of the energy/wavelength interval, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'.
:type unit: str
"""
self.mod_flux.elim(elow, ehigh, unit)
[docs] def flux_get(self, isect, icomp):
"""Get the fluxes and luminosities for component icomp in sector isect.
:param isect: Sector number of the component to obtain the flux from.
:type isect: int
:param icomp: Component number to obtain the flux from.
:type icomp: int
:return flux: A Flux object from pyspex.
:rtype flux: pyspex.model.Flux
"""
self.mod_flux.get(isect, icomp)
flux = copy.deepcopy(self.mod_flux)
return flux
# ======================================================
# Error
# ======================================================
[docs] def error(self, isect, icomp, name, dchi=None):
"""Calculate the error for a parameter (name) in sector (isect) and component (icomp).
:param isect: Sector number of the component to calculate the error for.
:type isect: int
:param icomp: Component number to calculate the error for.
:type icomp: int
:param name: Parameter name to calculate the error for.
:type name: str
:param dchi: (Optional) :math:`\Delta\chi^2` value to optimize for (Default: 1.0, 68% errors)
:type dchi: float
:return: An Error object from pyspex.
:rtype: pyspex.optimize.Error
"""
err = optimize.Error()
err.error(isect, icomp, name, dchi=dchi)
return err
# ======================================================
# Fit commands
# ======================================================
[docs] def fit(self, niter=100):
"""Fit command.
:param niter: (Optional) Maximum number of iterations.
:type niter: int
"""
self.opt_fit.fit(niter=niter)
self.mod.update()
[docs] def fit_print(self, status):
"""Set the fit output verbosity. Default is true.
:param status: Fit output verbosity setting.
:type status: bool
"""
self.opt_fit.print(status)
[docs] def fit_stat(self, stat):
"""Set the fit statistics.
:param stat: Fit statistics, for example: 'csta', 'chi2', 'wsta'
:type stat: str
"""
self.opt_fit.set_statistic(stat)
[docs] def fit_stat_inst(self, stat, inst, reg):
"""Set the fit statistics for a particular instrument and region.
:param stat: Fit statistics, for example: 'csta', 'chi2', 'wsta'
:type stat: str
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
"""
self.opt_fit.set_statistic_inst(stat, inst, reg)
[docs] def fit_method(self, method):
"""Set the fit statistics for a particular instrument and region.
:param method: Fit method: classical, simplex or anneal
:type stat: str
"""
self.opt_fit.set_method(method)
[docs] def fit_cstat(self):
"""Get the current C-statistics from the model.
:return: Tuple containing the C-stat value and the degrees of freedom.
:rtype: tuple
"""
self.opt_fit.update()
return float(self.opt_fit.cstat), int(self.opt_fit.nfree)
[docs] def fit_chisq(self):
"""Get the current Chi2-statistics from the model.
:return: Tuple containing the :math:`\chi^2` value and the degrees of freedom.
:rtype: tuple
"""
self.opt_fit.update()
return float(self.opt_fit.chisq), int(self.opt_fit.nfree)
[docs] def fit_set_ann(self, param, value):
"""Set the simulating annealing method parameters.
:param param: Type of annealing parameter (rt, t, eps, vm, ns, max, or print).
:type param: str
:param value: Value of the parameter (will be converted to the nearest int if necessary).
:type value: float
"""
self.opt_fit.set_ann(param, value)
# Help commands not implemented
# ======================================================
# Ibal commands
# ======================================================
[docs] def ibal(self, ref):
"""Set the ionisation balance.
:param ref: Abbreviation of the reference to the ionisation balance.
:type ref: str
"""
self.mod_ibal.set(ref)
[docs] def ibal_show(self):
"""Show the current ionisation balance.
:return: Reference to the used ionisation balance.
:rtype: str
"""
ibal = self.mod_ibal.get()
print(ibal)
return self.mod_ibal.ref
# ======================================================
# Ignore
# ======================================================
[docs] def ignore(self, inst, reg, elow, ehigh, unit=None):
"""Ignore the bins given by the energy/channel range.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.ignore(inst, reg, elow, ehigh, unit=unit)
# ======================================================
# Ions
# ======================================================
[docs] def ions_show(self):
"""Show the ions in use and with which parameters."""
self.mod_ions.show()
[docs] def ions_all(self, par, value):
"""Set parameter for all values.
:param par: Parameter name to set, e.g.: 'use', 'new', 'nmax', 'lmax'.
:type par: str
:param value: Value to set.
:type value: bool or int
"""
if par == 'new':
if value == 2:
self.mod_ions.new_all()
elif value == 1:
self.mod_ions.qc_all()
elif value == 0:
self.mod_ions.old_all()
else:
print("Value unknown. Please choose 0,1,2")
elif par == 'use':
if value == 1:
self.mod_ions.use_all()
else:
self.mod_ions.ignore_all()
elif par == 'nmax':
self.mod_ions.nmax_all(value)
elif par == 'lmax':
self.mod_ions.lmax_all(value)
[docs] def ions_iso(self, iso, par, value):
"""Set parameter for an iso-electronic sequence.
:param iso: Iso-electronic sequence number
:type iso: int
:param par: Parameter name to set, e.g.: 'use', 'new', 'nmax', 'lmax'.
:type par: str
:param value: Value to set.
:type value: bool or int
"""
if par == 'new':
if value == 2:
self.mod_ions.new_iso(iso)
elif value == 1:
self.mod_ions.qc_iso(iso)
elif value == 0:
self.mod_ions.old_iso(iso)
else:
print("Value unknown. Please choose 0,1,2")
elif par == 'use':
if value == 1:
self.mod_ions.use_iso(iso)
else:
self.mod_ions.ignore_iso(iso)
elif par == 'nmax':
self.mod_ions.nmax_iso(iso, value)
elif par == 'lmax':
self.mod_ions.lmax_iso(iso, value)
[docs] def ions_z(self, z, par, value):
"""Set parameter for an atom.
:param z: Atomic number of the element.
:type z: int
:param par: Parameter name to set, e.g.: 'use', 'new', 'nmax', 'lmax'.
:type par: str
:param value: Value to set.
:type value: bool or int
"""
if par == 'new':
if value == 2:
self.mod_ions.new_z(z)
elif value == 1:
self.mod_ions.qc_z(z)
elif value == 0:
self.mod_ions.old_z(z)
else:
print("Value unknown. Please choose 0,1,2")
elif par == 'use':
if value == 1:
self.mod_ions.use_z(z)
else:
self.mod_ions.ignore_z(z)
elif par == 'nmax':
self.mod_ions.nmax_z(z, value)
elif par == 'lmax':
self.mod_ions.lmax_z(z, value)
[docs] def ions_ion(self, z, i, par, value):
"""Set parameter for a specific ion.
:param z: Atomic number of the element.
:type z: int
:param i: Ionisation stage (in astrophysical notation, e.g. 1, 2, 3, etc.)
:type i: int
:param par: Parameter name to set, e.g.: 'use', 'new', 'nmax', 'lmax'.
:type par: str
:param value: Value to set.
:type value: bool or int
"""
if par == 'new':
if value == 2:
self.mod_ions.new_ion(z, i)
elif value == 1:
self.mod_ions.qc_ion(z, i)
elif value == 0:
self.mod_ions.old_ion(z, i)
else:
print("Value unknown. Please choose 0,1,2")
elif par == 'use':
if value == 1:
self.mod_ions.use_ion(z, i)
else:
self.mod_ions.ignore_ion(z, i)
elif par == 'nmax':
self.mod_ions.nmax_ion(z, i, value)
elif par == 'lmax':
self.mod_ions.lmax_ion(z, i, value)
[docs] def ions_line(self, z, i, lid, add):
"""Mutes or unmutes a spectral line from an emission spectrum.
:param z: Atomic number of the element.
:type z: int
:param i: Ionisation stage (in astrophysical notation, e.g. 1, 2, 3, etc.)
:type i: int
:param lid: Line ID number (see ascdump_line for the line ID number).
:type i: int
:param add: If False, the line is not added to the spectrum. True will add the line to the spectrum again.
:type mute: bool
"""
if add:
com = 'ions unmute line {:n} ion {:n} {:n}'.format(lid, z, i)
self.command(com)
print(com)
else:
com = 'ions mute line {:n} ion {:n} {:n}'.format(lid, z, i)
self.command(com)
print(com)
# ======================================================
# Log
# ======================================================
[docs] def log_exe(self, filename):
"""Execute log file.
:param filename: Filename of the ASCII file with the SPEX commands (including .com extension).
:type filename: str
"""
self.logs.execute(filename)
self.mod.update()
[docs] def log_out(self, outfile, status=None):
"""Save SPEX terminal output to file.
:param outfile: Filename of the output file (including .out extension).
:type outfile: str
:param status: (Optional) Overwrite or append an existing file name. String can be 'overwrite' or 'append'.
:type status: str
"""
self.logs.output(outfile, status=status)
[docs] def log_save(self, savefile, status=None):
"""Save the commands for a session in an ascii file. The status is optional and either overwrite or append.
:param savefile: Filename of the output command file (including .com extension).
:type savefile: str
:param status: (Optional) Overwrite or append an existing file name. String can be 'overwrite' or 'append'.
:type status: str
"""
self.logs.save(savefile, status=status)
[docs] def log_close(self, logtype):
"""Close the save or output file.
:param logtype: String indicating the type of log to close ('save' or 'out').
:type logtype: str
"""
self.logs.close(logtype)
# Menu (not implemented)
# ======================================================
# Model show
# ======================================================
# Multiply
# ======================================================
# Obin
# ======================================================
[docs] def obin(self, inst, reg, elow, ehigh, unit=None):
"""Bin the spectrum optimally given the instrument resolution and statistics.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.obin(inst, reg, elow, ehigh, unit=unit)
# ======================================================
# Optimally bin response
# ======================================================
[docs] def rbin(self, inst, reg, elow, ehigh, unit=None):
"""Bin the spectrum and response optimally given the instrument resolution and statistics.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.rbin(inst, reg, elow, ehigh, unit=unit)
# Optimize (not implemented)
# ======================================================
# Par
# ======================================================
[docs] def par(self, isect, icomp, name, value, thawn=False):
"""Set parameter value and status.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
:param value: New value of the parameter.
:type value: float
:param thawn: (Optional) Should the parameter be free (True) or frozen (False).
:type thawn: bool
"""
self.mod.par_set(isect, icomp, name, value, thawn=thawn)
[docs] def par_text(self, isect, icomp, name, value):
"""Set text type parameter value.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
:param value: New value of the parameter.
:type value: str
"""
self.mod.par_aset(isect, icomp, name, value)
def par_get(self, isect, icomp, name):
"""Return all the properties of a parameter with sector (isect), component (icomp) and name (name).
The output is a parameter object.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
:return par: Parameter properties
:rtype par: pyspex.model.Parameter
"""
par = self.mod.par_get(isect, icomp, name)
return par
[docs] def par_fix(self, isect, icomp, name):
"""Fix a model parameter during the fit.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
"""
self.mod.par_fix(isect, icomp, name)
[docs] def par_free(self, isect, icomp, name):
"""Free a model parameter during the fit.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
"""
self.mod.par_free(isect, icomp, name)
[docs] def par_range(self, isect, icomp, name, rlow, rupp):
"""Set parameter range.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
:param rlow: Lower limit of the parameter range.
:type rlow: float
:param rupp: Upper limit of the parameter range.
:type rupp: float
"""
self.mod.par_set_range(isect, icomp, name, rlow, rupp)
[docs] def par_couple(self, isect, icomp, iname, csect, ccomp, cname, factor):
"""Couple one parameter to another with a given multiplication factor.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param iname: Parameter name.
:type iname: str
:param csect: Sector number of the parameter to couple to.
:type csect: int
:param ccomp: Component number of the parameter to couple to.
:type ccomp: int
:param cname: Parameter name of the parameter to couple to.
:type cname: str
:param factor: Multiplication factor (value(name) = factor * value(cname)).
:type factor: float
"""
self.mod.par_couple(isect, icomp, iname, csect, ccomp, cname, factor)
[docs] def par_decouple(self, isect, icomp, name):
"""Decouple a parameter.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Parameter name.
:type name: str
"""
self.mod.par_decouple(isect, icomp, name)
[docs] def par_norm(self, ins, reg, value, thawn):
"""Set an instrument normalisation.
:param ins: Instrument number.
:type ins: int
:param reg: Region number.
:type reg: int
:param value: New value of the instrument normalisation.
:type value: float
:param thawn: (Optional) Should the parameter be free (True) or frozen (False).
:type thawn: bool
"""
self.mod.par_norm_set(ins, reg, value, thawn)
[docs] def par_show_classic(self, option=''):
"""Display parameter overview in terminal through the Fortran backend.
:param option: Select which information should be shown (free/couple/flux/stat/corr/all)
:type option: str
"""
self.mod.par_show_classic(option=option)
[docs] def par_show(self, option=''):
"""Display parameter overview in terminal. Options include:
- all: Show all information.
- free: Show only free parameters.
- couple: Show coupled parameters.
- flux: Show summary of the fluxes of each emission component.
:param option: Select which information should be shown (free/couple/flux/all)
:type option: str
"""
# Do the regular par show.
self.mod.update()
self.mod.par_show(option=option)
print("")
# Obtain the instrument normalisations if available and show.
self.dataset.update()
if self.dataset.ninst != 0:
for i in numpy.arange(self.dataset.ninst):
for j in numpy.arange(self.dataset.inst[i].nreg):
(norm, free) = self.mod.par_norm_get(self.dataset.inst[i].index, self.dataset.inst[i].reg[j].index)
if free:
status = 'thawn'
else:
status = 'frozen'
print("Instrument {0:4d} region {1:4d} has norm {2:#11.6G} and is {3:6s}."
.format(self.dataset.inst[i].index, self.dataset.inst[i].reg[j].index, norm, status))
print("")
# Show the fit statistics (if available)
if self.opt_fit.nfree != 0:
self.opt_fit.show()
def par_show_param(self):
"""Display the parameters through the Python interface and create Astropy table.
:return: tab
:rtype tab: astropy.table.QTable
"""
self.mod.update()
tab = self.mod.par_show_param()
return tab
def par_show_free(self):
"""Display the free parameters through the Python interface and create Astropy table.
:return: tab
:rtype tab: astropy.table.QTable
"""
self.mod.update()
tab = self.mod.par_show_free()
return tab
def par_show_couple(self):
"""Display the coupled parameters through the Python interface and create Astropy table.
:return: tab
:rtype tab: astropy.table.QTable
"""
self.mod.update()
tab = self.mod.par_show_couple()
return tab
def par_show_flux(self):
"""Display the fluxes per component through the Python interface and create Astropy table.
:return: tab
:rtype tab: astropy.table.QTable
"""
self.mod.update()
tab = self.mod.par_show_flux()
return tab
[docs] def par_write(self, comfile, overwrite=False):
"""Write parameters to a .com file.
:param comfile: Output file name (with .com extension!)
:type comfile: str
:param overwrite: (Optional) Overwrite existing files if necessary (True/False).
:type overwrite: bool
"""
self.mod.par_write(comfile, overwrite=overwrite)
# ======================================================
# Plot
# ======================================================
[docs] def plot_area(self, xlog=False, ylog=False, wave=False, title='Effective Area', show=True):
"""Plot the effective area of 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 effective area as 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: Plot object and optionally the matplotlib plt object.
:rtype: pyspex.plot.PlotArea, matplotlib.pyplot
"""
pl = plot.PlotArea()
if show:
pl.plot(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl
else:
plt = pl.plot(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl, plt
[docs] def plot_data(self, xlog=False, ylog=False, wave=False, title='SPEX', show=True):
"""Plot the observed spectrum with the convolved model 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 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: Plot object and optionally the matplotlib plt object.
:rtype: pyspex.plot.PlotData, matplotlib.pyplot
"""
pl = plot.PlotData()
if show:
pl.plot(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl
else:
plt = pl.plot(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl, plt
[docs] def plot_chi(self, xlog=False, ylog=False, chi='dchi', wave=False, title='SPEX', show=True):
"""Plot the observed spectrum with the convolved model and their residuals 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 chi: (Optional) Type of residual, either 'dchi' or 'rel'.
:type chi: str
:param wave: (Optional) Plot residuals as 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: Plot object and optionally the matplotlib plt object.
:rtype: pyspex.plot.PlotData, matplotlib.pyplot
"""
pl = plot.PlotData()
if show:
pl.chiplot(xlog=xlog, ylog=ylog, chi=chi, wave=wave, title=title, show=show)
return pl
else:
plt = pl.chiplot(xlog=xlog, ylog=ylog, chi=chi, wave=wave, title=title, show=show)
return pl, plt
[docs] def plot_comp(self, xlog=False, ylog=False, wave=False, title='SPEX', show=True):
"""Plot the individual 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 the components as 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: Plot object and optionally the matplotlib plt object.
:rtype: pyspex.plot.PlotData, matplotlib.pyplot
"""
pl = plot.PlotData()
if show:
pl.plot_comp(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl
else:
plt = pl.plot_comp(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl, plt
[docs] def plot_model(self, xlog=False, ylog=False, wave=False, title='SPEX Model Spectrum', show=True):
"""Plot the spectral 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 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: Plot object and optionally the matplotlib plt object.
:rtype: pyspex.plot.PlotModel, matplotlib.pyplot
"""
pl = plot.PlotModel()
if show:
pl.plot(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl
else:
plt = pl.plot(xlog=xlog, ylog=ylog, wave=wave, title=title, show=show)
return pl, plt
def plot_model_adum_xspec(self, isect, qdpfile, overwrite=False):
"""Save the model for sector isect into Xspec QDP format. This format is useful for use with
simulation software, 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
"""
pl = plot.PlotModel()
status = pl.adum_xspec_qdp(isect, qdpfile, overwrite=overwrite)
if status != 0:
print("Error: Provided sector number not found.")
return
return pl
# ======================================================
# Sector
# ======================================================
[docs] def sector(self):
"""Create a new sector."""
self.mod.sector_new()
[docs] def sector_copy(self, isect):
"""Copy sector with number isect.
:param isect: Sector number.
:type isect: int
"""
self.mod.sector_copy(isect)
[docs] def sector_del(self, isect):
"""Delete sector with number isect.
:param isect: Sector number.
:type isect: int
"""
self.mod.sector_delete(isect)
def sector_adump(self, isect):
"""Dump the spectrum of this sector to an object.
:param isect: Sector number.
:type isect: int
:return: Spectrum object.
:rtype: pyspex.model.Spectrum
"""
if (isect <= self.mod.nsector):
self.mod.sect[isect-1].update(isect)
spectrum = self.mod.sect[isect-1].spectrum
return spectrum
else:
return 1
# Shiftplot
# ======================================================
# Simulate
# ======================================================
[docs] def simulate(self, extime, inst=None, reg=None, ssys=None, bsys=None, noise=None, bnoise=None, seed=None):
"""Simulate a spectrum using a loaded response file.
:param extime: Exposure time to simulate.
:type extime: float
:param inst: (Optional) Instrument range to simulate (default all)
:type inst: str
:param reg: (Optional) Region range to simulate (default all)
:type reg: str
:param ssys: (Optional) Add a systematic error to the source spectrum (Default 0).
:type ssys: float
:param bsys: (Optional) Add a systematic error to the background spectrum (Default 0).
:type bsys: float
:param noise: (Optional) Add Poisson noise to the simulated source spectrum (Default True).
:type noise: bool
:param bnoise: (Optional) Add Poisson noise to the simulated background spectrum (Default False).
:type bnoise: bool
:param seed: (Optional) Set the random seed for the simulation (Default: system clock).
:type seed: int
"""
sim = data.Simulate()
sim.simulate(extime, inst=inst, reg=reg, ssys=ssys, bsys=bsys, noise=noise, bnoise=bnoise, seed=seed)
# Step
# ======================================================
# Syserr
# ======================================================
[docs] def syserr(self, inst, reg, elow, ehigh, src, bkg, unit=None):
"""Add an additional error to the source and background spectrum.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param src: Error value to be quadratically added to the current error bar of the source spectrum.
:type src: float
:param bkg: Error value to be quadratically added to the current error bar of the background spectrum.
:type bkg: float
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.syserr(inst, reg, elow, ehigh, src, bkg, unit=unit)
# System (not implemented)
# ======================================================
# Use
# ======================================================
[docs] def use(self, inst, reg, elow, ehigh, unit=None):
"""Use the bins given by the energy/channel range.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.use(inst, reg, elow, ehigh, unit=unit)
# ======================================================
# Var
# ======================================================
[docs] def var_gacc(self, value):
"""Set the free-bound emission accuracy.
:param value: Free-bound emission accuracy.
:type value: float
"""
self.mod_var.set_gacc(value)
def var_gacc_reset(self):
"""Reset the free-bound emission accuracy."""
self.mod_var.reset_gacc()
[docs] def var_line(self, ltype, status):
"""Set a line emission contribution on or off.
:param ltype: Line emission contribution ('ex', 'px', 'rr', 'dr', 'ds', 'ii', 'reset')
:type ltype: str
:param status: Boolean indicator whether contribution is on (True) or off (False).
:type status: bool
"""
self.mod_var.set_line(ltype, status)
[docs] def var_doppler(self, value):
"""Doppler broadening.
:param value: Doppler broadening type.
:type value: int
"""
self.mod_var.set_doppler(value)
[docs] def var_calc(self, status):
"""Perform SPEXACT 3 line calculations.
:param status: Use SPEXACT 2 (0), Quick SPEXACT 3 (1), or SPEXACT 3 (2)
:type status: int
"""
self.mod_var.set_calc(status)
[docs] def var_occstart(self, otype):
"""At which occupation level to start.
:param otype: Occupation level ('ground', 'boltz', 'last')
:type otype: str
"""
self.mod_var.set_occstart(otype)
[docs] def var_mekal(self, utype, status):
"""Switch old Mekal updates on/off.
:param utype: Update type ('wav', 'fe17', 'update', 'all')
:type utype: str
:param status: On (True) or off (False)
:type status: bool
"""
self.mod_var.set_mekal(utype, status)
[docs] def var_ibalmaxw(self, status):
"""Switch the Multi-Maxwellians for the ionisation balance on/off (True/False).
:param status: On (True) or off (False)
:type status: bool
"""
self.mod_var.set_ibalmaxw(status)
[docs] def var_newcoolexc(self, status):
"""Cooling by collisional excitation by Stofanova (SPEXACT 3) on/off (True/False).
:param status: On (True) or off (False)
:type status: bool
"""
self.mod_var.set_newcoolexc(status)
[docs] def var_newcooldr(self, status):
"""Cooling by dielectronic recombination (SPEXACT 3) on/off (True/False).
:param status: On (True) or off (False)
:type status: bool
"""
self.mod_var.set_newcooldr(status)
[docs] def var_cxcon(self, value):
"""Set charge exchange recombination and ionization rates according to either
1 = Arnaud & Rothenflug (1985) or 2 = Kingdon & Ferland (1996). Default is 2.
:param value: Recombination and ionisation rateset. 1 = Arnaud & Rothenflug, 2 = Kingdon & Ferland.
:type value: int
"""
self.mod_var.set_cxcon(value)
# ======================================================
# Vbin
# ======================================================
[docs] def vbin(self, inst, reg, elow, ehigh, factor, snr, unit=None):
"""Bin the spectrum using a variable bin size, given a minimum
bin factor and a minimum signal to noise ratio.
:param inst: Instrument number.
:type inst: int
:param reg: Region number.
:type reg: int
:param elow: Start point of energy/channel interval.
:type elow: float
:param ehigh: End point of energy/channel interval.
:type ehigh: float
:param factor: Minimal binning factor
:type factor: int
:param snr: Minimal signal to noise for a data point after binning.
:type snr: float
:param unit: Unit of the energy/channel range, for example: 'kev', 'ev', 'ryd', 'j', 'hz', 'ang', 'nm'
:type unit: str
"""
bins = data.Bins()
bins.vbin(inst, reg, elow, ehigh, factor, snr, unit=unit)
# Watch (not implemented)
# ======================================================
# Quit command (automatic)
# ======================================================
def __del__(self):
"""Automatic quit command."""
pyspex_f2py.spexapi.api_final()
