Toy model implementation of SE

Piet Hut and Jun Makino, June 30, 2002
Fortran implementation added July 1, 2002
Some bug fixes for Fortran implementation July 2, 2002

Table of contents

Description

One of the goals of our workshop was to specify ways to let existing computer codes for stellar dynamics (SD) and for stellar evolution (SE) communicate with minimum modification at either end. In fact, with a well-defined minimal interface between the two, SD should not care whether the information obtained results from running an actual SE code, or from a table look-up or fitting formula. The same holds for SE: each side should see the other as a black box. Soon, we will also discuss the role of stellar hydrodynamics (SH), but we'll leave that out for now, to get started.

We have just written a toy model version for the SD-SE interface, following the discussions during the workshop. In order to test our interface, we constructed a very simple implementation of both the SD and SE parts of a simulation code. For SD we simply took two unbound stars on a head-on collision course, to let them merge into a single star with an unusual composition. For SE we constructed the following stick-figure version. We would be grateful to get reactions and suggestions from those of you who have more experience with SE codes.

We approximate each physical quantity Q that describes the state of a star (mass, radius, Y, Z) with a piece-wise linear function with one discontinuity, specified by the following six numbers:

      t_endms   time at which the star leaves the main sequence
      t_endrg   time at which the star leaves the red giant branch.
      f_init    value of Q at birth of star (at time t=0)
      f_endms   value of Q at time t=t_endms
      f_endrg   value of Q just before time t=t_endrg
      f_remnant value of Q just after time t=t_endrg

                                         ....... f_endrg
                                        /|
                                       / |
                                      /  |
                                     /   |
     ^                              /    |
     |                             /     |
     Q                            /      |
                            _____/.......|...... f_endms
                  _____-----             |
        _____-----.......................|...... f_init
                                         |------ f_remnant
        +-----------------------+--------+------ 0
        0                       t_endms  t_endrg
                   t -->
Note that we stick-figure version of stellar evolution mimics only low-mass stars, while leaving out completely the horizontal branch part of the evolution. In addition, we will make the following simplified assumption: t_endrg = 1.1 t_endms. This leaves us with 17 independent values: one time scale and four Q values for each of the four Q choices. Here is what we choose, starting with the three free values mass M, Y_init Y0, Z_init Z0:

    
    t_endms = ( 10^4 / M^3 ) Myr

       Q (unit)     |  Q_init  |  Q_endms  |  Q_endrg  |  Q_remnant
    ----------------+----------+-----------+-----------+-------------
    mass (Msun)     |    M     |     M     |   M / 2   |    M / 4
    radius (Rsun)   |    M     |     M     |   100 M   |  0.01 / M
    He fraction Y   |    Y0    |     Y1    |     Y2    |     Y3
    metallicity Z   |    Z0    |     Z1    |     Z2    |     Z3

    where Y1 = minimum(Y0 + 0.1 , 1 - Z1) 
          Y2 = minimum(Y0 + 0.4 , 1 - Z2) 
          Y3 = minimum(Y0 + 0.8 , 1 - Z3) 
          Z1 = Z0
          Z2 = minimum(Z0 + 0.2 , 1) 
          Z3 = minimum(Z0 + 0.4 , 1) 

The idea is to let the star lose half of its mass on the red giant branch due to wind loss; and then again half of its mass in the instantaneous transition to a white dwarf. For the helium fraction Y we presume that it increases by 0.1 while on the main sequence, and by another amount 0.1 on the red giant branch; we then double that number to get Y = Y0 + 0.4, to account for the mass lost (the newly formed helium is not lost since it is in the inner regions). After losing again half the mass, the final Y is concentrated to Y = Y0 + 0.8, but of course in all cases the total Y+Z cannot exceed unity, hence the "minimum(,)" operators above. For the metallicity Z we assume that nothing is generated on the main sequence, and 0.1 on the red giant branch.
Does all this seem reasonable? Are there obvious improvements that we could make, without making the whole toy model much more complex? Is anybody aware of toy models resembling this in the literature?

