The HEASARC welcomes your participation in a brief survey to capture how users access and utilize HEASARC data, software, and services. The outcome(s) of this survey will be used to guide, prioritize, and plan our activities and development in the coming years. It contains 20 questions, generally takes just a few minutes to complete, and your answers will remain totally anonymous. The survey is open until Dec 18, 2023. We thank you in advance for your valuable feedback.

Search:

[Advanced Search]



XSPEC Home











XSPEC Internal Functions Guide



For XSPEC Version 12.13.1


Keith A. Arnaud


HEASARC
Code 662
Goddard Space Flight Center
Greenbelt, MD 20771
USA


August 2023

Contents

1  Introduction
2  General Utility Routines
    2.1  Numerics
        2.1.1  AdaptiveIntegrate
        2.1.2  AstroFunctions
        2.1.3  Special functions
        2.1.4  Random number generators
        2.1.5  BinarySearch
        2.1.6  Convolution
        2.1.7  CosmologyFunction
        2.1.8  EigensystemWrapper
        2.1.9  Histogram
        2.1.10  Integrate
        2.1.11  LinearInterp
        2.1.12  MathOperator
        2.1.13  ModularCounter
        2.1.14  SVDwrapper
        2.1.15  Useful constants
3  Model Function Utility Routines
    3.1  FunctionUtility
    3.2  xsFortran and xsCFortran
4  Model Function C++ Classes
    4.1  Aped
        4.1.1  ApedIonRecord class
        4.1.2  ApedElementRecord class
        4.1.3  ApedTemperatureRecord
        4.1.4  Aped class
        4.1.5  Examples
    4.2  IonBalNei
        4.2.1  IonBalElementRecord
        4.2.2  IonBalTemperatureRecord
        4.2.3  IonBalNei
        4.2.4  Examples
    4.3  NeutralOpacity
        4.3.1  Examples
    4.4  IonizedOpacity
        4.4.1  Examples
    4.5  MZCompRefl
        4.5.1  Examples

1 Introduction

This document describes XSPEC internal functions and C++ classes available for writers of models.
These fall under three broad classes
  • General utility routines (these are under XSUtil).
  • Utility routines for use in model functions (these are in XSFunctions/Utilities).
  • C++ classes used for some of the major XSPEC models (these are in XSFunctions).
While these are not controlled as rigorously as XSPEC external interfaces, we do try to keep them as stable as possible so we encourage their use in XSPEC models.

2 General Utility Routines

This chapter describes the C++ classes and functions available in the XSPEC utility library (XSUtil). The different sub-directories of the utility library are
  • Error - the error handler.
  • Numerics - a variety of numerical methods.
  • Parse - routines used to parse various expressions used in XSPEC.
  • Signals - the interrupt handler.
  • Utils - miscellaneous other classes and routines.
The Numerics sub-directory is likely to be most useful for anyone writing a model.

2.1  Numerics

The Numerics directory comprises general numerical routines which are used within XSPEC and may also be helpful for anyone who wants to write a model. These methods are all written in C++ although a limited number have C/Fortran wrappers in xsFortran (see below). All methods are in the Numerics namespace.

2.1.1  AdaptiveIntegrate

AdaptiveIntegrate uses the Gauss-Kronrod algorithm to integrate a function to a specified precision. For efficiency reasons the function is passed into the integrator as a template parameter. 
template<RealArray Integrand(const RealArray& x, void *p)> 
    int AdaptiveIntegrate(const Real LowerLimit, const Real UpperLimit, 
                          void *p, const Real Precision, Real& Integral, 
                          Real& IntegralError);
The integrand function should be written to take an input array of Real values on which the function is to be evaluated and return a Real array of values. Any parameters required can be passed through the pointer p. An example use of AdaptiveIntegrate can be found in gaussianAbsorptionLine.cxx in the XSFunctions directory.

2.1.2  AstroFunctions

This namespace is intended for astronomical functions and at present just comprises a sexagesimal converter and a test for whether a number is an integer. 
bool sexagesimalToDecimal(const string& input,  Real& decDegrees) 

bool isIntValue(Real inVal)

2.1.3  Special functions

The Numerics namespace includes routines to evaluate a number of special functions
betaI Incomplete beta function (in Beta.h)
E1::operator() Exponential integral (in ExpInt.h)
Faddeeva exp(-z2) erfc(-iz) (in Faddeeva.h)
GammaLN ln of the Gamma function (in Gamma.h)
GammaP Incomplete Gamma function (in Gamma.h)
GammaQ Complement of GammaP (in Gamma.h)
Erf Error function (in Gamma.h)
Erfc Complement of Erf (in Gamma.h)
IncGamma::operator() Incomplete Gamma function (in IncGamma.h)
LnBang::operator() ln n! (in LnBang.h)

2.1.4  Random number generators

XSPEC uses the Luxury Pseudorandom Number generator proposed by Marsaglia and Zaman then implemented by James and Luescher in Fortran 77. It was translated to C++ for XSPEC. On start-up XSPEC generates an initial seed based on clock time. This can then be used to make pseudorandom numbers drawn from the following distributions. All these routines take a RealArray& argument which contains the random numbers on output. For PoissonRand the array should contain the Poisson means on input. All the generators are defined in RandomGenerator.h.
GaussRand Random numbers on G(0,1)
PoissonRand Poisson random numbers based on the input means
UniformRand Random numbers on U[0,1]
CauchyRand Random numbers on Cauchy(0,1)

2.1.5  BinarySearch

The BinarySearch routines return the index to an array (x) of the element immediately before the input test value (y). They assume that the array is in ascending order. The version with multiple test values also assumes that the test values are in ascending order. If the y value is less than the minimum of the x array then -1 is returned while if the y value is greater than the maximum of the x array then -2 is returned. 
  int BinarySearch(const RealArray& x, const Real& y)
  IntegerVector BinarySearch(const RealArray& x, const RealArray& y)
If y is less than the first values of x then -1 is returned. If y is greater than the last value of x then -2 is returned.

2.1.6  Convolution

