Metisse integration test

src/cosmic/evolve.py

@@ -95,7 +95,7 @@
 
 INITIAL_CONDITIONS_MISC_COLUMN = ['bin_num']
 
-INITIAL_CONDITIONS_SSE_COLUMN = ['stellar_engine','tracks_path']
+INITIAL_CONDITIONS_SSE_COLUMN = ['stellar_engine','path_to_tracks','path_to_he_tracks']
 
 # Add the BSE COLUMSN and MISC COLUMN to the PASS_COLUMNS list
 INITIAL_CONDITIONS_PASS_COLUMNS.extend(INITIAL_CONDITIONS_BSE_COLUMNS)
@@ -478,11 +478,13 @@ def _evolve_single_system(f):
         if f["stellar_engine"] == "sse":
             _evolvebin.se_flags.using_sse = True
             _evolvebin.se_flags.using_metisse = False
-            metisse_path = ""
+            path_to_tracks = ""
+            path_to_he_tracks = ""
         elif f["stellar_engine"] == "metisse":
             _evolvebin.se_flags.using_metisse = True
             _evolvebin.se_flags.using_sse = False
-            metisse_path = f["tracks_path"]
+            path_to_tracks = f["path_to_tracks"]
+            path_to_he_tracks = f["path_to_he_tracks"]
 
         [bpp_index, bcm_index, kick_info] = _evolvebin.evolv2([f["kstar_1"], f["kstar_2"]],
                                                               [f["mass_1"], f["mass_2"]],
@@ -505,7 +507,8 @@ def _evolve_single_system(f):
                                                               np.zeros(20),
                                                               np.zeros(20),
                                                               f["kick_info"],
-                                                              metisse_path)
+                                                              path_to_tracks,
+                                                              path_to_he_tracks)
         bcm = _evolvebin.binary.bcm[:bcm_index].copy()
         bpp = _evolvebin.binary.bpp[:bpp_index].copy()
         _evolvebin.binary.bpp[:bpp_index] = np.zeros(bpp.shape)

src/cosmic/meson.build

