Fixes for tides, mean_anomaly, and metallicity error checking

changelog.md

@@ -112,4 +112,11 @@ This release contains _several_ fixes to how CO core masses/remnant masses are h
 
 - Documentation:
     - Start new settings gallery in the documentation
-    - Tag settings/options with the version they were added in the docs page and auto link them to release
+    - Tag settings/options with the version they were added in the docs page and auto link them to release
+
+## 3.7.5
+
+- Fixes:
+    - Prevent tides calculations based on NS/BHs
+    - Track mean_anomaly (save to kick_info) array even when a system is disrupted
+    - Throw an error if invalid metallicities are provided for evolution

debug/create_binary_in.py

@@ -39,7 +39,9 @@ def create_binary_in(mass0, tphysf, tb, kstar, Z, ecc, BSE_settings):
             ['wd_mass_lim', 'ecsn', 'ecsn_mlow', 'aic', 'ussn', 'sigmadiv', 'bhsigmafrac'],
             ['don_lim', 'acc_lim', 'bdecayfac', 'bconst', 'ck', 'qcflag', 'eddlimflag'],
             ['bhspinflag', 'bhspinmag', 'rejuv_fac', 'rejuvflag', 'htpmb', 'ST_cr'],
-            ['ST_tide', 'rembar_massloss', 'zsun']
+            ['ST_tide', 'rembar_massloss', 'zsun'],
+            ['natal_kick_1', 'theta_1', 'phi_1', 'mean_anomaly_1', 'randomseed_1'],
+            ['natal_kick_2', 'theta_2', 'phi_2', 'mean_anomaly_2', 'randomseed_2']
         ]
 
         for line in lines:
@@ -60,6 +62,14 @@ def convert_initC_row_to_binary_in(initC_file, bin_num):
     initC = pd.read_hdf(initC_file, key="initC")
     r = initC.loc[bin_num]
 
+    kick_cols = ['natal_kick_1', 'theta_1', 'phi_1', 'mean_anomaly_1', 'randomseed_1',
+                 'natal_kick_2', 'theta_2', 'phi_2', 'mean_anomaly_2', 'randomseed_2']
+    for col in kick_cols:
+        if col not in r:
+            BSE_settings[col] = -100.0 if 'randomseed' not in col else 0.0
+        else:
+            BSE_settings[col] = r[col].astype(float)
+
     # update BSE settings with those in the binary
     BSE_settings['idum'] = r['randomseed'].astype(int)
     for key in BSE_settings:

src/cosmic/_version.py

@@ -1 +1 @@
-__version__ = "3.7.4"
+__version__ = "3.7.5"

src/cosmic/evolve.py

@@ -417,6 +417,17 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No
         # is correct
         initial_conditions = initialbinarytable[INITIAL_CONDITIONS_PASS_COLUMNS].to_dict('records')
 
+        # ensure that metallicity is in the valid range (Z in [1e-4, 0.03])
+        low_met_mask = (initialbinarytable["metallicity"] < 1e-4)
+        high_met_mask = (initialbinarytable["metallicity"] > 0.03)
+        if any(low_met_mask | high_met_mask):
+            raise ValueError(
+                f"COSMIC only supports metallicities in the range [1e-4, 0.03]. You have {sum(low_met_mask)} "
+                f"systems with metallicity below 1e-4 and {sum(high_met_mask)} systems with metallicity "
+                "above 0.03. Some examples of problematic binaries have the following bin_nums: "
+                f"{initialbinarytable['bin_num'][low_met_mask | high_met_mask].values[:5]}."
+            )
+
         # we use different columns to save the BSE parameters because some
         # of the parameters are list/arrays which we instead save as
         # individual values because it makes saving to HDF5 easier/more efficient.

src/cosmic/src/evolv2.f

@@ -807,7 +807,8 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
 *
             dspint(k) = 0.d0
             if(((kstar(k).le.9.and.rad(k).ge.0.01d0*rol(k)).or.
-     &         (kstar(k).ge.10.and.k.eq.j1)).and.tflag.gt.0)then
+     &         (kstar(k).ge.10.and.k.eq.j1.and.kstar(k).le.12))
+     &         .and.tflag.gt.0)then
 *
                raa2 = (rad(k)/sep)**2
                raa6 = raa2**3

src/cosmic/src/kick.f

