From cadf30a51403e24f875accfd4eb5674b3edb520e Mon Sep 17 00:00:00 2001
From: David Hendriks <davidhendriks93@gmail.com>
Date: Mon, 7 Jun 2021 20:21:15 +0100
Subject: [PATCH] Working on fixing the segfault
---
binarycpython/utils/distribution_functions.py | 55 ++++++-------
binarycpython/utils/grid.py | 82 ++++++++++++++++++-
binarycpython/utils/grid_options_defaults.py | 2 +-
3 files changed, 107 insertions(+), 32 deletions(-)
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index a77c415a8..83aa6132d 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1803,9 +1803,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
- 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
+ # 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.
@@ -1898,9 +1898,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
- 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
+ # 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)
@@ -1941,9 +1941,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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
@@ -1971,9 +1971,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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.
@@ -2004,10 +2004,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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
@@ -2040,9 +2039,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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)
@@ -2081,9 +2080,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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:
@@ -2127,9 +2126,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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
@@ -2169,9 +2168,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
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
+ # 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
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 9ae0a1651..bda3c3527 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -3486,12 +3486,12 @@ class Population:
),
precode="""orbital_period = 10.0**log10per
qmin={}/M_1
-qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
- """.format(
+qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
+""".format(
options.get("Mmin", 0.07)
),
) # TODO: change the maximum_mass_ratio_for_RLOF
-
+
# binaries: mass ratio
self.add_grid_variable(
name="q",
@@ -3885,3 +3885,79 @@ eccentricity3=0
plt.ylabel("Probability")
plt.xlabel("Eccentricity")
plt.show()
+
+
+ def plot_MS_grid_ecc_q_distribution(self, options, mass, period):
+ """
+ Function to plot the ecc distribution along with its probability for a given mass.
+ """
+ import numpy as np
+ 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_ecc = const(
+ options["ranges"]["ecc"][0],
+ options["ranges"]["ecc"][1],
+ options["resolutions"]["ecc"][0])
+
+ sampling_range_q = const(
+ options["ranges"]["q"][0],
+ options["ranges"]["q"][1],
+ options["resolutions"]["M"][1])
+
+ total_prob = 0
+ probs_all = np.zeros(shape=(options["resolutions"]["ecc"][0], options["resolutions"]["M"][1]))
+ probs_ecc = np.zeros(shape=(options["resolutions"]["ecc"][0], options["resolutions"]["M"][1]))
+ probs_q = np.zeros(shape=(options["resolutions"]["ecc"][0], options["resolutions"]["M"][1]))
+
+ step_ecc = sampling_range_ecc[1]-sampling_range_ecc[0]
+ step_q = sampling_range_q[1]-sampling_range_q[0]
+
+ options['M1'] = mass
+ options['P'] = period
+
+ for i, ecc_value in enumerate(sampling_range_ecc):
+ for j, q_value in enumerate(sampling_range_q):
+ options['ecc'] = ecc_value
+ options['M2'] = q_value * options['M1']
+
+ prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])
+
+ probs_all[i][j] = prob.get('total_probdens', 0)
+ probs_ecc[i][j] = prob.get('ecc', 0)
+ probs_q[i][j] = prob.get('q', 0)
+
+ # Calculate total prob ecc
+ total_prob_ecc = 0
+ for ecc_value in sampling_range_ecc:
+ options['ecc'] = ecc_value
+ prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])
+ total_prob_ecc += prob.get('ecc', 0) * step_ecc
+
+ # Calculate total prob q:
+ total_prob_q = 0
+ for q_value in sampling_range_q:
+ options['M2'] = q_value * options['M1']
+ prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])
+ total_prob_q += prob.get('q', 0) * step_q
+
+ # Clean the interpolators
+ self._clean_interpolators()
+
+ # Set up plot
+ fig, axes = plt.subplots(nrows=1, ncols=1)
+ extent = [options["ranges"]["ecc"][0] , options["ranges"]["ecc"][1], options["ranges"]["q"][0], options["ranges"]["q"][1]]
+ im = axes.imshow(probs_all, origin='lower', extent=extent)
+ plt.title("ecc distribution for M={} P={}. Total prob ecc: {} Total prob q: {}".format(mass, period, total_prob_ecc, total_prob_q))
+ plt.ylabel("Eccentricity")
+ plt.xlabel("Mass ratio q")
+
+ fig.subplots_adjust(right=0.85)
+ cbar_ax = fig.add_axes([0.88, 0.15, 0.04, 0.7])
+ cb = fig.colorbar(im, cax=cbar_ax)
+ cb.set_label("Probability")
+ plt.show()
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index baf4b6130..3e53a23cb 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -508,7 +508,7 @@ moe_distefano_default_options = {
"ranges": {
# stellar masses (Msun)
"M": [
- 0.07, # maximum brown dwarf mass == minimum stellar mass
+ 0.08, # 0.08 is a tad bit above the minimum mass. Don't sample at 0.07, otherwise the first row of q values will have a phasevol of 0. Anything higher is fine.
80.0, # (rather arbitrary) upper mass cutoff
],
"q": [
--
GitLab