@@ -50,6 +50,8 @@ lib_source = [
     'src/SSE/SSE_star.f',
     'src/SSE/SSE_zcnsts.f',
     'src/SSE/SSE_zfuncs.f',
+    'src/SSE/SSE_gntage.f',
+    'src/METISSE/src/METISSE_gntage.f90',
     'src/METISSE/src/METISSE_deltat.f90',
     'src/METISSE/src/METISSE_mlwind.f90',
     'src/METISSE/src/METISSE_hrdiag.f90',

src/cosmic/src/SSE/SSE_gntage.f

@@ -0,0 +1,303 @@
+***
+      SUBROUTINE SSE_gntage(mc,mt,kw,zpars,m0,aj,k)
+      IMPLICIT NONE
+      INCLUDE '../const_bse.h'
+*
+* A routine to determine the age of a giant from its core mass and type.
+*
+*       Author : C. A. Tout
+*       Date   : 24th September 1996
+*       Revised: 21st February 1997 to include core-helium-burning stars
+*
+*       Rewritten: 2nd January 1998 by J. R. Hurley to be compatible with
+*                  the new evolution routines and to include new stellar
+*                  types.
+*
+*
+      integer kw, k
+      integer j,jmax
+      parameter(jmax=30)
+*
+      real*8 mc,mt,m0,aj,tm,tn,dtm
+      real*8 tscls(20),lums(10),GB(10),zpars(20)
+      real*8 mmin,mmax,mmid,dm,f,fmid,dell,derl,lum
+      real*8 macc,lacc,tiny
+      parameter(macc=0.00001d0,lacc=0.0001d0,tiny=1.0d-14)
+      real*8 mcx,mcy
+*
+      real*8 mcheif,mcagbf,mheif,mbagbf,mcgbf,lmcgbf,lbgbf,lbgbdf
+      external mcheif,mcagbf,mheif,mbagbf,mcgbf,lmcgbf,lbgbf,lbgbdf
+*
+* This should only be entered with KW = 3, 4, 5, 6 or 9
+*
+* First we check that we don't have a CheB star
+* with too small a core mass.
+*      write(UNIT=99,fmt=*) 'in gntage',kw,aj
+
+      if(kw.eq.4)then
+* Set the minimum CHeB core mass using M = Mflash
+         mcy = mcheif(zpars(2),zpars(2),zpars(10))
+         if(mc.le.mcy) kw = 3
+*        if(mc.le.mcy) WRITE(66,*)' GNTAGE4: changed to 3'
+      endif
+*
+* Next we check that we don't have a GB star for M => Mfgb
+      if(kw.eq.3)then
+* Set the maximum GB core mass using M = Mfgb
+         mcy = mcheif(zpars(3),zpars(2),zpars(9))
+         if(mc.ge.mcy)then
+            kw = 4
+            aj = 0.d0
+*           WRITE(66,*)' GNTAGE3: changed to 4'
+         endif
+      endif
+*
+      if(kw.eq.6)then
+*
+* We try to start the star from the start of the SAGB by
+* setting Mc = Mc,TP.
+*
+         mcy = 0.44d0*2.25d0 + 0.448d0
+         if(mc.gt.mcy)then
+* A type 6 with this sized core mass cannot exist as it should
+* already have become a NS or BH as a type 5.
+* We set it up so that it will.
+            mcx = (mc + 0.35d0)/0.773d0
+         elseif(mc.ge.0.8d0)then
+            mcx = (mc - 0.448d0)/0.44d0
+         else
+            mcx = mc
+         endif
+         m0 = mbagbf(mcx)
+         if(m0.lt.tiny)then
+* Carbon core mass is less then the minimum for the start of SAGB.
+* This must be the case of a low-mass C/O or O/Ne WD with only a
+* very small envelope added or possibly the merger of a helium star
+* with a main sequence star. We will set m0 = mt and then reset the
+* core mass to allow for some helium to be added to the C/O core.
+            kw = 14
+*           WRITE(66,*)' GNTAGE6: changed to 4'
+         else
+            CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+            aj = tscls(13) + 2.d0*tiny
+         endif
+      endif
+*
+      if(kw.eq.5)then
+*
+* We fit a Helium core mass at the base of the AGB.
+*
+         m0 = mbagbf(mc)
+         if(m0.lt.tiny)then
+* Helium core mass is less then the BAGB minimum.
+            kw = 14
+*           WRITE(66,*)' GNTAGE5: changed to 4'
+         else
+            CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+            aj = tscls(2) + tscls(3) + 2.d0*tiny
+         endif
+      endif
+*
+*
+      if(kw.eq.4)then
+*
+* The supplied age is actually the fractional age, fage, of CHeB lifetime
+* that has been completed, ie. 0 <= aj <= 1.
+*
+         if(aj.lt.0.d0.or.aj.gt.1.d0)then
+*           WRITE(99,*)' FATAL ERROR! GNTAGE4: fage out of bounds '
+*           WRITE(99,*)' FAGE ',aj
+*           WRITE(*,*)' STOP: FATAL ERROR '
+*           CALL exit(0)
+*           STOP
+            aj = 0.d0
+         endif
+* Get the minimum, fage=1, and maximum, fage=0, allowable masses
+         mcy = mcagbf(zpars(2))
+         if(mc.ge.mcy)then
+            mmin = mbagbf(mc)
+         else
+            mmin = zpars(2)
+         endif
+         mmax = mheif(mc,zpars(2),zpars(10))
+         if(aj.lt.tiny)then
+            m0 = mmax
+            goto 20
+         elseif(aj.ge.1.d0)then
+            m0 = mmin
+            goto 20
+         endif
+* Use the bisection method to find m0
+         fmid = (1.d0-aj)*mcheif(mmax,zpars(2),zpars(10)) +
+     &          aj*mcagbf(mmax) - mc
+         f = (1.d0-aj)*mcheif(mmin,zpars(2),zpars(10)) +
+     &       aj*mcagbf(mmin) - mc
+         if(f*fmid.ge.0.d0)then
+* This will probably occur if mc is just greater than the minimum
+* allowed mass for a CHeB star and fage > 0.
+            kw = 3
+*           WRITE(66,*)' GNTAGE4: changed to 3'
+            goto 90
+         endif
+         m0 = mmin
+         dm = mmax - mmin
+         do 10 , j = 1,jmax
+            dm = 0.5d0*dm
+            mmid = m0 + dm
+            fmid = (1.d0-aj)*mcheif(mmid,zpars(2),zpars(10)) +
+     &             aj*mcagbf(mmid) - mc
+            if(fmid.lt.0.d0) m0 = mmid
+            if(ABS(dm).lt.macc.or.ABS(fmid).lt.tiny) goto 20
+            if(j.eq.jmax)then
+*              WRITE(99,*)' FATAL ERROR! GNTAGE4: root not found '
+*              WRITE(*,*)' STOP: FATAL ERROR '
+*              CALL exit(0)
+*              STOP
+               m0 = mt
+               aj = 0.d0
+            endif
+ 10      continue
+ 20      continue
+*
+         CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+         aj = tscls(2) + aj*tscls(3)
+*
+      endif
+*
+ 90   continue
+*
+      if(kw.eq.3)then
+*
+* First we double check that we don't have a GB star for M => Mfgb
+         mcy = mcheif(zpars(3),zpars(2),zpars(9))
+         if(mc.ge.mcy)then
+*           WRITE(99,*)' GNTAGE3: star too big for GB '
+*           WRITE(*,*)' STOP: FATAL ERROR '
+*           CALL exit(0)
+*           STOP
+            mc = 0.99d0*mcy
+         endif
+* Next we find an m0 so as to place the star at the BGB
+         mcx = mcheif(zpars(2),zpars(2),zpars(9))
+         if(mc.gt.mcx)then
+            m0 = mheif(mc,zpars(2),zpars(9))
+         else
+* Use Newton-Raphson to find m0 from Lbgb
+            m0 = zpars(2)
+            CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+            lum = lmcgbf(mc,GB)
+            j = 0
+ 30         continue
+            dell = lbgbf(m0) - lum
+            if(ABS(dell/lum).le.lacc) goto 40
+            derl = lbgbdf(m0)
+            m0 = m0 - dell/derl
+            j = j + 1
+            if(j.eq.jmax)then
+*              WRITE(99,*)' FATAL ERROR! GNTAGE3: root not found '
+*              WRITE(*,*)' STOP: FATAL ERROR '
+*              CALL exit(0)
+*              STOP
+               m0 = zpars(2)
+               m0 = MAX(m0,mt)
+               goto 40
+            endif
+            goto 30
+ 40         continue
+         endif
+         CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+         aj = tscls(1) + 1.0d-06*(tscls(2) - tscls(1))
+*
+      endif
+*
+      if(kw.eq.8.or.kw.eq.9)then
+*
+* We make a post-MS naked helium star.
+* To make things easier we put the star at the TMS point
+* so it actually begins as type 8.
+*
+         kw = 8
+         mmin = mc
+         CALL star(kw,mmin,mc,tm,tn,tscls,lums,GB,zpars,dtm,k)
+         mcx = mcgbf(lums(2),GB,lums(6))
+         if(mcx.ge.mc)then
+*           WRITE(99,*)' FATAL ERROR! GNTAGE9: mmin too big '
+*           WRITE(*,*)' STOP: FATAL ERROR '
+*           CALL exit(0)
+*           STOP
+            m0 = mt
+            goto 80
+         endif
+         f = mcx - mc
+         mmax = mt
+         do 50 , j = 1,jmax
+            CALL star(kw,mmax,mc,tm,tn,tscls,lums,GB,zpars,dtm,k)
+            mcy = mcgbf(lums(2),GB,lums(6))
+            if(mcy.gt.mc) goto 60
+            mmax = 2.d0*mmax
+            if(j.eq.jmax)then
+*              WRITE(99,*)' FATAL ERROR! GNTAGE9: mmax not found '
+*              WRITE(*,*)' STOP: FATAL ERROR '
+*              CALL exit(0)
+*              STOP
+               m0 = mt
+               goto 80
+            endif
+ 50      continue
+ 60      continue
+         fmid = mcy - mc
+* Use the bisection method to find m0
+         if(f*fmid.ge.0.d0)then
+*           WRITE(99,*)' FATAL ERROR! GNTAGE9: root not bracketed '
+*           WRITE(*,*)' STOP: FATAL ERROR '
+*           CALL exit(0)
+*           STOP
+            m0 = mt
+            goto 80
+         endif
+         m0 = mmin
+         dm = mmax - mmin
+         do 70 , j = 1,jmax
+            dm = 0.5d0*dm
+            mmid = m0 + dm
+            CALL star(kw,mmid,mc,tm,tn,tscls,lums,GB,zpars,dtm,k)
+            mcy = mcgbf(lums(2),GB,lums(6))
+            fmid = mcy - mc
+            if(fmid.lt.0.d0) m0 = mmid
+            if(ABS(dm).lt.macc.or.ABS(fmid).lt.tiny) goto 80
+            if(j.eq.jmax)then
+*              WRITE(99,*)' FATAL ERROR! GNTAGE9: root not found '
+*              WRITE(*,*)' STOP: FATAL ERROR '
+*              CALL exit(0)
+*              STOP
+               m0 = mt
+               goto 80
+            endif
+ 70      continue
+ 80      continue
+*
+         CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+         aj = tm + 1.0d-10*tm
+*
+      endif
+*
+      if(kw.eq.14)then
+*
+         kw = 4
+         m0 = mt
+         mcy = mcagbf(m0)
+         aj = mc/mcy
+         CALL star(kw,m0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
+         if(m0.le.zpars(2))then
+            mcx = mcgbf(lums(4),GB,lums(6))
+         else
+            mcx = mcheif(m0,zpars(2),zpars(10))
+         end if
+         mc = mcx + (mcy - mcx)*aj
+         aj = tscls(2) + aj*tscls(3)
+      endif
+*
+*      print*, 'exiting gntage',kw,aj
+      RETURN
+      END
+***

src/cosmic/src/bpp_array.f

@@ -111,9 +111,14 @@ SUBROUTINE WRITEBCM(ip,tphys,kstar_1,mass0_1,mass_1,
         INTEGER ip
         REAL*8 yeardy,aursun,rsunau
         PARAMETER(yeardy=365.24d0,aursun=214.95d0)
-
+
+        bcm_err = .false.
+        if (ip+2>size(bcm,1)) then
+            bcm_err = .true.
+            return
+        endif
+
         rsunau = 1/aursun
-
         ip = ip + 1
         bcm(ip,1) = tphys
         bcm(ip,2) = float(kstar_1)

src/cosmic/src/comenv.f

@@ -896,6 +896,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
          if(output) write(*,*)'coel 2 1:',KW,KW1,KW2,M1,M2,MF,MC22,
      & TB,OORB
          IF(KW.EQ.2)THEN
+            if (using_METISSE) call set_star_type(star1)
             CALL star(KW,M1,M1,TM2,TN,TSCLS2,LUMS,GB,ZPARS,dtm,star1)
             IF(GB(9).GE.MC1)THEN
                M01 = M1
@@ -907,6 +908,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
      & TB,OORB
          ELSEIF(KW.EQ.7)THEN
             M01 = M1
+            if (using_METISSE) call set_star_type(star1)
             CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,dtm,star1)
             AJ1 = TM1*(FAGE1*MC1 + FAGE2*MC22)/(MC1 + MC22)
             if(output) write(*,*)'coel 2 3:',KW,KW1,KW2,M1,M01,MC22,
@@ -918,7 +920,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
 *
 * Obtain a new age for the giant.
 *
-            CALL gntage(MC1,M1,KW,ZPARS,M01,AJ1,dtm,star1)
+            CALL gntage(MC1,M1,KW,ZPARS,M01,AJ1,star1)
             CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,dtm,star1)
             if(output) write(*,*)'coel 2 4:',KW,KW1,KW2,M1,M01,MC22,
      & TB,OORB

src/cosmic/src/const_bse.h

@@ -61,4 +61,6 @@
 
       LOGICAL using_metisse, using_sse
       COMMON /SE_FLAGS/ using_metisse, using_sse
+      LOGICAL bcm_err
+      COMMON/ ER_FLAGS/ bcm_err
 *