The routines Convolution and ConvolutionInLnSpace use the fftw Fast Fourier Transform codes to convolve with a constant linear space kernel or a constant log space kernel, respectively. 
  template<void KernelFunction(const RealArray& kernelEnergyArray,
     const RealArray& kernelParams, int spectrumNumber, RealArray&
     kernelFlux, RealArray& kernelFluxErr, const string& kernelInitString)>
    void Convolution(const RealArray& energyArray, const RealArray&
     kernelParams, const int kernelFiducialEnergyIndex, const int
     spectrumNumber, const string& kernelInitString, RealArray&
     fluxArray, RealArray& fluxErrArray)
  template<void KernelFunction(const RealArray& kernelEnergyArray,
     const RealArray& kernelParams, int spectrumNumber, RealArray&
     kernelFlux, RealArray& kernelFluxErr, const string& kernelInitString)>
    void ConvolutionInLnSpace(const RealArray& energyArray, const
     RealArray& kernelParams, const int kernelFiducialEnergyIndex, const
     int spectrumNumber, const string& kernelInitString, RealArray&
     fluxArray, RealArray& fluxErrArray)
An example of the use of ConvolutionInLnSpace can be found in rdblur.cxx in the XSFunctions directory.

2.1.7  CosmologyFunction

The routine FZSQ takes as input Real variables specifying the redshift, q0 and λ0. Multiplying the output by c/H0 gives the luminosity distance using the approximation in Pen (ApJS 1990, 120, 49).

2.1.8  EigensystemWrapper

This is a wrapper for the gsl routines to calculate eigenvalues and eigenvectors for a symmetric matrix. 
  int calcEigensystem(double* matrix, const int N, double* eigenvals,
                double* eigenvects)
where matrix is the NxN matrix to be factored and which is destroyed on output, eigenvals is the eigenvalues vector, which must be allocated as size N, and eigenvects are the eigenvectors stored as columns in an array which must be allocated as NxN.

2.1.9  Histogram

Classes used to generate histograms from MCMC chains.

2.1.10  Integrate

The integrationKernel routine integrates a model flux array over a specified energy range and returns the photon and energy flux. Note that the energy array defines the ranges for each bin on which the input flux is given and hence has one more element. 
pair<Real,Real> integrationKernel (const RealArray& energy,
               const RealArray& fluxArray, const Real& eMin, const Real& eMax);

2.1.11  LinearInterp

The routines included in the Rebin namespace can be used to map XSPEC model flux arrays onto different energy bins. The first step is always to run the findFirstBins and initializeBins routines to set up the mapping between the current and target energy arrays. 
bool findFirstBins(const RealArray& currentBins, const RealArray& targetBins,
                   const Real FUZZY, size_t& currentStart, size_t& targetStart)
void initializeBins(const RealArray& currentBins, const RealArray& targetBins,
                    const Real FUZZY, size_t& currentStart, size_t& targetStart,
                    IntegerVector& startBin, IntegerVector& endBin, 
                    RealArray& startWeight, RealArray& endWeight)
FUZZY is a fractional fuzziness for deciding whether bin boundaries are equal. The output arrays from initializeBins define the mapping between the current and target energies.
For additive models the rebin function should be used and for multiplicative models the interpolate function should be used. Examples can be found in zashift.cxx and zmshift.cxx in the XSFunctions directory. The lowValue and highValue arguments can be used to define output values below and above the energies available from the input. If left out they default to 0. 
void rebin(const RealArray& inputArray, const IntegerVector& startBin, 
           const IntegerVector& endBin, const RealArray& startWeight, 
           const RealArray& endWeight, RealArray& outputArray, 
           const Real lowValue=0.0, const Real highValue=0.0);
void interpolate(const RealArray& inputArray, const IntegerVector& startBin, 
                 const IntegerVector& endBin, const RealArray& startWeight, 
                 const RealArray& endWeight, RealArray& outputArray, 
                 bool exponential, const Real lowValue=0.0,
                 const Real highValue=0.0)
The gainRebin function is a special case used within the gain command and linInterpInteg assumes that the input is in photons/cm2/s/keV at specific energy points and is integrated over the target bin sizes. 
void gainRebin(const RealArray& inputArray, const IntegerVector &startBin,
               const IntegerVector& endBin, const RealArray& startWeight,
               const RealArray& endWeight, RealArray& outputArray)
void linInterpInteg(const RealArray& currentPoints,
                    const RealArray& inputdValues,
                    const RealArray& targetBins, RealArray& outputValues,
                    const Real lowValue=0.0, const Real highValue=0.0)

2.1.12  MathOperator

MathOperator is used to evaluate models defined using the mdefine command. The following classes are all defined in MathOperator.h
PlusOp add two variables
MinusOp subtract two variables
MultOp multiply two variables
DivideOp divide two variables
PowOp take power of variable
MaxOp maximum of two variables
MinOp minimum of two variables
Atan2Op arctan using two variables
UnaryMinusOp multiply variable by minus one
ExpOp exponential of variable
SinOp sine of variable
SinDOp sine of variable in degrees
CosOp cosine of variable
CosDOp cosine of variable in degrees
TanOp tan of variable
TanDOp tan of variable in degrees
SinhOp sinh of variable
SinhDOp sinh of variable in degrees
CoshOp cosh of variable
CoshDOp cosh of variable in degrees
TanhOp tanh of variable
TanhDOp tanh of variable in degrees
LogOp log base 10 of variable
LnOp natural log of variable
SqrtOp sqrt of variable
AbsOp absolute value of variable
IntOp integer part of variable
SignOp -1 if negative, +1 if positive
HOp 0 if negative, +1 if positive
BoxcarOp +1 between 0 and 1, 0 otherwise
ASinOp inverse sine of variable
ACosOp inverse cosine of variable
ATanOp inverse tan of variable
ASinhOp inverse sinh of variable
ACoshOp inverse cosh of variable
ATanhOp inverse tanh of variable
ErfOp error function
ErfcOp complementary error function
GammaOp gamma function
Legendre2Op 2nd-order Legendre polynomial
Legendre3Op 3rd-order Legendre polynomial
Legendre4Op 4th-order Legendre polynomial
Legendre5Op 5th-order Legendre polynomial
MeanOp mean of a vector
DimOp length of a vector
SMinOp minimum value of a vector
SMaxOp maximum value of a vector

2.1.13  ModularCounter

The ModularCounter class is used to support the step and margin commands.

2.1.14  SVDwrapper

This is a wrapper for the gsl singular value decomposition code. 
int callSVD_square(double* matrix, double *eigenvals, double *eigenvects, 
                   const int N)
where matrix in input is the square matrix of size NxN to be factored into USVdT and on output is the matrix U. eigenvals is an array allocated to size N which on output contain S, the variances on the principal axes. eigenvects is an array allocated as NxN which on output contains V, the principal axes (not transposed).

2.1.15  Useful constants

The Numerics namespace defines a number of useful constants defined as static Real in Numerics.h.
KEVTOA 12.39841974 from CODATA 2014
KEVTOHZ 2.4179884076620228e17
KEVTOERG 1.60217733e-9
KEVTOJY 1.60217733e14
DEGTORAD 0.01745329252
LIGHTSPEED299792.458 defined in km/s
AMU 1.660539040e-24 unified atomic mass unit in g

3 Model Function Utility Routines

ComponentInfo, XSCCall, XSF77, and XSModelFunction are part of the internal XSPEC model interface and are not documented here since they should not be of general interest. FunctionUtility is the C++ class containing methods likely to be of use to people writing models or using the XSPEC model library. xsFortran contains C/Fortran wrappers for many FunctionUtility methods.

3.1  FunctionUtility

The FunctionUtility class stores information used by XSPEC models. There is a single instantiation. The private elements of FunctionUtility are all static and are as follows:
string s_XSECT Photoelectric cross-sections used
string s_ABUND Relative abundances used
const string CROSSSECTFILE File read for cross-sections
string s_abundanceFile File read for relative abundances
vector < string > s_elements Names of elements
string s_managerPath Directory path for model.dat
const size_t s_NELEMS Number of elements
string s_modelDataPath Directory path for files required by models
string s_NOT_A_KEY String returned if an xset variable is not set
FunctionUtility::Cosmology s_COSMO Cosmology used
string s_abundPath Directory path for the abundance file
string s_atomdbVersion version of AtomDB used
string s_neiVersion version of NEI code used
bool s_abundChanged Relative abundance choice has been changed
int s_xwriteChatter Chatter value
vector < double > s_tempsDEM Last CEI model temperatures calculated
vector < double > s_DEM Last CEI model DEMs calculated
map < int, map < string, Real > > s_XFLT XLT#### keyword values
map < string,double > s_valueDataBase Database of values calculated in model functions
map < string,vector < float > > s_abundanceVectors Stores sets of abundances
map < string,string > s_crossSections descriptions of cross-sections available
map < string,string > s_abundDoc descriptions of relative abundance sets
abundance sets
map < string,string > s_modelStringDataBase Keywords and values defined by the xset command
These elements can be accessed using either C++ or C/Fortran routines as follows. The following tables list the method name, the C/Fortran wrapper (where available) and a brief description. To see the calling sequences look in FunctionUtility.h, xsFortran.cxx, and xsCFortran.c for the C++, C and Fortran interfaces, respectively.
Any program which uses the XSPEC function library needs to call FNINIT() to initialize the FunctionUtility object. The managerPath is the directory which contains the model.dat file as well as information about the cross-sections and abundance tables available. Usually, it does not have to be changed from the default (Xspec/src/manager). The modelDataPath is the directory containing files used as input by models.
FNINIT FNINIT standard initialization
managerPath FGDATD get s_managerPath
managerPath FPDATD set s_managerPath
modelDataPath FGMODF get s_modelDataPath
modelDataPath set s_modelDataPath
NOT_A_KEY get s_NOT_A_KEY
XSPEC has several options for photo-electric cross-sections and relative abundances. A user-defined set of relative abundances can also be read from a file.
XSECT FGXSCT get s_XSECT
XSECT FPXSCT set s_XSECT
checkXsect check for valid cross-section table
crossSections get description of cross-section table
ABUND FGSOLR get s_ABUND
ABUND FPSOLR set s_ABUND
getAbundance FGABND get relative abundance for an element
getAbundance FGABNZ get relative abundance for an atomic number
getAbundance FGTABN get relative abundance for an element from a table
getAbundance FGTABZ get relative abundance for an atomic number from a table
readNewAbundances RFLABD read relative abundances from a file
checkAbund check for valid relative abundance table
abundanceVectors FPSLFL set the values of the file relative abundance table
abundDoc get description of relative abundance table
abundanceFile FGABFL get the name of the file used for abundances
abundanceFile FPABFL set the name of the file used for abundances
abundPath FGAPTH get the directory path for the file used for abundances
abundPath FPAPTH set the directory path for the file used for abundances
abundChanged get value of s_abundChanged flag
abundChanged set value of s_abundChanged flag
readInitializers set up cross-sections and relative abundances
elements FGELTI get name of element
NELEMS() FGNELT get number of elements
The AtomDB and NEI versions used can be set or retrieved
atomdbVersion FGATDV get AtomDB version in use
atomdbVersion FPATDV set AtomDB version in use
neiVersion get NEI version in use
neiVersion set NEI version in use
The (keyword, value) pairs defined using the xset command can be found using getModelString.
getModelString FGMSTR get entry from the model string database
setModelString FPMSTR set entry in the model string database
modelStringDataBase get database as a map<string,string>
eraseModelStringDataBase clear the model string database
XSPEC uses a cosmological model based on values of H0, q0, and λ0 to calculate distances from redshifts. A routine to calculate the luminosity distance is included in Numerics with a C/Fortran wrapper in xsFortran.
setFunctionCosmoParams CSMPALL set the cosmology H0, q0, and λ0
getq0 csmgq0 get cosmology q0
setq0 csmpq0 set cosmology q0
getH0 csmgh0 get cosmology H0
setH0 csmph0 set cosmology H0
getlambda0 csmgl0 get cosmology λ0
setlambda0 csmpl0 set cosmology λ0
XSPEC models have access to the values given by the keywords XFLT#### read from the SPECTRUM extension of the input file. These are stored internally as (key, value) pairs where the XFLT#### keywords are strings "key: value".
getNumberXFLT DGNFLT get number of XFLT#### keywords for spectrum
getXFLT DGFILT get value of XFLT#### keyword for spectrum
inXFLT DGQFLT check whether XFLT#### keyword exists for spectrum
loadXFLT DPFILT set XFLT value(s) for spectrum
clearXFLT DCLFLT clear all XFLT values
getAllXFLT get all XFLT values for spectrum
The calcMultiTempPlasma routine which underlies many of the collisional plasma models automatically saves its input temperatures and differential emission measures. This is used for the plot dem option but may have other uses.
tempsDEM get the vector of temperatures
tempsDEM set from a vector of temperatures
tempsDEM set from a RealArray of temperatures
DEM get the vector of DEMs
DEM set from a vector of DEMs
DEM set from a RealArray of DEMs
There is an internal database to which any model can add a (keyword,value) pair. This is useful for saving information during the calculation of the model which the user may need.
getDbValue GDBVAL get value for given keyword in internal database
loadDbValue PDBVAL set (keyword,value) pair in internal database
clearDb CDBASE clear all (keyword,value) pairs in internal database
getDbKeywords get all keywords in internal database as a string
getAllDbValues get all keywords in internal database as a map
The chattiness level set by the chatter command can be accessed from within in a model and strings can be output at the required level
xwriteChatter FGCHAT get current chatter level
xwriteChatter FPCHAT set current chatter level
xsWrite XWRITE write an output string
Usually table models are used in xspec through the atable or mtable model options but it also possible to write a model which reads a table and performs additional operations. There are two versions of tableInterpolate depending on whether additional xflt information. There is a C wrapper tabintxflt for the version which uses additional xflt information however it is not available from Fortran. There is also a helpful (C++ only) routine to return information about the table in a file whose name is input.
tableInterpolate TABINT get tabulated value
tableInfo return information about the table

