The fitsspec library: Writing XSPEC-ready FITS files
                   
Version 0.9

Introduction
------------
The routines in fitsspec are designed to open, write, and close
FITS files that contain spectral data, such as line emissivities and
continua.  There are 5 primary calls, one to open and initialize a
FITS file, three to write the four different types of files, and one
to close the files.   fitsspec uses the cfitsio library extensively to
write the files, and assumes that the cfitsio library is in the
include-file path during compilation (see the Makefile).  Currently,
all integer values are int's, and all real values are double although
when stored in the FITS file they are both 4 byte values.

The four file types are:

0: Line emissivity.  These files have one HDU for each set of
parameter values.  For example, a collisionally ionized plasma would
have one HDU for each pair of temperature and density (T,n).
Temperatures are given in Kelvin, while densities are in cm^-3.  Each HDU 
contains rows of the line wavelength lambda, and emissivity epsilon,
along with errors on each and the ion the line comes from.  The upper
and lower levels of the ion can also be included.

1: Binned continuum output.  This file contains the continuum
emissivities on a per process basis.  Again, each HDU is a (T,n) pair,
but this time the rows contain an evenly binned spectrum in units of photons
cm^3 / s / bin, along with each of major continua stored separately:
bremsstrahlung, radiative recombination, two-photon, and (optionally)
the pseudo-continuum.  The pseudo continuum consists of lines deemed too weak
to include in the individual line list, and so the are simply binned
and put here.

2: Linearly-compressed continuum output.  This is the primary
continuum output, and again each HDU contains one pair of (T,n).  Each
HDU then has a row for each atom, giving the total "true" continuum
from that atom (summed over all ions of the atom and over the
bremsstrahlung, two-photon, and radiative recombination continuum from
each ion) and the "pseudo" continuum.  When passed to the library, the
continuum is in a binned format.  However, in this case both the true
and pseudo continua are "linearly compressed" by the library
code. This means that instead of giving the emissivity in a per-bin
form, the spectrum is given only at a few points, in units of photons
cm^3 / s / keV.  The emissivity between any two points can then be
calculated by linearly interpolating between them.  The maximum error
induced by the interpolation is 1%; this can be set at the top of the
fits_spec_write_coco.c file.

3: Similar to #2 except the emissivity is done on a per-ion format
instead of per-atom.  

Installation
------------

In order to be compiled, the cfitsio v2.0.33 library or later is
required.  This can be obtained from
http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html
By default, fitsspec assumes that the cfitsio directory is in
../cfitsio/, but this can be changed in the Makefile.   

To create the library, edit the Makefile to ensure that CFITSIODIR
points to the correct location, and then type make.  This will compile
the library, and the test program fitsspec_test.  Running this program
will create two files, Test_line.fits and Test_coco.fits.  These files
can be compared to Compare_Test_line.fits and Compare_Test_coco.fits
to ensure that everything went well.

NOTE: To compile on a Solaris machine, the libraries -lnsl and
-lsocket are required.  For Linux or Dec Alphas running whatever Compaq
is calling Unix these days, these libraries must be removed from the
Makefile.

Procedure Calls
---------------

The procedure calls are stored in the file fitsspec.h, and are as 
follows:

void fits_spec_init(char *filehead, int filetype, int Nrows, fitsfile **fptr);

void fits_spec_write_lines(fitsfile *fptr, double T, double n_e, int nlines, 
  double lambda[], double lambda_err[], double epsilon[], double epsilon_err[],
  int el_Z[], int rmJ[], int UpperLev[], int LowerLev[]);

void fits_spec_write_continuum(fitsfile *fptr, double T, double n_e,
  int nbin, double Ebin_low[], double Ebin_high[],
  double pseudo[], double pseudo_err[], double brems[], double brems_err[],
  double twoph[], double twoph_err[], double radrec[], double radrec_err[]);

void fits_spec_write_coco(fitsfile *fptr, int N_Z, int list_Z[], 
  double T, double Ne, int nbin,  double E_low[], double E_high[], 
  double pseu_arr[], double cont_arr[]);

void fits_spec_close(fitsfile *fptr);


Sample Use
----------

Assume we want to write two files, a line emissivity file (type 0) and a
compressed continuum file (type 2).  Begin by initializing both
files as shown:


fits_spec_init("MySpec", 0, Nrows, &fptr_line);
fits_spec_init("MySpec", 2, Nrows, &fptr_coco);
/*
Here "MySpec" will be the first part of the filename; _line.fits and
_coco.fits will be automatically appended to the name to make the new
files MySpec_line.fits and MySpec_coco.fits.  The Nrows variable is
used if you know how many total HDUs you will be writing in advance.
For example, if you'll be calcuating 40 different temperatures at 2
densities each, Nrows = 40x2 = 80.  However, Nrows can be set to 0 if
the proper value is not known.  (Nrows is used to reserve space in the
initial binary table which lists all the parameter values used;
setting it to zero simply means some space will be wasted in the
file.)  The fptr variable is a fitsfile pointer, as defined by
cfitsio.
*/
/* These next two routines write the line and compressed continua files.
fptr_line and fptr_coco are the same fitsfile pointers from above. The
T and n_e variables are the temperature and electron density,
respectively.  
*/
while (!done) {
  fits_spec_write_lines(fptr_line, T, n_e, nlines, lambda, lambda_err, 
                        epsilon, epsilon_err, el_Z, rmJ, UpperLev, LowerLev);
/* 
In the line procedure, nlines gives the number of lines
in the lambda and epsilon arrays.   The lambda, lambda_err arrays
give the wavelengths (in Angstroms) and the errors on those
wavelengths, again in Angstroms.  The epsilon and epsilon_err arrays
give the emissivity (in photons cm^3 / s) and the error in the same
units.  If errors are not available, you can pass a NULL pointer
instead and the errors will be set to zero.  The el_Z, rmJ, UpperLev,
and LowerLev arrays give the elemental Z, the ionization state (in
standard roman numeral form, so 1 = I = neutral atom), and the upper
and lower levels of the ion, respectively.  The upper and lower levels
must be defined relative to some other energy level file, and so they
may also not be available.  If NULL is passed, their values will be
set to -1.
*/
  fits_spec_write_coco(fptr_coco, N_Z, list_Z, T, n_e, nbin,
                       E_low, E_high, pseu_arr, cont_arr);
/*
In the compressed continuum procedure, things get mildly complicated.
This is due to wanting to write routines that could be called either
from C or Fortran, and the difference in how 2d arrays are handled by
these languages.  So, here N_Z is the number of elements, and list_Z
is an array of the Z values.  T and n_e are the temperature and
density, and nbin is the number of energy bins.   The next four
variables, E_low, E_high, pseu_arr, and cont_arr are actually 2d
arrays of Z and Bin index.  Assume iZ is such that list_Z[iZ] = Z.
Then (iZ + iBin*N_Z) is the position of the (Z, BinIndex) element.  
So E_low[iZ + iBin*N_Z] is the lower edge of the energy bin (in keV) for
element Z and position iBin, and E_high[iZ + iBin*N_Z] is the upper
(also in keV).  The pseudo-continuum from this element and bin is
pseu_arr[iZ + iBin*N_Z] in phot cm^3 / s / bin, and the total
continuum emission is cont_arr[iZ + iBin*N_Z] in the same units. 
*/
}

fits_spec_close(fptr_line);
fits_spec_close(fptr_coco);
/* 
This closes each of the files, using the cfitsio close function. 
*/