Skip to main content

Come analyze HEASARC, IRSA, and MAST data in the cloud! The Fornax Initiative is now welcoming all interested beta users.

Xspec Home Page

next_inactive up previous












XSPEC Internal Functions Guide



For XSPEC Version 13.0.0


Keith A. Arnaud



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


August 2026


Contents

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.

The final two chapters describe the calling conventions that model functions themselves must follow: first the model function proper, then the optional analytic gradient a model may register alongside it.

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.

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.

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. The same header also provides GaussKronrodIntegrate, which performs a single (non-adaptive) Gauss-Kronrod integration of the templated integrand over the given limits.

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)

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(-z$^2$) 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)

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; the generator can be re-seeded and re-initialized at any point with the xset seed command. The generator 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)

BinarySearch

The BinarySearch routines return the index to an array (x) of the element immediately before the input test value (y). They will be more efficient if 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)

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.

CosmologyFunction

The routine FZSQ takes as input Real variables specifying the redshift, $q_0$ and $\lambda_0$. Multiplying the output by $c/H_0$ gives the luminosity distance using the approximation in Pen (ApJS 1990, 120, 49).

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.

Histogram

Classes used to generate histograms from MCMC chains.

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);

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/cm$^2$/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)

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

ModularCounter

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

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 USV$^{\rm T}$ 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).

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  
LIGHTSPEED 299792.458 defined in km/s
AMU 1.660539040e-24 unified atomic mass unit in g
EMASSINKEV 510.998950 electron mass in keV from CODATA 2018
THOMSON 6.6524587321E-25 Thomson x-section in cm$^{-2}$ from CODATA 2018
FINESTRUCT 7.2973525693E-3 Fine structure constant from CODATA 2018
PLANCK 6.62607015E-27 Planck's constant (in cgs : erg s) CODATA 2022
CLASSERAD 2.8179403262e-13 Classical electron radius in cm.

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.

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_spexVersion version of SPEX 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 CIE model temperatures calculated
vector$<$double$>$ s_DEM Last CIE model DEMs calculated
map$<$int, map$<$string, string$>$ $>$ s_XFLT XFLT#### keyword values
map$<$string,RealArray$>$ 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
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
modelDataPathGeneration   get a counter which is incremented every time the model data path changes; a model which caches loaded files can compare it against the value seen at load time to know when to re-read
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, SPEX, and NEI versions used can be set or retrieved

atomdbVersion FGATDV get AtomDB version in use
atomdbVersion FPATDV set AtomDB version in use
spexVersion FGSPXV get SPEX version in use
spexVersion FPSPXV set SPEX version in use
neiVersion FGNEIV get NEI version in use
neiVersion FPNEIV 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
getModelStringBool   get entry as a boolean (accepts yes/no/on/off/true/false/1/0, case-insensitive; used by the operational switches such as DISABLE_ANALYTIC_GRAD)
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 H$_0$, q$_0$, and $\lambda_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 H$_0$, q$_0$, and $\lambda_0$
getq0 csmgq0 get cosmology q$_0$
setq0 csmpq0 set cosmology q$_0$
getH0 csmgh0 get cosmology H$_0$
setH0 csmph0 set cosmology H$_0$
getlambda0 csmgl0 get cosmology $\lambda_0$
setlambda0 csmpl0 set cosmology $\lambda_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
getXFLTstr   get string 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
getAllXFLTstr   get all XFLT string 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
getDbArrayValues   get array of values for given keyword in internal database
loadDbArrayValues   set array of values for given keyword 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 is 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 is used. 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, tableInfo, to return information about the table in a file whose name is input.

tableInterpolate TABINT get tabulated value
tableInfo   return information about the table

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 wrapper 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

Model Function C++ Classes

The classes are

  • Aped - generates CIE 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.

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.

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;
  RealArray m_OscillatorStrength;

  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 CIE although for the latter the ElementDriver and IonDriver arrays are not required.

ApedElementRecord class 

class ApedElementRecord{
 public:

  int m_AtomicNumber;
  vector<ApedIonRecord> m_IonRecord;
  RealArray m_EquilibriumIonFraction;

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