3.2  xsFortran and xsCFortran

xsFortran and xsCFortran include the C and Fortran interfaces listed in the FunctionUtility section above as well as the following routines. xs_write should be used for all output.
xs_getChat XTGTCHT get terminal and log chatter levels
xs_getVersion XGVERS get XSPEC version
xs_write XWRITE writes to terminal and/or log file
xs_read XREAD interactive read from terminal
There are also the following wrappers to routines in the Numerics namespace.
xs_erf ERF wrapper for Numerics::Erf to get Error Function value
xs_erfc ERFC wrapper for Numerics::Erfc to get complementary Error Function value
gammap GAMMAP wrapper for Numerics::GammaP
gammaq GAMMQ wrapper for Numerics::GammaQ
fzsq FZSQ wrapper for Numerics::FZSQ to get luminosity distance is (c/H0)fzsq
findFirstBins FFBINS wrapper for Numerics::Rebin::findFirstBins
dfindFirstBins DFFBINS wrapper for Numerics::Rebin::findFirstBins
initBins INIBINS wrapper for Numerics::Rebin::initializeBins
dinitBins DINIBINS wrapper for Numerics::Rebin::initializeBins
rebinBins RBNBINS wrapper for Numerics::Rebin::rebin
drebinBins DRBNBINS wrapperm for Numerics::Rebin::rebin
interpBins INTBINS wrapper for Numerics::Rebin::interpolate
dinterpBins DINTBINS wrapper for Numerics::Rebin::interpolate
gainRebin GNREBIN wrapper for Numerics::Rebin::gainRebin
dgainRebin DGNREBIN wrapper for Numerics::Rebin::gainRebin
linInterpInteg LININTINTEG wrapper for Numerics::Rebin::linInterpInteg
dlinInterpInteg DLININTINTEG wrapper for Numerics::Rebin::linInterpInteg
The following handy functions return constants stored in Numerics.h
getkeVtoA KEVTOA The keV to Angstrom conversion factor
getkeVtoHz KEVTOHZ The keV to Hz conversion factor
getkeVtoErg KEVTOERG The keV to erg conversion factor
getkeVtoJy KEVTOJY The keV to Jy conversion factor
getdegtorad DEGTORAD The degrees to radians conversion factor
getlightspeed LIGHTSPEED The speed of light in km/s
getamu AMU The unified atomic mass unit in gm

Chapter 4
Model Function C++ Classes

The classes are
  • Aped - generates CEI and NEI spectra using AtomDB input files.
  • IonBalNei - calculates NEI ionization balances.
  • NeutralOpacity - calculates opacities for neutral material.
  • IonizedOpacity - calculates opacities for photo-ionized material.
  • MZCompRefl - calculates Compton reflection using the Magzdiarz and Zdziarski code.

4.1  Aped

The Aped classes store and use the data stored in the AtomDB files. The top-level class, Aped, comprises a vector of ApedTemperatureRecord objects, each of which comprises a vector of ApedElementRecord objects, each of which comprises a vector of ApedIonRecord objects. Thus, there is one ApedIonRecord for each temperature, element, and ion in the AtomDB files. The top-level Aped class provides the methods to load and use the atomic data.

4.1.1  ApedIonRecord class

class ApedIonRecord{
 public:

  int m_Ion;
  RealArray m_ContinuumEnergy;
  RealArray m_ContinuumFlux;
  RealArray m_ContinuumFluxError;
  RealArray m_PseudoContinuumEnergy;
  RealArray m_PseudoContinuumFlux;
  RealArray m_PseudoContinuumFluxError;
  RealArray m_LineEnergy;
  RealArray m_LineEnergyError;
  RealArray m_LineEmissivity;
  RealArray m_LineEmissivityError;
  IntegerVector m_ElementDriver;
  IntegerVector m_IonDriver;
  IntegerVector m_UpperLevel;
  IntegerVector m_LowerLevel;

  ApedIonRecord();    // default constructor
  ~ApedIonRecord();   // destructor

  ApedIonRecord& operator=(const ApedIonRecord&);    // deep copy

};
The class includes entries to store error estimates on the fluxes and emissivities although these are not used at the time of writing. Note that the same classes are used for NEI and CEI although for the latter the ElementDriver and IonDriver arrays are not required.

4.1.2  ApedElementRecord class

class ApedElementRecord{
 public:

  int m_AtomicNumber;
  vector<ApedIonRecord> m_IonRecord;

  ApedElementRecord();    // default constructor
  ~ApedElementRecord();   // destructor

  void LoadIonRecord(ApedIonRecord input);

