diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json
index 4ca0eec6..14286c98 100644
--- a/src/cosmic/data/cosmic-settings.json
+++ b/src/cosmic/data/cosmic-settings.json
@@ -123,7 +123,11 @@
},
{
"name": "martinez26",
- "description": "Uses sana12 for massive binaries (\\(m_1 ≥ 8 {\\rm M_{\\odot}}\\)) with upper limit at 3000 days and raghavan10 for low-mass binaries. Used in Martinez+2026."
+ "description": "Uses sana12 for binaries with a primary massive enough to undergo an FeCCSN with an upper limit at 3000 days and raghavan10 for low mass binaries. The mass threshold is metallicity dependent, and assumes zsun = 0.02, and standard winds. Used in Martinez+2026."
+ },
+ {
+ "name": "martinez26_ecsn",
+ "description": "Same as martinez26 but with the sana12 power law orbital period distribution for primaries massive enough to undergo an ECSN instead of an FeCCSN."
},
{
"name": "custom",
diff --git a/src/cosmic/sample/sampler/independent.py b/src/cosmic/sample/sampler/independent.py
index c135d5fe..80ec3006 100644
--- a/src/cosmic/sample/sampler/independent.py
+++ b/src/cosmic/sample/sampler/independent.py
@@ -75,8 +75,8 @@ def get_independent_sampler(
Model to sample eccentricity; choices include: thermal, uniform, sana12
porb_model : `str` or `dict`
- Model to sample orbital period; choices include: log_uniform, sana12, renzo19, raghavan10, moe19, martinez26
- or a custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
+ Model to sample orbital period; choices include: log_uniform, sana12, renzo19, raghavan10, moe19, martinez26,
+ martinez26_ecsn or a custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
(e.g. {"min": 0.15, "max": 0.55, "slope": -0.55}) would reproduce the Sana+2012 distribution
qmin : `float`
@@ -836,11 +836,14 @@ def sample_porb(self, mass1, mass2, rad1, rad2, porb_model, porb_max=None, size=
`Moe+2019 _`
martinez26 : piecewise model with a power law orbital period following
`Sana+2012 _`
- between 0.15 < log(P/day) < log(3000) for binaries with m1 >= 8Msun and following
+ between 0.15 < log(P/day) < log(3000) for binaries with primaries large enough to undergo an FeCCSN and following
`Raghavan+2010 _`
with a log normal orbital period in days with mean_logP = 4.9 and sigma_logP = 2.3 between
- 0 < log10(P/day) < 9 for binaries with m1 < 8Msun. Used in
- `Martinez+2026 _`.
+ 0 < log10(P/day) < 9 for binaries with low mass primaries. Used in
+ `Martinez+2026 _`. Assumes zsun = 0.02 for the
+ metallicity dependence to predict outcomes correctly.
+ martinez26_ecsn : same as martinez26 but with the Sana+2012 power law orbital period distribution
+ for primaries massive enough to undergo an ECSN instead of an FeCCSN.
Custom power law distribution defined with a dictionary with keys "min", "max", and "slope"
(e.g. porb_model={"min": 0.15, "max": 0.55, "slope": -0.55}) would reproduce the
Sana+2012 distribution.
@@ -1041,13 +1044,83 @@ def get_logP_dist(nsamp, norm_wide, norm_close, mu=4.4, sigma=2.1):
porb = 10**logP_dist
aRL_over_a = a_min / utils.a_from_p(porb,mass1,mass2)
- elif porb_model == "martinez26":
- # martinez+26 model: use sana12 for mass1 >= 8.0 and raghavan10 for mass1 < 8.0
- import scipy
+ elif porb_model == "martinez26" or porb_model == "martinez26_ecsn":
+ # martinez26 model: use sana12 for high mass primaries and raghavan10 for low mass primaries. We fit a metallicity dependent minimum mass for an FeCCSN for 'martinez26,
+ # and we fit for the minmum ecsn mass for 'martinez26_ecsn' using data from the population grid from Martinez+2026.
+ try:
+ met = kwargs.pop('met')
+ except:
+ raise ValueError(
+ "You have chosen martinez26 or martinez26_ecsn for the orbital period distribution which is a metallicity-dependent distribution. "
+ "Please specify a metallicity for the population."
+ )
+ import scipy
+
# Create mask for high-mass and low-mass systems
- (ind_massive,) = np.where(mass1 >= 8.0)
- (ind_lowmass,) = np.where(mass1 < 8.0)
+ def zams_threshold_mass(metallicity, include_ecsn=False, eps=1e-4):
+ z_by_zsun = np.array([
+ 0.005, 0.0061, 0.0074, 0.009, 0.011, 0.01335, 0.01625, 0.0198,
+ 0.0241, 0.02935, 0.03575, 0.0435, 0.05295, 0.0645, 0.0785,
+ 0.09555, 0.1163, 0.1416, 0.1724, 0.20985, 0.25545, 0.311,
+ 0.3786, 0.4609, 0.56105, 0.683, 0.83145, 1.0, 1.2322, 1.5
+ ])
+ fe_ccsne_low = np.array([
+ 6.800704688282865, 6.800897411951605, 6.800897411951605,
+ 6.801021970762962, 6.8025286354545775, 6.802769553556676,
+ 6.803089427396513, 6.803538527040527, 6.804243674222698,
+ 6.80467206993016, 6.80526072924524, 6.8059692527155,
+ 6.806721013425345, 6.81111108784108, 6.917788936421563,
+ 7.011947874916608, 7.08861989201708, 7.149902100414557,
+ 7.204919186852334, 7.238086587027953, 7.329104218495623,
+ 7.421043416504588, 7.514067557855885, 7.667766797129948,
+ 7.714185711362864, 7.866353261910183, 7.938786132196538,
+ 8.03205271544496, 8.222270760441472, 8.266114315589816
+ ]) - eps
+ incl_ecsne_low = np.array([
+ 6.400693921416194, 6.372129640761526, 6.3238994412467235,
+ 6.309580871018385, 6.318731679747097, 6.351048092438493,
+ 6.414738750162605, 6.44081895134046, 6.4754364194354626,
+ 6.520842448126225, 6.579606889898161, 6.626748726227675,
+ 6.697069204683146, 6.81111108784108, 6.917788936421563,
+ 7.011947874916608, 7.08861989201708, 7.149902100414557,
+ 7.204919186852334, 7.198121695356181, 7.221524796053249,
+ 7.350139621902266, 7.4732604008166295, 7.667766797129948,
+ 7.714185711362864, 7.866353261910183, 7.938786132196538,
+ 8.03205271544496, 8.222270760441472, 8.266114315589816
+ ]) - eps
+
+ masses = incl_ecsne_low if include_ecsn else fe_ccsne_low
+
+ #normalize all metallicities with zsun=0.02
+ metallicity_normalized = metallicity / 0.02
+
+ # take boundary value outside range
+ if metallicity_normalized <= z_by_zsun.min():
+ return masses[0]
+ if metallicity_normalized >= z_by_zsun.max():
+ return masses[-1]
+
+ # interpolate in log(Z)
+ logZ = np.log10(z_by_zsun)
+ logZ_target = np.log10(metallicity_normalized)
+
+ # find bracketing indices
+ idx = np.searchsorted(logZ, logZ_target) - 1
+
+ z1, z2 = logZ[idx], logZ[idx + 1]
+ m1, m2 = masses[idx], masses[idx + 1]
+
+ # linear interpolation in logZ
+ frac = (logZ_target - z1) / (z2 - z1)
+ mass_interp = m1 + frac * (m2 - m1)
+
+ return mass_interp
+
+ include_ecsn = porb_model == "martinez26_ecsn"
+ threshold_mass = zams_threshold_mass(met, include_ecsn=include_ecsn)
+ (ind_massive,) = np.where(mass1 >= threshold_mass)
+ (ind_lowmass,) = np.where(mass1 < threshold_mass)
# Initialize porb array
porb = np.zeros(size)
diff --git a/src/cosmic/tests/test_sample.py b/src/cosmic/tests/test_sample.py
index a782a840..d4ba8936 100644
--- a/src/cosmic/tests/test_sample.py
+++ b/src/cosmic/tests/test_sample.py
@@ -261,7 +261,7 @@ def test_msort(self):
def test_sample_porb(self):
# next do Sana12
np.random.seed(4)
- mass1, total_mass = SAMPLECLASS.sample_primary(primary_model='kroupa01', size=100000)
+ mass1, total_mass = SAMPLECLASS.sample_primary(primary_model='kroupa01', size=200000)
mass2 = SAMPLECLASS.sample_secondary(primary_mass = mass1, qmin=0.1)
rad1 = SAMPLECLASS.set_reff(mass=mass1, metallicity=0.02)
rad2 = SAMPLECLASS.set_reff(mass=mass2, metallicity=0.02)
@@ -319,23 +319,47 @@ def test_sample_porb(self):
self.assertTrue(np.round(log_porb_mean, 1) >= MEAN_RAGHAVAN-0.15)
self.assertEqual(np.round(log_porb_sigma, 0), np.round(SIGMA_RAGHAVAN, 0))
- # next check martinez26
+ # next check martinez26 for low mass
+ met = 2e-4 # 1/100 solar
+ feccsn_mass, ecsn_mass = 6.8, 6.4
porb,aRL_over_a = SAMPLECLASS.sample_porb(
- mass1, mass2, rad1, rad2, 'martinez26', size=mass1.size
+ mass1, mass2, rad1, rad2, 'martinez26', size=mass1.size, met=met
)
- # the part of the model with m1 < 8 M_sun should follow the Raghavan10 distribution
- porb_low_mass = porb[mass1 < 8]
+ # the part of the model with m1 < 6.8 M_sun should follow the Raghavan10 distribution
+ porb_low_mass = porb[mass1 < feccsn_mass]
log_porb_mean = np.mean(np.log10(porb_low_mass))
log_porb_sigma = np.std(np.log10(porb_low_mass))
self.assertTrue(np.round(log_porb_mean, 1) >= MEAN_RAGHAVAN-0.15)
self.assertEqual(np.round(log_porb_sigma, 0), np.round(SIGMA_RAGHAVAN, 0))
- # the part of the model with m1 >= 8 M_sun should follow the same power law as Sana12
- m1_high = mass1+8
- rad1_high = SAMPLECLASS.set_reff(mass=m1_high, metallicity=0.02)
+
+ # check martinez26_ecsn for low mass
+ porb,aRL_over_a = SAMPLECLASS.sample_porb(
+ mass1, mass2, rad1, rad2, 'martinez26_ecsn', size=mass1.size, met=met
+ )
+ # the part of the model with m1 < 5.3 M_sun should follow the Raghavan10 distribution
+ porb_low_mass = porb[mass1 < ecsn_mass]
+ log_porb_mean = np.mean(np.log10(porb_low_mass))
+ log_porb_sigma = np.std(np.log10(porb_low_mass))
+ self.assertTrue(np.round(log_porb_mean, 1) >= MEAN_RAGHAVAN-0.15)
+ self.assertEqual(np.round(log_porb_sigma, 0), np.round(SIGMA_RAGHAVAN, 0))
+
+ # check martinez_26 for high mass
+ m1_high = mass1+feccsn_mass
+ rad1_high = SAMPLECLASS.set_reff(mass=m1_high, metallicity=met)
+ porb,aRL_over_a = SAMPLECLASS.sample_porb(
+ m1_high, mass2, rad1_high, rad2, 'martinez26', size=m1_high.size, met=met
+ )
+ porb_high_mass = porb[(m1_high >= feccsn_mass) & (np.log10(porb) > 0.5)]
+ power_slope = power_law_fit(np.log10(porb_high_mass), n_bins=25)
+ self.assertEqual(np.round(power_slope, 2), SANA12_PORB_POWER_LAW)
+
+ # check martinez_26_ecsn for high mass
+ m1_high = mass1+ecsn_mass
+ rad1_high = SAMPLECLASS.set_reff(mass=m1_high, metallicity=met)
porb,aRL_over_a = SAMPLECLASS.sample_porb(
- m1_high, mass2, rad1_high, rad2, 'martinez26', size=m1_high.size
+ m1_high, mass2, rad1_high, rad2, 'martinez26_ecsn', size=m1_high.size, met=met
)
- porb_high_mass = porb[(m1_high >= 8) & (np.log10(porb) > 0.5)]
+ porb_high_mass = porb[(m1_high >= ecsn_mass) & (np.log10(porb) > 0.5)]
power_slope = power_law_fit(np.log10(porb_high_mass), n_bins=25)
self.assertEqual(np.round(power_slope, 2), SANA12_PORB_POWER_LAW)