The Input and Output of STARS

This page is almost finished, just has a couple of gaps.

Here, the input model and the output files of the code are discussed.

InputOutput
DATAMODOUT
MODINOUT
PLOT

DATA (fort.20)

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
NMSH2The 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)
ITER1The maximum number of iterations allowed on the first timestep. (10)
ITER2The maximum number of iterations allowed on later timesteps. (10)
JINThe number of independent variables of the H, DH arrays to be read in. (11)
JOUTThe number of independent variables written out at the end of a run. (11)
NCHThis 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)
JPThis 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)
ITHThe thermal energy generation rate is ITH*T*DS/Dt. Setting ITH to zero enables T*DS/Dt to be ignored. (1, sometimes 0)
IXDX/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)
IYas for IX, but relates to helium-4 instead. (1, sometimes 0)
IZas for IX, but relates to carbon-12 and oxygen-16. (1, sometimes 0)
ICLIf ICL=1, then PRESSI will include the effects of the Coloumb interaction. (1)
IONUsed by STATEF. This is the number of elements that that the routine will calculate ionizations for. (2)
IAMIf set to zero, the program will use integer atomic weights. (1)
IOPUsed 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)
INUCNo current function.
IBCNo current function.
ICNUsed 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)
NWRT1Prints the internal details of every NWRT1'th model to OUT. (e.g. 100)
NWRT2Prints the internal of every NWRT2'th meshpoint when the internal model is printed to OUT. (e.g. 1 or 2)
NWRT3Affects the number of 'pages' printed out for every NWRT1'th model. (1-3)
NWRT4Prints a short summary of every NWRT4'th model. (e.g. 1 or 4)
NWRT5Prints a one-line summary of each iteration of each model, excluding the first NWRT5 iterations of each model. (e.g. 0 or 2)
NSAVEAn 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)
EPSThe accuracy to which SOLVER is expected to solve the equations. (10-6 or 10-7)
DELThis 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)
DH0Affects 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)
DT3No current function.
DDDSets the modulus of the total increment that is desired in one timestep. (0.5 to 4)
NE1The number of 1st order equations the code will use. (5)
NE2The number of 2nd order equations the code will use. (5)
NE3This 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)
NBThe number of boundary conditions possessed by the 1st order equations at the stellar surface. (3)
NEVThe 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)
NFThis affects the number of variables being passed between the routines used by SOLVER. (30)
J1, J2, IH, JHThese 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:
  1. &Psi
  2. P
  3. &rho
  4. T
  5. &kappa
  6. a
  7. &nabla
  8. &nablar-&nablaa
  9. m
  10. 1H
  1. 4He
  2. 12C
  3. 14N
  4. 160
  5. 20Ne
  6. 24Mg
  7. r
  8. L
  9. &epsilonth
  10. &epsilonnuc
  1. &epsilon&nu
  2. &delta m
  3. k2
  4. n/n+1
  5. Uhom
  6. Vhom
  7. U
  8. S
  9. L/LEdd
  10. w.l
    23 is the square of the radius of gyration; 24-26 are homology invariants; 30 is the convective velocity times the mixing length.
    DT1Places a lower limit of DT1 * current timestep on the size of the next timestep. (0.8. 0.95 or 1.0)
    DT2Places 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.
    ZSThis 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)
    ALPHAThe mixing length. (2.00)
    CH-CFEThese 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)
    RCDThe diffusion coefficient for convective mixing is RCD*(&nablar-&nablaa)2/tnuc. (106)
    OSA convective overshoot parameter. Zero implies no overshoot. (0)
    RML, RMG and RMTAt 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)
    ECAUsed in the evolution of EC for pre-MS construction. (0)
    XFDefines the boundary of a core, for printout purposes only, to be at X(1H) or X(4He) = XF
    DRDefines the boundary between a convection and semiconvection zone, for printout purposes only, to be at &nablar-&nablaa = DR. (0.01)
    RHLThe 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)
    ACIf there is Roche-lobe overflow, the stellar wind is a fraction 1-AC of the mass being lost by star 1. (1.0)
    AK1The 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)
    AK2The 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)
    ECTA 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)
    TRBUsed in setting the surface conditions of binaries.

    MODIN (fort.25)

    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

    SMis the mass, in solar units of the star being evolved.
    DTYis the inital size of the timestep that will be used.
    AGEis 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.
    PERis 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).
    BMSis 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).
    ECis 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.
    NHis the number of mesh points in the current model.
    NPis the number of models that is desired for the current run.
    NMODis 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).
    IBis 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:

    1. ln f - a quantity closely related to electron degeneracy
    2. ln T - temperature in Kelvin
    3. X16 - the abundance of 16O
    4. ln m - mass in units of 1030kg
    5. X1 - the abundance of 1H
    6. C - the gradient of the mesh spacing function, Q(P,T,m,r). It does not depend on K, though it does vary with time. It is in effect an eigenvalue.
    7. ln r - the radius in units of 109m
    8. L - luminosity in units of 1026W
    9. X4 - the abundance of 4He
    10. X12 - the abundance of 12C
    11. X20 - the abundance of 20Ne

    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.

    MODOUT (fort.23)

    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.

    OUT (fort.22)

    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.

    PLOT (fort.10)

    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:

    1. Model number
    2. Age of model
    3. log-10 of the radius in solar units
    4. log-10 of the surface temperature
    5. log-10 of the surface luminosity in solar units
    6. mass in solar units
    7. location in mass of the boundary of the hydrogen burning shell
    8. location in mass of the boundary of the helium burning shell
    9. log-10 of the helium luminosity (or 1.01x10-10, whichever is larger)
    10. location in mass of the convective boundaries
    11. location in mass of point of maximum energy generation from H-burning
    12. location in mass of point of maximum energy generation from He-burning
    13. log-10 of the opacity
    14. timestep in years
    15. surface abundance of hydrogen
    16. surface abundance of helium
    17. surface abundance of carbon
    18. surface abundance of nitrogen
    19. surface abundance of oxygen
    20. mass location of burning shell boundaries
    21. d(log r)/d(log m)

    E-mail: stars@ast.cam.ac.uk
    Last modified: Mon May 31 11:42:47 2004
    by R. Stancliffe