\chapter{Theory of Operation}
\label{sec:internals}

\section{XSTAR}

\subsection{Introduction}

Although XSTAR is designed with the goal of maximizing the flexibility 
of parameter values and assumptions available to the user, there are 
likely to arise situations in which the standard set of input options are 
not sufficient.  Under these circumstances the user may want to attempt to 
modify the source code in order to effect a particular set of assumptions, 
geometry, etc.  
Problems for which customization is more likely to be necessary include 
models with additional heating or cooling processes (for example  adiabatic 
expansion cooling or cosmic ray heating) or different assumptions about line 
escape probability (e.g. Sobolev escape probability in a medium with a 
velocity gradient).  If so, the internal operation of the code must be 
confronted.  This appendix presents a summary of the code 
operation in order to aid in code customization.  

\subsection {Programming Philosophy}

XSTAR has been written and developed over a span of approximately 15 years, 
and much has changed during that time.  Changes include advances in the 
fortran language, changes in the speed and memory capacity of the available 
computers, new insights into the flexibility required of XSTAR, and insight 
into which physical processes are likely to have the greatest effect on the 
model results.  

XSTAR was written in standard fortran 77, and much of the 
code retains the programming style associated with the older versions 
of the fortran language.  Hard experience with moving the code from 
one type of machine to another has led to an attempt to avoid any machine 
dependent extensions to the fortran language, any word length dependent 
numerical constructs or any use of extended precision arithmetic 
(with one or two exceptions).

The code is structured in an attempt to be modular, 
and to separate the calculation of the atomic rates from the 
calculation of level populations, ion fractions, temperature, etc.
The goal is to make it relatively easy to add or change atomic data, 
requiring modification of just one subroutine if a new process is added
or a fitting formula is changed.  The data itself is all read in 
from an external file, so that it can be changed without any 
modification to the code itself if the existing fitting formulas 
are unchanged.

The remainder of this appendix is organized as follows.  First we present 
a list of the subroutines and a short description of their function.
This list is sorted by category:  subroutines used in input, output, 
defining ionizing spectrum, primary driving subroutines, and rate calculation.
Then we present a list of variable names and a description of their contents.
Finally we present a schematic flow chart of the operation of XSTAR.

\subsection {List of Subroutines}

\subsubsection {Primary Computational Subroutines}

\begin{description}

\item[xstar:] Main program.

\item[cndnse:] condense matrix for level population calculation

\item[dsec:]calculate thermal equilibrium and charge neutrality using secant 
method

\item[ener:] set up energy grid

\item[func:]calculate all ion fractions, level populations, heating, cooling.

\item[func1:]calculate rates affecting  ion fractions

\item[func2:]calculate rates affecting level populations

\item[func2a:]calculate rates affecting level populations due to Auger 
and inner shell photoionization

\item[func2i:]calculate number of bound levels for an ion

\item[func3:]calculate heating-cooling rates

\item[init:]initialization.  Zeroes most variables

\item[heatt:]total heating and cooling rate calculation

\item[hcor1:]Level population calculation

\item[ioneqm:]calculation ion fractions

\item[istruc:]calculate ion abundances

\item[invert:] used to prepare for next iteration of global calculation

\item[step:]calculates spatial step size

\item[stpcut:] update important quantities after each spatial step.

\item[trnfrc:]calculates dilution and attenuation of incident radiation

\item[trnfrn:]calculates transfer of diffusely emitted radiation.

\end{description}

\subsubsection {Rate Calculation Subroutines}

\begin{description}

\item[amcol:] calculate rate of collisional angular momentum changing
transitions in hydrogenic ions using method of Pengelly \& Seaton
(1964).

\item[amcrs:] calculate rate of collisional angular momentum changing 
transition in hydrogenic ions using either {\bf amcol} or {\bf impact}.

\item[anl1:] calculate allowed radiative transition rates for hydrogenic
ions.

\item[bkhsgo:]

