IFMR N20 Z

spisea/ifmr.py

@@ -45,7 +45,187 @@ def Kalirai_mass(self, MZAMS):
 
         return final
 
+
+class IFMR_N20_Sukhbold(IFMR):
+    """
+    BH/NS IFMR based on Sukhbold & Woosley 2014 for zero-Z models:
+        https://ui.adsabs.harvard.edu/abs/2014ApJ...783...10S/abstract
+    BH/NS IFMR based on Sukhbold et al. 2016 for solar-Z models::
+        https://ui.adsabs.harvard.edu/abs/2016ApJ...821...38S/abstract
+    PPISN based on Woosley 2017: 
+        https://ui.adsabs.harvard.edu/abs/2017ApJ...836..244W/abstract
+    WD IFMR from Kalirai et al. 2008:
+        https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract
+    """
+    # Linear fits to Sukhbold simulations.
+    zero_coeff = [0.46522639, -3.29170817]
+    solar_coeff = [-0.27079245, 24.74320755]
+
+    zero_BH_mass = np.poly1d(zero_coeff)
+    solar_BH_mass = np.poly1d(solar_coeff)
+
+    # Solar metallicity (what Sam is using)
+    Zsun = 0.014
+
+    def NS_mass(self, MZAMS):
+        """                                                                                                      
+        Paper: 9 < MZAMS 120                                                                                     
+        Drawing the mass from gaussian created using observational data
+         Done by Emily Ramey and Sergiy Vasylyev at University of Caifornia, Berkeley
+        sample oF NS from Ozel & Freire (2016) — J1811+2405 Ng et al. (2020), 
+        J2302+4442 Kirichenko et al. (2018), J2215+5135 Linares et al. (2018), 
+        J1913+1102 Ferdman & Collaboration (2017), J1411+2551 Martinez et al. (2017), 
+        J1757+1854 Cameron et al. (2018), J0030+0451 Riley et al. (2019), J1301+0833 Romani et al. (2016)
+        The gaussian distribution was created using this data and a Bayesian MCMC method adapted from
+        Kiziltan et al. (2010)
 
+        """
+        return np.random.normal(loc=1.36, scale=0.09, size=len(MZAMS))
+
+
+    def BH_mass_low(self, MZAMS):
+        """
+        9 < MZAMS < 40 Msun
+        """
+        mBH = zero_BH_mass(MZAMS)
+
+        return mBH
+
+
+    def BH_mass_high(self, MZAMS, Z):
+        """
+        39.6 Msun < MZAMS < 120 Msun
+        """
+        zfrac = Z/Zsun
+
+        # super-solar Z gives identical results as solar Z
+        if zfrac > 1:
+            zfrac = 1
+
+        if zfrac < 0:
+            raise ValueError('Z must be non-negative.')
+
+        # Linearly interpolate
+        mBH = (1 - zfrac) * zero_BH_mass(MZAMS) + zfrac*solar_BH_mass(MZAMS)
+
+        return mBH
+
+
+    def prob_BH_high(self, Z):
+        """
+        Probability of BH formation for 60 < Mzams < 120 Msun
+        """
+        zfrac = Z/Zsun
+
+        if Zfrac > 1:
+            Zfrac = 1
+
+        if zfrac < 0:
+            raise ValueError('Z must be non-negative.')
+
+        pBH = 1 - 0.8*zfrac
+
+        return pBH
+
+
+    def generate_death_mass(self, mass_array, metallicity_array):
+        """
+        The top-level function that assigns the remnant type 
+        and mass based on the stellar initial mass. 
+
+        Parameters
+        ----------
+        mass_array: array of floats
+            Array of initial stellar masses. Units are
+            M_sun.
+        metallicity_array: array of floats
+            Array of metallicities in terms of [Fe/H]
+        Notes
+        ------
+        The output typecode tells what compact object formed:
+
+        * WD: typecode = 101
+        * NS: typecode = 102
+        * BH: typecode = 103
+        A typecode of value -1 means you're outside the range of 
+        validity for applying the ifmr formula. 
+        A remnant mass of -99 means you're outside the range of 
+        validity for applying the ifmr formula.
+        Range of validity: MZAMS > 0.5 
+
+        Returns
+        -------
+        output_arr: 2-element array
+            output_array[0] contains the remnant mass, and 
+            output_array[1] contains the typecode
+        """
+        #output_array[0] holds the remnant mass
+        #output_array[1] holds the remnant type
+        output_array = np.zeros((2, len(mass_array)))
+
+        codes = {'WD': 101, 'NS': 102, 'BH': 103}
+
+        # Array to store the remnant masses
+        rem_mass_array = np.zeros(len(mass_array))
+
+        # Convert from [Fe/H] to Z
+        # FIXME: if have Fe/H = nan that makes Z = 0. Is that the behavior we want?
+        Z_array = np.zeros((len(metallicity_array)))
+        metal_idx = np.where(metallicity_array != np.nan)
+        Z_array[metal_idx] = self.get_Z(metallicity_array[metal_idx])
+
+        # Random array to get probabilities for what type of object will form
+        random_array = np.random.randint(1, 101, size = len(mass_array))
+
+        id_array0 = np.where((mass_array < 0.5) | (mass_array >= 120))
+        output_array[0][id_array0] = -99 * np.ones(len(id_array0))
+        output_array[1][id_array0]  = -1 * np.ones(len(id_array0))
+
+        id_array1 = np.where((mass_array >= 0.5) & (mass_array < 9))
+        output_array[0][id_array1] = self.Kalirai_mass(mass_array[id_array1])
+        output_array[1][id_array1]= codes['WD']
+
+        id_array2 = np.where((mass_array >= 9) & (mass_array < 15))
+        output_array[0][id_array2] = self.NS_mass(mass_array[id_array2])
+        output_array[1][id_array2] = codes['NS']
+
+        id_array3_BH = np.where((mass_array >= 15) & (mass_array < 21.8) & (random_array > 75))
+        output_array[0][id_array3_BH] = self.BH_mass_low(mass_array[id_array3_BH])
+        output_array[1][id_array3_BH] = codes['BH']
+
+        id_array3_NS = np.where((mass_array >= 15) & (mass_array < 21.8) & (random_array <= 75))
+        output_array[0][id_array3_NS] = self.NS_mass(mass_array[id_array3_NS])
+        output_array[1][id_array3_NS] = codes['NS']
+
+        id_array4 = np.where((mass_array >= 21.8) & (mass_array < 25.2))
+        output_array[0][id_array4] = self.BH_mass_low(mass_array[id_array4])
+        output_array[1][id_array4] = codes['BH']
+
+        id_array5 = np.where((mass_array >= 25.2) & (mass_array < 27.4))
+        output_array[0][id_array5] = self.NS_mass(mass_array[id_array5])
+        output_array[1][id_array5] = codes['NS']
+
+        id_array6 = np.where((mass_array >= 27.4) & (mass_array < 39.6))
+        output_array[0][id_array6] = self.BH_mass_low(mass_array[id_array6])
+        output_array[1][id_array6] = codes['BH']
+
+        id_array7 = np.where((mass_array >= 39.6) & (mass_array < 60))
+        output_array[0][id_array7] = self.BH_mass_high(mass_array[id_array7],
+                                                       Z_array[id_array7])
+        output_array[1][id_array7] = codes['BH']
+
+        id_array8 = np.where((mass_array >= 60) & (mass_array < 120))
+        for idx in id_array8:
+            pBH = prob_BH_high(Z_array[id_array8][idx])
+            if random_array[id_array8][idx] > 100*pBH:
+                output_array[0][id_array8][idx] = self.BH_mass_high(mass_array[id_array8][idx],
+                                                                    Z_array[id_array8][idx])
+                output_array[1][id_array8][idx] = codes['BH']
+            else:
+                output_array[0][id_array8][idx] = self.NS_mass(mass_array[id_array8][idx])
+                output_array[1][id_array8][idx] = codes['NS']
+
+        return(output_array)
 
 
 class IFMR_Spera15(IFMR):