This page is almost finished, just has a couple of gaps.
Here, the input model and the output files of the code are discussed.
This is file that effects how the model will be run. It consists of a block of 221 numbers, some of which relate to the model's evolution, while the others affect the format of the output. The structure of DATA is given below. I have tried to keep the variable names the same as in the code, at the expense of aesthetics in displaying the file's structure. The function of each variable is discussed below, with bracketed quantities representing the usual values they take.
NMSH2 ITER1 ITER2 JIN JOUT NCH JP ITH IX IY IZ ICL ION IAM IOP INUC IBC ICN NWRT1 NWRT2 NWRT3 NWRT4 NWRT5 NSAVE EPS DEL DH0 DT3 DDD ID(1-110) - 7 lines ISX(1-45) - 3 lines DT1 DT2 CT(1-10) ZS ALPHA CH CC CN CO CNE CMG CSI CFE RCD OS RML RMG ECA XF DR RMT RHL AC AK1 AK2 ECT TRB
|NMSH2||The desired number of mesh points. If this is different from that in modin, the code should interpolate the given model to give the new one, provided that NCH is greater than or equal to 1. (199)|
|ITER1||The maximum number of iterations allowed on the first timestep. (10)|
|ITER2||The maximum number of iterations allowed on later timesteps. (10)|
|JIN||The number of independent variables of the H, DH arrays to be read in. (11)|
|JOUT||The number of independent variables written out at the end of a run. (11)|
|NCH||This affects the remeshing of the model. For an NCH of 2 or 3, the model will alter the spacing of the mesh points so that they are optimally spaced according to the mesh spacing function. NCH=3 is used for changing the composition of a model. (1,2 or 3)|
|JP||This affects the guessed increments for DH for the next timestep. If it is set to 1, the values used are those of the previous timestep. If it is zero, then the increments are set to zero. (0, sometimes 1)|
|ITH||The thermal energy generation rate is ITH*T*DS/Dt. Setting ITH to zero enables T*DS/Dt to be ignored. (1, sometimes 0)|
|IX||DX/Dt = IX * (burning rate of Hydrogen). Setting this equal to 0 allows the model to ignore composition changes whilst still burning hydrogen. This is useful for creating intial models. (1,sometimes 0)|
|IY||as for IX, but relates to helium-4 instead. (1, sometimes 0)|
|IZ||as for IX, but relates to carbon-12 and oxygen-16. (1, sometimes 0)|
|ICL||If ICL=1, then PRESSI will include the effects of the Coloumb interaction. (1)|
|ION||Used by STATEF. This is the number of elements that that the routine will calculate ionizations for. (2)|
|IAM||If set to zero, the program will use integer atomic weights. (1)|
|IOP||Used to select spline interpolation in opacity. If set to 1, spline interpolation will be used. If set to zero, a linear interpolation in opacity will be used. (1)|
|INUC||No current function.|
|IBC||No current function.|
|ICN||Used for CNO equilibrium on the main sequence. If ICN=1, a baryon correction is applied to the hydrogen evolution equation. The timestep used is also affected by ICN=1. (0)|
|NWRT1||Prints the internal details of every NWRT1'th model to OUT. (e.g. 100)|
|NWRT2||Prints the internal of every NWRT2'th meshpoint when the internal model is printed to OUT. (e.g. 1 or 2)|
|NWRT3||Affects the number of 'pages' printed out for every NWRT1'th model. (1-3)|
|NWRT4||Prints a short summary of every NWRT4'th model. (e.g. 1 or 4)|
|NWRT5||Prints a one-line summary of each iteration of each model, excluding the first NWRT5 iterations of each model. (e.g. 0 or 2)|
|NSAVE||An output model is saved to MODOUT every NSAVE'th timestep, in the same format as the input model, so that it can be used for a further run. The final model is automatically saved. (300)|
|EPS||The accuracy to which SOLVER is expected to solve the equations. (10-6 or 10-7)|
|DEL||This is the maximum value in ERR for which the whole correction is applied in SOLVER. Above this limit, the correction applied is reduced by a factor of ERR/DEL. (0.01)|
|DH0||Affects the value of the increments of the variables during the numeric differentiation. It is no longer so important, now that everything is in double precision. (10-7)|
|DT3||No current function.|
|DDD||Sets the modulus of the total increment that is desired in one timestep. (0.5 to 4)|
|NE1||The number of 1st order equations the code will use. (5)|
|NE2||The number of 2nd order equations the code will use. (5)|
|NE3||This defines a subset of the 1st order equations that may be defined at 3, rather than two, adjacent meshpoints. At present, this is not implemented in the code. (0)|
|NB||The number of boundary conditions possessed by the 1st order equations at the stellar surface. (3)|
|NEV||The number of 'eigenvalues' (i.e. quantites that don't vary with the mesh) used by the model. At present there is only one, that being the gradient of the mesh spacing function. (1)|
|NF||This affects the number of variables being passed between the routines used by SOLVER. (30)|
|J1, J2, IH, JH||These variables are used for debugging. Suitable choice gives an output via PRINTC that can be used to see if FUNCS and EQUNS are setting up the difference equations correctly. (0,0,0,99 - this suppresses the debugging output.)|
|ID(1-110)||This is two blocks of numbers, one for FUNCS1/EQUNS1, the other for FUNCS2/EQUNS2. Each block consists of 1 lines of ten, and 3 lines of 15 integers. The line of 10 contains NE1,NE2,NB,NEV,NF,J1,J2,IH,JH (ie ID(1-10) is the previous line of the data file). The usual values are (5,5,0,3,1,30,0,0,0,99). Of the 3 lines of 15 integers, the first is a permutation of the independent variables (1,2,4,5,3,9,10,8,7,6,0,0,0,0,0), the second that of the the equations (6,7,8,9,4,2,1,3,5,0,0,0,0,0,0) and the last that of the boundary conditions (4,5,6,7,2,3,1,2,3,1,0,0,0,0,0).|
|ISX(1-45)||These values affect the variables of the internal structure that are printed to OUT. The first 15 values define what will be placed on the first 'page', with the next two sets of 15 defining the output to the extra pages. The options are:|
|23 is the square of the radius of gyration; 24-26 are homology invariants; 30 is the convective velocity times the mixing length.|
|DT1||Places a lower limit of DT1 * current timestep on the size of the next timestep. (0.8. 0.95 or 1.0)|
|DT2||Places an upper limit of DT2 * current timestep on the size of the next timestep. If both DT1 and DT2 are set to 1, then the timestep is constant, unless the model fails to converge (in which case it will become 0.8 times the current value for the next run). (1.2, 1.05 or 1)|
|CT(1-10)||These are coefficients used in the mesh spacing function, Q.|
|ZS||This is the star's metallicity. See the recipes section if you wish to change the metallicity of the model you currently have. (e.g. 0.02, 0.001)|
|ALPHA||The mixing length. (2.00)|
|CH-CFE||These are the values for initalising the abundances of hydrogen-1 through to iron-56 of a model. The metals are expressed as a fraction of the total metallicity. These are only used for ZAMS models, and with NCH set to 3. (0.7 or 0.73, 0.176, 0.502, 0.092, 0.034, 0.072, 0.072)|
|RCD||The diffusion coefficient for convective mixing is RCD*()2/tnuc. (106)|
|OS||A convective overshoot parameter. Zero implies no overshoot. (0)|
|RML, RMG and RMT||At the surface, the mass loss rate is equal to -RML*Lr/m + RMG*m - RMT*[rstar/rlobe-1]3. This enables inclusion of a Reimers' type mass-loss rate, a constant mass-loss/gain rate (used mainly for running up and down the main sequence - see the recipes), and/or Roche-lobe overflow. All the units are solar and years. (0 or 4x10-13, 0 or +/- 5x10-7, 103 respectively)|
|ECA||Used in the evolution of EC for pre-MS construction. (0)|
|XF||Defines the boundary of a core, for printout purposes only, to be at X(1H) or X(4He) = XF|
|DR||Defines the boundary between a convection and semiconvection zone, for printout purposes only, to be at = DR. (0.01)|
|RHL||The intrinsic rate of loss of angular momentum from the binary, excluding that from stellar wind by either star. The total angular momentum loss is a combination of the intrinsic loss (here just a constant logarithmic rate) combined with the wind from both stars. (0)|
|AC||If there is Roche-lobe overflow, the stellar wind is a fraction 1-AC of the mass being lost by star 1. (1.0)|
|AK1||The loss of angular momentum from star 1 as wind is AK1*the orbital angular momentum per unit mass as star 1. (0.0 or 1.0)|
|AK2||The loss of angular momentum from star 2 as wind is AK2*the orbital angular momentum per unit mass as star 2.(1.0 or 0.0)|
|ECT||A constant logarithmic increase or decrease for EC (see the section on MODIN), which can be used to push a ZAMS star back up it's Hayashi track to make a pre-MS star.(0)|
|TRB||Used in setting the surface conditions of binaries.|
This file contains the input model that the code will start to evolve when run. The first line of the code contains the following parameters:
SM DTY AGE PER BMS EC NH NP NMOD IB
|SM||is the mass, in solar units of the star being evolved.|
|DTY||is the inital size of the timestep that will be used.|
|AGE||is the current age of the model. Note that if you take a model from a previous run, you may wish to reset the age to zero. It has no effect on the calculations.|
|PER||is the period of the binary system in days. The code treats all models as if they are part of a binary system. If you only wish to evolve a single star, then set this variable to something large (1020 is my usual choice).|
|BMS||is the total mass of the binary; it must always be larger than the mass of the object star (2000 is a convenient choice for single star evolution).|
|EC||is the rate of energy generation used for evolving back up the Hayashi track to create pre-MS stars. For normal evolution, it is set to zero.|
|NH||is the number of mesh points in the current model.|
|NP||is the number of models that is desired for the current run.|
|NMOD||is the model number of the current model. If you are starting a new run, you may wish to set it to zero (this is not important, as it doesn't affect the calculation).|
|IB||is the number of the star in a binary system (i.e. it's either 1 or 2). It affects the boundary conditions being used at the surface - for stars losing mass, IB=1 and for stars accreting, IB=2. For single star evolution it is set to 1.|
The remaining numbers in this file are the actual details of the model. There are 11 independent variables recorded, and there should be NH rows of numbers. The variables, in the order of the columns that appear, are:
After these NH rows, come another NH rows. These are the changes made to the variables at the last timestep. In some cases, they will all be zero. It is important to copy all the numbers from the MODOUT file if you wish to start from exactly the same point as the last evolutionary run.
Output models from the code are written to here. The frequency at which this occurs is controlled by the value of NSAVE in DATA. The output format is the same as described above for MODIN. Any of the models written here can be placed into MODIN to be used in a subsequent run.
This file is useful as a reference for what is going on as a star evolves. The file starts with a copy of the data block used for the star's evolution. Below this is a print out of the first line of MODIN, so you know exactly where the model started from. The structure of the remaining output depends on the values that were chosen in DATA. Typically, the code will produce a short four line summary of the model (the frequency depends on the value of NWRT4), like the one below:
32 dt/age/MH/MHe tn/tKH/Mb P/rlf/dM LH/LHe/LC Lth/Lnu/m H1 He4 C12 N14 O16 Ne20 Mg24 psi dens temp 10.0000 2.942759062E+02 1.691E+07 1.006E+49 3.291E+03 2.833E+03 0.69942 0.28005 0.00027 0.00482 0.00963 0.00198 0.00076 -5.334 0.8548 7.4603 cntr 0.0000 3.441230970E+04 4.409E+04-7.522E+01 2.878E-22 2.029E+02 0.70000 0.28000 0.00346 0.00108 0.00964 0.00198 0.00076 -17.347 -8.9425 4.3838 srfc 0.0000 6.9984 2.8001 2000.0000 0.000E+00 0.000E+00 0.0001 0.69942 0.28005 0.00027 0.00482 0.00963 0.00198 0.00076 -5.334 0.8548 7.4603 Tmax
The first column contains the number of the model (in this case 579), then the mass of the star being evolved, in solar units (1.9953). The next number is the location in mass of the boundary of the hydrogen burning shell, and the last is the location in mass of the boundary of the He-burning zone. In the second column, we have the current timestep of the model and below it the current age of the model - both in years. The two numbers below these are the masses of H and He in the star, in solar units. The next column gives the nuclear and Kelvin-Helmholtz timescales in years, and the mass of the binary system in solar units. In column four, we find the period of the binary system in days, the difference between the stellar and Roche Lobe radii (in units of the stellar radius) and the rate of mass transfer. Next come two columns giving information about the luminosity of the system. First up are the luminosities from hydrogen, helium and carbon (actually, any nuclear burning that isn't H or He!), followed by the thermal and nuclear luminosities. The units are ? and the m at the base of this column is the mass location of the region of maximum temperature. The remaining entries give the abundances of nuclear species at the three points noted in the far right column.
After this block comes another 4 lines proving information about the model. It will take a form something like:
0.0000 -2.102 2.323 10.000 -10.000 10.000 -10.000 10.000 -10.000 -10.000 10.000 0.000 0.000 1.1233 0.0665 0.0197 0.6503 3.7891 32 199 199 -161 160 44 -44 39 -39 15 -15 -9 8 0 0 0.0001 0.0000 0.0000 0.0644 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.659E+04 4.743E-21 0.000E+00 199 0 0 186 0 0 0 0 0 0 0 0 0 0 0-Inf
On the first line, the first number is the mass interior to the point at which carbon burning occurs. The next twelve are the masses interior to the points at which the model changes from convective regions to radiative regions (and vice versa). The next 3 numbers are:(??). Then comes the base-10 log of the radius in solar units, followed by the base-10 log of the luminosity in solar units. The final number is the number of the model.
On the second line, the first number is the mesh point marking the H-shell burning boundary. The next number is the same quantity, but for the He-burning region. The last 12 numbers of the second line give information on where the boundaries between convective and radiative regions are. Boundaries between convective and semi-convective regions are marked by negative signs; convective/radiative boundaries have no sign.
On the third line, the first three numbers are the mass-locations of the points of maximum hydrogen, helium and carbon energy generation. The next 12 numbers are the mass-locations of the boundaries of burning regions. The final three numbers on this lines are the values of the energy generation rates for hydrogen, helium and carbon respectively, at their maximum points.
On the last line, the first three numbers are the the mesh point numbers of the points of maximum hydrogen, helium and carbon energy generation. The next twelve numbers are the mesh point numbers of the location of burning shell boundaries. The last value is d(log r)/d(log m).
Every NWRT1th model, the code will output details of the model. The printed details are selected via ISX in DATA.
The code writes select details of each model to this file, so that such things like evolutionary tracks can be plotted. The file contains the following information, in this order: