#!/usr/bin/env python
from . import pyspex_f2py
from .utils import Elements
import os
import numpy
from astropy.table import QTable
from astropy import units as u
[docs]class Model:
"""Top class containing the entire model, containing all sectors and components.
:ivar nsector: Number of sectors in this model
:vartype nsector: int
:ivar sect: List of sector objects
:vartype sect: list
"""
def __init__(self):
self.nsector = 0 # Number of sectors in this model
self.sect = [] # List of sector objects
self.update()
[docs] def update(self):
"""Obtain the Sectors from SPEX."""
pyspex_f2py.api_model_sector.api_sector_update()
self.nsector = int(pyspex_f2py.api_model_sector.sect_nsector)
self.sect = []
for i in numpy.arange(self.nsector):
s = Sector()
s.update(i+1)
self.sect.append(s)
# Sector commands
[docs] def sector_new(self):
"""Create a new sector."""
ier = pyspex_f2py.api_model_sector.api_sector_new()
self.update()
return ier
[docs] def sector_copy(self, isect):
"""Copy an existing sector to a new one.
:param isect: Sector number to copy.
:type isect: int
"""
ier = pyspex_f2py.api_model_sector.api_sector_copy(float(isect))
self.update()
return ier
[docs] def sector_delete(self, isect):
"""Delete an existing sector.
:param isect: Sector number to delete.
:type isect: int
"""
ier = pyspex_f2py.api_model_sector.api_sector_delete(float(isect))
self.update()
return ier
# Component commands
[docs] def comp_new(self, name, isect=1):
"""Add new component to the model. By default in sector 1.
: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
"""
ier = pyspex_f2py.api_model_component.api_comp_new(name, isect=float(isect))
if ier == 0:
self.update()
return ier
[docs] def comp_delete(self, isect, icomp):
"""Delete component from sector.
:param isect: Sector number of the component to delete.
:type isect: int
:param icomp: Component number to delete.
:type icomp: int
"""
ier = pyspex_f2py.api_model_component.api_comp_delete(float(isect), float(icomp))
if ier == 0:
self.update()
return ier
[docs] def comp_set_rel(self, isect, icomp, rel):
"""Set component relation.
: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
"""
nrel = rel.size
ier = pyspex_f2py.api_model_component.api_comp_set_rel(float(isect), float(icomp), rel, float(nrel))
if ier == 0:
self.update()
return ier
# Parameter commands
[docs] def par_set(self, isect, icomp, name, value, thawn=False):
"""Set parameter to value and optionally indicate thawn or frozen 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
"""
ier = pyspex_f2py.api_model_parameter.api_par_set(float(isect), float(icomp), name, value, thawn=thawn)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_aset(self, isect, icomp, name, value):
"""Set a 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
"""
ier = pyspex_f2py.api_model_parameter.api_par_aset(float(isect), float(icomp), name, value)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_fix(self, isect, icomp, name):
"""Fix a parameter in 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
"""
cmd = 'par {:d} {:d} {:s} status frozen'.format(isect, icomp, name)
ier = pyspex_f2py.api_command(cmd)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_free(self, isect, icomp, name):
"""Free a parameter in 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
"""
cmd = 'par {:d} {:d} {:s} status thawn'.format(isect, icomp, name)
ier = pyspex_f2py.api_command(cmd)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_set_range(self, isect, icomp, name, rlow, rupp):
"""Set the fit range for 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
:param rlow: Lower limit of the parameter range.
:type rlow: float
:param rupp: Upper limit of the parameter range.
:type rupp: float
"""
ier = pyspex_f2py.api_model_parameter.api_par_set_range(float(isect), float(icomp), name, rlow, rupp)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_get(self, isect, icomp, name):
"""Get the parameter object for sector isect, component icomp and parameter with name.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param name: Name of the parameter
:type name: str
:rtype: model.Parameter
"""
found = False
for k in numpy.arange(self.sect[isect-1].comp[icomp-1].npar):
if self.sect[isect-1].comp[icomp-1].par[k].name == name:
parout = self.sect[isect-1].comp[icomp-1].par[k]
found = True
if not found:
print("Name not found in component.")
return
return parout
[docs] def par_couple(self, isect, icomp, iname, csect, ccomp, cname, factor):
"""Couple parameter (iname) to another parameter with cname, with a coupling 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
"""
ier = pyspex_f2py.api_model_parameter.api_par_couple(float(isect), float(icomp), iname, float(csect), float(ccomp),
cname, factor)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_decouple(self, isect, icomp, iname):
"""Decouple parameter.
: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
"""
ier = pyspex_f2py.api_model_parameter.api_par_decouple(float(isect), float(icomp), iname)
if ier == 0:
self.sect[isect-1].comp[icomp-1].update(isect, icomp)
return ier
[docs] def par_norm_set(self, ins, reg, value, status):
"""Set the instrument normalisation for an instrument number and region number.
: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 status: (Optional) Should the parameter be free (True) or frozen (False).
:type status: bool
"""
ier = pyspex_f2py.api_model_instnorm.api_norm_set(ins, reg, value, status)
return ier
[docs] def par_norm_get(self, ins, reg):
"""Get the instrument normalisation for instrument ins and region reg.
:param ins: Instrument number.
:type ins: int
:param reg: Region number.
:type reg: int
"""
ier = pyspex_f2py.api_model_instnorm.api_norm_get(ins, reg)
return float(pyspex_f2py.api_model_instnorm.norm_renorm), bool(pyspex_f2py.api_model_instnorm.norm_free)
[docs] def par_show_classic(self, option=''):
"""Display parameter overview in terminal through the Fortran interface.
:param option: Select which information should be shown (free/couple/flux/stat/corr/all)
:type option: str
"""
if option == 'free':
pyspex_f2py.api_command('par show free')
elif option == 'couple':
pyspex_f2py.api_command('par show couple')
elif option == 'flux':
pyspex_f2py.api_command('par show flux')
elif option == 'stat':
pyspex_f2py.api_command('par show stat')
elif option == 'corr':
pyspex_f2py.api_command('par show corr')
elif option == 'all':
pyspex_f2py.api_command('par show all')
else:
pyspex_f2py.api_command('par show')
return
[docs] def par_show(self, option=''):
"""Display parameter overview through the Python interface.
:param option: Select which information should be shown (free/couple/flux/stat/corr/all)
:type option: str
"""
if option == 'free':
self.par_show_free()
self.par_show_flux()
elif option == 'couple':
self.par_show_couple()
elif option == 'all':
self.par_show_param()
self.par_show_couple()
self.par_show_flux()
else:
self.par_show_param()
self.par_show_flux()
return
[docs] def par_show_param(self):
"""Display the parameters through the Python interface and create Astropy table.
:return: par_show_table
:rtype par_show_table: astropy.table.QTable
"""
print("---------------------------------------------------------------------------------------------------")
print("sect comp mod acro parameter with unit value status minimum maximum lsec lcom lpar")
rows = []
for isect in numpy.arange(self.nsector):
csect = self.sect[isect]
print("")
for icomp in numpy.arange(csect.ncomp):
ccomp = csect.comp[icomp]
for ipar in numpy.arange(ccomp.npar):
cpar = ccomp.par[ipar]
# Format the thawn or frozen text
if cpar.free:
stat = 'thawn'
else:
stat = 'frozen'
# Set the link data if parameter is linked:
if cpar.linked:
link = "{0:4d} {1:4d} {2:4d}".format(cpar.link_sector, cpar.link_comp, cpar.link_par)
else:
link = ''
if cpar.type == 'text':
print("{0:4d} {1:4d} {2:4s} {3:4s} {4:20s} {5:64s}"
"".format(csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, cpar.text))
rows.append((csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, 0., cpar.text,
"", 0, 0, 0, 0, 0))
else:
print("{0:4d} {1:4d} {2:4s} {3:4s} {4:20s} {5:#11.6G} {6:6s} {7:#7.2G} {8:#7.2G} {9:10s}"
"".format(csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, cpar.value, stat,
cpar.low, cpar.upp, link))
rows.append((csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, cpar.value, "", stat,
cpar.low, cpar.upp, cpar.link_sector, cpar.link_comp, cpar.link_par))
print("")
print("")
par_show_table = QTable(rows=rows, names=['Sect', 'Comp', 'Model', 'Acro', 'Par', 'Value', 'AValue', 'Status',
'Min', 'Max', 'Link_sect', 'Link_comp', 'Link_par'],
meta={'name': 'Par show'})
return par_show_table
[docs] def par_show_free(self):
"""Display the free parameters through the Python interface.
:return: par_show_table
:rtype par_show_table: astropy.table.QTable
"""
print("---------------------------------------------------------------------------------------------------")
print("sect comp mod acro parameter with unit value status minimum maximum lsec lcom lpar")
rows = []
for isect in numpy.arange(self.nsector):
csect = self.sect[isect]
print("")
for icomp in numpy.arange(csect.ncomp):
ccomp = csect.comp[icomp]
for ipar in numpy.arange(ccomp.npar):
cpar = ccomp.par[ipar]
# Format the thawn or frozen text
if cpar.free:
stat = 'thawn'
# Set the link data if parameter is linked:
if cpar.linked:
link = "{0:4d} {1:4d} {2:4d}".format(cpar.link_sector, cpar.link_comp, cpar.link_par)
else:
link = ''
if cpar.type == 'text':
print("{0:4d} {1:4d} {2:4s} {3:4s} {4:20s} {5:64s}"
"".format(csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, cpar.text))
rows.append((csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, 0., cpar.text,
"", 0, 0, 0, 0, 0))
else:
print("{0:4d} {1:4d} {2:4s} {3:4s} {4:20s} {5:#11.6G} {6:6s} {7:#7.2G} {8:#7.2G} "
"{9:10s}".format(csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc,
cpar.value, stat, cpar.low, cpar.upp, link))
rows.append((csect.index, ccomp.index, ccomp.name, cpar.name, cpar.desc, cpar.value, "",
stat, cpar.low, cpar.upp, cpar.link_sector, cpar.link_comp, cpar.link_par))
print("")
print("")
par_show_table = QTable(rows=rows, names=['Sect', 'Comp', 'Model', 'Acro', 'Par', 'Value', 'AValue', 'Status',
'Min', 'Max', 'Link_sect', 'Link_comp', 'Link_par'],
meta={'name': 'Par show free'})
return par_show_table
[docs] def par_show_couple(self):
"""Display the coupled parameters through the Python interface and as Astropy table.
:return: par_show_table
:rtype par_show_table: astropy.table.QTable
"""
print("---------------------------------------------------------------------------------------------------")
print("Coupled parameter | coupling |Parameter on which it depends")
print("sect comp mod para | factor |sect comp mod para")
rows = []
for isect in numpy.arange(self.nsector):
csect = self.sect[isect]
print("")
for icomp in numpy.arange(csect.ncomp):
ccomp = csect.comp[icomp]
for ipar in numpy.arange(ccomp.npar):
cpar = ccomp.par[ipar]
# Set the link data if parameter is linked:
if cpar.linked:
print("{0:4d} {1:4d} {2:4s} {3:4s} | {4:#11.6G} |{5:4d} {6:4d} {7:4s} {8:4s}"
"".format(csect.index, ccomp.index, ccomp.name, cpar.name, cpar.coupling,
cpar.link_sector, cpar.link_comp, ccomp.name,
self.sect[cpar.link_sector-1].comp[cpar.link_comp-1].par[cpar.link_par-1].name))
rows.append((csect.index, ccomp.index, ccomp.name, cpar.name, cpar.coupling,
cpar.link_sector, cpar.link_comp, ccomp.name,
self.sect[cpar.link_sector-1].comp[cpar.link_comp-1].par[cpar.link_par-1].name))
print("")
par_show_table = QTable(rows=rows, names=['Sect', 'Comp', 'Model', 'Acro', 'Factor', 'CSect', 'CComp',
'CModel', 'CAcro'],
meta={'name': 'Par show couple'})
return par_show_table
[docs] def par_show_flux(self):
"""Display the flux information of the model components.
:return: par_show_table
:rtype par_show_table: astropy.table.QTable
"""
for isect in numpy.arange(self.nsector):
csect = self.sect[isect]
for icomp in numpy.arange(csect.ncomp):
if csect.comp[icomp].add:
elim = csect.comp[icomp].flux.elimflux
print("---------------------------------------------------------------------------------------------------")
print("Fluxes and restframe luminosities between {0:#11.6G} and {1:#11.6G}".format(elim[0], elim[1]))
print("")
print("sect comp mod photon flux energy flux nr of photons luminosity")
print(" (phot/m**2/s) (W/m**2) (photons/s) (W)")
rows = []
for isect in numpy.arange(self.nsector):
csect = self.sect[isect]
for icomp in numpy.arange(csect.ncomp):
cflux = csect.comp[icomp].flux
ccomp = csect.comp[icomp]
if ccomp.add:
print("{0:4d} {1:4d} {2:4s} {3:#11.6G} {4:#11.6G} {5:#11.6G} {6:#11.6G}"
"".format(csect.index, ccomp.index, ccomp.name, cflux.photflux.value,
cflux.enerflux.value, cflux.photlum.value, cflux.enerlum.value))
rows.append((csect.index, ccomp.index, ccomp.name, cflux.photflux,
cflux.enerflux, cflux.photlum, cflux.enerlum))
print("")
par_show_table = QTable(rows=rows, names=['Sect', 'Comp', 'Model', 'Phot_flux', 'Ener_flux', 'Nr_phot', 'Lum'],
meta={'name': 'Par show flux'})
return par_show_table
[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
"""
file = os.path.splitext(comfile)[0]
ext = os.path.splitext(comfile)[1]
if ext != '.com':
print("Error: File does not have the .com extension.")
return 1
if overwrite:
com = 'par write {:s} overwrite'.format(file)
else:
com = 'par write {:s}'.format(file)
ier = pyspex_f2py.api_command(com)
return ier
[docs]class Sector:
"""Class describing a Sector.
:ivar index: Sector number.
:vartype index: int
:ivar ncomp: Number of components in sector.
:vartype ncomp: int
:ivar comp: List of component objects.
:vartype comp: list
:ivar distance: Distance for this sector.
:vartype distance: float
:ivar spectrum: Model spectrum for this sector
:vartype spectrum: pyspex.model.Spectrum
"""
def __init__(self):
self.index = 0 # Sector number.
self.ncomp = 0 # Number of components in sector.
self.comp = [] # List of component objects.
self.distance = 0. # Distance for this sector.
self.spectrum = Spectrum() # The spectrum for this sector
[docs] def update(self, isect):
"""Obtain the information from SPEX for sector isect.
:param isect: Sector number.
:type isect: int
"""
self.index = int(pyspex_f2py.api_model_sector.sect_index[isect-1])
self.ncomp = int(pyspex_f2py.api_model_sector.sect_ncomp[isect-1])
self.distance = float(pyspex_f2py.api_model_sector.sect_distance[isect-1])
self.spectrum.get(isect)
self.comp = []
for i in numpy.arange(self.ncomp):
c = Component()
c.update(self.index, i+1)
self.comp.append(c)
[docs]class Component:
"""Class containing a model component.
:ivar index: Component number.
:vartype index: int
:ivar umodel: Model ID number.
:vartype umodel: int
:ivar add: Is this an additive component?
:vartype add: bool
:ivar mul: Is this a multiplicative component?
:vartype mul: bool
:ivar operator: True for complex components.
:vartype operator: bool
:ivar npar: Number of parameters in components.
:vartype npar: int
:ivar name: Name of model.
:vartype name: str
:ivar nmul: Number of elements in addmul (below).
:vartype nmul: int
:ivar addmul: Array with the numbers of the multiplicative components to multiply with.
:vartype addmul: numpy.ndarray
:ivar par: List of parameter objects.
:vartype par: list
"""
def __init__(self):
self.index = 0 # Component number.
self.umodel = 0 # Model ID number.
self.add = False # Is this an additive component?
self.mul = False # Is this a multiplicative component?
self.operator = False # True for complex components.
self.npar = 0 # Number of parameters in components.
self.name = '' # Name of model.
self.nmul = 0 # Number of elements in addmul (below).
self.addmul = numpy.array([]) # Array with the numbers of the multiplicative components to multiply with.
self.par = [] # List of parameter objects.
self.flux = Fluxes() # Flux for this component
[docs] def update(self, isect, icomp):
"""Obtain the information from SPEX for sector isect and component icomp.
:param isect: Sector number of the component.
:type isect: int
:param icomp: Component number of the component.
:type icomp: int
"""
pyspex_f2py.api_model_component.api_comp_update(float(isect))
# Update the name and relational information for this component
pyspex_f2py.api_model_component.api_comp_get_rel(float(isect), float(icomp))
pyspex_f2py.api_model_component.api_comp_get_name(float(isect), float(icomp))
self.index = int(pyspex_f2py.api_model_component.comp_index[icomp-1])
self.umodel = int(pyspex_f2py.api_model_component.comp_umodel[icomp-1])
self.add = bool(pyspex_f2py.api_model_component.comp_add[icomp-1])
self.mul = bool(pyspex_f2py.api_model_component.comp_mul[icomp-1])
self.operator = bool(pyspex_f2py.api_model_component.comp_operator[icomp-1])
self.npar = int(pyspex_f2py.api_model_component.comp_npar[icomp-1])
self.name = pyspex_f2py.api_model_component.comp_name.view('S4').tobytes().decode('utf-8')
self.nmul = int(pyspex_f2py.api_model_component.comp_nmul)
self.addmul = numpy.zeros(self.nmul, dtype=int)
for i in numpy.arange(self.nmul):
self.addmul[i] = int(pyspex_f2py.api_model_component.comp_addmul[i])
self.par = []
pyspex_f2py.api_model_component.api_comp_final()
for i in numpy.arange(self.npar):
p = Parameter()
p.update(isect, icomp, i+1)
self.par.append(p)
if self.add:
self.flux.get(isect, icomp)
[docs]class Parameter:
"""Class for model parameters.
:ivar isect: Sector number
:vartype index: int
:ivar icomp: Component number
:vartype index: int
:ivar index: Parameter number
:vartype index: int
:ivar name: Name of parameter
:vartype name: str
:ivar type: Type of parameter (norm, abun, fitp, cons or text)
:vartype type: str
:ivar desc: Description of parameter
:vartype desc: str
:ivar value: Value
:vartype value: float
:ivar low: Lower boundary of parameter range
:vartype low: float
:ivar upp: Upper boundary of parameter range
:vartype upp: float
:ivar err_low: Lower error boundary
:vartype err_low: float
:ivar err_upp: Upper error boundary
:vartype err_upp: float
:ivar free: True if parameter is free
:vartype free: bool
:ivar linked: True if parameter is linked
:vartype linked: bool
:ivar link_sector: Sector to which parameter is linked
:vartype link_sector: int
:ivar link_comp: Component to which parameter is linked
:vartype link_comp: int
:ivar link_par: Parameter to which parameter is linked
:vartype link_par: int
:ivar coupling: Coupling factor
:vartype coupling: float
"""
def __init__(self):
self.isect = 0 # Sector number
self.icomp = 0 # Component number
self.index = 0 # Parameter number
self.name = '' # Name of parameter
self.type = '' # Type of parameter (norm, abun, fitp, cons or text)
self.desc = '' # Description of parameter
self.value = 0. # Value
self.text = '' # Text value
self.low = 0. # Lower boundary of parameter range
self.upp = 0. # Upper boundary of parameter range
self.err_low = 0. # Lower error boundary
self.err_upp = 0. # Upper error boundary
self.free = False # True if parameter is free
self.linked = False # True if parameter is linked
self.link_sector = 0 # Sector to which parameter is linked
self.link_comp = 0 # Component to which parameter is linked
self.link_par = 0 # Parameter to which parameter is linked
self.coupling = 1.0 # Coupling factor
[docs] def update(self, isect, icomp, ipar):
"""Obtain the parameter values from SPEX memory.
:param isect: Sector number of the parameter.
:type isect: int
:param icomp: Component number of the parameter.
:type icomp: int
:param ipar: Parameter number.
:type ipar: int
"""
pyspex_f2py.api_model_parameter.api_par_update(float(isect), float(icomp), float(ipar))
self.isect = isect
self.icomp = icomp
self.index = int(pyspex_f2py.api_model_parameter.par_index)
self.name = pyspex_f2py.api_model_parameter.par_acronym.view('S4').tobytes().decode('utf-8')
self.name = self.name.strip()
self.type = pyspex_f2py.api_model_parameter.par_type.view('S4').tobytes().decode('utf-8')
self.type = self.type.strip()
self.desc = pyspex_f2py.api_model_parameter.par_name.view('S20').tobytes().decode('utf-8')
self.desc = self.desc.strip()
self.text = pyspex_f2py.api_model_parameter.par_text.view('S128').tobytes().decode('utf-8')
self.text = self.text.strip()
self.value = float(pyspex_f2py.api_model_parameter.par_value)
self.low = float(pyspex_f2py.api_model_parameter.par_low)
self.upp = float(pyspex_f2py.api_model_parameter.par_upp)
self.err_low = float(pyspex_f2py.api_model_parameter.par_err_low)
self.err_upp = float(pyspex_f2py.api_model_parameter.par_err_upp)
self.free = bool(pyspex_f2py.api_model_parameter.par_free)
self.linked = bool(pyspex_f2py.api_model_parameter.par_linked)
self.link_sector = int(pyspex_f2py.api_model_parameter.par_link_sector)
self.link_comp = int(pyspex_f2py.api_model_parameter.par_link_comp)
self.link_par = int(pyspex_f2py.api_model_parameter.par_link_par)
self.coupling = float(pyspex_f2py.api_model_parameter.par_coupling)
[docs]class Abundance:
"""This class manages the SPEX abundance setting.
:ivar index: Abundance index in list
:vartype index: int
:ivar ref: Reference to abundance set (string)
:vartype ref: str
:ivar list: The available abundance sets.
:vartype list: tuple
"""
def __init__(self):
self.index = 0 # Abundance index in list
self.ref = "" # Reference to abundance set (string)
self.list = ("Lodders et al. (2009)", "Allen (1973)", "Ross & Aller (1976)",
"Grevesse et al. (1992)", "Grevesse & Sauval (1998)", "Lodders proto-Solar (2003)",
"Lodders Solar Photospheric (2003)", "Anders & Grevesse (1989)",
"Asplund et al. (2009)") # The available abundance sets.
[docs] def set(self, abun):
"""Set the abundance in SPEX to another set.
:param abun: Abbreviation of the abundance set.
:type abun: str
"""
ier = pyspex_f2py.api_model_abundance.api_abun_set(abun)
self.update()
return ier
[docs] def get(self):
"""Get the current Abundance setting (reference)."""
self.update()
return self.list[self.index]
[docs] def update(self):
"""Update the abundance setting in pyspex."""
pyspex_f2py.api_model_abundance.api_abun_update()
self.index = int(pyspex_f2py.api_model_abundance.abun_index)
self.ref = pyspex_f2py.api_model_abundance.abun_ref.view('S48').tobytes().decode('utf-8')
[docs]class Distance:
"""This class is used to set and get the distances in SPEX.
:ivar m: Distance in meter
:vartype m: astropy.units.quantity.Quantity
:ivar au: Distance in Astronomical Units
:vartype au: astropy.units.quantity.Quantity
:ivar ly: Distance in light year
:vartype ly: astropy.units.quantity.Quantity
:ivar pc: Distance in parsec
:vartype pc: astropy.units.quantity.Quantity
:ivar kpc: Distance in kiloparsec
:vartype kpc: astropy.units.quantity.Quantity
:ivar mpc: Distance in megaparsec
:vartype mpc: astropy.units.quantity.Quantity
:ivar z: Distance in redshift
:vartype z: astropy.units.quantity.Quantity
:ivar cz: Distance in cz
:vartype cz: astropy.units.quantity.Quantity
:ivar age: Lookback time (yr)
:vartype age: astropy.units.quantity.Quantity
:ivar h0: Hubble constant (km/s/Mpc)
:vartype h0: astropy.units.quantity.Quantity
:ivar omega_m: Omega Matter
:vartype omega_m: astropy.units.quantity.Quantity
:ivar omega_l: Omega Lambda
:vartype omega_l: astropy.units.quantity.Quantity
:ivar omega_r: Omega R
:vartype omega_r: astropy.units.quantity.Quantity
"""
def __init__(self):
self.m = 0. * u.m # Distance in meter
self.au = 0. * u.AU # Distance in Astronomical Units
self.ly = 0. * u.lyr # Distance in light year
self.pc = 0. * u.pc # Distance in parsec
self.kpc = 0. * u.kpc # Distance in kiloparsec
self.mpc = 0. * u.Mpc # Distance in megaparsec
self.z = 0. * u.dimensionless_unscaled # Distance in redshift
self.cz = 0. * u.dimensionless_unscaled # Distance in cz
self.age = 0. * u.yr # Lookback time (yr)
# Cosmology
self.h0 = 70.0 * u.km / (u.s * u.Mpc) # Hubble constant (km/s/Mpc)
self.omega_m = 0.3 * u.dimensionless_unscaled # Omega Matter
self.omega_l = 0.7 * u.dimensionless_unscaled # Omega Lambda
self.omega_r = 0.0 * u.dimensionless_unscaled # Omega R
[docs] def set(self, isect, dist, unit):
"""Set the distance in units for a particular sector (isect).
:param isect: Sector number.
:type isect: int
:param dist: Distance value.
:type dist: float
:param unit: Unit of the distance, for example: 'm', 'au', 'ly', 'pc', 'kpc', 'mpc', 'z', 'cz'.
:type unit: str
"""
ier = 0
ier = pyspex_f2py.api_model_distance.api_distance_set(float(isect), dist, unit)
if ier == 0:
self.m = pyspex_f2py.api_model_distance.dist_m * u.m
self.au = pyspex_f2py.api_model_distance.dist_au * u.AU
self.ly = pyspex_f2py.api_model_distance.dist_ly * u.lyr
self.pc = pyspex_f2py.api_model_distance.dist_pc * u.pc
self.kpc = pyspex_f2py.api_model_distance.dist_kpc * u.kpc
self.mpc = pyspex_f2py.api_model_distance.dist_mpc * u.Mpc
self.z = pyspex_f2py.api_model_distance.dist_z * u.dimensionless_unscaled
self.cz = pyspex_f2py.api_model_distance.dist_cz * u.dimensionless_unscaled
self.age = pyspex_f2py.api_model_distance.dist_age * u.yr
else:
print("Error: Distances not set correctly.")
[docs] def get(self, isect):
"""Get the distances for a sector from SPEX.
:param isect: Sector number.
:type isect: int
"""
ier = pyspex_f2py.api_model_distance.api_distance_calc(int(isect))
if ier == 0:
self.m = pyspex_f2py.api_model_distance.dist_m * u.m
self.au = pyspex_f2py.api_model_distance.dist_au * u.AU
self.ly = pyspex_f2py.api_model_distance.dist_ly * u.lyr
self.pc = pyspex_f2py.api_model_distance.dist_pc * u.pc
self.kpc = pyspex_f2py.api_model_distance.dist_kpc * u.kpc
self.mpc = pyspex_f2py.api_model_distance.dist_mpc * u.Mpc
self.z = pyspex_f2py.api_model_distance.dist_z * u.dimensionless_unscaled
self.cz = pyspex_f2py.api_model_distance.dist_cz * u.dimensionless_unscaled
self.age = pyspex_f2py.api_model_distance.dist_age * u.yr
else:
print("Error: Distances could not be retrieved.")
[docs] def set_cosmo(self, h0, omega_m, omega_l, omega_r):
"""Set the cosmological constants for the distance calculation.
:param h0: Hubble constant (km/s/Mpc).
:type h0: astropy.units.quantity.Quantity
:param omega_m: Omega matter.
:type omega_m: astropy.units.quantity.Quantity
:param omega_l: Omega lambda.
:type omega_l: astropy.units.quantity.Quantity
:param omega_r: Omega R.
:type omega_r: astropy.units.quantity.Quantity
"""
ier = pyspex_f2py.api_model_distance.api_distance_set_cosmo(h0, omega_m, omega_l, omega_r)
if ier == 0:
self.h0 = h0 * u.km / (u.s * u.Mpc)
self.omega_m = omega_m * u.dimensionless_unscaled
self.omega_l = omega_l * u.dimensionless_unscaled
self.omega_r = omega_r * u.dimensionless_unscaled
else:
print("Error: Cosmology not set.")
[docs]class Egrid:
"""This class is used to specify the model energy grid.
:ivar nbins: Number of bins
:vartype nbins: int
:ivar energy: Centroids of bins (Energy, keV)
:vartype energy: astropy.units.quantity.Quantity
:ivar energy_upper: Upper boundaries of bins (Energy, keV)
:vartype energy_upper: astropy.units.quantity.Quantity
:ivar energy_width: Widths of bins (Energy, keV)
:vartype energy_width: astropy.units.quantity.Quantity
"""
def __init__(self):
self.nbins = 0 # Number of bins
self.energy = numpy.array([0.], dtype=float) * u.keV # Centroids of bins
self.energy_upper = numpy.array([], dtype=float) * u.keV # Upper boundaries of bins
self.energy_width = numpy.array([], dtype=float) * u.keV # Widths of bins
def allocate(self):
self.energy = numpy.zeros(self.nbins, dtype=float) * u.keV
self.energy_upper = numpy.zeros(self.nbins, dtype=float) * u.keV
self.energy_width = numpy.zeros(self.nbins, dtype=float) * u.keV
return 0
[docs] def set(self, elow, ehigh, nbins, unit, log):
"""Set egrid using limits and number of bins.
: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
"""
ier = pyspex_f2py.api_model_egrid.api_egrid_set(elow, ehigh, float(nbins), unit, log)
return ier
[docs] def set_step(self, elow, ehigh, step, unit, log):
"""Set egrid using limits and step size.
: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
"""
ier = pyspex_f2py.api_model_egrid.api_egrid_set_step(elow, ehigh, step, unit, log)
return ier
[docs] def get(self):
"""Get the current energy grid from SPEX."""
pyspex_f2py.api_model_egrid.api_egrid_get()
self.nbins = pyspex_f2py.api_model_egrid.nbins
self.allocate()
for i in numpy.arange(self.nbins):
self.energy[i] = pyspex_f2py.api_model_egrid.energy_cent[i] * u.keV
self.energy_upper[i] = pyspex_f2py.api_model_egrid.energy_bupp[i] * u.keV
self.energy_width[i] = pyspex_f2py.api_model_egrid.energy_width[i] * u.keV
pyspex_f2py.api_model_egrid.api_egrid_final()
return 0
[docs] def save(self, savefile):
"""Save the energy grid to a file named savefile (extension: .egr).
:param savefile: Filename to save the energy grid to (including .egr extension)
:type savefile: str
"""
file = os.path.splitext(savefile)[0]
ext = os.path.splitext(savefile)[1]
if ext != '.egr':
print("Error: File does not have the .egr extension.")
return 1
ier = pyspex_f2py.api_model_egrid.api_egrid_save(file)
return ier
[docs] def read(self, readfile):
"""Read the energy grid from a file named readfile (extension: .egr).
:param readfile: Filename to read the energy grid from (including .egr extension)
:type readfile: str
"""
file = os.path.splitext(readfile)[0]
ext = os.path.splitext(readfile)[1]
if ext != '.egr':
print("Error: File does not have the .egr extension.")
return 1
ier = pyspex_f2py.api_model_egrid.api_egrid_read(file)
return ier
[docs] def grid(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 (keV).
:type ebounds: numpy.ndarray
"""
if type(ebounds) is numpy.ndarray:
ier = pyspex_f2py.api_model_egrid.api_egrid_grid(ebounds)
else:
print("Error: Input object is not a numpy array.")
return 1
return ier
[docs]class Fluxes:
"""This class is used to calculate fluxes and luminosities of spectra in a spectral band.
:ivar sector: Sector number of flux calculation
:vartype sector: int
:ivar component: Component number of flux calculation
:vartype component: int
:ivar photflux: Photon flux (phot/m**2/s)
:vartype photflux: astropy.units.quantity.Quantity
:ivar enerflux: Energy flux (W/m**2)
:vartype enerflux: astropy.units.quantity.Quantity
:ivar photlum: Photon luminosity (photons/s)
:vartype photlum: astropy.units.quantity.Quantity
:ivar enerlum: Energy luminosity (W)
:vartype enerlum: astropy.units.quantity.Quantity
:ivar elimflux: Flux energy limits (keV)
:vartype elimflux: astropy.units.quantity.Quantity
"""
def __init__(self):
self.sector = 0 # Sector number of flux calculation
self.component = 0 # Component number of flux calculation
self.photflux = 0. * u.ph / (u.m**2 * u.s) # Photon flux (phot/m**2/s)
self.enerflux = 0. * u.W / u.m**2 # Energy flux (W/m**2)
self.photlum = 0. * u.ph / u.s # Photon luminosity (photons/s)
self.enerlum = 0. * u.W # Energy luminosity (W)
self.elimflux = numpy.zeros(2, dtype=float) * u.keV #: Flux energy limits (keV)
[docs] def elim(self, elow, ehigh, unit):
"""Set the energy limits for the flux and luminosity calculation.
: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
"""
ier = pyspex_f2py.api_model_fluxes.api_flux_elim(elow, ehigh, unit)
if ier != 0:
print("Error: Setting of energy limits failed.")
return ier
[docs] def get(self, isect, icomp):
"""Get the flux and luminosity from SPEX for a given sector and component number.
: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
"""
pyspex_f2py.api_model_fluxes.api_flux_get(float(isect), float(icomp))
self.sector = isect
self.component = icomp
self.photflux = pyspex_f2py.api_model_fluxes.flux_photflux * u.ph / (u.m**2 * u.s)
self.enerflux = pyspex_f2py.api_model_fluxes.flux_enerflux * u.W / u.m**2
self.photlum = pyspex_f2py.api_model_fluxes.flux_photlum * u.ph / u.s
self.enerlum = pyspex_f2py.api_model_fluxes.flux_enerlum * u.W
self.elimflux = pyspex_f2py.api_model_fluxes.flux_elim * u.keV
return
[docs] def calc(self, isect, icomp):
"""Calculate and get flux and luminosity from SPEX for a given sector and component number.
:param isect: Sector number of the component to calculate.
:type isect: int
:param icomp: Component number to calculate.
:type icomp: int
"""
pyspex_f2py.api_model.api_model_calculate()
self.get(isect, icomp)
return
[docs]class Ibal:
"""This class manages the SPEX ionisation balance setting.
:ivar index: Index number for list of ionisation balance sets.
:vartype index: int
:ivar ref: Reference to the current ionisation balance.
:vartype ref: str
:ivar list: List of available ionisation balance data.
:vartype list: tuple
"""
def __init__(self):
self.index = 0 # Index number for list.
self.ref = "" # Reference to ibal set (string)
self.list = ("", "Arnaud & Raymond (1992) for Fe, Arnaud & Rothenflug (1985) for other elements",
"Arnaud & Rothenflug (1985)", "Arnaud & Raymond (1992) for Fe, Arnaud & Rothenflug (1985) "
"for other elements", "Old Bryans et al. data (NOT recommended)", "Bryans et al. (2009)",
"Urdampilleta et al. (2017)") # List of available ibal data.
[docs] def set(self, ibal):
"""Set the abundance in SPEX to another set.
:param ibal: Abbreviation of the reference to the ionisation balance.
:type ibal: str
"""
ier = pyspex_f2py.api_model_ibal.api_ibal_set(ibal)
self.update()
return ier
[docs] def get(self):
"""Get the current Abundance setting (reference)."""
self.update()
return self.list[self.index]
[docs] def update(self):
"""Update the abundance setting in pyspex."""
pyspex_f2py.api_model_ibal.api_ibal_update()
self.index = int(pyspex_f2py.api_model_ibal.ibal_index)
self.ref = pyspex_f2py.api_model_ibal.ibal_ref.view('S80').tobytes().decode('utf-8')
# Classes for ions in the model calculation
class Ion:
"""Class containing the settings of an ion.
:ivar z: Atomic number.
:vartype z: int
:ivar i: Ionisation stage.
:vartype i: int
:ivar use: Ion used in calculation.
:vartype use: bool
:ivar new: Use ion in the new calculation.
:vartype new: bool
:ivar nmax: Maximum principle quantum number.
:vartype nmax: int
:ivar lmax: Maximum angular momentum number.
:vartype lmax: int
"""
def __init__(self, z, i):
self.z = z # Atomic number
self.i = i # Ionisation stage
self.use = True # Ion used in calculation
self.new = True # Use ion in the new calculation
self.nmax = 10 # Maximum principle quantum number
self.lmax = 10 # Maximum angular momentum number
self.update()
def update(self):
"""Update the information for an ion with atomic number z and ionisation stage i."""
(self.nmax, self.lmax, use, new) = pyspex_f2py.api_model_ions.api_ion_get(float(self.z), float(self.i))
if use == 1:
self.use = True
else:
self.use = False
if new == 1:
self.new = True
else:
self.new = False
[docs]class Aerror:
"""Class for calculating the uncertainty in a parameter due to atomic data uncertainties.
:ivar aerr: Error value for the requested parameter.
:vartype aerr: float
"""
def __init__(self):
"""Initialize Aerror."""
self.aerr = 0.
[docs] def get(self, isect, icomp, iname, shell=0):
"""Obtain the atomic data error for parameter `iname` in sector `isect` and component number `icomp`.
: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 shell: Shell number (L-shell = 1, K-shell = 2)
:type shell: int
"""
ier = 0
if shell == 1 or shell == 2:
ier = pyspex_f2py.api_model_aerror.api_aerror_get_shell(isect, icomp, iname, shell)
else:
ier = pyspex_f2py.api_model_aerror.api_aerror_get(isect, icomp, iname)
if ier == 0:
self.aerr = float(pyspex_f2py.api_model_aerror.aerr_error)
else:
return
return self.aerr
class Atom:
"""Class containing all ions for an atom.
:ivar nion: Number of ions for this atom.
:vartype nion: int
:ivar z: Atomic number for this atom.
:vartype z: int
:ivar ions: List of ion objects with ion information.
:vartype ions: list
"""
def __init__(self, z):
"""Initialize the information for an atom.
:param z: Atomic number
:type z: int
"""
self.nion = z + 1
self.z = z
self.ions = []
self.update()
def update(self):
"""Update the properties for this atom."""
for ii in numpy.arange(self.nion):
ion = Ion(self.z, ii+1)
self.ions.append(ion)
[docs]class Ions:
"""Class to manage the ions taken into account in the model calculation.
:ivar nz: Total number of atoms considered.
:vartype nz: int
:ivar atoms: List of atoms (with information about each ion).
:vartype atoms: list
"""
def __init__(self):
self.nz = 30
self.atoms = []
self.update()
[docs] def update(self):
"""Update the properties of all atoms."""
self.atoms = []
for iz in numpy.arange(self.nz):
atom = Atom(iz+1)
self.atoms.append(atom)
[docs] def show(self):
"""Show the settings for the ions."""
names = Elements()
for iz in numpy.arange(self.nz):
for jz in numpy.arange(self.atoms[iz].nion):
if self.atoms[iz].ions[jz].new:
new = 'v3'
else:
new = 'v2'
nmax = self.atoms[iz].ions[jz].nmax
lmax = self.atoms[iz].ions[jz].lmax
if self.atoms[iz].ions[jz].use:
print("{:9} {:2} {:3d} {:3d}".format(names.ion(iz+1, jz+1), new, nmax, lmax))
# Set nmax
[docs] def nmax_all(self, nmax):
"""Set all ions to maximum quantum number nmax.
:param nmax: Maximum principle quantum number to use.
:type nmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_nmax_all(float(nmax))
self.update()
return ier
[docs] def nmax_iso(self, iso, nmax):
"""Set this iso-electronic sequence to maximum quantum number nmax.
:param iso: Iso-electronic sequence.
:type iso: int
:param nmax: Maximum principle quantum number to use.
:type nmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_nmax_iso(float(nmax), float(iso))
self.update()
return ier
[docs] def nmax_z(self, z, nmax):
"""Set all ions of this element to maximum quantum number nmax.
:param z: Atomic number
:type z: int
:param nmax: Maximum principle quantum number to use.
:type nmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_nmax_z(float(nmax), float(z))
self.update()
return ier
[docs] def nmax_ion(self, z, i, nmax):
"""Set this ion to maximum quantum number nmax.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
:param nmax: Maximum principle quantum number to use.
:type nmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_nmax_ion(float(nmax), float(z), float(i))
self.update()
return ier
# Set lmax
[docs] def lmax_all(self, lmax):
"""Set all ions to maximum angular momentum lmax.
:param lmax: Maximum angular momentum to use.
:type lmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_lmax_all(float(lmax))
self.update()
return ier
[docs] def lmax_iso(self, iso, lmax):
"""Set this iso-electronic sequence to maximum angular momentum lmax.
:param iso: Iso-electronic sequence.
:type iso: int
:param lmax: Maximum angular momentum to use.
:type lmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_lmax_iso(float(lmax), float(iso))
self.update()
return ier
[docs] def lmax_z(self, z, lmax):
"""Set all ions of this element to maximum angular momentum lmax.
:param z: Atomic number
:type z: int
:param lmax: Maximum angular momentum to use.
:type lmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_lmax_z(float(lmax), float(z))
self.update()
return ier
[docs] def lmax_ion(self, z, i, lmax):
"""Set this ion to maximum angular momentum lmax.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
:param lmax: Maximum angular momentum to use.
:type lmax: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_lmax_ion(float(lmax), float(z), float(i))
self.update()
return ier
# Set use
[docs] def use_all(self):
"""Use all ions."""
ier = pyspex_f2py.api_model_ions.api_ion_use_all()
self.update()
return ier
[docs] def use_iso(self, iso):
"""Use this iso-electronic sequence.
:param iso: Iso-electronic sequence.
:type iso: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_use_iso(float(iso))
self.update()
return ier
[docs] def use_z(self, z):
"""Use all ions of this element.
:param z: Atomic number
:type z: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_use_z(float(z))
self.update()
return ier
[docs] def use_ion(self, z, i):
"""Use this ion.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_use_ion(float(z), float(i))
self.update()
return ier
# Set ignore
[docs] def ignore_all(self):
"""Ignore all ions."""
ier = pyspex_f2py.api_model_ions.api_ion_ignore_all()
self.update()
return ier
[docs] def ignore_iso(self, iso):
"""Ignore this iso-electronic sequence.
:param iso: Iso-electronic sequence.
:type iso: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_ignore_iso(float(iso))
self.update()
return ier
[docs] def ignore_z(self, z):
"""Ignore all ions of this element.
:param z: Atomic number
:type z: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_ignore_z(float(z))
self.update()
return ier
[docs] def ignore_ion(self, z, i):
"""Ignore this ion.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_ignore_ion(float(z), float(i))
self.update()
return ier
# Set new
[docs] def new_all(self):
"""Use new calculations for all ions."""
ier = pyspex_f2py.api_model_ions.api_ion_new_all()
self.update()
return ier
[docs] def new_iso(self, iso):
"""Use new calculations for this iso-electronic sequence.
:param iso: Iso-electronic sequence.
:type iso: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_new_iso(float(iso))
self.update()
return ier
[docs] def new_z(self, z):
"""Use new calculations for all ions of this element.
:param z: Atomic number
:type z: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_new_z(float(z))
self.update()
return ier
[docs] def new_ion(self, z, i):
"""Use new calculations for this ion.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_new_ion(float(z), float(i))
self.update()
return ier
# Set old
[docs] def old_all(self):
"""Use old calculations for all ions."""
ier = pyspex_f2py.api_model_ions.api_ion_old_all()
self.update()
return ier
[docs] def old_iso(self, iso):
"""Use old calculations for this iso-electronic sequence.
:param iso: Iso-electronic sequence.
:type iso: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_old_iso(float(iso))
self.update()
return ier
[docs] def old_z(self, z):
"""Use old calculations for all ions of this element.
:param z: Atomic number
:type z: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_old_z(float(z))
self.update()
return ier
[docs] def old_ion(self, z, i):
"""Use old calculations for this ion.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_qc_ion(float(z), float(i))
self.update()
return ier
# Set qc
[docs] def qc_all(self):
"""Use qc calculations for all ions."""
ier = pyspex_f2py.api_model_ions.api_ion_qc_all()
self.update()
return ier
[docs] def qc_iso(self, iso):
"""Use qc calculations for this iso-electronic sequence.
:param iso: Iso-electronic sequence.
:type iso: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_qc_iso(float(iso))
self.update()
return ier
[docs] def qc_z(self, z):
"""Use qc calculations for all ions of this element.
:param z: Atomic number
:type z: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_qc_z(float(z))
self.update()
return ier
[docs] def qc_ion(self, z, i):
"""Use qc calculations for this ion.
:param z: Atomic number
:type z: int
:param i: Ion number
:type i: int
"""
ier = pyspex_f2py.api_model_ions.api_ion_qc_ion(float(z), float(i))
self.update()
return ier
[docs]class Spectrum:
"""This class obtains and stores the model spectrum from SPEX.
:ivar nbins: Number of bins
:vartype nbins: int
:ivar energy: Centroids of bins (Energy, keV)
:vartype energy: astropy.units.quantity.Quantity
:ivar energy_upper: Upper boundaries of bins (Energy, keV)
:vartype energy_upper: astropy.units.quantity.Quantity
:ivar energy_width: Widths of bins (Energy, keV)
:vartype energy_width: astropy.units.quantity.Quantity
:ivar spectrum: Spectrum of bins (in ph/s/m**2/bin at observatory)
:vartype spectrum: astropy.units.quantity.Quantity
:ivar luminosity: Spectrum of bins (in 10^44 ph/s/keV at source distance)
:vartype luminosity: astropy.units.quantity.Quantity
:ivar table: Astropy QTable containing spectrum.
:vartype table: astropy.table.QTable
"""
def __init__(self):
self.nbins = 0 # Total number of bins
self.energy = numpy.array([0.], dtype=float) * u.keV # Centroids of bins
self.energy_upper = numpy.array([0.], dtype=float) * u.keV # Upper boundaries of bins
self.energy_width = numpy.array([0.], dtype=float) * u.keV # Widths of bins
self.spectrum = numpy.array([0.], dtype=float) * u.ph / (u.m**2 * u.s * u.bin) # Spectrum of bins (in ph/s/m**2/bin at observatory)
self.luminosity = numpy.array([0.], dtype=float) * u.ph * 1E+44 / (u.s * u.keV) # Spectrum of bins (in 10^44 ph/s/keV at source distance)
self.table = QTable()
def allocate(self):
self.energy = numpy.zeros(self.nbins, dtype=float) * u.keV
self.energy_upper = numpy.zeros(self.nbins, dtype=float) * u.keV
self.energy_width = numpy.zeros(self.nbins, dtype=float) * u.keV
self.spectrum = numpy.zeros(self.nbins, dtype=float) * u.ph / (u.m**2 * u.s * u.bin)
self.luminosity = numpy.zeros(self.nbins, dtype=float) * u.ph * 1E+44 / (u.s * u.keV)
[docs] def get(self, isect):
"""Get the model spectra from SPEX for sector number isect."""
pyspex_f2py.api_model_spectrum.api_model_spectrum_total(isect)
self.nbins = pyspex_f2py.api_model_spectrum.nbins
self.allocate()
for i in numpy.arange(self.nbins):
self.energy[i] = pyspex_f2py.api_model_spectrum.energy_cent[i] * u.keV
self.energy_upper[i] = pyspex_f2py.api_model_spectrum.energy_bupp[i] * u.keV
self.energy_width[i] = pyspex_f2py.api_model_spectrum.energy_width[i] * u.keV
self.spectrum[i] = pyspex_f2py.api_model_spectrum.spectrum[i] * u.ph / (u.m**2 * u.s * u.bin)
self.luminosity[i] = pyspex_f2py.api_model_spectrum.luminosity[i] * u.ph * 1E+44 / (u.s * u.keV)
self.table = QTable([self.energy, self.energy_upper, self.energy_width, self.spectrum, self.luminosity],
names=('X_ctr', 'X_upp', 'X_wid', 'Spectrum', 'Luminosity'),
meta={'name': 'Sector Model spectrum'})
pyspex_f2py.api_model_spectrum.api_model_spectrum_final()
[docs]class Var:
"""Various settings for the plasma models.
:ivar gacc: Free-bound accuracy
:vartype gacc: float
:ivar line_ex: Electron excitation included
:vartype line_ex: bool
:ivar line_px: Proton excitation included
:vartype line_px: bool
:ivar line_rr: Radiative recombination included
:vartype line_rr: bool
:ivar line_dr: Di-electronic recombination included
:vartype line_dr: bool
:ivar line_ds: Di-electronic satellites included
:vartype line_ds: bool
:ivar line_ii: Inner shell ionisation included
:vartype line_ii: bool
:ivar doppler: Doppler broadening
:vartype doppler: int
:ivar newcalc: SPEXACT 3 calculations (False: SPEXACT 2)
:vartype newcalc: bool
:ivar occstart: Occupation calculations
:vartype occstart: int
:ivar mekal_wav: Wavelength corrections according to the work of Phillips et al. (1999)
:vartype mekal_wav: bool
:ivar mekal_fe17: The strongest Fe XVII lines by Doron & Behar (2002).
:vartype mekal_fe17: bool
:ivar mekal_update: Several minor corrections
:vartype mekal_update: book
:ivar ibalmaxw: Multi-Maxwellians for the ionisation balance
:vartype ibalmaxw: bool
:ivar newcoolexc: Cooling by collisional excitation by Stofanova (SPEXACT 3)
:vartype newcoolexc: bool
:ivar newcooldr: Cooling by dielectronic recombination (SPEXACT 3)
:vartype newcooldr: bool
"""
def __init__(self):
self.gacc = 0. # Free-bound accuracy
self.line_ex = True # Electron excitation included
self.line_px = True # Proton excitation included
self.line_rr = True # Radiative recombination included
self.line_dr = True # Dielectronic recombination included
self.line_ds = True # Dielectronic satellites included
self.line_ii = True # Inner shell ionisation included
self.doppler = 2 # Doppler broadening
# 1 = No broadening at all
# 2 = Only Doppler broadening (default)
# 3 = Only natural broadening (works only for var calc new)
# 4 = Doppler and natural broadening, Voigt profiles (best physical representation,
# works only for var calc new)
self.newcalc = 0 # SPEXACT 3 calculations (0 = SPEXACT 2; 1 = Quick CIE; 2 = SPEXACT 3)
self.occstart = 0 # Occupation calculations
# 0 = from ground state
# 1 = from Boltzmann distribution
# 2 = from last calculation
self.mekal_wav = True # Wavelength corrections according to the work of Phillips et al. (1999)
self.mekal_fe17 = True # The strongest Fe XVII lines by Doron & Behar (2002).
self.mekal_update =True # Several minor corrections
self.ibalmaxw = True # Multi-Maxwellians for the ionisation balance
self.newcoolexc = False # Cooling by collisional excitation by Stofanova (SPEXACT 3)
self.newcooldr = False # Cooling by dielectronic recombination (SPEXACT 3)
self.cxcon = 2 # Charge exchange rates (1 = Arnaud & Rothenflug, 2 = Kingdon & Ferland)
[docs] def update(self):
"""Obtain the current plasma model settings from SPEX."""
pyspex_f2py.api_model_var.api_var_update()
self.gacc = float(pyspex_f2py.api_model_var.var_gacc)
self.line_ex = bool(pyspex_f2py.api_model_var.var_line_ex)
self.line_px = bool(pyspex_f2py.api_model_var.var_line_px)
self.line_rr = bool(pyspex_f2py.api_model_var.var_line_rr)
self.line_dr = bool(pyspex_f2py.api_model_var.var_line_dr)
self.line_ds = bool(pyspex_f2py.api_model_var.var_line_ds)
self.line_ii = bool(pyspex_f2py.api_model_var.var_line_ii)
self.doppler = int(pyspex_f2py.api_model_var.var_line_doppler)
self.newcalc = int(pyspex_f2py.api_model_var.var_calc)
self.occstart = int(pyspex_f2py.api_model_var.var_occstart)
self.mekal_wav = bool(pyspex_f2py.api_model_var.var_mekal_wav)
self.mekal_fe17 = bool(pyspex_f2py.api_model_var.var_mekal_fe17)
self.mekal_update = bool(pyspex_f2py.api_model_var.var_mekal_update)
self.ibalmaxw = bool(pyspex_f2py.api_model_var.var_ibalmaxw)
self.newcoolexc = bool(pyspex_f2py.api_model_var.var_newcoolexc)
self.newcooldr = bool(pyspex_f2py.api_model_var.var_newcooldr)
self.cxcon = int(pyspex_f2py.api_model_var.var_cxcon)
[docs] def set_gacc(self, value):
"""Set the free-bound emission accuracy.
:param value: Free-bound emission accuracy.
:type value: float
"""
com = 'var gacc {:e}'.format(value)
ier = pyspex_f2py.api_command(com)
return ier
[docs] def reset_gacc(self):
"""Reset the free-bound emission accuracy."""
ier = pyspex_f2py.api_command('var gacc reset')
self.update()
return ier
[docs] def set_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
"""
if ltype == 'reset':
com = 'var line reset'
else:
if status:
com = 'var line {:s} true'.format(ltype)
else:
com = 'var line {:s} false'.format(ltype)
ier = pyspex_f2py.api_command(com)
self.update()
[docs] def set_doppler(self, value):
"""Doppler broadening.
:param value: Doppler broadening type.
:type value: int
"""
com = "var doppler {:d}".format(value)
ier = pyspex_f2py.api_command(com)
self.update()
return ier
[docs] def set_calc(self, status):
"""Perform SPEXACT 3 line calculations.
:param status: Use SPEXACT 2 (0), Quick SPEXACT 3 (1), or SPEXACT 3 (1)
:type status: int
"""
if (status==0):
ier = pyspex_f2py.api_command('var calc old')
elif (status==1):
ier = pyspex_f2py.api_command('var calc qc')
elif (status==2):
ier = pyspex_f2py.api_command('var calc new')
else:
print("Option "+str(status)+" is unknown.")
self.update()
return ier
[docs] def set_occstart(self, otype):
"""At which occupation level to start.
:param otype: Occupation level ('ground', 'boltz', 'last')
:type otype: str
"""
com = 'var occstart '+otype
ier = pyspex_f2py.api_command(com)
self.update()
return ier
[docs] def set_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
"""
if status:
com = 'var newmekal {:s} true'.format(utype)
else:
com = 'var newmekal {:s} false'.format(utype)
ier = pyspex_f2py.api_command(com)
self.update()
return ier
[docs] def set_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
"""
if status:
com = 'var ibalmaxw true'
else:
com = 'var ibalmaxw false'
ier = pyspex_f2py.api_command(com)
self.update()
return ier
[docs] def set_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
"""
if status:
com = 'var newcoolexc true'
else:
com = 'var newcoolexc true'
ier = pyspex_f2py.api_command(com)
self.update()
return ier
[docs] def set_newcooldr(self, status):
"""Cooling by dielectronic recombination (SPEXACT 3) on/off (True/False).
:param status: On (True) or off (False)
:type status: bool
"""
if status:
com = 'var newcooldr true'
else:
com = 'var newcooldr false'
ier = pyspex_f2py.api_command(com)
self.update()
return ier
[docs] def set_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
"""
ier = pyspex_f2py.api_command('var cxcon {0}'.format(value))
self.update()
return ier
[docs]class Dem:
"""SPEX DEM modeling interface.
:ivar nr: Number of temperature bins.
:vartype nr: int
:ivar tw: Temperature bins in keV.
:vartype tw: astropy.quantity.Quantity
:ivar yw: Differential emission measure as a function of temperature (in 10^64 m**3 / keV).
:vartype yw: astropy.quantity.Quantity
:ivar ywerr: Error on the differential emission measure as a function of temperature (in 10^64 m**3 / keV).
:vartype ywerr: astropy.quantity.Quantity
:ivar chisq: Chi^2 value for the DEM fit.
:vartype chisq: float
:ivar dempen: Dem penalty, number of bins that are less or equal to 0.
:vartype dempen: float
:ivar table: Astropy QTable containing DEM model.
:vartype table: astropy.table.QTable
"""
def __init__(self):
self.nr = 0
self.tw = numpy.array([0], dtype=float) * u.keV
self.yw = numpy.array([0], dtype=float) * 10**64 * u.m**3 / u.keV
self.ywerr = numpy.array([0], dtype=float) * 10**64 * u.m**3 / u.keV
self.mult_n = 0
self.mult_t = numpy.array([0], dtype=float) * u.keV
self.mult_dt = numpy.array([0], dtype=float) * u.keV
self.mult_y = numpy.array([0], dtype=float) * 10**64 * u.m**3
self.chi_n = 0
self.chi_reg = numpy.array([0], dtype=float) * u.dimensionless_unscaled
self.chi_chi = numpy.array([0], dtype=float) * u.dimensionless_unscaled
self.chi_pen = numpy.array([0], dtype=float) * u.dimensionless_unscaled
self.chisq = 0.
self.dempen = 0.
self.table = QTable()
self.mult_table = QTable()
self.chi_table = QTable()
def update(self):
pyspex_f2py.api_model_dem.api_model_dem_update()
self.nr = pyspex_f2py.api_model_dem.dem_nr
self.tw = numpy.zeros(self.nr, dtype=float) * u.keV
self.yw = numpy.zeros(self.nr, dtype=float) * u.m**3 / u.keV
self.ywerr = numpy.zeros(self.nr, dtype=float) * u.m**3 / u.keV
self.mult_n = pyspex_f2py.api_model_dem.dem_mul_n
self.mult_t = numpy.zeros(self.mult_n, dtype=float) * u.keV
self.mult_dt = numpy.zeros(self.mult_n, dtype=float) * u.keV
self.mult_y = numpy.zeros(self.mult_n, dtype=float) * u.m**3
self.chi_n = pyspex_f2py.api_model_dem.dem_nreg
self.chi_reg = numpy.zeros(self.chi_n, dtype=float) * u.dimensionless_unscaled
self.chi_chi = numpy.zeros(self.chi_n, dtype=float) * u.dimensionless_unscaled
self.chi_pen = numpy.zeros(self.chi_n, dtype=float) * u.dimensionless_unscaled
self.chisq = float(pyspex_f2py.api_model_dem.dem_chisq)
self.dempen = float(pyspex_f2py.api_model_dem.dem_dempen)
for i in numpy.arange(self.nr):
self.tw[i] = pyspex_f2py.api_model_dem.dem_tw[i] * u.keV
self.yw[i] = pyspex_f2py.api_model_dem.dem_yw[i] * 1E+64 * u.m**3 / u.keV
self.ywerr[i] = pyspex_f2py.api_model_dem.dem_ywerr[i] * 1E+64 * u.m**3 / u.keV
for i in numpy.arange(self.mult_n):
self.mult_t[i] = pyspex_f2py.api_model_dem.dem_mul_t[i] * u.keV
self.mult_dt[i] = pyspex_f2py.api_model_dem.dem_mul_dt[i] * u.keV
self.mult_y[i] = pyspex_f2py.api_model_dem.dem_mul_y[i] * 1E+64 * u.m**3
for i in numpy.arange(self.chi_n):
self.chi_reg[i] = pyspex_f2py.api_model_dem.dem_reg[i] * u.dimensionless_unscaled
self.chi_chi[i] = pyspex_f2py.api_model_dem.dem_chi[i] * u.dimensionless_unscaled
self.chi_pen[i] = pyspex_f2py.api_model_dem.dem_pen[i] * u.dimensionless_unscaled
self.table = QTable([self.tw, self.yw, self.ywerr], names=('kT', 'DY', 'DY_Err'),
meta={'name': 'DEM Model'})
self.mult_table = QTable([self.mult_t, self.mult_dt, self.mult_y],
names=('kT', 'dkT', 'DY'), meta={'name': 'DEM Mult Components'})
self.chi_table = QTable([self.chi_reg, self.chi_chi, self.chi_pen],
names=('Reg', 'Chi2', 'Penalty'), meta={'name': 'DEM Chireg Results'})
