Here is a brief overview of the STARS code. It is based on the comments in the code, P.P. Eggleton's instructions (July 1995) and my own understanding from working through it. If I seem to have given an unfair emphasis to REMESH and EQUNS1, it is because these are the areas I have primarily been involved with. Any errors are of my own making.
The program consists of a main routine, MAIN, and 20 subroutines. It requires input from two files, modin and data (referred to as fort.25 and fort.20 in the code), and it generates three output files: modout, out and plot (fort.23, fort.22 and fort.10 in the code). There is one more possible output (fort.24), used for debugging. A schematic of the program is shown below, and the routines and their functions are listed:
This is the main routine, controlling the number of timesteps the model is evolved for. It calls PRINTA to obtain the input details, and send out the newly calculated model. SOLVER is called to solve the simultaneous difference equations for the independent variables by a Newton-Raphson method. It will alter the timestep if the model fails to converge.
The input/output workhorse. It reads in model information and evolutionary parameters from MODIN and DATA, and between timesteps, updates the independent variables. It also chooses the length of the next timestep, though this is not always done effectively. Calls to COMPOS, PRINTB and REMESH (if desired) are made from this subroutine.
This routine spaces the meshpoints out at equal intervals of the mesh spacing function, if NCH is greater than or equal to 2. This is acheived my calculating a gradient from the first and last mesh points, and redistributing them as necessary. I am engaged in modifying this routine so that the points may be locked into place during rapid evolutionary phases.
This evaluates the functions of the independent variables at each mesh point that are required for the difference equations. It also holds the equation for the mesh spacing function.
This sets up the boundary conditions and the difference equations, which may be a mix of first and second order. It takes it's input from the calculations of FUNCS1.
As FUNCS1, but only deals with minor compositon variables. The functions are analytically differentiated, which saves computer time, but is harder to alter.
As EQUNS1, but relating to FUNCS2. Again, it involves analytic integration.
This handles most of the output, giving summaries or details of the stellar models. It also writes an output file that can be used for plotting. The rountine will also locate the boundaries between different zones in the model. It is also the place to evaluate characteristics (e.g. moment of inertia, gravitational binding energy) of a model that has just been calculated.
This routine will set compositions to zero if they get lower than 1e-12. It is supposed to stop the compositions from being given negative values.
Evaluates nuclear reaction rates, neutrino loss rates, etc. for use by the FUNCSs and EQUNSs.
Acts as a storage buffer between FUNCS1 and STATEF. This saves recomputing the equation of state when only a variable that doesn't affect the state has been changed since the last call to STATEF.
Routine for calculating the equation of state. It calculates all the required quantities (save the reaction rates) which are functions of only the state of the gas. It calls both PRESSI and OPACTY.
Computes the effects of the free energy contribution on the equation of state. It also contains an approximation to the Coloumb interaction. It evaluates changes to ionization potential, pressure and entropy and their derivatives.
Along with DIFRNS, ELIMN8 and DIVIDE, this is a general solution package for systems of simultaneous difference equations (some first and some second order), with appropriate boundary conditions at each end. It does not depend on the particular character of the stellar structure equations. Its purpose is to acheive, by Newton-Raphson iteration, a solution to the equations evaluated in the EQUNSs routines. SOLVER may write to OUT, giving information on the meshpoints where it had trouble acheiving the desired accuracy.
Sets up the required difference equations. It evaluates derivatives of these equations numerically, by varying each independent variable in turn and passes the derivatives back to SOLVER for implementation of the scheme for solving the difference equations.
Does some necessary matrix multiplication needed for successive elimination.
A custom-built matrix inverter.
Both DIVIDE and ELIMN8 call it. It writes an output file that may be used in debugging, if the user has sufficient grasp of the code.
Calculates a bicubic spline interpolation fit for the temperature and density opacity fit.
Used to create the spline interpolation data. Uses a bicubic spline.
A routine which initializes the physical constants used by the code.