As for the interface specification, we will comment on that later. Below is a sample output of the output from the above model, starting with a time step of 1 Myr, and letting the time step double each time, unless the changes in mass, radius, Y, or Z became too large (this is something that SD could limit through the interface to SE, as we discussed at the workshop). We imposed the following requirements:

  | delta mass |    < 0.1 * mass
  | delta radius |  < radius
  | delta Y |       < 0.1
  | delta Z |       < 0.1

Sample output from Fortran code

Output for a one-solar-mass star:

Age =1.000000000     M = 1.000000000     R = 1.000000000     Y = .2500100000     Z = .2000000000E-01
Age =3.000000000     M = 1.000000000     R = 1.000000000     Y = .2500300000     Z = .2000000000E-01
Age =7.000000000     M = 1.000000000     R = 1.000000000     Y = .2500700000     Z = .2000000000E-01
Age =15.00000000     M = 1.000000000     R = 1.000000000     Y = .2501500000     Z = .2000000000E-01
Age =31.00000000     M = 1.000000000     R = 1.000000000     Y = .2503100000     Z = .2000000000E-01
Age =63.00000000     M = 1.000000000     R = 1.000000000     Y = .2506300000     Z = .2000000000E-01
Age =127.0000000     M = 1.000000000     R = 1.000000000     Y = .2512700000     Z = .2000000000E-01
Age =255.0000000     M = 1.000000000     R = 1.000000000     Y = .2525500000     Z = .2000000000E-01
Age =511.0000000     M = 1.000000000     R = 1.000000000     Y = .2551100000     Z = .2000000000E-01
Age =1023.000000     M = 1.000000000     R = 1.000000000     Y = .2602300000     Z = .2000000000E-01
Age =2047.000000     M = 1.000000000     R = 1.000000000     Y = .2704700000     Z = .2000000000E-01
Age =4095.000000     M = 1.000000000     R = 1.000000000     Y = .2909500000     Z = .2000000000E-01
Age =8191.000000     M = 1.000000000     R = 1.000000000     Y = .3319100000     Z = .2000000000E-01
Age =10010.10098     M = .9949495117     R = 1.999996680     Y = .3530302930     Z = .2202019531E-01
Age =10030.30293     M = .9848485352     R = 3.999990040     Y = .3590908789     Z = .2606058594E-01
Age =10070.70684     M = .9646465820     R = 7.999976758     Y = .3712120508     Z = .3414136719E-01
Age =10151.51465     M = .9242426758     R = 15.99995020     Y = .3954543945     Z = .5030292969E-01
Age =10313.13027     M = .8434348633     R = 31.99989707     Y = .4439390820     Z = .8262605469E-01
Age =10481.81719     M = .7590914062     R = 48.69990156     Y = .4945451563     Z = .1163634375    
Age =10633.63545     M = .6831822754     R = 63.72990947     Y = .5400906348     Z = .1467270898    
Age =10770.27188     M = .6148640625     R = 77.25691563     Y = .5810815625     Z = .1740543750    
Age =10893.24463     M = .5533776855     R = 89.43121826     Y = .6179733887     Z = .1986489258    
Age =10999.99990     M = .5000000488     R = 99.99999033     Y = .6499999707     Z = .2199999805    
Age =11000.00010     M = .2500000000     R = .1000000000E-01 Y = .5800000000     Z = .4200000000    
Age =19192.00010     M = .2500000000     R = .1000000000E-01 Y = .5800000000     Z = .4200000000    

Note the at first the time step doubling was the limiting factor. Then during the first five steps on the red giant branch, the time step was constrained by the requirement to at most double the radius. During the next five steps the constraint was to allow a change in mass of not more than ten percent. The final step was again constrained by time step doubling, while the two before were forced to occur at both sides of the discontinuity at the end of the life of the star.

Driver source code in Fortran