  void LoadIonRecord(ApedIonRecord input);
  void LoadEquilibriumIonFractions(const RealArray& 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.

ApedTemperatureRecord 

class ApedTemperatureRecord{
 public:

  Real m_Temperature;
  vector<ApedElementRecord> m_ElementRecord;

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

  void LoadElementRecord(ApedElementRecord input);
  ApedElementRecord& ElementRecord(const int Z);

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

};

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

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;
  string m_ionbalname;

  bool m_noLines;
  bool m_thermalBroadening;
  Real m_velocityBroadening;
  Real m_minimumLinefluxForBroadening;
  int m_broadenPseudoContinuum;
  bool m_logTempInterpolation;
  bool m_multiThread;
  bool m_useEEbremss;
  bool m_doElectronDensityCorrection;
  Real m_rsColumn;

  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, line, and ion balance  files and stores
  // the names and temperatures. Does not read the actual data.

  int Read(string cocofilename, string linefilename, string ionbalfilename);

  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. When type is set to APEC these are APECNOLINES, APECTHERMAL, APECVELOCITY, APECMINFLUX, APECBROADPSEUDO, APECLOGTINTERP, APECMULTITHREAD, APECEEBREMSS, APECDOECOR respectively. When type is set to SPEX these are SPEXNOLINES, SPEXTHERMAL, SPEXVELOCITY, SPEXMINFLUX, SPEXBROADPSEUDO, SPEXLOGTINTERP, SPEXMULTITHREAD, SPEXEEBREMSS, SPEXDOECOR respectively. For SetVelocityBroadening a warning will be issued if the input value is non-zero and differs from APECVELOCITY/SPEXVELOCITY. If m_noLines is set then the spectrum will only be continuum. 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. m_broadenPseudoContinuum selects whether the pseudo-continuum is also broadened: 0 means no, 1 the fast method (a single broadening using the abundance-averaged thermal width), and 2 the slow method (broadening each element separately). 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. If m_useEEbremss is set then the e-e bremsstrahlung contribution will be included. If m_doElectronDensityCorrection is set then the electron to ion ratio will be calculated (instead of assuming it is 1.2) and the interpretation of the normalization changed appropriately. m_rsColumn (see SetRSColumn) sets the hydrogen column used in the resonance-scattering factor applied to the line fluxes. 

  void SetNoLines(const bool qno, const string type);
  void SetThermalBroadening(const bool qtherm, const string type);
  void SetVelocityBroadening(const Real velocity, const string type);
  void SetMinimumLinefluxForBroadening(const Real flux, const string type);
  void SetBroadenPseudoContinuum(const int pbroad, const string type);
  void SetLogTempInterpolation(const bool qltemp, const string type);
  void SetMultiThread(const bool qmulti, const string type);
  void SetUseEEbremss(const bool qeebremss, const string type);
  void SetDoElectronDensityCorrection(const bool qdoecor, const string type);
  void SetRSColumn(const Real rsColumn);

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 CIE 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 RealArray& Dem,
                  const vector<vector<RealArray> >& IonFrac, const bool isCIE,
                  RealArray& fluxArray, RealArray& fluxErrArray);

Underneath this there are routines to evaluate the continuum and line spectra for a given temperature and when including resonance scattering routines to evaluate the line spectrum for a given element (over all temperatures) and sort the lines into increasing order. 

  int calcContinuumSpectrumForTemperature(const size_t TRecordIndex,
               const IntegerVector& Zinput, const RealArray& abunZ,
               const Real Dem, const vector<RealArray>& IonFrac, 
               const RealArray& sourceFrameEnergy, const bool isCIE,
               RealArray& fluxArray);

  int calcLineSpectrumForTemperature(const size_t TRecordIndex,
               const IntegerVector& Zinput, const RealArray& abunZ,
               const Real Dem, const vector<RealArray>& IonFrac,
               const RealArray& sourceFrameEnergy, const bool isCIE,
               const Real minLineflux, RealArray& fluxArray,
               Real& maxLineflux);

  int calcLineSpecForElt(const int eltZ, const Real abund,
               const vector<int> TRecordIndexArr,
               const vector<Real> DemArr,
               const vector<RealArray>& IonFracForElt,
               const RealArray& sourceFrameEnergy, const bool isCIE,
               const Real minLineflux, RealArray& fluxArray,
               Real& maxLineflux, RealArray& sumTauArray);