\item[calt60\_62:] returns Callaway Upsilons for H-like
ions.

\item[calt66:] returns Kato \& Nakazaki (1989) Upsilons for He-like
ions.

\item[calt67:] returns Keenan et al.(1987) Upsilons for He-like ions.

\item[calt68:] returns Sampson \& Zhang Upsilons for He-like ions.

\item[calt69:] returns Kato \& Nakazaki (1989) Upsilons for He-like
ions.

\item[calt70:] takes type 70 data and calculates recombination rates and
phot. cross sections of superlevels. 

\item[calt71:] returns radiative trans. rates from superlevels to
spectroscopic levels.

\item[calt72:] returns capture rates for DR through satellite levels.

\item[calt73:] returns Upsilons for satellite levels of He-like ions.

\item[calt74:] returns levels specific DR and calculates photoionization
rates over DR associated resonances. 

\item[calt77:] returns collisional rates between superlevels and
spectroscopic levels.

\item[erc:] calculate collisional excitation rates for hydrogenic ions.

\item[hgf:] calculates hypergeometric functions.

\item[impact:] impact parameter cross sections. 

\item[impactn:] calculate collisional excitation rates for hydrogenic
ions using the impact parameter. 

\item[impcfn:] calculate the functions used in the impact parameter.

\item[intin:] 

\item[intin2:]

\item[intino:]

\item[irc:] calculate collisional ionization rates using the Johnson
(1972) method for hydrogenic ions or routine {\bf szirc} for other 
ions. 

\item[levwk:] calculates partition functions

\item[milne:] calculates milne integral

\item[phintf:] photoionization rate and milne rate integrator.

\item[szcoll:] calculates electron impact excitation rates from semiempirical 
 formula (eq.35) from Sampson $\&$ Zhang (1988, apj 335, 516)

\item[szirc:]  calculates electron impact ionization rates from semiempirical
  formula (eq.35) from Sampson $\&$ hang (1988, apj 335, 516)


\item[ucalc:] main rate calculation routine

\item[upsil:] returns Upsilons from Burgess \& Tully fits.

\item[velimp:] calculate collisional excitation rate from the impact
parameter. 

\end{description}


\subsubsection {Input Subroutines}

\begin{description}

\item[readtabl1:] reads atomic pointer fits file

\item[readtabl2:] reads atomic data fits file

\item[rread1:] reads in commands.  calls xpi routines


\item[unsavd:] the opposite of savd

\end{description}


\subsubsection {Output Subroutines}

\begin{description}

\item[bnchmrk:] comparison with Lexington benchmarks

\item[deletefile:] deletes unneeded fits files

\item[fheader:] writes fits header

\item[fitsclose:] close fits files

\item[fparmlist:] writes fits parameter list

\item[fstepr:]

\item[fwrtascii:] writes ascii fits  files

\item[savd:]  Saves temporary data after each spatial zone for 
iterative calculation

\item[pprint:] main printout routine.

\item[printerror:] Prints errors encountered during fits i/o.

\item[writespectra:] writes out fits file with spectrum at end of calculation.


\end{description}


\subsubsection {Atomic Database Subroutines}

\begin{description}

\item[dbwk:] main database manipulation routine. 
sets up atomic database pointers among other things

\item[dprint:] write out a record to the database

\item[dprinto:] write out a record in ascii

\item[dprints:] write a one line summary of record in ascii

\item[dread:] read a record

\item[dreado:] read a record in ascii


\end{description}

\subsubsection {Miscellaneous Subroutines}

\begin{description}

\item[bremem:]calculate bremsstrahlung emissivity


\item[comp:]calculate non-relativistic Compton heating and cooling

\item[dexpo:] real*8 expo

\item[dfact:] real*8 fact

\item[ee1exp:]first exponential integral*exponential

\item[ee1expo:]first exponential integral*exponential

\item[eint:]first exponential integral

\item[enxt:] find next energy for trapezoid quadrature

\item[expo:]exponential function with boundaries

\item[exint1:]