Here is the fortran source code for testing the toy model implementation described above, with the header file for interface functions. This is not meant for distribution; it is still very rough, quickly cobbled together to get some results; but if you're interested to look under the hood at what we've done, here it is.

      subroutine print_star(idstar)
      integer idstar
#include "modest_star.h"
      write(6,600) getTime(idstar),getMass(idstar),getRadius(idstar),
     $     getY(idstar),getZ(idstar)
 600  format('Age =',g15.10,' M = ',g15.10,' R = ',g15.10,
     $     ' Y = ',g15.10,' Z = ',g15.10)
      end

      program one_star
#include "modest_star.h"      
      integer star
      real*8 dtmax, t
      integer discon_reached, countdown
      real*8 dMmax, dRmax, dYmax, dZmax, new_time
      star =  CreateStar(1.0d0,0.25d0,0.02d0)
      dtmax = 1d0
      discon_reached = 0
      t = 0
      countdown = 3
 100  if (countdown .gt. 0) then
         dMmax = 0.1d0*getMass(star)
         dRmax = 1d0*getRadius(star)
         dYmax = 0.1d0
         dZmax = 0.1d0
         new_time = EvolveStar(star,dtmax,dMmax,dRmax,dYmax,dZmax)
         call print_star(star)
         if ((new_time -t).gt. dtmax - 1d-10) dtmax = dtmax*2
         t = new_time
         if (get_discontinuity_flag(star) .ne. 0) discon_reached = 1
         if (discon_reached .ne. 0 ) then
            countdown  = countdown - 1
         endif
         goto 100
      endif
      end

      
        


modest_star.h
      real*8 getMass, getRadius, getY, getZ, EvolveStar
      real*8 getTime
      integer CreateStar, get_discontinuity_flag

Toy model source code in Fortran

Here is the Fortran source code for toy model implementation along with the include file for the common block. This too is not meant for distribution; it is still very rough, quickly cobbled together to get some results; but if you're interested to look under the hood at what we've done, here it is.