  void sortLines(const vector<vector<Real> >& lineEnergy,
               const vector<vector<RealArray> >& lineWidth,
               const vector<vector<Real> >& lineFlux,
               const vector<vector<Real> >& ilneTau,
               vector<Real>& allLineEnergy,
               vector<RealArray>& allLineWidth,
               vector<Real>& allLineFlux,
               vector<Real>& allLineTau);

The number of electrons per ion for this temperature and ionization fractions are calculated by 

  void calcElectronsAndIonsForTemperature(const IntegerVector& Zinput,
               const RealArray& abunZ, const vector<RealArray>& Ionfrac,
               Real& electronNumbers, Real& ionNumbers,
               Real& HionNumbers);

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 CIE 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 calcCIESpectrum(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 string Type,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCIESpectrum(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 string Type,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCIESpectrum(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 string Type,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCIESpectrum(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 string Type, RealArray& fluxArray,
                    RealArray& fluxErrArray);

int calcCIESpectrum(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,
                    const string Type, RealArray& fluxArray,
                    RealArray& fluxErrArray);

int calcCIESpectrum(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,
                    const string Type, RealArray& fluxArray,
                    RealArray& fluxErrArray);

int calcCIESpectrum(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,
                    const string Type, RealArray& fluxArray,
                    RealArray& fluxErrArray);

int calcCIESpectrum(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,
                    const string Type, RealArray& fluxArray,
                    RealArray& fluxErrArray);

int calcCIESpectrum(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, const string Type,
                    RealArray& fluxArray, RealArray& fluxErrArray);

int calcCIESpectrum(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, const Real rsColumn,
                    const string Type, 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);

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,
                    const Real rsColumn, RealArray& fluxArray,
                    RealArray& fluxErrArray);

A few other useful routines also live in Aped.h. getAtomicMass returns the atomic mass for the given atomic number. updateNeiVersion returns the NEI version. getCIEApedFileNames returns the version number string, and the continuum, line, and ionization balance filenames for CIE models. getNEIApedFileNames returns the NEI version string and the continuum and line filenames. constructApedFilename constructs the requested filename. apedInterpFlux is used to interpolate the continuum and pseudo-continuum arrays. calcEEbrems calculates the e-e bremsstrahlung continuum. calcResScatFactor calculates the resonance scattering multiplication factor for the input optical depth (see Chakraborty et al 2023) 

Real getAtomicMass(const int& AtomicNumber);
bool updateNeiVersion(string& neiVersion);
bool getCIEApedFileNames(string& version, string& continuumFile,
                         string& lineFile, string& ionbalFile,
                         const string type);
bool getNEIApedFileNames(string& version, string& continuumFile,
                         string& lineFile);
std::string constructApedFilename(const string type, const string version,
                                  const string fileSuffix);

void apedInterpFlux(const RealArray& inputEnergy, const RealArray& inputFlux, 
                    const Real& z15, const Real& coeff,
                    const RealArray& energyArray, RealArray& fluxArray);

RealArray calcEEbrems(const RealArray& energyArray, const Real Temperature,
                      const Real ElectronRatio);

Real calcResScatFactor(const Real tau);

Examples

An example use of the Aped class to calculate a CIE 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.

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 eigenvector data files. The top-level IonBalNei class provides methods to load and use the eigenvector data.

IonBalElementRecord 

class IonBalElementRecord{
 public:

  int m_AtomicNumber;
  RealArray m_EquilibriumPopulation;
  RealArray m_Eigenvalues;
  vector<RealArray> m_LeftEigenvectors;
  vector<RealArray> m_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 CIE ion fractions.

IonBalTemperatureRecord 

class IonBalTemperatureRecord{
 public:

  Real m_Temperature;
  vector<IonBalElementRecord> m_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.

IonBalNei 

class IonBalNei{
 public:

  vector<IonBalTemperatureRecord> m_TemperatureRecord;
  RealArray m_Temperature;

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

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

  void Clear();  // clear out the contents

The filename containing the eigenvector data is found by checking the NEI version and reading the SEIGEN keyword in the NEI continuum file for v3.x otherwise it uses hardcoded names for v1.0, v1.1, and v2.0. This ensures consistency between the eigenvector data and the spectrum data. 

  string getEigenFilename();

The eigenvector data is read using ReadElements. 

  int ReadElements(const int Z);
  int ReadElements(const vector<int> Z);

  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 CIE method returns the ionization fractions for the specified element for a collisional ionization equilibrium plasma with the given electron temperature. 

  RealArray CIE(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 NEIVERS 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 calcCIEfractions(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. 

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);

Examples

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

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 

};

Examples

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

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 

};

Examples

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

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 $m_e c^2$ and the input spectrum is $E F_E$. The output spectrum Sptot is also $E F_E$. 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.

Examples

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

Model Functions

A model function is the routine that calculates a model component's spectrum, given an array of energy-bin boundaries and an array of parameter values. This chapter is the reference for the calling conventions such a function must follow; the user-level how-to for building and loading a local model package (initpackage, lmod, and the full field-by-field description of the model.dat format) is in the “Adding models to XSPEC” appendix of the XSPEC manual and is not repeated here.

Registration

A model is declared to XSPEC by an entry in a model.dat-format file. The first line names the model, gives the number of parameters, the valid energy range, the function to call, and the model type (add, mul, con, mix, acn, or amx); subsequent lines describe the parameters: 

modelentry        4  0.    1.e20      C_funcName    add  0 0
lowT    keV     0.1   0.0808  0.0808 79.9      79.9       0.001
...

The prefix on the function-name field selects the calling convention: no prefix for single-precision Fortran, F_ for double-precision Fortran, c_ for C-style arguments, C_ for C++-style arguments, and Py_ for a Python function. The prefix appears only in model.dat; the function in the source code is the bare funcName. The cxsetup and initpackage tools read this file and generate the dispatch code that binds the name to the function, so a new function never has to be registered by hand.

Calling conventions

The four compiled calling conventions are defined as typedefs in funcType.h: 

// single-precision Fortran (no prefix)
typedef void (xsf77Call) (const float* energyArray,
                          const int& Nenergy,
                          const float* parameterValues,
                          const int& spectrumNumber,
                          float* flux,
                          float* fluxError);

// double-precision Fortran (F_)
typedef void (xsF77Call) (const double* energyArray,
                          const int& Nenergy,
                          const double* parameterValues,
                          const int& spectrumNumber,
                          double* flux,
                          double* fluxError);

// C-style (c_)
typedef void (xsccCall)  (const Real* energyArray,
                          int Nenergy,
                          const Real* parameterValues,
                          int spectrumNumber,
                          Real* flux,
                          Real* fluxError,
                          const char* initString);

// C++-style (C_)
typedef void (XSCCall) (const RealArray& energyArray,
                        const RealArray& parameterValues,
                        int spectrumNumber,
                        RealArray& flux,
                        RealArray& fluxError,
                        const std::string& initString);

We recommend the C++ convention for new models written in C++ (it avoids copying between C arrays and the internal data structures) and, to prevent unresolved-symbol linkage errors, prefacing the function definition with extern "C". Conventions common to all four forms:

  • energyArray holds the boundaries of contiguous, ascending energy bins, so it has one more entry than the flux array: nE+1 boundaries for nE flux bins (in the Fortran forms Nenergy is nE and the array is dimensioned ear(0:ne)). The energies are set by the responses in use; the function must make no assumptions about the range or bin sizes.
  • parameterValues holds the parameter values in model.dat order. The norm of an additive model is not included; it is applied by the dispatcher after the call.
  • For an additive model, flux bin i is the spectrum integrated over the bin, in photons cm$^{-2}$ s$^{-1}$ (not per keV). For a multiplicative model it is the dimensionless multiplicative factor for the bin. On input the flux array contents are undefined except for convolution models (below).
  • fluxError allows a model to return variances on the flux. Almost all models do not: a C++ function then leaves the array at size 0, a C or Fortran function simply does not write to it. A model which does fill it must also set the variance flag (the first of the two trailing flags) in its model.dat line.
  • spectrumNumber identifies which spectrum's response generated these energies. Most models ignore it; it exists for multi-dimensional models whose output depends on the spectrum through more than the energy grid (such models must also set the second trailing model.dat flag to force a calculation per spectrum).
  • initString (C and C++ forms only) receives the optional initialization string appended to the model's model.dat line, read once at initialization – typically the name of a data file, so one function can serve several models that differ only in their input data.
  • Memory: a C++ function must resize flux to energyArray.size()-1; C and Fortran functions must not allocate the output arrays, which the wrapper has already sized.

Real is a typedef for double (in xsTypes.h) and RealArray is a resizeable array of Real, currently implemented as std::valarray<Real>. Avoid valarray-specific features such as slicing in case the typedef changes.

Convolution models

A convolution (con or acn) model uses the same signatures as above, but its flux argument is both input and output: on entry it contains the calculated flux of the model components the convolution operates on, and the function overwrites it with the convolved result. A convolution model must therefore cope with whatever spectrum it is handed and, like every model, with whatever energy grid the responses impose – in particular a kernel that redistributes flux beyond the grid edges will lose flux there unless the user extends the grid with the energies command.

Mixing and response models

Mixing models (mix and amx), which transform fluxes across data groups rather than computing a spectrum, use separate conventions (also in funcType.h): XSMixCCall receives per-data-group energy arrays and a GroupFluxContainer of fluxes to transform in place, together with a MixUtility object that carries the model's per-session state; xsmixcall is the C-style equivalent. Response-modifying functions use XSRespCCall/xsrespcall, which receive only the parameter values and a RespUtility object. Writing a mixing model also requires implementing the MixUtility initialization interface; see the “Writing new mixing models” section of the manual appendix.

Model functions in Python

A Py_<module>.<func> entry names a function in a Python module distributed beside the compiled package. The function takes the same arguments as the compiled forms – energies, parameters, flux, and optionally a flux-error array and the spectrum number – passed as Python tuples with the flux as a pre-sized list to fill in place, or as zero-copy numpy arrays if the function is decorated with @xspec.vectorized. The full interface, including the flux-error convention, is in the “Local Models in Python” section of the PyXspec manual.

Analytic Gradients

As of version 13.0 most built-in models carry an analytic gradient: a function which returns the derivative of the model flux with respect to each of the model's parameters. When every component of the active models has one, the default Levenberg-Marquardt fitter and the Minuit methods compute the fit derivatives analytically instead of by finite differences of the folded model, and the hmc sampler uses the same machinery for its posterior gradient (in reverse mode). The user can revert the fitters to finite differences with xset DISABLE_ANALYTIC_GRAD yes. This chapter describes how a model registers a gradient and the conventions its implementation must follow. The user-level description of the pipeline is in the Algorithms appendix of the XSPEC manual.

Registration

A model declares its gradient in model.dat by appending a grad= flag to the model definition line: 

vlorentz  2  0.  1.e20  C_vlorentzianLine  add  0
    grad=gv:vlorentzianLineGradient,vlorentzianLineVJP

(shown wrapped; in model.dat the flag is on the definition line). The flag takes the forms grad=g (forward gradient only) or grad=gv (gradient plus a reverse-mode vector-Jacobian product, see below), optionally followed by a colon and explicit function names. Without explicit names the functions are assumed to be called nameGradient and nameVJP for model name. The cxsetup and initpackage tools read these flags and auto-generate gradientFunctionMap.cxx, which registers the functions through 

XSModelFunction::addGradientFunctionPointer(name, func);
XSModelFunction::addVJPFunctionPointer(name, func);
XSModelFunction::addConvGradientFunctionPointer(name, func);
XSModelFunction::addConvVJPFunctionPointer(name, func);

Do not edit gradientFunctionMap.cxx by hand; it is generated.

The forward gradient (XSCCGrad)

All forward gradients share the XSCCGrad signature defined in funcType.h: 

typedef void (XSCCGrad) (const RealArray& energyArray,
                         const RealArray& parameterValues,
                         int spectrumNumber,
                         std::vector<RealArray>& dFlux_dParam,
                         const std::vector<bool>& parThawed,
                         const std::string& initString);

The implementation must resize dFlux_dParam to (nModelParams, nE-1) and fill row k with the derivative of the flux in each energy bin with respect to parameter k. Conventions:

  • The additive-model norm is not included: its derivative is supplied generically by the dispatcher via the chain rule. A model with no parameters other than the norm still registers (its gradient is empty) so that it does not force the whole fit onto the finite-difference fallback.
  • parThawed (when non-empty) flags which columns the fit will actually read. The implementation may zero the column and skip the work for entries that are false; an empty vector means treat all parameters as thawed.
  • An “analytic” gradient need not be closed-form. Many registered gradients evaluate a five-point central finite difference of the model's own forward function with a step chosen to suit that model (clearing table-cell quantization noise, respecting hard limits by switching to one-sided stencils, treating integer switch parameters as frozen by convention). This still beats the generic fit-level finite difference because the stencil is applied to the unfolded model with a model-appropriate step, once, rather than to the folded statistic with a generic step.
  • A parameter whose value is a discrete selector (an integer cast inside the forward, or a $-prefixed switch) has no meaningful derivative; its column should be zero.

The reverse-mode product (XSCCGradVJP)

For high-dimensional or expensive models the full Jacobian is wasteful when only the scalar gradient of the statistic is needed (as in HMC). A model may therefore also register a vector-Jacobian product: 

typedef void (XSCCGradVJP) (const RealArray& energyArray,
                            const RealArray& parameterValues,
                            int spectrumNumber,
                            const RealArray& dStat_dFlux,
                            RealArray& dStat_dParam,
                            const std::vector<bool>& parThawed,
                            const std::string& initString);

Given the per-bin adjoint dStat/dFlux it returns, in one backward pass, dStat_dParam[p] $= \sum_i$ dStat_dFlux[i] $\cdot$ dF$_i$/d$\theta_p$ without materializing the Jacobian. Models without a VJP transparently fall back to the forward gradient plus a dot product, so registering one is an optimization, never a requirement.

Convolution models

Convolution components use separate signatures, XSCCConvGrad and XSCCConvGradVJP (also in funcType.h), because they must propagate derivatives through the convolution. XSCCConvGrad receives the upstream flux and the upstream Jacobian and returns the output flux together with the output Jacobian over [own parameters] ++ [upstream parameters]; the chain rule through the convolution operator is the convolution's responsibility, because only it knows the operator's structure (a diagonal scaling for cflux, a sparse rebin for zashift, and so on). A convolution whose output is not linear in its input flux (for example one that renormalizes its output) must not apply the linear shortcut Conv(dF$_{\rm in}$/dp) for the upstream rows; the standard treatment is a finite difference over a scalar $\epsilon$ in Conv(F $_{\rm in} + \epsilon\,$dF$_{\rm in}$/dp). XSCCConvGradVJP returns the gradient over the convolution's own parameters together with the upstream adjoint dStat/dF$_{\rm in}$ which the backward walker propagates further up the chain; where the kernel transpose is not available in closed form the dispatcher falls back to the forward Jacobian path, which is correct but slower.

Validating a new gradient

Two validation layers exist and both should be used when adding a gradient:

  • xset HYBRID_GRAD_CROSS yes makes the Levenberg-Marquardt fitter compute every iteration's derivatives by both the analytic and the finite-difference path and print a comparison, with per-parameter detail at higher chatter levels.
  • The gradient-audit regression suite (testCases/extended/gradient-audit in the test-suite tree) checks every registered gradient against a five-point finite difference of the same forward function and freezes the result in a baseline, so a change to any gradient or forward function is caught. A new model with a registered gradient is picked up automatically.

When a gradient and the finite-difference reference disagree, be aware that the disagreement is not always the gradient's fault: a parameter sitting on a hard limit makes a central stencil invalid, a tabulated forward function has plateaus whose finite difference is noise, and a breakpoint-energy parameter cannot be finite-differenced at any fixed step. Scanning the finite difference against the step size at the worst bin usually identifies which side is wrong.

About this document ...

This document was generated using the LaTeX2HTML translator Version 2024 (Released January 1, 2024)

The command line arguments were:
latex2html -split 0 –local_icons -no_auto_link XspecInternalFunctionsGuide

The translation was initiated on 2026-08-13

next_inactive up previous