Making Pre-supernova mass quantities accessible in the BPP file

docs/pages/config/config_insert_bse.html

@@ -404,12 +404,12 @@ <h3 class="name"><code>kickflag</code></h3>
 <p class="description">Sets the particular natal kick prescription to use. Note that <code class="docutils literal notranslate"><span class="pre">sigmadiv</span></code>, <code class="docutils literal notranslate"><span class="pre">bhflag</span></code>, <code class="docutils literal notranslate"><span class="pre">bhsigmafrac</span></code>, <code class="docutils literal notranslate"><span class="pre">aic</span></code>, and <code class="docutils literal notranslate"><span class="pre">ussn</span></code>, which are described below, are only used when <code class="docutils literal notranslate"><span class="pre">abs(kickflag)=1</span></code>. Positive values use the Pfahl+2002 prescription for handling natal kicks.</p>
 <p class="default">Default: 1</p>
 </div>
-<div class="col-3"><select class="form-control"><option selected="true" value="1">1</option><option value="2">2</option><option value="3">3</option><option value="4">4</option><option value="negative values">negative values</option></select></div>
+<div class="col-3"><select class="form-control"><option selected="true" value="1">1</option><option value="2">2</option><option value="3">3</option><option value="4">4</option><option value="5">5</option><option value="negative values">negative values</option></select></div>
 </div>
 <div class="options-expander" style="color: lightcoral;">Option details <i class="fa fa-chevron-down"></i></div>
 <div class="row options hide">
 <p class="options-preface"></p>
