diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 6731570134d11930f287092cce8a185cf53b43e9..e872326b5eea8260a24123be5cf034df4d285b4f 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -181,8 +181,8 @@ def powerlaw(
# Handle faulty value
if k == -1:
- print("wrong value for k")
- raise ValueError
+ msg = "wrong value for k"
+ raise ValueError(msg)
if (x < min_val) or (x > max_val):
print("input value is out of bounds!")
@@ -590,8 +590,8 @@ def imf_chabrier2003(m: Union[int, float]) -> Union[int, float]:
chabrier_a2 = 4.43e-2
chabrier_x = -1.3
if m <= 0:
- print("below bounds")
- raise ValueError
+ msg = "below bounds"
+ raise ValueError(msg)
if 0 < m < 1.0:
A = 0.158
dm = math.log10(m) - chabrier_logmc
@@ -1010,8 +1010,14 @@ def _poisson(lambda_val, n):
def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
"""
Function that returns a list of multiplicity fractions for a given list
- TODO: make an extrapolation functionality in this. log10(1.6e1) is low, we can probably go a bit further
+
+ TODO: make an extrapolation functionality in this. log10(1.6e1)
+ is low, we can probably go a bit further
+
+ The default result that is returned when sampling the mass outside
+ of the mass range is now the last known value
"""
+
# Returns a list of multiplicity fractions for a given input list of masses
# Use the global Moecache
@@ -1019,7 +1025,6 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
result = {}
- # TODO: decide what to do with this.
# if(0)
# {
# # we have these in a table, so could interpolate, but...
@@ -1058,6 +1063,12 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
else 0
)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: using model {}".format("Poisson"),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
elif options["multiplicity_model"] == "data":
# use the fractions calculated from Moe's data directly
#
@@ -1072,10 +1083,23 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
* options["multiplicity_modulator"][n]
)
result[3] = 0.0 # no quadruples
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: using model {}".format("data"),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
# Normalisation:
if options["normalize_multiplicities"] == "raw":
+ # TODO: Not sure what actually happens here. It eems the same as above
# do nothing : use raw Poisson predictions
+
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Not normalizing (using raw results)",
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
pass
elif options["normalize_multiplicities"] == "norm":
@@ -1084,6 +1108,12 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
for key in result.keys():
result[key] = result[key] / sum_result
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Normalizing with {}".format("norm"),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
elif options["normalize_multiplicities"] == "merge":
# if multiplicity == 1, merge binaries, triples and quads into singles
# (de facto this is the same as "norm")
@@ -1093,24 +1123,18 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
# TODO: ask rob about this part. Not sure if i understand it.
- multiplicity = 0
+ max_multiplicity = 0
for n in range(4):
if options["multiplicity_modulator"][n] > 0:
- multiplicity = n + 1
+ max_multiplicity = n + 1
- verbose_print(
- "\tM&S: Multiplicity: {}: {}".format(multiplicity, str(result)),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
-
- if multiplicity == 1:
+ if max_multiplicity == 1:
result[0] = 1.0
for n in range(1, 4):
result[0] += result[n]
result[n] = 0
- elif multiplicity == 2:
+ elif max_multiplicity == 2:
# we have only singles and binaries:
# triples and quads are treated as binaries
@@ -1118,19 +1142,19 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
result[1] += result[n]
result[n] = 0
- elif multiplicity == 3:
+ elif max_multiplicity == 3:
# we have singles, binaries and triples:
# quads are treated as triples
result[2] += result[3]
result[3] = 0
- elif multiplicity == 4:
+ elif max_multiplicity == 4:
# we have singles, binaries, triples and quads,
# so do nothing
pass
else:
msg = "\tM&S: Error: in the Moe distribution, we seem to want no stars at all (multiplicity == {})".format(
- multiplicity
+ max_multiplicity
)
verbose_print(
msg,
@@ -1139,10 +1163,17 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
)
raise ValueError(msg)
+ # Normalising results again
sum_result = sum([result[key] for key in result.keys()])
for key in result.keys():
result[key] = result[key] / sum_result
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Normalizing with {}, max_multiplicity={}".format("merge", max_multiplicity),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
verbose_print(
"\tM&S: Multiplicity array: {}".format(str(result)),
verbosity,
@@ -1306,8 +1337,8 @@ def build_q_table(options, m, p, verbosity=0):
qs = sorted(qdata.keys())
if len(qs) == 0:
- print("No qs found error")
- raise ValueError
+ msg = "No qs found error"
+ raise ValueError(msg)
elif len(qs) == 1:
# only one q value : pretend there are two
@@ -1400,14 +1431,15 @@ def build_q_table(options, m, p, verbosity=0):
_MS_VERBOSITY_LEVEL,
)
else:
- verbose_print(
- "\tM&S: build_q_table: Error no other methods available. The chosen method ({}) does not exist!".format(
+ msg = "\tM&S: build_q_table: Error no other methods available. The chosen method ({}) does not exist!".format(
method
- ),
+ )
+ verbose_print(
+ msg,
verbosity,
_MS_VERBOSITY_LEVEL,
)
- raise ValueError
+ raise ValueError(msg)
# TODO: consider implementing the poly method (see perl version)
@@ -1536,14 +1568,15 @@ def get_integration_constant_q(q_interpolator, verbosity=0):
for q in np.arange(0, 1 + 2e-6, dq):
x = q_interpolator.interpolate([q])
if len(x) == 0:
- verbose_print(
- "\tM&S: build_q_table: Q interpolator table interpolation failed.\n\t\ttmp_table = {}\n\t\tq_data = {}".format(
+ msg = "\tM&S: build_q_table: Q interpolator table interpolation failed.\n\t\ttmp_table = {}\n\t\tq_data = {}".format(
str(tmp_table), str(qdata)
- ),
+ )
+ verbose_print(
+ msg,
verbosity,
_MS_VERBOSITY_LEVEL,
)
- raise ValueError
+ raise ValueError(msg)
else:
I += x[0] * dq
# verbose_print(
@@ -1580,6 +1613,57 @@ def fill_data(sample_values, data_dict):
return data
+def calc_e_integral(
+ options, integrals_string, interpolator_name, mass_string, period_string, verbosity=0
+):
+ """
+ Function to calculate the P integral
+
+ We need to renormalize this because min_per > 0, and not all periods should be included
+ """
+
+ global Moecache
+ min_ecc = 0
+ max_ecc = 0.9999
+
+ mass_period_string = "{}_{}".format(options[mass_string], options[period_string])
+
+ # Check if the dict exists
+ if not Moecache.get(integrals_string, None):
+ Moecache[integrals_string] = {}
+
+ # Check for cached value. If it doesn't exist: calculate
+ if not Moecache[integrals_string].get(mass_period_string, None):
+ I = 0
+ decc = 1e-3
+
+ for ecc in np.arange(min_ecc, max_ecc, decc):
+ # Loop over all the values in the table, between the min and max P
+ dp_decc = Moecache[interpolator_name].interpolate(
+ [np.log10(options[mass_string]), np.log10(options[period_string]), ecc]
+ )[0]
+
+ I += dp_decc * decc
+
+ # Set the integral value in the dict
+ Moecache[integrals_string][mass_period_string] = I
+ verbose_print(
+ "\tM&S: calc_ecc_integral: min_ecc: {} max ecc: {} integrals_string: {} interpolator_name: {} mass_string: {} period_string: {} mass: {} period: {} I: {}".format(
+ min_ecc, max_ecc, integrals_string, interpolator_name, mass_string, period_string, options[mass_string], options[period_string], I
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ else:
+ verbose_print(
+ "\tM&S: calc_ecc_integral: Found cached value for min_ecc: {} max ecc: {} integrals_string: {} interpolator_name: {} mass_string: {} period_string: {} mass: {} period: {} I: {}".format(
+ min_ecc, max_ecc, integrals_string, interpolator_name, mass_string, period_string, options[mass_string], options[period_string], Moecache[integrals_string][mass_period_string]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+
def calc_P_integral(
options, min_P, integrals_string, interpolator_name, mass_string, verbosity=0
):
@@ -1592,32 +1676,58 @@ def calc_P_integral(
global Moecache
max_logP = 10
- I = 0
- dlogP = 1e-3
+ # Check if the dict exists
+ if not Moecache.get(integrals_string, None):
+ Moecache[integrals_string] = {}
- for logP in np.arange(np.log10(min_P), max_logP, dlogP):
- # Loop over all the values in the table, between the min and max P
- dp_dlogP = Moecache[interpolator_name].interpolate(
- [np.log10(options[mass_string]), logP]
- )[0]
+ # Check for cached value. If it doesn't exist: calculate
+ if not Moecache[integrals_string].get(options[mass_string], None):
+ I = 0
+ dlogP = 1e-3
- I += dp_dlogP * dlogP
- Moecache[integrals_string][options[mass_string]] = I
- verbose_print(
- "\tM&S: calc_P_integral: min_P: {} integrals_string: {} interpolator_name: {} mass_string: {} I: {}".format(
- min_P, integrals_string, interpolator_name, mass_string, I
- ),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
+ for logP in np.arange(np.log10(min_P), max_logP, dlogP):
+ # Loop over all the values in the table, between the min and max P
+ dp_dlogP = Moecache[interpolator_name].interpolate(
+ [np.log10(options[mass_string]), logP]
+ )[0]
+ I += dp_dlogP * dlogP
+
+ # Set the integral value in the dict
+ Moecache[integrals_string][options[mass_string]] = I
+ verbose_print(
+ "\tM&S: calc_P_integral: min_P: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
+ min_P, integrals_string, interpolator_name, mass_string, options[mass_string], I
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ else:
+ verbose_print(
+ "\tM&S: calc_P_integral: Found cached value for min_P: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
+ min_P, integrals_string, interpolator_name, mass_string, options[mass_string], Moecache[integrals_string][options[mass_string]]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+def calc_total_probdens(prob_dict):
+ """
+ Function to calculate the total probability density
+ """
+
+ total_probdens = 1
+ for key in prob_dict:
+ total_probdens *= prob_dict[key]
+ prob_dict['total_probdens'] = total_probdens
+
+ return prob_dict
def Moe_de_Stefano_2017_pdf(options, verbosity=0):
"""
Moe & distefano function to calculate the probability density.
-
- takes a dctionary as input (in options) with options:
+ takes a dictionary as input (in options) with options:
M1, M2, M3, M4 => masses (Msun) [M1 required, rest optional]
P, P2, P3 => periods (days) [number: none=binary, 2=triple, 3=quadruple]
@@ -1635,18 +1745,18 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
_MS_VERBOSITY_LEVEL,
)
- prob = [] # Value that we will return
prob_dict = {} # Dictionary containing all the pdf values for the different parameters
# Get the multiplicity from the options, and if its not there, calculate it based on the
# TODO: the function below makes no sense. We NEED to pass the multiplicity in the
if not options.get("multiplicity", None):
+ msg = "\tM&S: Moe_de_Stefano_2017_pdf: Did not find a multiplicity value in the options dictionary"
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Did not find a multiplicity value in the options dictionary",
+ msg,
verbosity,
_MS_VERBOSITY_LEVEL,
)
- raise ValueError
+ raise ValueError(msg)
# multiplicity = 1
# for n in range(2, 5):
# multiplicity += 1 if options.get("M{}".format(n), None) else 0
@@ -1665,17 +1775,16 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
############################################################
# multiplicity fraction
# Calculate the probabilty, or rather, fraction, of stars that belong to this mass
- prob.append(
- Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity)[multiplicity - 1]
- )
+
+ multiplicity_probability = Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity)[multiplicity - 1]
+ prob_dict['multiplicity'] = multiplicity_probability
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob list ({})".format(
- options["M1"], prob[-1], prob
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob dict ({})".format(
+ options["M1"], prob_dict['multiplicity'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
- prob_dict['multiplicity'] = prob[-1]
############################################################
# always require an IMF for the primary star
@@ -1685,15 +1794,18 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# TODO: Create an n-part-powerlaw method that can have breakpoints and slopes. I'm using a three-part powerlaw now.
# TODO: is this actually the correct way? putting the M1 in there? Do we sample in logspace?
- prob.append(Kroupa2001(options["M1"]) * options["M1"])
+ M1_probability = Kroupa2001(options["M1"]) * options["M1"]
+ prob_dict['M1'] = M1_probability
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob list ({})".format(
- options["M1"], prob[-1], prob
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob dict ({})".format(
+ options["M1"], prob_dict['M1'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
- prob_dict['M1'] = prob[-1]
+ if M1_probability == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
"""
From here we go through the multiplicities.
@@ -1714,137 +1826,82 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Created new period interpolator: {}".format(
- Moecache["rinterpolator_q"]
- ),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
-
- # Make a table storing Moe's data for q distributions
- if not Moecache.get("rinterpolator_q", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache["q_distributions"], # Contains the table of data
- nparams=3, # log10M, log10P, q
- ndata=1, #
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
- Moecache["rinterpolator_q"]
+ Moecache["rinterpolator_log10P"]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# Make a table storing Moe's data for q distributions
- if not Moecache.get("rinterpolator_e", None):
- Moecache["rinterpolator_e"] = py_rinterpolate.Rinterpolate(
- table=Moecache["ecc_distributions"], # Contains the table of data
- nparams=3, # log10M, log10P, e
- ndata=1, #
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new e interpolator: {}".format(
- Moecache["rinterpolator_e"]
- ),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
-
- if not Moecache.get("rinterpolator_log10P2", None):
- Moecache["rinterpolator_log10P2"] = py_rinterpolate.Rinterpolate(
- table=Moecache[
- "period_distributions"
- ], # Contains the table of data
- nparams=2, # log10(M1+M2), log10P # TODO: Really? Is the first column of that table M1+M2?
- ndata=2, # binary, triple
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
-
- if not Moecache.get("rinterpolator_q", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache["q_distributions"], # Contains the table of data
- nparams=3, # log10(M1+M2), log10P2, q
- ndata=1, #
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
- Moecache["rinterpolator_q"]
- ),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
-
- if not Moecache.get("rinterpolator_log10P2", None):
- Moecache[
- "rinterpolator_log10P2"
- ] = py_rinterpolate.Rinterpolate(
- table=Moecache[
- "period_distributions"
- ], # Contains the table of data
- nparams=2, # log10(M1+M2), log10P3
- ndata=2, # binary, triple
- verbosity=verbosity
- - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
+
+ if options.get("M2", None) or options.get("M3", None) or options.get("M4", None):
+ if not Moecache.get("rinterpolator_q", None):
+ Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["q_distributions"], # Contains the table of data
+ nparams=3, # log10M, log10P, q
+ ndata=1, #
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
+ Moecache["rinterpolator_q"]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
- if not Moecache.get("rinterpolator_q", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache[
- "q_distributions"
- ], # Contains the table of data
- nparams=3, # log10(M1+M2), log10P2, q
- ndata=1, #
- verbosity=verbosity
- - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
+ # Make a table storing Moe's data for q distributions, but only if the ecc is actually sampled
+ if "ecc" in options:
+ if not options['ecc']==None:
+ if not Moecache.get("rinterpolator_e", None):
+ Moecache["rinterpolator_e"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["ecc_distributions"], # Contains the table of data
+ nparams=3, # log10M, log10P, e
+ ndata=1, #
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new e interpolator: {}".format(
+ Moecache["rinterpolator_e"]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
###############
# Calculation for period of the binary
- # Separation of the inner binary
- options["sep"] = calc_sep_from_period(
- options["M1"], options["M2"], options["P"]
- )
- # TODO: add check for min_P with instant RLOF?
-
- # Total mass inner binary:
- options["M1+M2"] = options["M1"] + options["M2"]
-
- # Set up the interpolator for the periods
- if not Moecache.get("rinterpolator_log10P", None):
- Moecache["rinterpolator_log10P"] = py_rinterpolate.Rinterpolate(
- table=Moecache["period_distributions"], # Contains the table of data
- nparams=2, # log10M, log10P
- ndata=2, # binary, triple
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
-
- if not Moecache.get("P_integrals", None):
- Moecache["P_integrals"] = {}
-
- # Calculate normalisation contant for P1
- if not Moecache["P_integrals"].get(options["M1"], None):
- calc_P_integral(
- options, 1, "P_integrals", "rinterpolator_log10P", "M1", verbosity
+ if options.get("M2", None):
+ # Separation of the inner binary
+ options["sep"] = calc_sep_from_period(
+ options["M1"], options["M2"], options["P"]
)
+ # TODO: add check for min_P with instant RLOF?
+ # TODO: Actually use the value above.
+ # Total mass inner binary:
+ options["M1+M2"] = options["M1"] + options["M2"]
+
+ # Calculate P integral or use cached value
+ calc_P_integral(
+ options, 1, "P_integrals", "rinterpolator_log10P", "M1", verbosity
+ )
# Set probabilty for P1
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M1"]), np.log10(options["P"])]
)[0]
p_val = p_val / Moecache["P_integrals"][options["M1"]]
- prob.append(p_val)
+ prob_dict['P'] = p_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended period (m={}, P={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["P"], prob[-1], prob
+ options["M1"], options["P"], prob_dict['P'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
- prob_dict['P'] = prob[-1]
+ if prob_dict['P'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
############################################################
# mass ratio (0 < q = M2/M1 < qmax)
@@ -1852,75 +1909,72 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# we need to construct the q table for the given M1
# subject to qmin = Mmin/M1
- # Build the table for q
- primary_mass = options["M1"]
- secondary_mass = options["M2"]
- m_label = "M1"
- p_label = "P"
+ if options.get("M2", None):
+ # Build the table for q
+ primary_mass = options["M1"]
+ secondary_mass = options["M2"]
+ m_label = "M1"
+ p_label = "P"
- # Make a table storing Moe's data for q distributions
- if not Moecache.get("rinterpolator_q", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache["q_distributions"], # Contains the table of data
- nparams=3, # log10M, log10P, q
- ndata=1, #
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
+ # Construct the q table
+ build_q_table(options, m_label, p_label, verbosity=verbosity)
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
- Moecache["rinterpolator_q"]
- ),
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created q_table ({}) for m={} p={}".format(Moecache[
+ "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
+ ], options[m_label], options[p_label]),
verbosity,
_MS_VERBOSITY_LEVEL,
)
- # Construct the q table
- build_q_table(options, m_label, p_label, verbosity=verbosity)
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created q_table ({}) for m={} p={}".format(Moecache[
- "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
- ], options[m_label], options[p_label]),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
-
- # Add probability for the mass ratio
- prob.append(
- Moecache[
- "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
- ].interpolate([secondary_mass / primary_mass])[0]
- )
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M={} P={} q={}) probability ({}) to the prob list ({}) ".format(
- options["M1"],
- options["P"],
- options["M2"] / options["M1"],
- prob[-1],
- prob,
- ),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
-
- ############################################################
- # Eccentricity
- # TODO add eccentricity calculations
-
- # Make a table storing Moe's data for q distributions
- if not Moecache.get("rinterpolator_e", None):
- Moecache["rinterpolator_e"] = py_rinterpolate.Rinterpolate(
- table=Moecache["ecc_distributions"], # Contains the table of data
- nparams=3, # log10M, log10P, e
- ndata=1, #
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
+ # Add probability for the mass ratio
+ q_prob = Moecache[
+ "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
+ ].interpolate([secondary_mass / primary_mass])[0]
+ prob_dict['q'] = q_prob
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new e interpolator: {}".format(
- Moecache["rinterpolator_e"]
+ "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+ options["M1"],
+ options["P"],
+ options["M2"] / options["M1"],
+ prob_dict['q'],
+ prob_dict,
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ if prob_dict['q'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
+
+ ############################################################
+ # Eccentricity
+ # TODO: ask rob if the eccentricity requires an extrapolation as well.
+
+ # Only do this if the eccentricity is sampled
+ if "ecc" in options:
+ if not options['ecc']==None:
+ # Calculate ecc integral or use cached value
+ calc_e_integral(
+ options, "ecc_integrals", "rinterpolator_e", "M1", "P", verbosity
+ )
+ mass_period_string = "{}_{}".format(options["M1"], options["P"])
+
+ # Set probabilty for ecc
+ ecc_val = Moecache["rinterpolator_e"].interpolate(
+ [np.log10(options["M1"]), np.log10(options["P"]), options["ecc"]]
+ )[0]
+ ecc_val = ecc_val / Moecache["ecc_integrals"][mass_period_string]
+ prob_dict['ecc'] = ecc_val
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended eccentricity (m={}, P={}, ecc={}) probability ({}) to the prob list ({})".format(
+ options["M1"], options["P"], options["ecc"], prob_dict['ecc'], prob_dict
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ if prob_dict['ecc'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
# Calculations for when multiplicity is bigger than 3
# BEWARE: binary_c does not evolve these systems actually and the code below should be revised for when binary_c actually evolves triples.
@@ -1943,7 +1997,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
if options["P2"] < min_P2:
# period is too short : system is not hierarchical
- prob.append(0)
+ prob_dict['P2'] = 0
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: period2 is too short: {} < {}, system is not hierarchichal. Added 0 to probability list".format(
options["P1"], min_P2
@@ -1951,7 +2005,10 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
- # TODO: consider directly returning here
+
+ if prob_dict['P2'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
else:
# period is long enough that the system is hierarchical
@@ -1961,68 +2018,38 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
options["M3"], options["M1+M2"], options["P2"]
)
- # TODO: Consider actually removing these double interpolators.
- # Not necessary to have these around, its the same dataset thats loaded
- if not Moecache.get("rinterpolator_log10P2", None):
- Moecache["rinterpolator_log10P2"] = py_rinterpolate.Rinterpolate(
- table=Moecache[
- "period_distributions"
- ], # Contains the table of data
- nparams=2, # log10(M1+M2), log10P # TODO: Really? Is the first column of that table M1+M2?
- ndata=2, # binary, triple
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
-
- # Create the dictionary where we store the area of the P2 integrals
- if not Moecache.get("P2_integrals", None):
- Moecache["P2_integrals"] = {}
-
- # Calculate the P2 integrals if they don't exist
- if not Moecache["P2_integrals"].get(options["M1+M2"], None):
- calc_P_integral(
- options,
- min_P2,
- "P2_integrals",
- "rinterpolator_log10P2",
- "M1+M2",
- verbosity,
- )
+ # Check for cached value of P integral or calculate
+ calc_P_integral(
+ options,
+ min_P2,
+ "P2_integrals",
+ "rinterpolator_log10P",
+ "M1+M2",
+ verbosity,
+ )
# Add the probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M1+M2"]), np.log10(options["P2"])]
)[0]
p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
- prob.append(p_val)
+ prob_dict['P2'] = p_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["M2"], options["P2"], prob[-1], prob
+ options["M1"], options["M2"], options["P2"], prob_dict['P2'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ if prob_dict['P2'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
############################################################
# mass ratio 2 = q2 = M3 / (M1+M2)
#
# we need to construct the q table for the given M1
# subject to qmin = Mmin/(M1+M2)
- #
- # Make a table storing Moe's data for q distributions
- if not Moecache.get("rinterpolator_q", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache["q_distributions"], # Contains the table of data
- nparams=3, # log10(M1+M2), log10P2, q
- ndata=1, #
- verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
- Moecache["rinterpolator_q"]
- ),
- verbosity,
- _MS_VERBOSITY_LEVEL,
- )
# Set the variables for the masses and their names
primary_mass = options["M1+M2"]
@@ -2039,26 +2066,27 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
# Add the probability
- prob.append(
- Moecache[
- "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
- ].interpolate([secondary_mass / primary_mass])[0]
- )
+ q2_val = Moecache[
+ "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
+ ].interpolate([secondary_mass / primary_mass])[0]
+ prob_dict['q2'] = q2_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
options["M1+M2"],
options["M3"],
options["P"],
secondary_mass / primary_mass,
- prob[-1],
- prob,
+ prob_dict['q2'],
+ prob_dict,
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ if prob_dict['q2'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
# TODO: Implement ecc2 calculation
-
if multiplicity == 4:
# quadruple system.
# TODO: Ask Rob about the strructure of the quadruple. Is this only double binary quadrupples?
@@ -2076,46 +2104,33 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
options["M1+M2"], options["mmin"], max_sep3
)
- if not Moecache.get("rinterpolator_log10P2", None):
- Moecache[
- "rinterpolator_log10P2"
- ] = py_rinterpolate.Rinterpolate(
- table=Moecache[
- "period_distributions"
- ], # Contains the table of data
- nparams=2, # log10(M1+M2), log10P3
- ndata=2, # binary, triple
- verbosity=verbosity
- - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
-
- if not Moecache.get("P2_integrals", None):
- Moecache["P2_integrals"] = {}
-
- if not Moecache["P2_integrals"].get(options["M1+M2"], None):
- calc_P_integral(
- options,
- min_P2,
- "P2_integrals",
- "rinterpolator_log10P2",
- "M1+M2",
- verbosity,
- )
+ # Calculate P integral or use the cached value
+ # TODO: Make sure we use the correct period idea here.
+ calc_P_integral(
+ options,
+ min_P2,
+ "P2_integrals",
+ "rinterpolator_log10P",
+ "M1+M2",
+ verbosity,
+ )
- # TODO: should this really be the log10P interpolator?
- #
+ # Set probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M1+M2"]), np.log10(options["P2"])]
)[0]
p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
- prob.append(p_val)
+ prob_dict['P3'] = p_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (M=4) (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["M2"], options["P2"], prob[-1], prob
+ options["M1"], options["M2"], options["P2"], prob_dict['P3'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ if prob_dict['P3'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
############################################################
# mass ratio 2
@@ -2123,17 +2138,6 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# we need to construct the q table for the given M1
# subject to qmin = Mmin/(M1+M2)
# Make a table storing Moe's data for q distributions
- # TODO: Check if this is correct.
- if not Moecache.get("rinterpolator_q", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache[
- "q_distributions"
- ], # Contains the table of data
- nparams=3, # log10(M1+M2), log10P2, q
- ndata=1, #
- verbosity=verbosity
- - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
- )
# Build the table for q2
primary_mass = options["M1+M2"]
@@ -2141,6 +2145,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
m_label = "M1+M2"
p_label = "P2"
+ # Calculate new q table
build_q_table(options, m_label, p_label, verbosity=verbosity)
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Created q_table ".format(),
@@ -2149,93 +2154,102 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
# Add the probability
- prob.append(
- Moecache[
+ q3_prob = Moecache[
"rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
- ].interpolate([secondary_mass / primary_mass])[0]
- )
+ ].interpolate([secondary_mass / primary_mass])[0]
+ prob_dict['q3'] = q3_prob
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
options["M1+M2"],
options["M3"],
options["P"],
secondary_mass / primary_mass,
- prob[-1],
- prob,
+ prob_dict['q3'],
+ prob_dict,
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ if prob_dict['q3'] == 0: # If the probability is 0 then we don't have to calculate more
+ calc_total_probdens(prob_dict)
+ return prob_dict
- # todo ecc 3
+ # TODO ecc 3
+ # check for input of multiplicity
elif multiplicity not in range(1, 5):
- verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Unknown multiplicity {}".format(
+ msg = "\tM&S: Moe_de_Stefano_2017_pdf: Unknown multiplicity {}".format(
multiplicity
- ),
+ )
+ verbose_print(
+ msg,
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ raise ValueError(msg)
- # the final probability density is the product of all the
- # probability density functions
- prob_dens = 1 * np.prod(prob)
+ # Calculate total probdens:
+ prob_dict = calc_total_probdens(prob_dict)
+ # Some info
if multiplicity == 1:
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: M1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
- options["M1"], len(prob), str(prob), prob_dens
+ options["M1"],
+ len(prob_dict),
+ str(prob_dict),
+ prob_dict['total_probdens']
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif multiplicity == 2:
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ({}): {} -> {}\n".format(
+ "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} ({}): {} -> {}\n".format(
options["M1"],
- options["M2"] / options["M1"],
+ options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
np.log10(options["P"]),
- len(prob),
- str(prob),
- prob_dens,
+ options['ecc'] if options.get("ecc", None) else "N/A",
+ len(prob_dict),
+ str(prob_dict),
+ prob_dict['total_probdens']
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif multiplicity == 3:
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ({}): M3={}} P2={} ecc2={} : {} - {}".format(
+ "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} M3={} log10P2={} ecc2={} ({}): {} -> {}".format(
options["M1"],
- options["M2"] / options["M1"],
+ options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
np.log10(options["P"]),
- len(prob),
+ options["ecc"] if options.get("ecc", None) else "N/A",
options["M3"],
- np.log10("P2"),
- np.log10("ecc2"),
- str(prob),
- prob_dens,
+ np.log10(options["P2"]),
+ options["ecc2"] if options.get("ecc2", None) else "N/A",
+ len(prob_dict),
+ str(prob_dict),
+ prob_dict['total_probdens']
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif multiplicity == 4:
verbose_print(
- "M&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ({}) : M3={} P2={}} ecc2={} : M4={} P3={} ecc3={} : {} -> {}".format(
+ "M&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} M3={} log10P2={} ecc2={} M4={} log10P3={} ecc3={} ({}) : {} -> {}".format(
options["M1"],
- options["M2"] / options["M1"],
+ options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
np.log10(options["P"]),
- len(prob),
+ options['ecc'] if options.get("ecc", None) else "N/A",
options["M3"],
- np.log10("P2"),
- np.log10(options["ecc2"]),
+ np.log10(options["P2"]),
+ options["ecc2"] if options.get("ecc2", None) else "N/A",
options["M4"],
- np.log10("P3"),
- np.log10(options["ecc3"]),
- str(prob),
- prob_dens,
+ np.log10(options["P3"]),
+ options["ecc3"] if options.get("ecc3", None) else "N/A",
+ len(prob_dict), str(prob_dict), prob_dict['total_probdens']
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
- return prob_dens
+ return prob_dict
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 26c52d937d9a1a3939fefee0475b82ee319d6078..9ae0a1651bf0e3264124acd606032455c0d2cb89 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -79,6 +79,7 @@ from binarycpython.utils.distribution_functions import (
Moecache,
LOG_LN_CONVERTER,
fill_data,
+ Moe_de_Stefano_2017_pdf,
)
from binarycpython.utils.spacing_functions import (
@@ -3459,7 +3460,7 @@ class Population:
gridtype="centred",
dphasevol="dlnm1",
precode='M_1 = np.exp(lnm1); options["M1"]=M_1',
- probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity']) if multiplicity == 1 else 1".format(
+ probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens'] if multiplicity == 1 else 1".format(
str(options)[1:-1], "'multiplicity': multiplicity", "'M1': M_1"
),
)
@@ -3598,7 +3599,7 @@ q2max=maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
M_3 = q2 * (M_1 + M_2)
sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
eccentricity2=0
- """,
+""",
spacingfunc="const({}, {}, {})".format(
options["ranges"]["q"][0]
if options.get("ranges", {}).get("q", None)
@@ -3728,7 +3729,7 @@ eccentricity3=0
# the real probability. The trick we use is to strip the options_dict as a string
# and add some keys to it:
- probdist_addition = "Moe_de_Stefano_2017_pdf({{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}, verbosity=self.grid_options['verbosity'])".format(
+ probdist_addition = "Moe_de_Stefano_2017_pdf({{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens']".format(
str(options)[1:-1],
'"multiplicity": multiplicity',
'"M1": M_1',
@@ -3776,20 +3777,111 @@ eccentricity3=0
Moecache[key].destroy()
del Moecache[key]
- def plot_MS_grid(self, options, mass):
+ def plot_MS_grid_period_distribution(self, options, mass):
"""
Function to plot the period distribution along with its probability for a given mass.
"""
import matplotlib.pyplot as plt
- plt.title("Period distribution for M={}".format(mass))
- plt.plot()
- plt.show()
+ default_options = copy.deepcopy(moe_distefano_default_options)
+ # Take the option dictionary that was given and override.
+ options = update_dicts(default_options, options)
+ options['multiplicity'] = 2
sampling_range = const(
- options["ranges"]["logP"][0],
- options["ranges"]["logP"][1],
+ options["ranges"]["logP"][0],
+ options["ranges"]["logP"][1],
options["resolutions"]["logP"][0])
+ total_prob = 0
+ probs = []
+ step = sampling_range[1]-sampling_range[0]
+ for value in sampling_range:
+ options['M1'] = mass
+ options['P'] = 10**value
+ prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])['P']
+ probs.append(prob)
+ total_prob += step * prob
+
+ self._clean_interpolators()
+
+ plt.title("Period distribution for M={}. Total prob: {}".format(mass, total_prob))
+ plt.plot(sampling_range, probs)
+ plt.ylabel("Probability")
+ plt.xlabel("log10(period)")
+ plt.show()
+
+ def plot_MS_grid_q_distribution(self, options, mass, period):
+ """
+ Function to plot the q distribution along with its probability for a given mass.
+ """
+
+ import matplotlib.pyplot as plt
+
+ default_options = copy.deepcopy(moe_distefano_default_options)
+
+ # Take the option dictionary that was given and override.
+ options = update_dicts(default_options, options)
+ options['multiplicity'] = 2
+ sampling_range = const(
+ options["ranges"]["q"][0],
+ options["ranges"]["q"][1],
+ options["resolutions"]["M"][1])
+
+ total_prob = 0
+ probs = []
+ step = sampling_range[1]-sampling_range[0]
+ print(sampling_range)
+ options['M1'] = mass
+ options['P'] = period
+ for value in sampling_range:
+ options['M2'] = value * options['M1']
+ prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity']).get('q', 0)
+ probs.append(prob)
+ total_prob += step * prob
+
+ self._clean_interpolators()
+
+ plt.title("q distribution for M={} P={}. Total prob: {}".format(mass, period, total_prob))
+ plt.plot(sampling_range, probs)
+ plt.ylabel("Probability")
+ plt.xlabel("mass ratio q")
+ plt.show()
+
+ def plot_MS_grid_ecc_distribution(self, options, mass, period):
+ """
+ Function to plot the ecc distribution along with its probability for a given mass.
+ """
+
+ import matplotlib.pyplot as plt
+ default_options = copy.deepcopy(moe_distefano_default_options)
+
+ # Take the option dictionary that was given and override.
+ options = update_dicts(default_options, options)
+ options['multiplicity'] = 2
+ sampling_range = const(
+ options["ranges"]["ecc"][0],
+ options["ranges"]["ecc"][1],
+ options["resolutions"]["ecc"][0])
+
+ total_prob = 0
+ probs = []
+ step = sampling_range[1]-sampling_range[0]
+ print(sampling_range)
+ options['M1'] = mass
+ options['P'] = period
+ for value in sampling_range:
+ options['ecc'] = value
+ prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity']).get('ecc', 0)
+ probs.append(prob)
+ total_prob += step * prob
+
+ self._clean_interpolators()
+
+ plt.title("ecc distribution for M={} P={}. Total prob: {}".format(mass, period, total_prob))
+ plt.plot(sampling_range, probs)
+ plt.ylabel("Probability")
+ plt.xlabel("Eccentricity")
+ plt.show()
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index ae1f87736c032741edb209130bf50c5f86ed2851..baf4b613003a167e5a1020d2c134e650c845700f 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -516,10 +516,11 @@ moe_distefano_default_options = {
None, # artificial qmax : set to None to use default
],
"logP": [0.0, 8.0], # 0 = log10(1 day) # 8 = log10(10^8 days)
- "ecc": [0.0, 1.0],
+ "ecc": [0.0, 0.99],
},
# minimum stellar mass
"Mmin": 0.07,
+
# multiplicity model (as a function of log10M1)
#
# You can use 'Poisson' which uses the system multiplicty
@@ -530,19 +531,20 @@ moe_distefano_default_options = {
# from the data, but then triples and quadruples are
# combined (and there are NO quadruples).
"multiplicity_model": "Poisson",
+
# multiplicity modulator:
# [single, binary, triple, quadruple]
#
# e.g. [1,0,0,0] for single stars only
# [0,1,0,0] for binary stars only
#
- # defaults to [1,1,1,1] i.e. all types
+ # defaults to [1,1,0,0] i.e. all types
#
"multiplicity_modulator": [
1, # single
1, # binary
- 1, # triple
- 1, # quadruple
+ 0, # triple
+ 0, # quadruple
],
# given a mix of multiplities, you can either (noting that
# here (S,B,T,Q) = appropriate modulator * model(S,B,T,Q) )
@@ -575,17 +577,7 @@ moe_distefano_default_options = {
# This is probably not useful except for
# testing purposes or comparing to old grids.
"normalize_multiplicities": "merge",
- # q extrapolation (below 0.15) method
- # none # q extrapolation (below 0.15) method
- # none
- # flat
- # linear2
- # plaw2
- # nolowq
- # flat
- # linear2
- # plaw2
- # nolowq
+ # q extrapolation (below 0.15 and above 0.9) method. We can choose from ['flat', 'linear', 'plaw2', 'nolowq']
"q_low_extrapolation_method": "linear",
"q_high_extrapolation_method": "linear",
}