\item[exp10:] $10^x$ with boundaries

\item[expint:]first exponential integral*exponential

\item[fact:] factorial

\item[fact8:] real*8 factorial

\item[fbg:]used in calculating bremsstrahlung emissivity (\cite{Raymond1976})

\item[hlike:]hydrogenic photoionization cross section

\item[hunt:]table search (\cite{Press})

\item[huntf:]table search assuming logarithmic spacing

\item[ispec:]input spectrum from thermal brems with unit gaunt factor

\item[ispcg2:]

\item[ispec4:]input spectrum, single power law

\item[ispecg:]input spectrum, from atable

\item[ispecgg:]input spectrum, generic renormalization

\item[leqt2f:]solves linear system

\item[lubksb:]used in linear system solution.  from \cite{Press}.

\item[ludcmp:]used in linear system solution.  from \cite{Press}.

\item[mprove:]used in linear system solution.  from \cite{Press}.

\item[nbinc:] finds continuum bin for given energy

\item[remtms:]fake ibm routine to calculate remaining cpu time in msec.

\item[spline:] Spline fit

\item[splinem:] Spline fit

\item[splint:] Spline fit

\item[pescl:] calculates line escape probability

\item[pescv:] calculates continuum escape probability

\item[starf:]calculates blackbody spectrum


\end{description}


\subsection {Units}  

An attempt has been made to retain the same units for 
a given physical quantity throughout the code.

\begin{description}
\item[Temperatures:] 10$^4$ K
\item[Distances:] cm
\item[Total Luminosities:] 10$^{38}$ erg/s
\item[Photon energies:] eV
\item[Continuum Emissivities (Specific):] erg cm$^{-3}$s$^{-1}$ erg$^{-1}$
\item[Continuum Luminosities (Specific):] 10$^{38}$ erg s$^{-1}$ erg$^{-1}$
\item[Line Emissivities:] erg cm$^{-3}$ s$^{-1}$
\item[Line Luminosities:] 10$^{38}$ erg s$^{-1}$
\item[Heating and Cooling rates:] erg cm$^{-3}$ s$^{-1}$
\item[Photoionization Rate Coefficients:] s$^{-1}$
\item[Collisional Rate Coefficients:] cm$^3$ s$^{-1}$
\end{description}

\subsection {List of Variable Names}

\begin{description}


\item[tp] is the radiation temperature, in units of KeV
\item[xlum] is the x-ray source luminosity, in units of 10**38 erg/s
\item[ecut] is the low energy cutoff of the input spectrum
\item[lpri] is the print switch ,1=on, 0=off.
\item[lwri] is the write switch. 
\item[nel] is the number of elements (max=13)
\item[nnnl] is the maximum  number of lines (1530)
\item[nni] is the number of ions, excluding fully stripped (max=168)
\item[ktitle] is the model title
\item[r] is the radius (cm)
\item[delr]=r-rl
\item[rstp]=delr/r
\item[rdel] is the distance from the current r to the first r of
        the calculation (cm)
