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

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 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"
\$     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)
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)
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)
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)
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

integer idstar
#include "modest_common.F"
real*8 stellar_quantity
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 standard_timestep(nmax)
real*8 discontinuity_time(nmax)
integer  discontinuity_flag(nmax)
integer  is_used(nmax)
real*8 age(nmax)

\$   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 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:
double f0=0,
double f1=0,
double f2=0,
double f3=0)
{
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:
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(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;
double r = m;// simple linear mass-radius relation
//  in Msol-Rsol units
double z0 = Z_init;
double z1 = Z_init;
double z2 = fmin(1,z0+0.2);
double z3 = fmin(1,z0+0.4);
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);
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;
return age;
}
double old_age = age;
if (discontinuity_flag == 1)discontinuity_flag =0;
m0 = mass.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;
double newt = (t1+t0)/2;
if(dmax_exceeded(newt)){
t1 = newt;
}else{
t0 = newt;
}
}
if(fabs(mass.at(t1)-mass.at(t0))>dMmax ||
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 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);
}
};
```