updated comprad for METISSE

Author: poojanagrawal

.gitignore

@@ -246,7 +246,7 @@ def evolve(cls, initialbinarytable, pool=None, **kwargs):
         # check the initial conditions of the system and warn user if
         # anything is weird about them, such as the star starts
         # in Roche Lobe overflow
-        utils.check_initial_conditions(initialbinarytable)
+        utils.check_initial_conditions(initialbinarytable,SSEDict)
 
         # assign some columns based on keyword arguments but that
         # can be overwritten by the params or BSEDict
@@ -486,6 +486,8 @@ def _evolve_single_system(f):
             _evolvebin.se_flags.using_sse = False
             path_to_tracks = f["path_to_tracks"]
             path_to_he_tracks = f["path_to_he_tracks"]
+        else:
+            raise ValueError("Use either 'sse' or 'metisse' as stellar engine")
 
         [bpp_index, bcm_index, kick_info] = _evolvebin.evolv2([f["kstar_1"], f["kstar_2"]],
                                                               [f["mass_1"], f["mass_2"]],

meson.build

@@ -45,6 +45,7 @@ def get_independent_sampler(
     SF_duration,
     binfrac_model,
     met,
+    SSEDict=None,
     size=None,
     total_mass=np.inf,
     sampling_target="size",
@@ -159,6 +160,8 @@ def get_independent_sampler(
     n_binaries : `int`
         Number of binaries needed to generate a population
     """
+    
+    
     if sampling_target == "total_mass" and (total_mass is None or total_mass == np.inf):
         raise ValueError("If `sampling_target == 'total mass'` then `total_mass` must be supplied")
     if size is None and (total_mass is None or total_mass == np.inf):
@@ -284,8 +287,8 @@ def get_independent_sampler(
 
     zsun = kwargs.pop("zsun", 0.02)
 
-    rad1 = initconditions.set_reff(mass1_binary, metallicity=met, zsun=zsun)
-    rad2 = initconditions.set_reff(mass2_binary, metallicity=met, zsun=zsun)
+    rad1 = initconditions.set_reff(SSEDict, mass1_binary, metallicity=met, zsun=zsun)
+    rad2 = initconditions.set_reff(SSEDict, mass2_binary, metallicity=met, zsun=zsun)
 
     # sample periods and eccentricities
     # if the porb_model is moe19, the metallicity needs to be supplied
@@ -1072,7 +1075,7 @@ def set_kstar(self, mass):
 
         return kstar
 
-    def set_reff(self, mass, metallicity, zsun=0.02):
+    def set_reff(self, SSEDict ,mass, metallicity, zsun=0.02):
         """
         Better way to set the radii from BSE, by calling it directly
 
@@ -1092,12 +1095,25 @@ def set_reff(self, mass, metallicity, zsun=0.02):
 
         _evolvebin.metvars.zsun = zsun
 
+        if (SSEDict == None) or (SSEDict["stellar_engine"] == "sse"):
+            _evolvebin.se_flags.using_sse = True
+            _evolvebin.se_flags.using_metisse = False
+            path_to_tracks = ""
+            path_to_he_tracks = ""
+        elif SSEDict["stellar_engine"] == "metisse":
+            _evolvebin.se_flags.using_metisse = True
+            _evolvebin.se_flags.using_sse = False
+            path_to_tracks = SSEDict["path_to_tracks"]
+            path_to_he_tracks = SSEDict["path_to_he_tracks"]
+        else:
+            raise ValueError("Use either 'sse' or 'metisse' as stellar engine")
+            
         idx = 0
         while total_length > max_array_size:
             ## cycle through the masses max_array_size number at a time
             temp_mass = mass[idx*max_array_size:(idx+1)*max_array_size]
 
-            temp_radii = _evolvebin.compute_r(temp_mass,metallicity,max_array_size)
+            temp_radii = _evolvebin.compute_r(temp_mass,metallicity,max_array_size,path_to_tracks,path_to_he_tracks)
 
             ## put these in the radii array
             radii[idx*max_array_size:(idx+1)*max_array_size] = temp_radii
@@ -1111,7 +1127,7 @@ def set_reff(self, mass, metallicity, zsun=0.02):
         temp_mass = np.zeros(max_array_size)
         temp_mass[:length_remaining] = mass[-length_remaining:]
 
-        temp_radii = _evolvebin.compute_r(temp_mass,metallicity,length_remaining)
+        temp_radii = _evolvebin.compute_r(temp_mass,metallicity,length_remaining,path_to_tracks,path_to_he_tracks)
 
         #finish up the array
         radii[-length_remaining:] = temp_radii[:length_remaining]

src/cosmic/evolve.py

@@ -1,5 +1,6 @@
 ***
-      SUBROUTINE compute_r(mass,z,num,rad)
+      SUBROUTINE compute_r(mass,z,num,rad,
+     &           path_to_tracks,path_to_he_tracks)
       IMPLICIT NONE
       INCLUDE 'const_bse.h'
 
@@ -17,8 +18,9 @@ SUBROUTINE compute_r(mass,z,num,rad)
       PARAMETER(loop=100000)
       real*8 mass(loop),rad(loop),z
       integer k,kstar,num
-      real*8 mt,tm,tn,mass0,age,lum,mc,rc,me,re
+      real*8 mt,tm,tn,mass0,age,lum,mc,rc,me,re,dtm
       REAL*8 tscls(20),lums(10),GB(10),zpars(20),k2,bhspin
+      CHARACTER*256 path_to_tracks,path_to_he_tracks
 
 ***
 * f2py directives go here; we'll return the radii as a 10^5 array
@@ -29,7 +31,11 @@ SUBROUTINE compute_r(mass,z,num,rad)
 Cf2py intent(in) num
 Cf2py intent(out) rad
 
-      CALL zcnsts(z,zpars)
+      if(using_METISSE) CALL initialize_front_end('cosmic')
+      CALL zcnsts(z,zpars,path_to_tracks,path_to_he_tracks)
+      
+      if(using_METISSE) call allocate_track(num,mass)
+
 
 ***
 * Then just loop through everything 
@@ -46,11 +52,15 @@ SUBROUTINE compute_r(mass,z,num,rad)
          kstar = 0
          if(mt.ge.0.7) kstar = 1
          bhspin = 0.d0
-         rc = 0.d0 
-         CALL star(kstar,mass0,mt,tm,tn,tscls,lums,GB,zpars)
+         rc = 0.d0
+         dtm = 0.d0
+         CALL star(kstar,mass0,mt,tm,tn,tscls,lums,GB,zpars,dtm,k)
          CALL hrdiag(mass0,age,mt,tm,tn,tscls,lums,GB,zpars,
      &               rad(k),lum,kstar,mc,rc,me,re,k2,bhspin,k)
-
+         
   10  continue
+  
+        
+        if (using_METISSE) call dealloc_track()
 
       END SUBROUTINE compute_r

src/cosmic/sample/sampler/independent.py

@@ -1602,13 +1602,14 @@ def error_check(BSEDict, SSEDict, filters=None, convergence=None, sampling=None)
     return
 
 
-def check_initial_conditions(full_initial_binary_table):
+def check_initial_conditions(full_initial_binary_table, SSEDict):
     """Checks initial conditions and reports warnings
 
     Only warning provided right now is if star begins in Roche lobe
     overflow
     """
-
+    from cosmic import _evolvebin
+    
     def rzamsf(m):
         """A function to evaluate Rzams
         ( from Tout et al., 1996, MNRAS, 281, 257 ).
@@ -1635,10 +1636,20 @@ def rzamsf(m):
 
     if np.all(mass2 == 0.0):
         return
-    else:
+    elif (SSEDict["stellar_engine"] == "sse"):
         rzams1 = rzamsf(mass1)
         rzams2 = rzamsf(mass2)
-
+        
+#        TODO: the following section with METISSE is incomplete
+#    elif (SSEDict["stellar_engine"] == "metisse"):
+#        _evolvebin.se_flags.using_sse = False
+#        _evolvebin.se_flags.using_metisse = True
+#        path_to_tracks = SSEDict["path_to_tracks"]
+#        path_to_he_tracks = SSEDict["path_to_he_tracks"]
+        
+#        rzams1 = _evolvebin.compute_r(mass1,z,len(mass1),path_to_tracks,path_to_he_tracks)
+#        rzams1 = _evolvebin.compute_r(mass2,z,len(mass2),path_to_tracks,path_to_he_tracks)
+        
         # assume some time step in order to calculate sep
         yeardy = 365.24
         aursun = 214.95