  ApedElementRecord& operator=(const ApedElementRecord&);    // deep copy

};
Each ApedElementRecord object stores the atomic number of the element and the records for each of its ions.

4.1.3  ApedTemperatureRecord

class ApedTemperatureRecord{
 public:

  Real m_Temperature;
  vector<ApedElementRecord> m_ElementRecord;

  ApedTemperatureRecord();     //default constructor
  ~ApedTemperatureRecord();    //destructor

  void LoadElementRecord(ApedElementRecord input);

  ApedTemperatureRecord& operator=(const ApedTemperatureRecord&);    // deep copy

};
Each ApedTemperatureRecord object stores the temperature and the records for each element.

4.1.4  Aped class

class Aped{
 public:

  vector<ApedTemperatureRecord> m_TemperatureRecord;

  // These store the information in the initial PARAMETERS extension

  RealArray m_Temperatures;
  IntegerVector m_NelementsLine;
  IntegerVector m_NelementsCoco;
  IntegerVector m_Nline;
  IntegerVector m_Ncont;

  string m_coconame;
  string m_linename;

  bool m_noLines;
  bool m_noResonanceLines;
  bool m_thermalBroadening;
  Real m_velocityBroadening;
  Real m_minimumLineFluxForBroadening;
  bool m_broadenPseudoContinuum;
  bool m_logTempInterpolation;
  bool m_multiThread;

  Aped();     //default constructor
  ~Aped();    //destructor

  Aped& operator=(const Aped&);    // deep copy
The class provides methods to load the data from the AtomDB continuum and line files. The Read method just gets the information from the PARAMETERS extension, stores the filenames used and sizes the subsidiary objects correctly. To actually load the data use ReadTemperature. 
  // Reads the continuum and line files and stores the names and
  // temperatures. Does not read the actual continuum and line data.

  int Read(string cocofilename, string linefilename);

  int ReadTemperature(const int TemperatureIndex);
  int ReadTemperature(const vector<int>& TemperatureIndex);
The following methods set the switches which determine which options are used during the calculation. Each of these methods also checks for the corresponding xset variables: APECNOLINES, APECNORES, APECTHERMAL, APECVELOCITY, APECMINFLUX, APECBROADPSEUDO, APECLOGINTERP, APECMULTITHREAD respectively. For SetVelocityBroadening a warning will be issued if the input value is non-zero and differs from APECVELOCITY. If m_noLines is set then the spectrum will only be continuum. At the moment m_noResonanceLines has no effect but is there for a possible approach to resonance scattering. If m_thermalBroadening is set then any lines will be thermally broadened while if m_velocityBroadening is set then any lines will be broadened by the appropriate velocity. The line shape is assumed Gaussian and thermal and velocity broadening are added in quadrature. If m_minimumLinefluxForBroadening is set then all lines below that flux will not be broadened. If m_broadenPseudoContinuum is set then the pseudo-continuum will also be broadened. Note that line broadened spectra take longer to calculate. If m_multiThread is set then the calculation of the line spectrum is multi-threaded over temperatures. 
  void SetNoLines(const bool qno);
  void SetNoResonanceLines(const bool qnores);
  void SetThermalBroadening(const bool qtherm);
  void SetVelocityBroadening(const Real velocity);
  void SetMinimumLinefluxForBroadening(const Real flux);
  void SetBroadenPseudoContinuum(const bool qpbroad);
  void SetLogTempInterpolation(const bool qltemp);
  void SetMultiThread(const bool qmulti);
  void SetLineSpecs(const bool qno, const bool qnores, const bool qtherm,
                    const Real velocity, const Real flux,
                    const bool qpbroad, const bool qltinterp,
                    const bool qmulti);
The following provide methods to extract some useful numbers. 
  int NumberTemperatures();     // return number of tabulated temperatures
  int NumberElements();         // return number of tabulated elements
  int NumberIons(int Z);        // return number of ions across all temperatures
                                // for element Z.
Whether the data for a given temperature index has been loaded can be checked using IsTemperatureLoaded. 
 bool IsTemperatureLoaded(const int TemperatureIndex);
The method to generate spectra for CEI plasmas is SumEqSpectra. This comes in several overloaded options depending on whether a single temperature of distribution of temperatures are required. It is also possible to use a different temperature for thermal broadening of lines than that used for the ionization balance. The energy bins on which the spectrum is to be calculated are specified by standard XSPEC energyArray. The elements required and their abundances are given by Zinput and abundance. Tinput and Dem give the temperature(s) and emission measure(s) required. 
  void SumEqSpectra(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const Real& Dem, RealArray& fluxArray, 
                    RealArray& fluxErrArray);

  void SumEqSpectra(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const RealArray& Dem, RealArray& fluxArray, 
                    RealArray& fluxErrArray);

  // case where the temperature used for the thermal broadening differs from that
  // for the ionization

  void SumEqSpectra(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput, 
                    const Real& Tbinput, const Real& Dem,
                    RealArray& fluxArray, RealArray& fluxErrArray);

  void SumEqSpectra(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput, 
                    const RealArray& Tbinput, const RealArray& Dem, 
                    RealArray& fluxArray, RealArray& fluxErrArray);
For NEI plasmas, the method to calculate the spectrum is SumNeqSpectra. The inputs differ from SumEqSpectra only in including IonFrac, which specifies the ionization fractions required for each element, for each temperature. 
  void SumNeqSpectra(const RealArray& energyArray, 
                     const IntegerVector& Zinput, const RealArray& abundance,
                     const Real Redshift, const Real& Tinput,
                     const vector<RealArray>& IonFrac, 
                     RealArray& fluxArray, RealArray& fluxErrArray);

  void SumNeqSpectra(const RealArray& energyArray, 
                     const IntegerVector& Zinput, const RealArray& abundance,
                     const Real Redshift, const RealArray& Tinput,
                     const vector<vector<RealArray> >& IonFrac, 
                     RealArray& fluxArray, RealArray& fluxErrArray);

  // case where the temperature used for the thermal broadening differs from that
  // for the ionization

  void SumNeqSpectra(const RealArray& energyArray, 
                     const IntegerVector& Zinput, const RealArray& abundance,
                     const Real Redshift, const Real& Tinput,
                     const Real& Tbinput, 
                     const vector<RealArray>& IonFrac, 
                     RealArray& fluxArray, RealArray& fluxErrArray);

