Here is a brief overview of the 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.
The program consists of a main routine, MAIN, and 26 subroutines. It requires input from the files data (fort.1 -- which essentially tells the code how to evolve the model), modin and nucmodin (fort.30 and fort.31 -- the input structure model and nucleosynthesis model). A second pair of input files, modin2 and nucmodin2 (fort.60 and fort.61) can be used if you wish to evolve a binary system. It generates three output files: modout, out and plot (fort.32, fort.33 and fort.34) of running in single star mode without nucleosynthesis. Additional output files are produced if you evolve a binary and/or include nucleosynthesis. There is one more possible output (fort.24), used for debugging (though extremely rarely!). 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. This routine is also responsible for reseting the composition of a model (NCH=3), as well as resetting the angular momentum of the orbit of a binary and the spin angular momenta of the two stars.
This evaluates the functions of the independent variables at each mesh point that are required for the difference equations. It contains most of the large-scale physics used by the code, with STATEF dealing with the equation of state of the stellar material. 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 its input from the calculations of FUNCS1.
As FUNCS1, but only deals with minor compositon variables. This routine is used, together with EQUNS2, to compute the nucleosynthesis of 40 energetically unimportant isotopes. The functions are analytically differentiated, which saves computer time, but is harder to alter. Temperatures and densities (plus a few other important quantities) as computed and stored during the pass through FUNCS1 are used here to save us having to call STATEF again.
As EQUNS1, but relating to FUNCS2. Again, it involves analytic integration.
This routine computes the diffusion coefficients need if you wish to evolve a model with gravitational settling. Currently this is written assuming that hydrogen is the dominant element. EQUNS1 also assumes this. EQUNS2 has not yet been adapted to include this physics.
Works out the mass loss from the stars in the binary you are evolving. This routine also works out any stellar interactions, such as tides modifying the orbit, mass transfer by Roche Lobe overflow or wind accretion, etc.
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. The viscous mesh control is in this routine, as is the AGB timestep control.
This routine will set compositions (and angular momenta) 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 and energy generation for use by FUNCS1.
Evaluates nuclear reaction rates for the isotopes in the nucleosynthesis network, for use in FUNCS2.
Calculates the nucleosynthesis associated with neutron captures.
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. This includes the neutrino loss rates and the opacities.
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.
Computes Fermi-Dirac integrals need for computing the equation of state.
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.
A more recent opacity routine, which includes opacities which vary with the carbon and oxygen content of the stellar material.
Used to create the spline interpolation data. Uses a bicubic spline.
A routine which initializes the physical constants used by the code.