C------------------------------------------------------------
C        MODEST_STAR.F
C
C       Piet Hut and Jun Makino
C       Version 0.0 July 1, 2001
C------------------------------------------------------------
C
C  Each physical quantity Q that describes the state of a star
C  can be encapsulated as an object of class `stellar_quantity',
C  where Q can stand for mass, radius, etc.  Such an object
C  contains the full information about the evolution of Q over
C  the life time of a star.  The value of a quantity Q at time t
C  can be obtained by invoking its member function Q.at(t) .
C 
C  In our current toy-model implementation, Q(t) is described by
C  a piece-wise linear function with one discontinuity, specified
C  by the following six numbers:
C 
C    t_endms   time at which the star leaves the main sequence
C    t_endrg   time at which the star leaves the red giant branch.
C    f_init    value of Q at birth of star (at time t=0)
C    f_endms   value of Q at time t=t_endms
C    f_endrg   value of Q just before time t=t_endrg
C    f_remnant value of Q just after time t=t_endrg
C 
C                                       ....... f_endrg
C                                      /|
C                                     / |
C                                    /  |
C                                   /   |
C   ^                              /    |
C   |                             /     |
C   Q                            /      |
C                          _____/.......|...... f_endms
C                _____-----             |
C      _____-----.......................|...... f_init
C                                       |------ f_remnant
C      +-----------------------+--------+------ 0
C      0                       t_endms  t_endrg
C                 t -->
C 
C  Note that we stick-figure version of stellar evolution mimics
C  only low-mass stars, while leaving out completely the horizontal
C  branch part of the evolution.  In addition, we will make the
C  following simplified assumption: t_endrg = 1.1 t_endms.
C  With no need to specify t_endrg separately, we thus are left
C  with only five independent values.
C
        
      function interpolate(t0,t1,t,f0,f1)
      real*8 interpolate, t0,t1,t,f0,f1
      interpolate= f0+(f1-f0)*(t-t0)/(t1-t0)
      end

      subroutine setup_stellar_quantity(idstar,idfunc,
     $     mstime,f0,f1,f2,f3)
      integer idstar, idfunc
      real*8 mstime, f0,f1,f2,f3
#include "modest_common.F"
      ttable(1,idfunc,idstar)= mstime
      ttable(2,idfunc,idstar)= mstime*1.1d0
      ftable(1,idfunc,idstar)= f0
      ftable(2,idfunc,idstar)= f1
      ftable(3,idfunc,idstar)= f2
      ftable(4,idfunc,idstar)= f3
      end

      function stellar_quantity(idstar,idfunc,t)
      real*8 stellar_quantity,t, value, interpolate
      integer idstar, idfunc
#include "modest_common.F"
      if (t .lt. ttable(1,idfunc,idstar)) then
         value = interpolate(0d0, ttable(1,idfunc,idstar),
     $        t, ftable(1,idfunc,idstar), ftable(2,idfunc,idstar))
      elseif(t .lt. ttable(2,idfunc,idstar)) then
         value = interpolate(ttable(1,idfunc,idstar),
     $        ttable(2,idfunc,idstar),
     $        t, ftable(2,idfunc,idstar), ftable(3,idfunc,idstar))
      else
         value = ftable(4,idfunc,idstar)
      endif
      stellar_quantity = value
      end


      
      block data block
#include "modest_common.F"
      data is_used/nmax*0/
      end

      function first_unused_star()
      integer first_unused_star
#include "modest_common.F"
      integer i
      do i=1,nmax
         if (is_used(i) .eq. 0) then
            first_unused_star = i
            return
         endif
      enddo
      first_unused_star = -1
      end
      

      

C     FORTRAN:integer CreateStar(M, Y, Z)
C     should be the constructor here...
      function CreateStar(m_init, Y_init, Z_init)
      integer CreateStar
      integer idstar
      integer first_unused_star
      real*8 m_init, Y_init, Z_init
      real*8 m,r,z0,z1,z2,z3,y0,y1,y2,y3,mstime
#include "modest_common.F"
      idstar = first_unused_star()
      if (idstar .lt. 0) then
         CreateStar= -1
         return
      endif
      CreateStar=idstar
      m = m_init
      mstime =1d4 /(m*m*m)
      ms_lifetime(idstar) = mstime
      standard_timestep(idstar) = mstime*1d-5
      call setup_stellar_quantity(idstar,1,
     $     mstime,m,m, m*0.5D0,m*0.25d0)
      r = m
      call setup_stellar_quantity(idstar,2,mstime,
     $     r,r,r*100d0,0.01d0/m)
      z0 = Z_init
      z1 = Z_init
      z2 = min(1.0d0,z0+0.2d0)
      z3 = min(1.0d0,z0+0.4d0)
      call setup_stellar_quantity(idstar,4,mstime,z0,z1,z2,z3)
      y0 = Y_init
      y1 = min(y0+0.1d0,1d0-z1)
      y2 = min(y0+0.4d0,1d0-z2)
      y3 = min(y0+0.8d0,1d0-z3)
      call setup_stellar_quantity(idstar,3,mstime,y0,y1,y2,y3)
      age(idstar) = 0
      discontinuity_flag(idstar) = 0
      end
      
      function EvolveStar(idstar,dtmax,dMmax, dRmax, dYmax, dZmax)
      real*8 EvolveStar
      integer idstar
      real*8 dtmax,dMmax, dRmax, dYmax, dZmax
      real*8 old_age, last_age, max_age
      real*8 m0,r0,y0,z0
      integer dmax_exceeded
      real*8 t0,t1,newt
      real*8 stellar_quantity
#include "modest_common.F"
      if (discontinuity_flag(idstar) .eq. -1)then
         discontinuity_flag(idstar) = 1
         age(idstar) = age(idstar)
     $        +2*discon_dt*ms_lifetime(idstar)
         EvolveStar = age(idstar)
         return
      endif
      old_age = age(idstar)
      if (discontinuity_flag(idstar) .eq. 1) then
         discontinuity_flag(idstar) =0
      endif
      m0 = stellar_quantity(idstar,1,age)
      r0 = stellar_quantity(idstar,2,age)
      y0 = stellar_quantity(idstar,3,age)
      z0 = stellar_quantity(idstar,4,age)
      last_age = age(idstar)
      max_age = old_age + dtmax
 100  if (dmax_exceeded(idstar,age(idstar),old_age,
     $     dMmax,dRmax,dYmax,dZmax)
     $     .eq. 0 .and. age(idstar) .lt. max_age) then
         last_age = age(idstar)
         age(idstar) = age(idstar)+standard_timestep(idstar)
         if (age(idstar) .gt. max_age ) age(idstar) = max_age
         goto 100
      endif
      if (dmax_exceeded(idstar,age(idstar),old_age,
     $     dMmax,dRmax,dYmax,dZmax)
     $     .ne. 0) then
         t0 = last_age
         t1 = age(idstar)
 200     if ((t1-t0) .gt. discon_dt*ms_lifetime(idstar))then
            newt = (t1+t0)/2.0d0
            if(dmax_exceeded(idstar,newt,old_age,
     $           dMmax,dRmax,dYmax,dZmax) .ne. 0) then
               t1 = newt
            else
               t0 = newt
            endif
            goto 200
         endif
         if(dmax_exceeded(idstar,t1,t0,
     $        dMmax,dRmax,dYmax,dZmax) .ne. 0) then
            discontinuity_flag(idstar) = -1
         endif
         age(idstar) = t0
      endif
      EvolveStar = age(idstar)
      end
      
      function dmax_exceeded(idstar, new_age,old_age,
     $     dMmax,dRmax,dYmax,dZmax)
      integer dmax_exceeded,idstar
      real*8 stellar_quantity
      real*8 new_age,old_age,dMmax,dRmax,dYmax,dZmax
      if (abs(stellar_quantity(idstar,1,new_age)
     $     -stellar_quantity(idstar,1,old_age)) .ge. dMmax .or.
     $     abs(stellar_quantity(idstar,2,new_age)
     $     -stellar_quantity(idstar,2,old_age)) .ge. dRmax .or.
     $     abs(stellar_quantity(idstar,3,new_age)
     $     -stellar_quantity(idstar,3,old_age)) .ge. dYmax .or.
     $     abs(stellar_quantity(idstar,4,new_age)
     $     -stellar_quantity(idstar,4,old_age)) .ge. dZmax)then
          dmax_exceeded = 1
       else
          dmax_exceeded = 0
       endif
       end
      function get_discontinuity_flag(idstar)
      integer get_discontinuity_flag, idstar
#include "modest_common.F"
      get_discontinuity_flag = discontinuity_flag(idstar)
      end
      function DestroyStar(idstar)
      integer DestroyStar
#include "modest_common.F"
      integer idstar
      if (idstar .lt. 0 .or.
     $     idstar .gt. nmax) then
         DestroyStar = -1
      else
         if (is_used(idstar) .eq. 0) then
            DestroyStar = -1
         else
            DestroyStar = idstar
            is_used(idstar) = 0
         endif
      endif
      end

      function getMass(idstar)
      real*8 getMass
      integer idstar
#include "modest_common.F"
      real*8 stellar_quantity
      getMass = stellar_quantity(idstar,1,age(idstar))
      end
      
      function getRadius(idstar)
      real*8 getRadius
      integer idstar
#include "modest_common.F"
      real*8 stellar_quantity
      getRadius = stellar_quantity(idstar,2,age(idstar))
      end
      
      function getTime(idstar)
      real*8 getTime
      integer idstar
#include "modest_common.F"
      getTime = age(idstar)
      end
      
      function getY(idstar)
      real*8 getY
      integer idstar
#include "modest_common.F"
      real*8 stellar_quantity
      getY = stellar_quantity(idstar,3,age(idstar))
      end
      
      
      function getZ(idstar)
      real*8 getZ
      integer idstar
#include "modest_common.F"
      real*8 stellar_quantity
      getZ = stellar_quantity(idstar,4,age(idstar))
      end
      

modest_common.F

C------------------------------------------------------------
C	MODEST_INC.F
C
C       Piet Hut and Jun Makino
C       Version 0.0 July 1, 2001
C------------------------------------------------------------
C      
      integer nmax
      parameter (nmax = 100000)
      integer nfmax
      parameter (nfmax = 4)
      real*8 discon_dt
      parameter (discon_dt = 1d-8)
      real*8 ftable(4,nfmax,nmax), ttable(2,nfmax,nmax)
      real*8 ms_lifetime(nmax)
      real*8 standard_timestep(nmax)
      real*8 discontinuity_time(nmax)
      integer  discontinuity_flag(nmax)
      integer  is_used(nmax)
      real*8 age(nmax)

      common /modest/ftable,ttable,ms_lifetime,standard_timestep,
     $   discontinuity_time, discontinuity_flag, age

      

Sample output from C++ code

Output for a one-solar-mass star:

 Age = 1            M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 3            M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 7            M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 15           M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 31           M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 63           M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 127          M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 255          M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 511          M = 1            R = 1            Y = 0.25   Z = 0.02
 Age = 1023         M = 1            R = 1            Y = 0.26   Z = 0.02
 Age = 2047         M = 1            R = 1            Y = 0.27   Z = 0.02
 Age = 4095         M = 1            R = 1            Y = 0.29   Z = 0.02
 Age = 8191         M = 1            R = 1            Y = 0.33   Z = 0.02
 Age = 10010.1010   M = 0.99494951   R = 1.99999668   Y = 0.35   Z = 0.02
 Age = 10030.3029   M = 0.98484854   R = 3.99999004   Y = 0.35   Z = 0.02
 Age = 10070.7068   M = 0.96464658   R = 7.99997676   Y = 0.37   Z = 0.03
 Age = 10151.5146   M = 0.92424268   R = 15.9999502   Y = 0.39   Z = 0.05
 Age = 10313.1303   M = 0.84343486   R = 31.9998971   Y = 0.44   Z = 0.08
 Age = 10481.8172   M = 0.75909141   R = 48.6999016   Y = 0.49   Z = 0.11
 Age = 10633.6354   M = 0.68318228   R = 63.7299095   Y = 0.54   Z = 0.14
 Age = 10770.2719   M = 0.61486406   R = 77.2569156   Y = 0.58   Z = 0.17
 Age = 10893.2446   M = 0.55337769   R = 89.4312183   Y = 0.61   Z = 0.19
 Age = 10999.9999   M = 0.50000005   R = 99.9999903   Y = 0.64   Z = 0.21
 Age = 11000.0001   M = 0.25         R = 0.01         Y = 0.58   Z = 0.42
 Age = 19192.0001   M = 0.25         R = 0.01         Y = 0.58   Z = 0.42
Note the at first the time step doubling was the limiting factor. Then during the first five steps on the red giant branch, the time step was constrained by the requirement to at most double the radius. During the next five steps the constraint was to allow a change in mass of not more than ten percent. The final step was again constrained by time step doubling, while the two before were forced to occur at both sides of the discontinuity at the end of the life of the star.

Driver source code in C++

Here is the C++ source code for testing the toy model implementation described above. This is not meant for distribution; it is still very rough, quickly cobbled together to get some results; but if you're interested to look under the hood at what we've done, here it is.


#include  <iostream>
#include  <cmath>
#include  <cstdlib>
#include  "modest_star.H"
using namespace std;

void print_star(modest_star star)
{
    cout.precision(10);
    cout << " Age = " << star.getTime()
	 << " M = " << star.getMass()
	 << " R = " << star.getRadius()
	 << " Y = " << star.getY()
	 << " Z = " << star.getZ() << endl;
}



int main()
{
    modest_star star = modest_star(1,0.25,0.02);
    double dtmax = 1;
    int discon_reached = 0;
    double t = 0;
    int countdown = 3;
    while(countdown>0){
	double dMmax = 0.1*star.getMass();
	double dRmax = 1*star.getRadius();
	double dYmax = 0.1;
	double dZmax = 0.1;
	double new_time = star.EvolveStar(dtmax,dMmax,dRmax,dYmax,dZmax);
	print_star(star);
	if ((new_time -t)> dtmax - 1e-10) dtmax *= 2;
	t = new_time;
	if (star.get_discontinuity_flag()) discon_reached = 1;
	if (discon_reached)countdown --;
    }
}

      
	

Toy model source code in C++(modest_star.H)

Here is the C++ source code for toy model implementation (as a header file containing class definitions). This too is not meant for distribution; it is still very rough, quickly cobbled together to get some results; but if you're interested to look under the hood at what we've done, here it is.

//
// Each physical quantity Q that describes the state of a star
// can be encapsulated as an object of class `stellar_quantity',
// where Q can stand for mass, radius, etc.  Such an object
// contains the full information about the evolution of Q over
// the life time of a star.  The value of a quantity Q at time t
// can be obtained by invoking its member function Q.at(t) .
//
// In our current toy-model implementation, Q(t) is described by
// a piece-wise linear function with one discontinuity, specified
// by the following six numbers:
//
//   t_endms   time at which the star leaves the main sequence
//   t_endrg   time at which the star leaves the red giant branch.
//   f_init    value of Q at birth of star (at time t=0)
//   f_endms   value of Q at time t=t_endms
//   f_endrg   value of Q just before time t=t_endrg
//   f_remnant value of Q just after time t=t_endrg
//
//                                      ....... f_endrg
//                                     /|
//                                    / |
//                                   /  |
//                                  /   |
//  ^                              /    |
//  |                             /     |
//  Q                            /      |
//                         _____/.......|...... f_endms
//               _____-----             |
//     _____-----.......................|...... f_init
//                                      |------ f_remnant
//     +-----------------------+--------+------ 0
//     0                       t_endms  t_endrg
//                t -->
//
// Note that we stick-figure version of stellar evolution mimics
// only low-mass stars, while leaving out completely the horizontal
// branch part of the evolution.  In addition, we will make the
// following simplified assumption: t_endrg = 1.1 t_endms.
// With no need to specify t_endrg separately, we thus are left
// with only five independent values.