  void SumNeqSpectra(const RealArray& energyArray, 
                     const IntegerVector& Zinput, const RealArray& abundance,
                     const Real Redshift, const RealArray& Tinput,
                     const RealArray& Tbinput, 
                     const vector<vector<RealArray> >& IonFrac, 
                     RealArray& fluxArray, RealArray& fluxErrArray);
All versions of SumEqSpectra and SumNeqSpectra operate through 
  void SumSpectra(const RealArray& energyArray, 
                  const IntegerVector& Zinput, const RealArray& abundance,
                  const Real Redshift, const RealArray& Tinput,
                  const RealArray& Tbinput, 
                  const vector<vector<RealArray> >& IonFrac, const bool isCEI,
                  RealArray& fluxArray, RealArray& fluxErrArray);
Underneath this there are routines to evaluate the continuum and line spectra for a given temperature 
  int calcContinuumSpectrumForTemperature(const size_t TRecordIndex,
               const IntegerVector& Zinput, const RealArray& abunZ,
               const vector<RealArray>& IonFrac, 
               const RealArray& sourceFrameEnergy, const bool isCEI,
               RealArray& fluxArray);

  int calcLineSpectrumForTemperature(const size_t TRecordIndex,
               const IntegerVector& Zinput, const RealArray& abunZ,
               const vector<RealArray>& IonFrac,
               const RealArray& sourceFrameEnergy, const bool isCEI,
               const Real minLineflux, RealArray& fluxArray,
               Real& maxLineflux);
Outside the class there are wrap-up routines which create an internal, static Aped object then call SumSpectra to calculate the output spectrum. There are overloaded versions corresponding to the options for SumEqSpectra and SumNeqSpectra. The int returned will be non-zero in the event of an error on reading the AtomDB files. For the CEI case these routines check the APECROOT variable to find the files to read. For the NEI case they check the NEIVERS and NEIAPECROOT variables. The calcRSSpectrum routines are for the old Raymond-Smith model and read the RS files in the AtomDB format. 
int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const Real& Dem, const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const RealArray& Dem, const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const Real& Tbinput, const Real& Dem, 
                    const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const RealArray& Tbinput, const RealArray& Dem, 
                    const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const Real& Dem, const bool qtherm, const Real velocity,
                    const bool noLines,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const RealArray& Dem, const bool qtherm, const Real velocity,
                    const bool noLines,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const Real& Tbinput, const Real& Dem, 
                    const bool qtherm, const Real velocity,
                    const bool noLines,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const RealArray& Tbinput, const RealArray& Dem, 
                    const bool qtherm, const Real velocity,
                    const bool noLines,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcRSSpectrum(const RealArray& energyArray, 
                   const IntegerVector& Zinput, const RealArray& abundance,
                   const Real Redshift, const Real& Tinput,
                   const Real& Dem, const bool qtherm, const Real velocity,
                   RealArray& fluxArray, RealArray& fluxErrArray);

int calcRSSpectrum(const RealArray& energyArray, 
                   const IntegerVector& Zinput, const RealArray& abundance,
                   const Real Redshift, const RealArray& Tinput,
                   const RealArray& Dem, const bool qtherm, const Real velocity,
                   RealArray& fluxArray, RealArray& fluxErrArray);

int calcRSSpectrum(const RealArray& energyArray, 
                   const IntegerVector& Zinput, const RealArray& abundance,
                   const Real Redshift, const Real& Tinput,
                   const Real& Tbinput, const Real& Dem, 
                   const bool qtherm, const Real velocity,
                   RealArray& fluxArray, RealArray& fluxErrArray);

int calcRSSpectrum(const RealArray& energyArray, 
                   const IntegerVector& Zinput, const RealArray& abundance,
                   const Real Redshift, const RealArray& Tinput,
                   const RealArray& Tbinput, const RealArray& Dem, 
                   const bool qtherm, const Real velocity,
                   RealArray& fluxArray, RealArray& fluxErrArray);
int calcNEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const vector<RealArray>& IonFrac, 
                    const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcNEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const vector<vector<RealArray> >& IonFrac, 
                    const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcNEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const Real& Tinput,
                    const Real& Tbinput,
                    const vector<RealArray>& IonFrac, 
                    const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcNEISpectrum(const RealArray& energyArray, 
                    const IntegerVector& Zinput, const RealArray& abundance,
                    const Real Redshift, const RealArray& Tinput,
                    const RealArray& Tbinput,
                    const vector<vector<RealArray> >& IonFrac, 
                    const bool qtherm, const Real velocity,
                    RealArray& fluxArray, RealArray& fluxErrArray);
A few other useful routines also live in Aped.h. getAtomicMass returns the atomic mass for the given atomic number. getApedFileNames returns the AtomDB version number string, and the continuum and line filenames. getNEIApedFileNames returns the AtomDB version number string, the NEI eigenvectors version number string, and the continuum and line filenames. apedInterpFlux is used to interpolate the continuum and pseudo-continuum arrays. 
Real getAtomicMass(const int& AtomicNumber);
bool getApedFileNames(string& version, string& continuumFile, string& lineFile);
bool getNEIApedFileNames(string& version, string& eigenVersion,
                         string& continuumFile, string& lineFile);
void apedInterpFlux(const RealArray& inputEnergy, const RealArray& inputFlux, 
                    const Real& z15, const Real& coeff,
                    const RealArray& energyArray, RealArray& fluxArray);

4.1.5  Examples

An example use of the Aped class to calculate a CEI spectrum can be found in the file calcMultiTempPlasma.cxx in Xspec/src/XSFunctions. An example use to calculate an NEI spectrum can be found in vvgnei.cxx in the same directory.

4.2  IonBalNei

The IonBalNei classes store the eigenvector data and calculate ionization fractions for an NEI collisional plasma. The top-level class IonBalNei comprises a vector of IonBalTemperatureRecord classes, each of which comprises a vector of IonBalElementRecord classes. Thus there is one IonBalElementRecord for each temperature and element in the eigevector data files. The top-level IonBalNei class provides methods to load and use the eigenvector data.

4.2.1  IonBalElementRecord

class IonBalElementRecord{
 public:

  int AtomicNumber;
  RealArray EquilibriumPopulation;
  RealArray Eigenvalues;
  vector<RealArray> LeftEigenvectors;
  vector<RealArray> RightEigenvectors;
  
  IonBalElementRecord();     // default constructor
  ~IonBalElementRecord();    // destructor

  void Clear();  // clear out the contents

  IonBalElementRecord& operator=(const IonBalElementRecord&);  // deep copy

};
This class stores the eigenvector data for the element with specified AtomicNumber. EquilibriumPopulation stores the CEI ion fractions.

4.2.2  IonBalTemperatureRecord

class IonBalTemperatureRecord{
 public:

  Real Temperature;
  vector<IonBalElementRecord> ElementRecord;

  IonBalTemperatureRecord();        // default constructor
  ~IonBalTemperatureRecord();       // destructor

  void LoadElementRecord(IonBalElementRecord input);

  void Clear();   // clear out the contents

  IonBalTemperatureRecord& operator=(const IonBalTemperatureRecord&);

};
Each IonBalTemperatureRecord object stores the temperature and the eigenvector data for each element.

4.2.3  IonBalNei

class IonBalNei{
 public:

  vector<IonBalTemperatureRecord> TemperatureRecord;
  RealArray Temperature;
  string Version;

  IonBalNei();       // default constructor
  ~IonBalNei();      // destructor

  IonBalNei& operator=(const IonBalNei&);    // deep copy

  void Clear();  // clear out the contents
The Version string is used to control which version of the eigenvector data files are to be used. If the version is changed then the object is reset using Clear and setVersion returns true. The setVersion method with no input checks NEIVERS and uses that if it has been set. 
  string getVersion();
  bool setVersion();
  bool setVersion(string version);
This class provides a Read method to load data from the eigenvector files. This can either be done using a vector of atomic numbers and the directory name and version string, so that the relevant filenames can be constructed, or using a single atomic number and filename. 
  //  Read loads eigenvector data for a set of elements

  int Read(vector<int> Z, string dirname, string versionname);

  //  ReadElement loads eigenvector data for that element

  int ReadElement(int Z, string filename);

  void LoadTemperatureRecord(IonBalTemperatureRecord input);
The following provide methods to extract some useful numbers. 
  RealArray Temperatures();     // return tabulated temperatures
  int NumberTemperatures();     // return number of tabulated temperatures
  int NumberElements();         // return number of tabulated elements
  bool ContainsElement(const int& Z); // returns true if data for Z
The CEI method returns the ionization fractions for the specified element for a collisional ionization equilibrium plasma with the given electron temperature. 
  RealArray CEI(const Real& Te, const int& Z);
The method to calculate ion fractions is Calc. It is overloaded to give a number of different options. If the initIonFrac array is used then this gives the initial ion fractions for the element Z, if not the element is assumed to start un-ionized. 
  // calculate the NEI ion fractions for electron temperature Te, ionization
  // parameter tau and element Z.

  RealArray Calc(const Real& Te, const Real& tau, const int& Z);
  RealArray Calc(const Real& Te, const Real& tau, const int& Z, 
                 const RealArray& initIonFrac);

  // Calculates ionization fractions at electron temperatures
  // Te and a set of ionization parameters tau(i), i=1,..,n,
  // where each tau is given weight(i). Electron temperature is
  // assumed to be linear function of tau.
  // Based on the old noneq.f.

  RealArray Calc(const RealArray& Te, const RealArray& tau, 
                 const RealArray& weight, const int& Z);
  RealArray Calc(const RealArray& Te, const RealArray& tau, 
                 const RealArray& weight, const int& Z, 
                 const RealArray& initIonFrac);

  // Calculates ionization fractions at electron temperature
  // Te(n) and ionization parameter tau(n), for electron 
  // temperatures Te given in a tabular form as a function of 
  // ionization parameter tau.
  // Based on the old noneqr.f.
  // Does initIonFrac make sense in this case.

  RealArray Calc(const RealArray& Te, const RealArray& tau, 
                 const int& Z);
  RealArray Calc(const RealArray& Te, const RealArray& tau, 
                 const int& Z, const RealArray& initIonFrac);

};
Outside the class there are wrap-up functions to read (if necessary) the eigenvector files and calculate ion fractions. They use the NEISPEC xset variable to choose which version of the files to use. The various overloaded versions match to the Calc methods. 
void calcNEIfractions(const Real& Te, const Real& tau, const int& Z, 
                      RealArray& IonFrac);
void calcNEIfractions(const Real& Te, const Real& tau, const IntegerVector& Z, 
                      vector<RealArray>& IonFrac);

void calcNEIfractions(const Real& Te, const Real& tau, const int& Z, 
                      const RealArray& initIonFrac, RealArray& IonFrac);
void calcNEIfractions(const Real& Te, const Real& tau, const IntegerVector& Z, 
                      const vector<RealArray>& initIonFrac, 
                      vector<RealArray>& IonFrac);

void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const RealArray& weight, const int& Z, RealArray& IonFrac);
void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const RealArray& weight, const IntegerVector& Z, 
                      vector<RealArray>& IonFrac);

void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const RealArray& weight, const int& Z, 
                      const RealArray& initIonFrac, RealArray& IonFrac);
void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const RealArray& weight, const IntegerVector& Z, 
                      const vector<RealArray>& initIonFrac, 
                      vector<RealArray>& IonFrac);

void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const int& Z, RealArray& IonFrac);
void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const IntegerVector& Z, vector<RealArray>& IonFrac);

void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const int& Z, const RealArray& initIonFrac, 
                      RealArray& IonFrac);
void calcNEIfractions(const RealArray& Te, const RealArray& tau, 
                      const IntegerVector& Z, const RealArray& initIonFrac, 
                      vector<RealArray>& IonFrac);
The wrapper routine to return the collisional ionization equilibrium fractions is 
void calcCEIfractions(const Real Te, const IntegerVector& Z, 
                      vector<RealArray>& IonFrac);
There are also handy routines to return and index into the temperatures and the number of temperatures 
int getNEItempIndex(const Real& tkeV);
int getNEInumbTemp();
to providing debugging information 
string writeIonFrac(const IntegerVector& Zarray, 
                    const vector<RealArray>& IonFrac);
