Analytic Model Components
Additive, Multiplicative, Convolution/Pileup Components
If you have a model that will be used frequently, you may want to include it among the standard models so that it will be listed after the model ? command. To do so, you first must create a subroutine that calculates the model spectrum given an input array of energy bins and an array of parameter values. The input array of energy bins gives the boundaries of the energy bins and hence has one more entry than the output flux arrays. The energy bins are assumed to be contiguous and will be determined by the response matrix in use. Your model should thus make no assumptions about the energy range and bin sizes. The output flux array for an additive model should be in terms of photons cm2 s-1 (not photons cm2 s-1 keV-1) so is the model spectrum integrated over the energy bin. For a multiplicative model the output array is the multiplicative factor for that bin. Convolution models are operators on additive model fluxes (perhaps modified by multiplicative factors).
Except on the Cygwin platform, user model components are added to XSPEC12 using the two commands initpackage, which prepares and compiles a library module containing them, and the lmod command which actually loads them into the program. These commands are described in the XSPEC commands section of the manual, to which the user is referred. Any number of different user model packages may be added to XSPEC from the user prompt, and the user has control over the directory from which models are loaded. On Cygwin, local model libraries must be built and linked to XSPEC statically prior to runtime, and therefore cannot be loaded with lmod. Cygwin users should consult the static_initpackage section of the initpackage command documentation.
Note that the lmod command requires that the user has write-access to the particular directory they specify. This is because lmod utilizes the Tcl ‘make package’ and ‘package require’ mechanisms for doing automatic library loads, which then require that Tcl be able to write an index file to the directory (pkgIndex.tcl). Therefore for those who wish to load a local model library for many users across a system (by placing the command in the global_customize.tcl script, see the section “Customizing XSPEC”), we recommend directly calling the Tcl load command instead.
XSPEC12 explicitly supports model components written in fortran, either in single or double precision, in C++ using either C++-style arguments or C style arguments, and in C, as described below. The language and argument types of the implementation are specified in the initialization “model.dat” format file. In the spirit of the revisions, XSPEC11’s single implementation of model components in single precision fortran requires no change in the syntax previously used.
A sample model.dat entry has the following form:
modelentry 5 0. 1.e20 modelfunc add 0
lowT keV 0.1 0.0808 0.0808 79.9 79.9 0.001
highT keV 4. 0.0808 0.0808 79.9 79.9 0.001
Abundanc " " 1. 0. 0. 5. 5. 0.01
*redshift " " 0.0
$switch 1
The first line for each model gives the model name, the number of parameters, the low and high energies for which the model is valid, the name of the subroutine to be called and the type of model (add, mul, mix, or con, or acn). The final argument is a flag, which should be set to 1 if model variances are calculated by modelfunc.
The remaining lines give the specifications for each parameter in the model. For regular model parameters the first two fields are the parameter name followed by an optional units label. If there is no units label then there must be a quoted blank (“ “) placeholder. The remaining 6 numerical entries are the default parameter value, hard min, soft min, soft max, hard max, and fit delta, which are described in the newpar command section. Note that scale parameters only require the parameter value since they are always frozen. Switch parameters only have 2 fields: the parameter name and an integer value.
Two new features of XSPEC12 are:
Designation of model function call type. The function call type is indicated by a prefix, if different from the previous XSPEC type, as in the following table. Note that in all cases, the function/subroutine name in the code is the same (i.e., in the example below it is always modelfunc, not the modified name that must appear in the file).
|
Call Type |
Specification |
Arguments and Type |
Meaning |
|
Single precision fortran |
modelfunc |
real*4 ear(0:ne) |
Energy array |
|
integer ne |
Size of flux array |
||
|
real*4 param(*) |
Parameter values. (Dimension must be specified inside the function) |
||
|
integer ifl |
The spectrum number being calculated |
||
|
real*4 photar(ne) |
Output flux array |
||
|
real*4 photer(ne) |
Output flux error array (optional) |
||
|
Double precision fortran |
F_modelfunc |
real*8 ear(0:ne) |
As above |
|
integer ne |
“ |
||
|
real*8 param() |
“ |
||
|
integer ifl |
“ |
||
|
real*8 photar(ne) |
“ |
||
|
real*8 photer(ne) |
“ |
||
|
C/C++, C-style |
c_modelfunc |
const Real* energy |
Energy array (size Nflux+1) |
|
int Nflux |
Size of flux array |
||
|
const Real* parameter |
Parameter values |
||
|
int spectrum |
Spectrum number of model component being calculated |
||
|
Real* flux |
Output flux array |
||
|
Real* fluxError |
Output flux error array (optional) |
||
|
const char* init |
Initialization string (see below) |
||
|
C++, C++ style |
C_modelfunc |
const RealArray& energy |
Energy array |
|
const RealArray& parameter |
Parameter values |
||
|
int spectrum |
Spectrum number of model component being calculated |
||
|
RealArray& flux |
Output flux array |
||
|
RealArray& fluxError |
Output flux error array (optional) |
||
|
const string& init |
Initialization string (see below) |
For example, if one were to implement a model component in double precision fortran, the top line of this description would be:
modelentry 5 0. 1.e20 F_modelfunc add 0
XSPEC12 would then pick out the right function definition, and call the function modelfunc which expects double precision arguments, and similarly for the C and C++ call types. The C-style call can clearly be compiled and implemented by either a C or a C++ compiler: however we recommend using the C++ call if the model is written in C++ as it will reduce overhead in copying C arrays in and out XSPEC’s internal data structures. To prevent unresolved symbol linkage errors, we also recommend prefacing your C++ local model function definitions with the extern “C” directive.
Example C/C++ function definitions:
/* C -style */
extern “C” void modelfunc(const Real* energy, int Nflux, const Real* parameter, int spectrum, Real* flux, Real* fluxError, const char* init)
{
/* Your model code: Do not allocate memory for flux and fluxError arrays. XSPEC’s
C-function wrapper will allocate arrays prior to calling your function (and will free them
afterwards). */
}
// C++
extern “C” void modelfunc(const RealArray& energy, const RealArray& parameter, int spectrum, RealArray& flux, RealArray& fluxError, const string& init)
{
// Your model code: Should resize flux RealArray to energy.size() - 1. Do the same for
// fluxError array if calculating errors, otherwise leave it at size 0.
}
Note on type definitions for (C and C++): XSPEC12 provides a typedef for Real, in the xsTypes.h header file. The distributed code has
typedef Real double;
i.e. all calculations are performed in double precision. This is used for C models and C++ models with C-style arguments.
The type RealArray is a dynamic (resizeable) array of Real. XSPEC12 uses the std::valarray template class to implement RealArray. The internal details of XSPEC12 require that the RealArray typedef supports vectorized assignments and mathematical operations, and indirect addressing (see C++ documentation for details). However, we do not recommend using specific features of the std::valarray class, such as array slicing, in case there is a future need to change the typedef.
The input energies are set by the response matrices of the detectors in use. IFL is an integer which specifies to which response (and therefore which spectrum) these energies correspond. It exists to allow multi-dimensional models where the function might also depend on eg pulse-phase in a variable source. The output flux array should not be assumed to have any particular values on input. It is assumed to contain previously calculated values only by convolution/pileup models, which have the nature of operators. The output flux error array allows the function to return model variances.
The C and C++ call types allow one extra argument, which is a character string that can be appended to the top line of the model component description. This string is read on initialization and available to the model during execution. An example of its use might be the name of a file with specific data used in the model calculation: this allows different models to be implemented the same way except for different input data by specifying different names and input strings.
Parameter Types: The second new feature of XSPEC12 is the ability to designate parameters of different type, as described earlier. Using these different parameter types is optional, and designed to allow a greater level of robustness.
· If the name of the parameter is prefixed with a “*” the parameter is a “scale” parameter and cannot be made variable or linked to an other kind of parameter other than another scale parameter. This allows the designation of parameters that should never be made variable.
· If the name of the parameter is prefixed with a “$” the parameter is a “switch” parameter which is not used directly as part of the calculation, but switches the model component function’s mode of operation (i.e. calculate or interpolate). Switch parameters are of integer type.
· If neither * or $ begins the parameter name, the parameter is allowed to be variable.
Mixing models are fundamentally different from the other kinds of models since they apply a transformation to a set of modeled fluxes (as enumerated by the spectra in the fit), rather than modify the flux designed to fit a single spectrum. The need to store temporary results, as well as the requirements of the model calculation, uses a lot of workspace arrays within the model: further, the transformations applied are often fixed during a fit, or can be split into parts that are fixed and parts that change during iteration in order to avoid redundant calculations. Since XSPEC12’s internal organization (data structures) can be mapped straightforwardly to the requirements of these models, and therefore to implement them efficiently and handle memory allocation, we recommend that these models should be written in C++ or C. In the initial release only a C++ implementation for mixing models will be available. Users considering adding new mixing model types should contact the developers of XSPEC at xanprob@athena.gsfc.nasa.gov
A very simple way of fitting with user-defined models is available for a particular class of models. These are models that can be defined by a grid of spectra, with the elements of the grid covering the range of values of the parameters of the model. For instance, for a one-parameter model, a set of model spectra can be tabulated for different values of the parameter (P1, P2, P3, etc.) The correct model spectra for a value P then is calculated by interpolation on the grid. The generalization to more parameters works in the obvious way.
Table model components can be much slower than most analytical models if there are significant numbers of parameters. The memory requirements increase as 2n where n is the number of parameters in the model. A table model with more than 3 or 4 fitting parameters is not recommended. Additionally, the interpolation is linear, which implies that the second derivatives used by the default Levenberg-Marquadt algorithm may not be accurate. If difficult in fitting ensues, the user may try the migrad (minuit library) algorithm which makes no assumptions about the second derivative.
XSPEC12 permits any number of table model components to be used simultaneously.
As with standard models, the spectra should be in terms of flux-per-bin and not flux-per-keV. Any set of energy bins can be used, and XSPEC will interpolate the model spectra onto the appropriate energy bins for the detectors in use. It is therefore a good idea to choose energy bins such that the spectrum is well-sampled over the range of interest. The file structure for these models is a FITS format and is described by OGIP memo OGIP/92-009 also available by anonymous ftp from ftp://legacy.gsfc.nasa.gov/fits_info/fits_formats/docs/general/ogip_92_009