\item[t] is the temperature (10**4 k)
\item[xee]: electron relative abundance, 
\item[xpx]: total particle density (/cm**3)
\item[xnx]: electron number density (/cm**3).
\item[abel]: the relative element abundances.
\item[xii]: the relative ion abundances.
\item[epi]: the continuum energy bins and the difference between
\item[bremsa]: the local ionizing flux ( erg/erg/cm**2/sec)
\item[zrems]: the specific luminosity.  (10**38 erg/erg/sec).
\item[zeta]: ionization parameter log(l/(n*r**2))  (
\item[opakc]: continuum opacities.  units are /cm.
\item[tau0]: line center optical depths from the cloud center.

\end{description}

\section {Flow Chart}

\begin{figure}
\epsfxsize=5.6in  % narrow the plot
\epsfysize=7.0in  % shorten the plot
\epsffile{xstarflow.eps}
\caption{The basic program flow of XSTAR.}
\label{fig:xstarflow}
\end{figure}

%Main Program:

% 1) Initialize (Define default parameter values, call collrd, sigmao, sigsum, init1
% 2) Get an Input Parameters (call rread1), calculate input spectrum
% 3) Read in atomic data
% 4) initialize the database (call dbwk)
% 5) Set up (call init)
% 6) Step through iterations
% A) Step through spatial zones
% a) Calculate step size 
%   if niter=1, then calculate (call step)
%            else use previous value (call unsavd)
% b) calculate ionizing flux (call trnfrc)
% c) Calculate state variables (call dsec)
%  i) calculate ionization, heating-cooling (call func)
%    I) calculate total ionization and recombination rates (call func1)
%    II) calculate ionization balance (call istruc)
%    III) calculate level rates (call func2)
%    III) calculate level populations (call hcor)
%    IV) calculate ionization  balance (call istruc)
%    V) calculate heating, cooling, emission, absorption (call func3)
%  ii) if heating=cooling done
%    else update temperature, go to (i)
% d) print out (call pprint)
% e) update radiation spectrum (call stpcut, call trnfrn)
% e) Done with spatial zones? n=go to (a)
% B) done with iterations? n=go to (A)
% 7) final printout



\section{XSTAR2XSPEC}

XSTAR2XSPEC consists of three major components (in addition to XSTAR 
itself).

\begin{description}
	\item[{\bf xstar2xspec:}] xstar2xspec is a Perl script which 
	manages the overall program flow of XSTAR2XSPEC.

	\item[{\bf xstinitable:}] xstinitable is an FTOOL used in the 
	initialization phase of XSTAR2XSPEC.  It builds the FITS PARAMETER 
	table from the input data (xstinitable.fits) and a text file 
	(xstinitable.lis) which is basically of list of all the calls to 
	XSTAR to generate the required spectra.

	\item[{\bf xstar2table:}] xstar2table is an FTOOL called after xstar in each 
	iteration of the xstar2xspec loop.
\end{description}

\subsection{xstar2xspec (script)}

This is a Perl script.

One current limitation is that the file tagging in the {\tt -save} 
option file names are limited to 9999 calls of XSTAR.

\begin{figure}
\epsfxsize=5.6in  % narrow the plot
\epsfysize=7.0in  % shorten the plot
\epsffile{xstar2xspec.eps}
\caption{The basic program flow of XSTAR2XSPEC.}
\label{fig:xstar2xspec}
\end{figure}

\subsection{xstinitable}

This FTOOL is written in C.  It generates an initial FITS file 
(xstinitable.fits) with 
the appropriate PRIMARY and PARAMETERS extension from the 
atables and mtable file are build.  It also generates a text file 
(xstinitable.lis) which contains a complete XSTAR calling command on 
each line.  This file is processed by the XSTAR2XSPEC Perl script

Examples of changes in XSTAR that would require changes in this FTOOL 
(just meant as a sample, not necessarily an exhaustive list).

\begin{enumerate}
	\item  Changing the number of physical parameters in XSTAR.

	\item  Changing or altering the control parameters in XSTAR.

\end{enumerate}

\subsection{xstar2table}

This FTOOL is also written in C.
Full use is made of dynamic memory allocation.  If the number of 
spectral channels is changed in XSTAR, this program should adapt 
appropriately, automatically.

Examples of changes in XSTAR that would require changes in this FTOOL 
(just meant as a sample, not necessarily an exhaustive list).

\begin{enumerate}
	\item  Changing the number of physical parameters in XSTAR.

	\item  At present, this program is written assuming that the XSTAR runs are 
performed sequentially and in a particular order.  In fact, it checks the 
sequence by comparing the LASTSPEC keyword to the LOOPCONTROL variable.
If, at some time in the future, the XPI interface is modified in such 
a way that it becomes possible to simultaneously submit multiple XSTAR 
runs on multiple processors, xstar2table must be modified to ignore 
this comparison.

\end{enumerate}