string writeIonFrac(const int& Z, const IntegerVector& Zarray, 
                    const vector<RealArray>& IonFrac);
to do a binary search on a RealArray and return the index of the element immediately less than the input target 
int locateIndex(const RealArray& xx, const Real x);
and to check whether arrays are identical (note that the C++11 standard includes this as part of the valarray class but this is not yet implemented in all compilers. 
bool identicalArrays(const vector<RealArray>& a, const vector<RealArray>& b);
bool identicalArrays(const RealArray& a, const RealArray& b);
bool identicalArrays(const IntegerVector& a, const IntegerVector& b);

4.2.4  Examples

A calcNEIfractions use can be found in the file vvgnei.cxx in the directory Xspec/src/XSFunctions.

4.3  NeutralOpacity

This class serves as a limited C++ interface to the Fortran gphoto and photo routines to parallel the IonizedOpacity class. It will use the cross-sections set using the xsect command. The method Setup should be used to initialize the object then GetValue or Get to return opacities for a single or multiple energies, respectively. IronAbundance is used to set the iron abundance (relative to the defined Solar) and Abundance the abundances of all other elements. IncludeHHe specifies whether to include the contributions of hydrogen and helium in the total opacity. 
class NeutralOpacity{
 public:

  IntegerVector AtomicNumber;
  vector<string> ElementName;

  string CrossSectionSource;

  NeutralOpacity();     // default constructor
  ~NeutralOpacity();    // destructor

  void Setup();   // set up opacities
  void Get(RealArray inputEnergy, Real Abundance, Real IronAbundance, 
           bool IncludeHHe, RealArray& Opacity);  // return opacities 
  void GetValue(Real inputEnergy, Real Abundance, Real IronAbundance, 
                bool IncludeHHe, Real& Opacity);  // return single opacity 

};

4.3.1  Examples

An example use of NeutralOpacity can be found in the routine calcCompReflTotalFlux in the file MZCompRefl.cxx in the directory Xspec/src/XSFunctions.

4.4  IonizedOpacity

This class calculates opacity of a photo-ionized material. At present it uses the Reilman & Manson (1979) opacities although this could be generalized in future. The Setup method reads the input files if necessary and calculate ion fractions for the requested ionization parameter, temperature and spectrum. GetValue or Get then return opacities for a single or multiple energies, respectively. IronAbundance is used to set the iron abundance (relative to the defined Solar) and Abundance the abundances of all other elements. IncludeHHe specifies whether to include the contributions of hydrogen and helium in the total opacity. 
class IonizedOpacity{
 public:

  IntegerVector AtomicNumber;
  vector<string> ElementName;

  RealArray Energy;

  RealArray **ion;
  RealArray **sigma;
  RealArray *num;

  IonizedOpacity();     // default constructor
  ~IonizedOpacity();    // destructor

  void LoadFiles();     // internal routine to load model data files
  void Setup(Real Xi, Real Temp, RealArray inputEnergy, 
             RealArray inputSpectrum);   // set up opacities
  void Get(RealArray inputEnergy, Real Abundance, Real IronAbundance, 
           bool IncludeHHe, RealArray& Opacity);  // return opacities 
  void GetValue(Real inputEnergy, Real Abundance, Real IronAbundance, 
                bool IncludeHHe, Real& Opacity);  // return single opacity 

};

4.4.1  Examples

An example use of IonizedOpacity can be found in the routine calcCompReflTotalFlux in the file MZCompRefl.cxx in the directory Xspec/src/XSFunctions.

4.5  MZCompRefl

This class calculates Compton reflection using the Magzdiarz and Zdziarski code. 
class MZCompRefl{
 public:

  MZCompRefl();         // default constructor
  ~MZCompRefl();        // default destructor

  void CalcReflection(string RootName, Real cosIncl, Real xnor, 
                      Real Xmax, RealArray& InputX, RealArray& InputSpec, 
                      RealArray& Spref);
The easiest way to use the MZCompRefl class is through the associated function calcCompReflTotalFlux. 
void calcCompReflTotalFlux(string ModelName, Real Scale, Real cosIncl, 
                           Real Abund, Real FeAbund, Real Xi, Real Temp, 
                           Real inXmax, RealArray& X, RealArray& Spinc, 
                           RealArray& Sptot);
ModelName is the name of model and is used to check the ModelName_PRECISION xset variable to determine the precision to which internal integrals should be calculated.
Scale is the fraction of reflected emission to include. If Scale is zero then no reflection component is included, a value of one corresponds to an isotropic source above an infinite disk. The special case of minus one will return only the reflected component.
cosIncl is the cosine of the inclination angle of the disk to the line of sight. The iron abundance is specified by FeAbund and the abundances of all other elements by Abund.
For an ionized disk Xi and Temp give the ionization parameter and temperature for the IonizedOpacity class. If Xi is zero then the NeutralOpacity class is used instead. In both these cases the boolean IncludeHHe is set to false.
The input energies and spectrum are given by the X and Spinc arrays. The energies are in units of me c2 and the input spectrum is E FE. The output spectrum Sptot is also E FE. The input variable inXmax is the maximum value of X for which the reflected spectrum is calculated. This is useful because X and Spinc should be specified to a higher energy than required in the output because the energy downscattering in Compton reflection. Setting inXmax to the highest output energy required will save computation time.

4.5.1  Examples

An example use of calcCompReflTotalFlux can be found in the routine doreflect in the file ireflct.cxx in the directory Xspec/src/XSFunctions.


File translated from TEX by TTH, version 4.16.
On 22 Aug 2023, 15:06.
HEASARC Home | Observatories | Archive | Calibration | Software | Tools | Students/Teachers/Public
Last modified: Tuesday, 22-Aug-2023 15:44:29 EDT

The HEASARC welcomes your participation in a brief survey to capture how users access and utilize HEASARC data, software, and services. The outcome(s) of this survey will be used to guide, prioritize, and plan our activities and development in the coming years. It contains 20 questions, generally takes just a few minutes to complete, and your answers will remain totally anonymous. The survey is open until Dec 18, 2023. We thank you in advance for your valuable feedback.