class stellar_quantity
{
private:
    double t_endms, t_endrg;
    double f_init, f_endms, f_endrg, f_remnant;
    double interpolate(double t0,
		       double t1,
		       double t,
		       double f0,
		       double f1)
    {
	return f0+(f1-f0)*(t-t0)/(t1-t0);
    }
public:
    stellar_quantity(double ms_lifetime=0,
			      double f0=0,
				   double f1=0,
				   double f2=0,
				   double f3=0)
    {
	t_endms = ms_lifetime;
	t_endrg = ms_lifetime*1.1;
	f_init = f0;
	f_endms = f1;
	f_endrg = f2;
	f_remnant = f3;
    }
    double at(double t)
    {
	if (t < t_endms)
	    return interpolate(0.0,    t_endms, t, f_init, f_endms);
	if (t < t_endrg)
	    return interpolate(t_endms,t_endrg, t, f_endms,f_endrg);
	return f_remnant;
    }
    double discontinuity_time()
    {
	return t_endrg;
    }
};

double discon_dt = 1e-8;

class modest_star
{
private:
    stellar_quantity mass, radius, Y, Z;
    double ms_lifetime;
    double discontinuity_time;
    int  discontinuity_flag; // Should be enum...
                             // for now, -1= just before discontinuity
                             //          +1= just after  discontinuity
                             //           0= no discontinuity
    double age;
    double standard_timestep;
    double m0,r0,y0,z0,dmmax,drmax,dymax,dzmax;
    int dmax_exceeded(double age)
    {
	return ( (fabs(mass.at(age)-m0)>=dmmax)||
		 (fabs(radius.at(age)-r0)>=drmax)||
		 (fabs(Y.at(age)-y0)>=dymax)||
		 (fabs(Z.at(age)-z0)>=dzmax));
    }

public:
    //      FORTRAN:integer CreateStar(M, Y, Z)
    //      should be the constructor here...
    modest_star(double m_init, double Y_init, double Z_init){
	double m = m_init;
	ms_lifetime = 1e4 /(m*m*m);
	standard_timestep = ms_lifetime*1e-5;
	mass = stellar_quantity(ms_lifetime,m,m, m*0.5,m*0.25);
	double r = m;// simple linear mass-radius relation
                             //  in Msol-Rsol units
	radius = stellar_quantity(ms_lifetime,r,r,r*100,0.01/m);
	double z0 = Z_init;
	double z1 = Z_init;
	double z2 = fmin(1,z0+0.2);
	double z3 = fmin(1,z0+0.4);
	Z = stellar_quantity(ms_lifetime,z0,z1,z2,z3);
	double y0 = Y_init;
	double y1 = fmin(y0+0.1,1-z1);
	double y2 = fmin(y0+0.4,1-z2);
	double y3 = fmin(y0+0.8,1-z3);
	Y = stellar_quantity(ms_lifetime,y0,y1,y2,y3);
	age = 0;
	discontinuity_flag = 0;
    }
    //      FORTRAN:real*8 EvolveStar(id, dtmax, dMmax, dRmax, dYmax, dZmax)
    //
    // Sample implementation for checking dMmax etc.
    // return time to the accuracy of discon_dt

