From 3b2d9aa1ad84a9eaed2fe2ff443ac1e1f6774d82 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Thu, 29 Oct 2020 09:49:42 -0700 Subject: [PATCH 1/8] start of new ifmr class --- spisea/ifmr.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index b8a03ee0..f8451e8b 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -45,8 +45,18 @@ def Kalirai_mass(self, MZAMS): return final - +class IFMR_N20_Sukhbold(IFMR): + """ + BH/NS IFMR comes from Tuguldur Sukhbold + Based on Sukhbold & Woosley 2014: + https://ui.adsabs.harvard.edu/abs/2014ApJ...783...10S/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 + """ + pass class IFMR_Spera15(IFMR): """ From caebbeedf729d204b12d78e5031b4e8820b82e19 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Thu, 29 Oct 2020 19:31:33 -0700 Subject: [PATCH 2/8] work in progress --- spisea/ifmr.py | 113 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 112 insertions(+), 1 deletion(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index f8451e8b..4f720c99 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -56,7 +56,118 @@ class IFMR_N20_Sukhbold(IFMR): WD IFMR from Kalirai et al. 2008: https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract """ - pass + # Comes from fit (needs a function) + zero_coeff = [ 1.73877791e-06, -2.29152875e-04, 1.20555396e-02, 1.35831935e-01, 2.31381928e-01] + solar_coeff = [ 7.22032444e-03, -1.07320209e+00, 4.61793972e+01] + zfrac_cut = 0.3 + Zsun = 0.014 # Following what Sam has + + def BH_mass_1(self, MZAMS): + """ + 9 < MZAMS < 40 Msun + """ + BH_mass = np.poly1d(solar_coeff) + mBH = BH_mass(MZAMS) + + return mBH + + + def BH_mass_2(self, MZAMS, Z): + """ + 40 Msun < MZAMS < ... + """ + solar_BH_mass = np.poly1d(solar_coeff) + zero_BH_mass = np.poly1d(zero_coeff) + + zfrac = Z/Zsun + + # super-solar Z gives identical results as solar Z + if zfrac > 1: + zfrac = 1 + + # Linearly interpolate + mBH = solar_BH_mass(MZAMS) - zfrac * zero_BH_mass(MZAMS) + + return mBH + + + def BH_mass_PPISN_1(self, MZAMS): + """ + For 70 < MZAMS < 120 at Z = 0 + """ + return 0.4 * MZAMS + 4 + + + def BH_mass_PPISN_2(self, MZAMS): + """ + For 120 < MZAMS < 140 at Z = 0 + """ + return -0.75 * MZAMS + 142 + + + 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, 1001, size = len(mass_array)) + + id_array0 = np.where(...) # FIXME + 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'] + + # Need probabilities... + class IFMR_Spera15(IFMR): """ From 3f440050b64449bdc1bb805540e7c3671380b7c3 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Mon, 2 Nov 2020 16:11:26 -0800 Subject: [PATCH 3/8] more work on ifmr --- spisea/ifmr.py | 101 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 79 insertions(+), 22 deletions(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index 4f720c99..d8b2575d 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -56,13 +56,32 @@ class IFMR_N20_Sukhbold(IFMR): WD IFMR from Kalirai et al. 2008: https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract """ - # Comes from fit (needs a function) - zero_coeff = [ 1.73877791e-06, -2.29152875e-04, 1.20555396e-02, 1.35831935e-01, 2.31381928e-01] - solar_coeff = [ 7.22032444e-03, -1.07320209e+00, 4.61793972e+01] - zfrac_cut = 0.3 - Zsun = 0.014 # Following what Sam has + # Comes from fit (needs a function here) + 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) + + Zsun = 0.014 # What sam is using - def BH_mass_1(self, MZAMS): + 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 """ @@ -72,9 +91,9 @@ def BH_mass_1(self, MZAMS): return mBH - def BH_mass_2(self, MZAMS, Z): + def BH_mass_high(self, MZAMS, Z): """ - 40 Msun < MZAMS < ... + 39.6 Msun < MZAMS < 120 Msun """ solar_BH_mass = np.poly1d(solar_coeff) zero_BH_mass = np.poly1d(zero_coeff) @@ -85,25 +104,30 @@ def BH_mass_2(self, MZAMS, Z): if zfrac > 1: zfrac = 1 + if zfrac < 0: + raise ValueError('Z must be non-negative.') + # Linearly interpolate - mBH = solar_BH_mass(MZAMS) - zfrac * zero_BH_mass(MZAMS) + # CHECK + mBH = zfrac*solar_BH_mass(MZAMS) + (1 - zfrac) * zero_BH_mass(MZAMS) return mBH + def prob_BH_high(self, Z): + zfrac = Z/Zsun - def BH_mass_PPISN_1(self, MZAMS): - """ - For 70 < MZAMS < 120 at Z = 0 - """ - return 0.4 * MZAMS + 4 + if Zfrac > 1: + Zfrac = 1 + if zfrac < 0: + raise ValueError('Z must be non-negative.') - def BH_mass_PPISN_2(self, MZAMS): - """ - For 120 < MZAMS < 140 at Z = 0 - """ - return -0.75 * MZAMS + 142 + pBH_sun = 0.8 + pBH_zero = 1.0 + pBH = zfrac*pBH_sun + (1-zfrac)*pBH_zero + + return pBH def generate_death_mass(self, mass_array, metallicity array): """ @@ -154,7 +178,7 @@ def generate_death_mass(self, mass_array, metallicity array): # Random array to get probabilities for what type of object will form random_array = np.random.randint(1, 1001, size = len(mass_array)) - id_array0 = np.where(...) # FIXME + id_array0 = np.where((mass_array < 0.5) | (mass_array >= 120)) # FIXME output_array[0][id_array0] = -99 * np.ones(len(id_array0)) output_array[1][id_array0] = -1 * np.ones(len(id_array0)) @@ -166,8 +190,41 @@ def generate_death_mass(self, mass_array, metallicity array): output_array[0][id_array2] = self.NS_mass(mass_array[id_array2]) output_array[1][id_array2] = codes['NS'] - # Need probabilities... - + id_array3_BH = np.where((mass_array >= 15) & (mass_array < 21.8) & (random_array > 750)) + 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 <= 750)) + 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] > 1000*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'] class IFMR_Spera15(IFMR): """ From ccc8d84f3590efdfac2932b7a38f8b0e907f548f Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Mon, 2 Nov 2020 16:18:19 -0800 Subject: [PATCH 4/8] parentheses --- spisea/ifmr.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index d8b2575d..7eac5054 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -113,6 +113,7 @@ def BH_mass_high(self, MZAMS, Z): return mBH + def prob_BH_high(self, Z): zfrac = Z/Zsun @@ -129,6 +130,7 @@ def prob_BH_high(self, Z): return pBH + def generate_death_mass(self, mass_array, metallicity array): """ The top-level function that assigns the remnant type @@ -226,6 +228,9 @@ def generate_death_mass(self, mass_array, metallicity array): 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): """ This IFMR comes from Spera et. al. 2015 Appendix C (Used for all MZAMS>= 7 M_sun) From a78ebf55d7e8bac22a2439d5dae316bd5756d3b4 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Mon, 2 Nov 2020 19:20:40 -0800 Subject: [PATCH 5/8] some changes to be consistent with paper draft --- spisea/ifmr.py | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index 7eac5054..80448db4 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -56,14 +56,15 @@ class IFMR_N20_Sukhbold(IFMR): WD IFMR from Kalirai et al. 2008: https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract """ - # Comes from fit (needs a function here) + # 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) - Zsun = 0.014 # What sam is using + # Solar metallicity (what Sam is using) + Zsun = 0.014 def NS_mass(self, MZAMS): """ @@ -85,8 +86,7 @@ def BH_mass_low(self, MZAMS): """ 9 < MZAMS < 40 Msun """ - BH_mass = np.poly1d(solar_coeff) - mBH = BH_mass(MZAMS) + mBH = zero_BH_mass(MZAMS) return mBH @@ -95,9 +95,6 @@ def BH_mass_high(self, MZAMS, Z): """ 39.6 Msun < MZAMS < 120 Msun """ - solar_BH_mass = np.poly1d(solar_coeff) - zero_BH_mass = np.poly1d(zero_coeff) - zfrac = Z/Zsun # super-solar Z gives identical results as solar Z @@ -108,13 +105,15 @@ def BH_mass_high(self, MZAMS, Z): raise ValueError('Z must be non-negative.') # Linearly interpolate - # CHECK - mBH = zfrac*solar_BH_mass(MZAMS) + (1 - zfrac) * zero_BH_mass(MZAMS) + 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: @@ -123,10 +122,7 @@ def prob_BH_high(self, Z): if zfrac < 0: raise ValueError('Z must be non-negative.') - pBH_sun = 0.8 - pBH_zero = 1.0 - - pBH = zfrac*pBH_sun + (1-zfrac)*pBH_zero + pBH = 1 - 0.8*zfrac return pBH @@ -178,9 +174,9 @@ def generate_death_mass(self, mass_array, metallicity array): 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, 1001, size = len(mass_array)) + random_array = np.random.randint(1, 101, size = len(mass_array)) - id_array0 = np.where((mass_array < 0.5) | (mass_array >= 120)) # FIXME + 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)) @@ -192,11 +188,11 @@ def generate_death_mass(self, mass_array, metallicity array): 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 > 750)) + 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 <= 750)) + 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'] @@ -220,7 +216,7 @@ def generate_death_mass(self, mass_array, metallicity array): 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] > 1000*pBH: + 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'] From d11de9ce2a5d51c569c6e01ca70ed95aa464d2b6 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Mon, 2 Nov 2020 19:30:41 -0800 Subject: [PATCH 6/8] references --- spisea/ifmr.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index 80448db4..81c17d0a 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -49,11 +49,13 @@ def Kalirai_mass(self, MZAMS): class IFMR_N20_Sukhbold(IFMR): """ BH/NS IFMR comes from Tuguldur Sukhbold - Based on Sukhbold & Woosley 2014: + Based on Sukhbold & Woosley 2014 for zero-Z models: https://ui.adsabs.harvard.edu/abs/2014ApJ...783...10S/abstract - PPISN based on Woosley 2017: + 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 for PPISN: https://ui.adsabs.harvard.edu/abs/2017ApJ...836..244W/abstract - WD IFMR from Kalirai et al. 2008: + WD IFMR from Kalirai et al. 2008 for white dwarfs: https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract """ # Linear fits to Sukhbold simulations. From de63ac3492c9f0193dc3e259059f088f3cbcc249 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Mon, 2 Nov 2020 19:33:48 -0800 Subject: [PATCH 7/8] typos --- spisea/ifmr.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index 81c17d0a..0663d171 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -48,15 +48,14 @@ def Kalirai_mass(self, MZAMS): class IFMR_N20_Sukhbold(IFMR): """ - BH/NS IFMR comes from Tuguldur Sukhbold - Based on Sukhbold & Woosley 2014 for zero-Z models: - https://ui.adsabs.harvard.edu/abs/2014ApJ...783...10S/abstract - 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 for PPISN: - https://ui.adsabs.harvard.edu/abs/2017ApJ...836..244W/abstract - WD IFMR from Kalirai et al. 2008 for white dwarfs: - https://ui.adsabs.harvard.edu/abs/2008ApJ...676..594K/abstract + 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] @@ -124,9 +123,9 @@ def prob_BH_high(self, Z): if zfrac < 0: raise ValueError('Z must be non-negative.') - pBH = 1 - 0.8*zfrac + pBH = 1 - 0.8*zfrac - return pBH + return pBH def generate_death_mass(self, mass_array, metallicity array): From 0ca5619493e22c699d3e20cbf89c63e950bb85d6 Mon Sep 17 00:00:00 2001 From: Casey Lam Date: Mon, 9 Nov 2020 19:22:13 -0800 Subject: [PATCH 8/8] add missing underscore --- spisea/ifmr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spisea/ifmr.py b/spisea/ifmr.py index 0663d171..25ec8086 100755 --- a/spisea/ifmr.py +++ b/spisea/ifmr.py @@ -128,7 +128,7 @@ def prob_BH_high(self, Z): return pBH - def generate_death_mass(self, mass_array, metallicity array): + 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.