@@ -483,6 +483,12 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2,
       call VectorMagnitude(h, h_mag)
       sep = h_mag * h_mag / (G_const * mtot * (1 - ecc_2)) / rsunkm
 
+* Record the mean anomaly in the arrays
+      kick_info(sn,6) = mean_anom * 180 / pi
+      if (using_cmc.eq.0) then
+         natal_kick_array(sn,4) = mean_anom * 180 / pi
+      endif
+
 * ----------------------------------------------------------------------
 * -------- Split based on whether this kick disrupts the system --------
 * ----------------------------------------------------------------------
@@ -585,11 +591,6 @@ SUBROUTINE kick_pfahl(kw,m1,m1c,m1n,m2,ecc,sep,jorb,vk,sn,r2,
 * ----------------------------------------------------------------------
 * The system is still bound
       else
-* Record the mean anomaly in the arrays
-         kick_info(sn,6) = mean_anom * 180 / pi
-         if (using_cmc.eq.0) then
-            natal_kick_array(sn,4) = mean_anom * 180 / pi
-         endif
 
 * Update the total orbital angular momentum (in Msun Rsun^2/yr)
          jorb = m1n * m2 / mtot * h_mag / rsunkm / rsunkm * yearsc

src/cosmic/src/test_bse.f

@@ -106,6 +106,11 @@ PROGRAM bse
       READ(22,*)don_lim,acc_lim,bdecayfac,bconst,ck,qcflag,eddlimflag
       READ(22,*)bhspinflag,bhspinmag,rejuv_fac,rejuvflag,htpmb,st_cr
       READ(22,*)st_tide,rembar_massloss,zsun
+      READ(22,*)natal_kick_array(1,1),natal_kick_array(1,2),
+     &natal_kick_array(1,3),natal_kick_array(1,4),natal_kick_array(1,5)
+      READ(22,*)natal_kick_array(2,1),natal_kick_array(2,2),
+     &natal_kick_array(2,3),natal_kick_array(2,4),natal_kick_array(2,5)
+
 
       if(kstar(1).lt.0.or.kstar(2).lt.0)then
          READ(22,*)tphys
@@ -150,13 +155,6 @@ PROGRAM bse
 * file. This is a temporary measure until the input file is updated.
 *
         polar_kick_angle = 90
-        do i = 1,4
-            do j = 1,2
-               natal_kick_array(j, i) = -100.d0
-            enddo
-        enddo
-        natal_kick_array(1,5) = 0.0
-        natal_kick_array(2,5) = 0.0
         do i = 1,8
             qcrit_array(i) = 0.0
         enddo
@@ -273,22 +271,22 @@ PROGRAM bse
 * of evolution stage.
 *
  50   j = 0
-      WRITE(*,*)'     TIME      M1       M2   K1 K2        SEP    ECC',  
-     &          '  R1/ROL1 R2/ROL2  TYPE'
- 52   j = j + 1
-      if(bpp(j,1).lt.0.0) goto 60
-      kstar(1) = INT(bpp(j,4))
-      kstar(2) = INT(bpp(j,5))
-      kw = INT(bpp(j,11))
-      WRITE(*,100)(bpp(j,k),k=1,3),kstar,(bpp(j,k),k=6,9),label(kw)
-      goto 52
- 60   continue
- 100  FORMAT(g30.18,2g30.18,2i3,g30.18,g30.18,2g30.18,2x,a8)
-      WRITE(*,*)
+!       WRITE(*,*)'     TIME      M1       M2   K1 K2        SEP    ECC',  
+!      &          '  R1/ROL1 R2/ROL2  TYPE'
+!  52   j = j + 1
+!       if(bpp(j,1).lt.0.0) goto 60
+!       kstar(1) = INT(bpp(j,4))
+!       kstar(2) = INT(bpp(j,5))
+!       kw = INT(bpp(j,11))
+!       WRITE(*,100)(bpp(j,k),k=1,3),kstar,(bpp(j,k),k=6,9),label(kw)
+!       goto 52
+!  60   continue
+!  100  FORMAT(g30.18,2g30.18,2i3,g30.18,g30.18,2g30.18,2x,a8)
+!       WRITE(*,*)
 *
 ************************************************************************
 *
-      WRITE(*,*) bcm(2,31), bpp(j, 6)
+      ! WRITE(*,*) bcm(2,31), bpp(j, 6)
 
       STOP
       END