    double  EvolveStar(double dtmax,
		       double dMmax,
		       double dRmax,
		       double dYmax,
		       double dZmax )
    {
	if (discontinuity_flag == -1){
	    discontinuity_flag = 1;
	    age += 2*discon_dt*ms_lifetime;
	    return age;
	}
	double old_age = age;
	if (discontinuity_flag == 1)discontinuity_flag =0;
	m0 = mass.at(age);
	r0 = radius.at(age);
	y0 = Y.at(age);
	z0 = Z.at(age);
	dmmax = dMmax;
	drmax = dRmax;
	dymax = dYmax;
	dzmax = dZmax;
	double last_age = age;
        double max_age = old_age + dtmax;
	while((!dmax_exceeded(age))
	      &&(age < max_age)) {
	    last_age = age;
	    age += standard_timestep;
	    if (age  > max_age ) age = max_age;
	}
	if (dmax_exceeded(age)){
	    double t0 = last_age;
	    double t1 = age;
	    while ((t1-t0)>discon_dt*ms_lifetime){
		double newt = (t1+t0)/2;
		if(dmax_exceeded(newt)){
		    t1 = newt;
		}else{
		    t0 = newt;
		}
	    }
	    if(fabs(mass.at(t1)-mass.at(t0))>dMmax ||
	       fabs(radius.at(t1)-radius.at(t0))>dRmax ||
	       fabs(Y.at(t1)-Y.at(t0))>dYmax ||
	       fabs(Z.at(t1)-Z.at(t0))>dZmax ){
		discontinuity_flag = -1;
	    }
	    age = t0;
	    
	}
	return age;
    }
    int get_discontinuity_flag()
    {
	return discontinuity_flag;
    }
    
    //      integer DestroyStar(id)   --- replaced by default destructor
    
    //     real*8 getMass(id)
    double getMass()
    {
	return mass.at(age);
    }
	
    //      real*8 getRadius(id)
    double getRadius()
    {
	return radius.at(age);
    }
    //      real*8 getTime(id)
    double getTime()
    {
	return age;
    }
    //      real*8 getY(id)
    double getY()
    {
	return Y.at(age);
    }
    //      real*8 getZ(id)
    double getZ()
    {
	return Z.at(age);
    }
};