-<ul style="margin-left: 2rem; max-width: calc(100% - 2rem)"><li><code class="docutils literal notranslate"><span class="pre opt-val">1</span></code>: <span class="opt-desc">The standard COSMIC kick prescription, where kicks are drawn from a bimodal distribution with standard FeCCSN getting a kick drawn from a Maxwellian distribution with dispersion parameter <code class="docutils literal notranslate"><span class="pre">sigma</span></code> and ECSN/USSN are drawn according to <code class="docutils literal notranslate"><span class="pre">sigmadiv</span></code>. This setting has additional possible options for <code class="docutils literal notranslate"><span class="pre">bhflag</span></code>, <code class="docutils literal notranslate"><span class="pre">bhsigmafrac</span></code>, <code class="docutils literal notranslate"><span class="pre">aic</span></code> and <code class="docutils literal notranslate"><span class="pre">ussn</span></code>.</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">2</span></code>: <span class="opt-desc">Natal kicks are drawn according to <code class="docutils literal notranslate"><span class="pre">sigma</span></code> and scaled by the ejecta mass and remnant mass following Eq. 1 of <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2020ApJ...891..141G/abstract">Giacobbo &amp; Mapelli 2020</a> with their default parameters (\(m_{\rm NS = 1.2 {\rm M_\odot}\), \(m_{\rm ej = 9 {\rm M_\odot}\))</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">3</span></code>: <span class="opt-desc">Natal kicks are drawn according to <code class="docutils literal notranslate"><span class="pre">sigma</span></code> and scaled by just the ejecta mass following Eq. 2 of <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2020ApJ...891..141G/abstract">Giacobbo &amp; Mapelli 2020</a>, which does not scale the kick by (\(m_{\rm NS\)</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">4</span></code>: <span class="opt-desc">Natal kicks are drawn according to Eq. 1 of <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2016MNRAS.461.3747B/abstract">Bray &amp; Eldridge 2016</a>, with their default parameters (\(\alpha=70 \, {\rm km/s}, \beta = 120 \, {\rm km/s)}</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">negative values</span></code>: <span class="opt-desc">Same as above settings but using the old Kiel &amp; Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work</span></li></ul>
+<ul style="margin-left: 2rem; max-width: calc(100% - 2rem)"><li><code class="docutils literal notranslate"><span class="pre opt-val">1</span></code>: <span class="opt-desc">The standard COSMIC kick prescription, where kicks are drawn from a bimodal distribution with standard FeCCSN getting a kick drawn from a Maxwellian distribution with dispersion parameter <code class="docutils literal notranslate"><span class="pre">sigma</span></code> and ECSN/USSN are drawn according to <code class="docutils literal notranslate"><span class="pre">sigmadiv</span></code>. This setting has additional possible options for <code class="docutils literal notranslate"><span class="pre">bhflag</span></code>, <code class="docutils literal notranslate"><span class="pre">bhsigmafrac</span></code>, <code class="docutils literal notranslate"><span class="pre">aic</span></code> and <code class="docutils literal notranslate"><span class="pre">ussn</span></code>.</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">2</span></code>: <span class="opt-desc">Natal kicks are drawn according to <code class="docutils literal notranslate"><span class="pre">sigma</span></code> and scaled by the ejecta mass and remnant mass following Eq. 1 of <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2020ApJ...891..141G/abstract">Giacobbo &amp; Mapelli 2020</a> with their default parameters (\(m_{\rm NS = 1.2 {\rm M_\odot}\), \(m_{\rm ej = 9 {\rm M_\odot}\))</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">3</span></code>: <span class="opt-desc">Natal kicks are drawn according to <code class="docutils literal notranslate"><span class="pre">sigma</span></code> and scaled by just the ejecta mass following Eq. 2 of <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2020ApJ...891..141G/abstract">Giacobbo &amp; Mapelli 2020</a>, which does not scale the kick by (\(m_{\rm NS\)</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">4</span></code>: <span class="opt-desc">Natal kicks are drawn according to Eq. 1 of <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2016MNRAS.461.3747B/abstract">Bray &amp; Eldridge 2016</a>, with their default parameters (\(\alpha=70 \, {\rm km/s}, \beta = 120 \, {\rm km/s)}</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">5</span></code>: <span class="opt-desc">Follows the same prescription as 1, but uses the kick prescription described in <a class="reference external" href="https://ui.adsabs.harvard.edu/abs/2025arXiv250522102D/abstract">Disberg &amp; Mandel 2025</a> for CCSN.</span></li><li><code class="docutils literal notranslate"><span class="pre opt-val">negative values</span></code>: <span class="opt-desc">Same as above settings but using the old Kiel &amp; Hurley 2009 prescription for changing the orbital configuration of the binary, available for reproducibility purposes but not recommended for new work</span></li></ul>
 </div>
 </div><div class="setting">
 <div class="row align-items-center setting-chooser">

docs/pages/inifile.rst

@@ -17,7 +17,7 @@ The buttons below link to the most recent stable and unstable default inifiles f
 .. raw:: html
 
     <div class="toms-nav-container" style="margin-bottom: 3rem; height: 90px; grid-template-rows: 90px">
-        <div class="box" data-href="https://github.com/COSMIC-PopSynth/COSMIC/blob/master/examples/Params.ini">Latest stable INIFILE</div>
+        <div class="box" data-href="https://github.com/COSMIC-PopSynth/COSMIC/blob/v3.6.1/examples/Params.ini">Latest stable INIFILE</div>
         <div class="box" data-href="https://github.com/COSMIC-PopSynth/COSMIC/blob/develop/examples/Params.ini">Latest development INIFILE</div>
     </div>
 

meson.build

@@ -1,7 +1,7 @@
 project('cosmic', 
     'c',
     'fortran',
-    version : '3.6.1',
+    version : '3.6.2',
     default_options: ['warning_level=0', 'optimization=3'],
 )
 
@@ -37,20 +37,16 @@ lib_source = [
     'src/cosmic/src/mrenv.f',
     'src/cosmic/src/ran3.f',
     'src/cosmic/src/rl.f',
+    'src/cosmic/src/hrdiag.f',
+    'src/cosmic/src/star.f',
     'src/cosmic/src/concatkstars.f',
     'src/cosmic/src/comprad.f',
     'src/cosmic/src/bpp_array.f',
     'src/cosmic/src/checkstate.f',
     'src/cosmic/src/deltat.f',
     'src/cosmic/src/mlwind.f',
-    'src/cosmic/src/hrdiag.f',
-    'src/cosmic/src/star.f',
     'src/cosmic/src/zcnsts.f',
     'src/cosmic/src/deltat.f',
-    'src/cosmic/src/mlwind.f',
-    'src/cosmic/src/hrdiag.f',
-    'src/cosmic/src/star.f',
-    'src/cosmic/src/zcnsts.f',
     'src/cosmic/src/zfuncs.f']
 
 # Detect operating system and set appropriate linker flags

requirements.txt

@@ -1,5 +1,5 @@
 scipy >= 0.12.1
-numpy == 1.23.5
+numpy >= 1.26.0
 astropy >= 1.1.1
 configparser
 tqdm >= 4.0

src/cosmic/evolve.py

@@ -49,30 +49,31 @@
 # Make this match the ordering of all_cols in bpp_array.f
 ALL_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb',
                'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1',
-               'tms_2', 'massc_1', 'massc_2', 'rad_1', 'rad_2', 'mass0_1',
+               'tms_2', 'massc_he_layer_1', 'massc_he_layer_2', 'massc_co_layer_1', 'massc_co_layer_2',
+               'rad_1', 'rad_2', 'mass0_1',
                'mass0_2', 'lum_1', 'lum_2', 'teff_1', 'teff_2', 'radc_1',
                'radc_2', 'menv_1', 'menv_2', 'renv_1', 'renv_2', 'omega_spin_1',
                'omega_spin_2', 'B_1', 'B_2', 'bacc_1', 'bacc_2', 'tacc_1',
                'tacc_2', 'epoch_1', 'epoch_2', 'bhspin_1', 'bhspin_2',
-               'deltam_1', 'deltam_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type']
+               'deltam_1', 'deltam_2', 'SN_1', 'SN_2', 'bin_state', 'merger_type', 'metallicity']
 
 INTEGER_COLUMNS = ["bin_state", "bin_num", "kstar_1", "kstar_2", "SN_1", "SN_2", "evol_type"]
 
 
 BPP_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2',
                'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type',
                'aj_1', 'aj_2', 'tms_1', 'tms_2',
-               'massc_1', 'massc_2', 'rad_1', 'rad_2',
+               'massc_he_layer_1', 'massc_he_layer_2', 'massc_co_layer_1', 'massc_co_layer_2', 'rad_1', 'rad_2',
                'mass0_1', 'mass0_2', 'lum_1', 'lum_2', 'teff_1', 'teff_2',
                'radc_1', 'radc_2', 'menv_1', 'menv_2', 'renv_1', 'renv_2',
                'omega_spin_1', 'omega_spin_2', 'B_1', 'B_2', 'bacc_1', 'bacc_2',
                'tacc_1', 'tacc_2', 'epoch_1', 'epoch_2',
                'bhspin_1', 'bhspin_2']
 
 BCM_COLUMNS = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lum_1', 'rad_1',
-               'teff_1', 'massc_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1',
+               'teff_1', 'massc_he_layer_1', 'massc_co_layer_1', 'radc_1', 'menv_1', 'renv_1', 'epoch_1',
                'omega_spin_1', 'deltam_1', 'RRLO_1', 'kstar_2', 'mass0_2', 'mass_2',
-               'lum_2', 'rad_2', 'teff_2', 'massc_2', 'radc_2', 'menv_2',
+               'lum_2', 'rad_2', 'teff_2', 'massc_he_layer_2', 'massc_co_layer_2', 'radc_2', 'menv_2',
                'renv_2', 'epoch_2', 'omega_spin_2', 'deltam_2', 'RRLO_2',
                'porb', 'sep', 'ecc', 'B_1', 'B_2',
                'SN_1', 'SN_2', 'bin_state', 'merger_type']
@@ -368,14 +369,14 @@ def evolve(self, initialbinarytable, pool=None, bpp_columns=None, bcm_columns=No
         for i in range(len(initial_conditions)):
             initial_conditions[i]["n_col_bpp"] = len(bpp_columns)
             initial_conditions[i]["col_inds_bpp"] = col_inds_bpp
-
+        
         # same for bcm
         col_inds_bcm = np.zeros(len(ALL_COLUMNS), dtype=int)
         col_inds_bcm[:len(bcm_columns)] = [ALL_COLUMNS.index(col) + 1 for col in bcm_columns]
         for i in range(len(initial_conditions)):
             initial_conditions[i]["n_col_bcm"] = len(bcm_columns)
             initial_conditions[i]["col_inds_bcm"] = col_inds_bcm
-
+        
         # check if a pool was passed
         if pool is None:
             with MultiPool(processes=nproc) as pool:

src/cosmic/sample/sampler/independent.py

@@ -1158,10 +1158,8 @@ def set_reff(self, mass, metallicity, zsun=0.02):
         of length 10^5.  If your masses are more than that, you'll
         need to divide it into chunks
         """
-
         from cosmic import _evolvebin
 
-
         max_array_size = 100000
         total_length = len(mass)
         radii = np.zeros(total_length)
@@ -1183,7 +1181,7 @@ def set_reff(self, mass, metallicity, zsun=0.02):
 
         length_remaining = total_length
 
-        ## if smaller than 10^5, need to pad out the array
+        # if smaller than 10^5, need to pad out the array
         temp_mass = np.zeros(max_array_size)
         temp_mass[:length_remaining] = mass[-length_remaining:]
 

src/cosmic/src/assign_remnant.f

@@ -4,6 +4,8 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,mt,kw,bhspin,kidx)
       INCLUDE 'const_bse.h'
 
       common /fall/fallback
+*      REAL*8 mass_preSN, mHe_preSN, massc_preSN
+*      COMMON mass_preSN, mHe_preSN, massc_preSN
       REAL*8 fallback
       REAL ran3
       EXTERNAL ran3
@@ -70,6 +72,10 @@ SUBROUTINE assign_remnant(zpars,mc,mcbagb,mass,mt,kw,bhspin,kidx)
                mass = mt
 *
             else
+* Store values in common block
+*               mass_preSN = mt
+*               mHe_preSN = mc_he(kidx)
+*               massc_preSN = mc_co(kidx)
                if(ecsn.gt.0.d0.and.mcbagb.lt.ecsn_mlow)then
 *
 * Star is not massive enough to ignite C burning.

src/cosmic/src/bpp_array.f

@@ -3,26 +3,30 @@ SUBROUTINE WRITETAB(jp,tphys,evolve_type,
      &                      mass1,mass2,kstar1,kstar2,sep,
      &                      tb,ecc,rrl1,rrl2,
      &                      aj1,aj2,tms1,tms2,
-     &                      massc1,massc2,rad1,rad2,
+     &                      massc_he_1,massc_he_2,
+     &                      massc_co_1,massc_co_2, rad1,rad2,
      &                      mass0_1,mass0_2,lumin1,lumin2,
      &                      teff1,teff2,radc1,radc2,menv1,
      &                      menv2,renv1,renv2,ospin1,ospin2,
      &                      b_0_1,b_0_2,bacc1,bacc2,tacc1,tacc2,
      &                      epoch1,epoch2,bhspin1,bhspin2,
      &                      deltam_1,deltam_2,SN_1,SN_2,
-     &                      bin_state,merger_type,tabname)
+     &                      bin_state,merger_type,metallicity,
+     &                      tabname)
         IMPLICIT NONE
         INCLUDE 'const_bse.h'
 
 *
 * Write results to bpp or bcm array.
 *
-*     Author : Scott Coughlin, Tom Wagg
-*     Date :   12th March 2019, September 2024
+*     Author : Scott Coughlin, Tom Wagg, Katie Breivik
+*     Date :   12th March 2019, September 2024, June 2025
 *
-        REAL*8 mass1,mass2
+        REAL*8 mass1,mass2,metallicity
         REAL*8 evolve_type,sep,tb,ecc,tphys,rrl1,rrl2
-        REAL*8 aj1,aj2,tms1,tms2,massc1,massc2,rad1,rad2
+        REAL*8 aj1,aj2,tms1,tms2,rad1,rad2
+        REAL*8 massc_he_1,massc_he_2
+        REAL*8 massc_co_1,massc_co_2
         REAL*8 mass0_1,mass0_2,lumin1,lumin2,radc1,radc2
         REAL*8 menv1,menv2,renv1,renv2,ospin1,ospin2
         REAL*8 b_0_1,b_0_2,bacc1,bacc2,tacc1,tacc2,epoch1,epoch2
@@ -33,7 +37,7 @@ SUBROUTINE WRITETAB(jp,tphys,evolve_type,
         INTEGER jp, col_ind
         INTEGER kstar1,kstar2
         REAL*8 yeardy,aursun,rsunau
-        REAL*8 all_cols(49)
+        REAL*8 all_cols(52)
         CHARACTER*3 tabname
         PARAMETER(yeardy=365.24d0,aursun=214.95d0)
 
@@ -60,40 +64,43 @@ SUBROUTINE WRITETAB(jp,tphys,evolve_type,
         all_cols(13) = aj2
         all_cols(14) = tms1
         all_cols(15) = tms2
-        all_cols(16) = massc1
-        all_cols(17) = massc2
-        all_cols(18) = rad1
-        all_cols(19) = rad2
-        all_cols(20) = mass0_1
-        all_cols(21) = mass0_2
-        all_cols(22) = lumin1
-        all_cols(23) = lumin2
-        all_cols(24) = teff1
-        all_cols(25) = teff2
-        all_cols(26) = radc1
-        all_cols(27) = radc2
-        all_cols(28) = menv1
-        all_cols(29) = menv2
-        all_cols(30) = renv1
-        all_cols(31) = renv2
-        all_cols(32) = ospin1
-        all_cols(33) = ospin2
-        all_cols(34) = b_0_1
-        all_cols(35) = b_0_2
-        all_cols(36) = bacc1
-        all_cols(37) = bacc2
-        all_cols(38) = tacc1
-        all_cols(39) = tacc2
-        all_cols(40) = epoch1
-        all_cols(41) = epoch2
-        all_cols(42) = bhspin1
-        all_cols(43) = bhspin2
-        all_cols(44) = deltam_1
-        all_cols(45) = deltam_2
-        all_cols(46) = float(SN_1)
-        all_cols(47) = float(SN_2)
-        all_cols(48) = bin_state
-        all_cols(49) = merger_type
+        all_cols(16) = massc_he_1
+        all_cols(17) = massc_he_2
+        all_cols(18) = massc_co_1
+        all_cols(19) = massc_co_2
+        all_cols(20) = rad1
+        all_cols(21) = rad2
+        all_cols(22) = mass0_1
+        all_cols(23) = mass0_2
+        all_cols(24) = lumin1
+        all_cols(25) = lumin2
+        all_cols(26) = teff1
+        all_cols(27) = teff2
+        all_cols(28) = radc1
+        all_cols(29) = radc2
+        all_cols(30) = menv1
+        all_cols(31) = menv2
+        all_cols(32) = renv1
+        all_cols(33) = renv2
+        all_cols(34) = ospin1
+        all_cols(35) = ospin2
+        all_cols(36) = b_0_1
+        all_cols(37) = b_0_2
+        all_cols(38) = bacc1
+        all_cols(39) = bacc2
+        all_cols(40) = tacc1
+        all_cols(41) = tacc2
+        all_cols(42) = epoch1
+        all_cols(43) = epoch2
+        all_cols(44) = bhspin1
+        all_cols(45) = bhspin2
+        all_cols(46) = deltam_1
+        all_cols(47) = deltam_2
+        all_cols(48) = float(SN_1)
+        all_cols(49) = float(SN_2)
+        all_cols(50) = bin_state
+        all_cols(51) = merger_type
+        all_cols(52) = metallicity
 
 * check which table we are writing to and write the appropriate columns
         if (tabname .eq. 'bpp') then