"""
Module containing the predefined distribution functions
The user can use any of these distribution functions to
generate probability distributions for sampling populations
There are distributions for the following parameters:
- mass
- period
- mass ratio
- binary fraction
Tasks:
- TODO: make some things globally present? rob does this in his module..i guess it saves
calculations but not sure if I'm gonna do that now
- TODO: add eccentricity distribution: thermal
- TODO: Add SFH distributions depending on redshift
- TODO: Add metallicity distributions depending on redshift
- TODO: Add initial rotational velocity distributions
- TODO: make an n-part power law that's general enough to fix the three part and the 4 part
"""
import gc
import math
import json
from typing import Union
import numpy as np
from binarycpython.utils.useful_funcs import calc_period_from_sep, calc_sep_from_period
from binarycpython.utils.functions import verbose_print
from binarycpython.utils.grid_options_defaults import (
_MS_VERBOSITY_LEVEL,
_MS_VERBOSITY_INTERPOLATOR_LEVEL,
)
###
# File containing probability distributions
# Mostly copied from the Perl modules
LOG_LN_CONVERTER = 1.0 / math.log(10.0)
distribution_constants = {} # To store the constants in
def prepare_dict(global_dict: dict, list_of_sub_keys: list) -> None:
"""
Function that makes sure that the global dict is prepared to have a value set there.
This dictionary will store values and factors for the distribution functions,
so that they don't have to be calculated each time.
Args:
global_dict: globally accessible dictionary where factors are stored in
list_of_sub_keys: List of keys that must become be(come) present in the global_dict
"""
internal_dict_value = global_dict
# This loop almost mimics a recursive loop into the dictionary.
# It checks whether the first key of the list is present, if not; set it with an empty dict.
# Then it overrides itself to be that (new) item, and goes on to do that again, until the list
# exhausted
for k in list_of_sub_keys:
# If the sub key doesnt exist then make an empty dict
if not internal_dict_value.get(k, None):
internal_dict_value[k] = {}
internal_dict_value = internal_dict_value[k]
def set_opts(opts: dict, newopts: dict) -> dict:
"""
Function to take a default dict and override it with newer values.
# TODO: consider changing this to just a dict.update
Args:
opts: dictionary with default values
newopts: dictionary with new values
Returns:
returns an updated dictionary
"""
if newopts:
for opt in newopts.keys():
if opt in opts.keys():
opts[opt] = newopts[opt]
return opts
def flat() -> float:
"""
Dummy distribution function that returns 1
Returns:
a flat uniform distribution: 1
"""
return 1.0
def number(value: Union[int, float]) -> Union[int, float]:
"""
Dummy distribution function that returns the input
Args:
value: the value that will be returned by this function.
Returns:
the value that was provided
"""
return value
def const(
min_bound: Union[int, float], max_bound: Union[int, float], val: float = None
) -> Union[int, float]:
"""
a constant distribution function between min=min_bound and max=max_bound.
Args:
min_bound: lower bound of the range
max_bound: upper bound of the range
Returns:
returns the value of 1/(max_bound-min_bound). If val is provided, it will check whether min_bound < val <= max_bound. if not: returns 0
"""
if val:
if not min_bound < val <= max_bound:
print("out of bounds")
prob = 0
return prob
prob = 1.0 / (max_bound - min_bound)
return prob
def powerlaw_constant(
min_val: Union[int, float], max_val: Union[int, float], k: Union[int, float]
) -> Union[int, float]:
"""
Function that returns the constant to normalise a power law
TODO: what if k is -1?
Args:
min_val: lower bound of the range
max_val: upper bound of the range
k: power law slope
Returns:
constant to normalise the given power law between the min_val and max_val range
"""
k1 = k + 1.0
# print(
# "Powerlaw consts from {} to {}, k={} where k1={}".format(
# min_val, max_val, k, k1
# )
# )
powerlaw_const = k1 / (max_val ** k1 - min_val ** k1)
return powerlaw_const
def powerlaw(
min_val: Union[int, float],
max_val: Union[int, float],
k: Union[int, float],
x: Union[int, float],
) -> Union[int, float]:
"""
Single power law with index k at x from min to max
Args:
min_val: lower bound of the power law
max_val: upper bound of the power law
k: slope of the power law
x: position at which we want to evaluate
Returns:
`probability` at the given position(x)
"""
# Handle faulty value
if k == -1:
msg = "wrong value for k"
raise ValueError(msg)
if (x < min_val) or (x > max_val):
print("input value is out of bounds!")
return 0
powerlaw_const = powerlaw_constant(min_val, max_val, k)
# power law
prob = powerlaw_const * (x ** k)
# print(
# "Power law from {} to {}: const = {}, y = {}".format(
# min_val, max_val, const, y
# )
# )
return prob
def calculate_constants_three_part_powerlaw(
m0: Union[int, float],
m1: Union[int, float],
m2: Union[int, float],
m_max: Union[int, float],
p1: Union[int, float],
p2: Union[int, float],
p3: Union[int, float],
) -> Union[int, float]:
"""
Function to calculate the constants for a three-part power law
TODO: use the power law_constant function to calculate all these values
Args:
m0: lower bound mass
m1: second boundary, between the first slope and the second slope
m2: third boundary, between the second slope and the third slope
m_max: upper bound mass
p1: first slope
p2: second slope
p3: third slope
Returns:
array of normalisation constants
"""
# print("Initialising constants for the three-part powerlaw: m0={} m1={} m2={}\
# m_max={} p1={} p2={} p3={}\n".format(m0, m1, m2, m_max, p1, p2, p3))
array_constants_three_part_powerlaw = [0, 0, 0]
array_constants_three_part_powerlaw[1] = (
((m1 ** p2) * (m1 ** (-p1)))
* (1.0 / (1.0 + p1))
* (m1 ** (1.0 + p1) - m0 ** (1.0 + p1))
)
array_constants_three_part_powerlaw[1] += (
(m2 ** (1.0 + p2) - m1 ** (1.0 + p2))
) * (1.0 / (1.0 + p2))
array_constants_three_part_powerlaw[1] += (
((m2 ** p2) * (m2 ** (-p3)))
* (1.0 / (1.0 + p3))
* (m_max ** (1.0 + p3) - m2 ** (1.0 + p3))
)
array_constants_three_part_powerlaw[1] = 1.0 / (
array_constants_three_part_powerlaw[1] + 1e-50
)
array_constants_three_part_powerlaw[0] = array_constants_three_part_powerlaw[1] * (
(m1 ** p2) * (m1 ** (-p1))
)
array_constants_three_part_powerlaw[2] = array_constants_three_part_powerlaw[1] * (
(m2 ** p2) * (m2 ** (-p3))
)
return array_constants_three_part_powerlaw
# $$array[1]=(($m1**$p2)*($m1**(-$p1)))*
# (1.0/(1.0+$p1))*
# ($m1**(1.0+$p1)-$m0**(1.0+$p1))+
# (($m2**(1.0+$p2)-$m1**(1.0+$p2)))*
# (1.0/(1.0+$p2))+
# (($m2**$p2)*($m2**(-$p3)))*
# (1.0/(1.0+$p3))*
# ($mmax**(1.0+$p3)-$m2**(1.0+$p3));
# $$array[1]=1.0/($$array[1]+1e-50);
# $$array[0]=$$array[1]*$m1**$p2*$m1**(-$p1);
# $$array[2]=$$array[1]*$m2**$p2*$m2**(-$p3);
# #print "ARRAY SET @_ => @$array\n";
# $threepart_powerlaw_consts{"@_"}=[@$array];
def three_part_powerlaw(
m: Union[int, float],
m0: Union[int, float],
m1: Union[int, float],
m2: Union[int, float],
m_max: Union[int, float],
p1: Union[int, float],
p2: Union[int, float],
p3: Union[int, float],
) -> Union[int, float]:
"""
Generalised three-part power law, usually used for mass distributions
Args:
m: mass at which we want to evaluate the distribution.
m0: lower bound mass
m1: second boundary, between the first slope and the second slope
m2: third boundary, between the second slope and the third slope
m_max: upper bound mass
p1: first slope
p2: second slope
p3: third slope
Returns:
'probability' at given mass m
"""
# TODO: add check on whether the values exist
three_part_powerlaw_constants = calculate_constants_three_part_powerlaw(
m0, m1, m2, m_max, p1, p2, p3
)
#
if m < m0:
prob = 0 # Below lower bound TODO: make this clear.
elif m0 < m <= m1:
prob = three_part_powerlaw_constants[0] * (m ** p1) # Between M0 and M1
elif m1 < m <= m2:
prob = three_part_powerlaw_constants[1] * (m ** p2) # Between M1 and M2
elif m2 < m <= m_max:
prob = three_part_powerlaw_constants[2] * (m ** p3) # Between M2 and M_MAX
else:
prob = 0 # Above M_MAX
return prob
def gaussian_normalizing_const(
mean: Union[int, float],
sigma: Union[int, float],
gmin: Union[int, float],
gmax: Union[int, float],
) -> Union[int, float]:
"""
Function to calculate the normalisation constant for the Gaussian
Args:
mean: mean of the Gaussian
sigma: standard deviation of the Gaussian
gmin: lower bound of the range to calculate the probabilities in
gmax: upper bound of the range to calculate the probabilities in
Returns:
normalisation constant for the Gaussian distribution(mean, sigma) between gmin and gmax
"""
# First time; calculate multiplier for given mean and sigma
ptot = 0
resolution = 1000
d = (gmax - gmin) / resolution
for i in range(resolution):
y = gmin + i * d
ptot += d * gaussian_func(y, mean, sigma)
# TODO: Set value in global
return ptot
def gaussian_func(
x: Union[int, float], mean: Union[int, float], sigma: Union[int, float]
) -> Union[int, float]:
"""
Function to evaluate a Gaussian at a given point, but this time without any boundaries.
Args:
x: location at which to evaluate the distribution
mean: mean of the Gaussian
sigma: standard deviation of the Gaussian
Returns:
value of the Gaussian at x
"""
gaussian_prefactor = 1.0 / math.sqrt(2.0 * math.pi)
r = 1.0 / (sigma)
y = (x - mean) * r
return gaussian_prefactor * r * math.exp(-0.5 * y ** 2)
def gaussian(
x: Union[int, float],
mean: Union[int, float],
sigma: Union[int, float],
gmin: Union[int, float],
gmax: Union[int, float],
) -> Union[int, float]:
"""
Gaussian distribution function. used for e..g Duquennoy + Mayor 1991
Args:
x: location at which to evaluate the distribution
mean: mean of the Gaussian
sigma: standard deviation of the Gaussian
gmin: lower bound of the range to calculate the probabilities in
gmax: upper bound of the range to calculate the probabilities in
Returns:
'probability' of the Gaussian distribution between the boundaries, evaluated at x
"""
# # location (X value), mean and sigma, min and max range
# my ($x,$mean,$sigma,$gmin,$gmax) = @_;
if (x < gmin) or (x > gmax):
prob = 0
else:
# normalise over given range
# TODO: add loading into global var
normalisation = gaussian_normalizing_const(mean, sigma, gmin, gmax)
prob = gaussian_func(x, mean, sigma) / normalisation
return prob
#####
# Mass distributions
#####
def Kroupa2001(m: Union[int, float], newopts: dict = None) -> Union[int, float]:
"""
Probability distribution function for Kroupa 2001 IMF, where the default values to the
three_part_powerlaw are: default = {"m0": 0.1, "m1": 0.5, "m2": 1, "mmax": 100, "p1": -1.3, "p2": -2.3,"p3": -2.3}
Args:
m: mass to evaluate the distribution at
newopts: optional dict to override the default values.
Returns:
'probability' of distribution function evaluated at m
"""
# Default parameters and override them
default = {
"m0": 0.1,
"m1": 0.5,
"m2": 1,
"mmax": 100,
"p1": -1.3,
"p2": -2.3,
"p3": -2.3,
}
value_dict = default.copy()
if newopts:
value_dict.update(newopts)
return three_part_powerlaw(
m,
value_dict["m0"],
value_dict["m1"],
value_dict["m2"],
value_dict["mmax"],
value_dict["p1"],
value_dict["p2"],
value_dict["p3"],
)
def ktg93(m: Union[int, float], newopts: dict = None) -> Union[int, float]:
"""
Probability distribution function for KTG93 IMF, where the default values to the three_part_powerlaw are: default = {"m0": 0.1, "m1": 0.5, "m2": 1, "mmax": 80, "p1": -1.3, "p2": -2.2,"p3": -2.7}
Args:
m: mass to evaluate the distribution at
newopts: optional dict to override the default values.
Returns:
'probability' of distribution function evaluated at m
"""
# TODO: ask rob what this means
# if($m eq 'uncertainties')
# {
# # return (pointer to) the uncertainties hash
# return {
# m0=>{default=>0.1,
# fixed=>1},
# m1=>{default=>0.5,
# fixed=>1},
# m2=>{default=>1.0,
# fixed=>1},
# mmax=>{default=>80.0,
# fixed=>1},
# p1=>{default=>-1.3,
# low=>-1.3,
# high=>-1.3},
# p2=>{default=>-2.2,
# low=>-2.2,
# high=>-2.2},
# p3=>{default=>-2.7,
# low=>-2.7,
# high=>-2.7}
# };
# }
# set options
# opts = set_opts({'m0':0.1, 'm1':0.5, 'm2':1.0, 'mmax':80, 'p1':-1.3, 'p2':-2.2, 'p3':-2.7},
# newopts)
defaults = {
"m0": 0.1,
"m1": 0.5,
"m2": 1.0,
"mmax": 80,
"p1": -1.3,
"p2": -2.2,
"p3": -2.7,
}
value_dict = defaults.copy()
if newopts:
value_dict.update(newopts)
return three_part_powerlaw(
m,
value_dict["m0"],
value_dict["m1"],
value_dict["m2"],
value_dict["mmax"],
value_dict["p1"],
value_dict["p2"],
value_dict["p3"],
)
# sub ktg93_lnspace
# {
# # wrapper for KTG93 on a ln(m) grid
# my $m=$_[0];
# return ktg93(@_) * $m;
# }
def imf_tinsley1980(m: Union[int, float]) -> Union[int, float]:
"""
Probability distribution function for Tinsley 1980 IMF (defined up until 80Msol): three_part_powerlaw(m, 0.1, 2.0, 10.0, 80.0, -2.0, -2.3, -3.3)
Args:
m: mass to evaluate the distribution at
Returns:
'probability' of distribution function evaluated at m
"""
return three_part_powerlaw(m, 0.1, 2.0, 10.0, 80.0, -2.0, -2.3, -3.3)
def imf_scalo1986(m: Union[int, float]) -> Union[int, float]:
"""
Probability distribution function for Scalo 1986 IMF (defined up until 80Msol): three_part_powerlaw(m, 0.1, 1.0, 2.0, 80.0, -2.35, -2.35, -2.70)
Args:
m: mass to evaluate the distribution at
Returns:
'probability' of distribution function evaluated at m
"""
return three_part_powerlaw(m, 0.1, 1.0, 2.0, 80.0, -2.35, -2.35, -2.70)
def imf_scalo1998(m: Union[int, float]) -> Union[int, float]:
"""
From Scalo 1998
Probability distribution function for Scalo 1998 IMF (defined up until 80Msol): three_part_powerlaw(m, 0.1, 1.0, 10.0, 80.0, -1.2, -2.7, -2.3)
Args:
m: mass to evaluate the distribution at
Returns:
'probability' of distribution function evaluated at m
"""
return three_part_powerlaw(m, 0.1, 1.0, 10.0, 80.0, -1.2, -2.7, -2.3)
def imf_chabrier2003(m: Union[int, float]) -> Union[int, float]:
"""
Probability distribution function for IMF of Chabrier 2003 PASP 115:763-795
Args:
m: mass to evaluate the distribution at
Returns:
'probability' of distribution function evaluated at m
"""
chabrier_logmc = math.log10(0.079)
chabrier_sigma2 = 0.69 * 0.69
chabrier_a1 = 0.158
chabrier_a2 = 4.43e-2
chabrier_x = -1.3
if m <= 0:
msg = "below bounds"
raise ValueError(msg)
if 0 < m < 1.0:
A = 0.158
dm = math.log10(m) - chabrier_logmc
prob = chabrier_a1 * math.exp(-(dm ** 2) / (2.0 * chabrier_sigma2))
else:
prob = chabrier_a2 * (m ** chabrier_x)
prob = prob / (0.1202462 * m * math.log(10))
return prob
########################################################################
# Binary fractions
########################################################################
def Arenou2010_binary_fraction(m: Union[int, float]) -> Union[int, float]:
"""
Arenou 2010 function for the binary fraction as f(M1)
GAIA-C2-SP-OPM-FA-054
www.rssd.esa.int/doc_fetch.php?id=2969346
Args:
m: mass to evaluate the distribution at
Returns:
binary fraction at m
"""
return 0.8388 * math.tanh(0.688 * m + 0.079)
# print(Arenou2010_binary_fraction(0.4))
def raghavan2010_binary_fraction(m: Union[int, float]) -> Union[int, float]:
"""
Fit to the Raghavan 2010 binary fraction as a function of
spectral type (Fig 12). Valid for local stars (Z=Zsolar).
The spectral type is converted mass by use of the ZAMS
effective temperatures from binary_c/BSE (at Z=0.02)
and the new "long_spectral_type" function of binary_c
(based on Jaschek+Jaschek's Teff-spectral type table).
Rob then fitted the result
Args:
m: mass to evaluate the distribution at
Returns:
binary fraction at m
"""
return min(
1.0,
max(
(m ** 0.1) * (5.12310e-01) + (-1.02070e-01),
(1.10450e00) * (m ** (4.93670e-01)) + (-6.95630e-01),
),
)
# print(raghavan2010_binary_fraction(2))
########################################################################
# Period distributions
########################################################################
def duquennoy1991(logper: Union[int, float]) -> Union[int, float]:
"""
Period distribution from Duquennoy + Mayor 1991. Evaluated the function gaussian(logper, 4.8, 2.3, -2, 12)
Args:
logper: logarithm of period to evaluate the distribution at
Returns:
'probability' at gaussian(logper, 4.8, 2.3, -2, 12)
"""
return gaussian(logper, 4.8, 2.3, -2, 12)
def sana12(
M1: Union[int, float],
M2: Union[int, float],
a: Union[int, float],
P: Union[int, float],
amin: Union[int, float],
amax: Union[int, float],
x0: Union[int, float],
x1: Union[int, float],
p: Union[int, float],
) -> Union[int, float]:
"""
distribution of initial orbital periods as found by Sana et al. (2012)
which is a flat distribution in ln(a) and ln(P) respectively for stars
* less massive than 15Msun (no O-stars)
* mass ratio q=M2/M1<0.1
* log(P)<0.15=x0 and log(P)>3.5=x1
and is be given by dp/dlogP ~ (logP)^p for all other binary configurations (default p=-0.55)
arguments are M1, M2, a, Period P, amin, amax, x0=log P0, x1=log P1, p
example args: 10, 5, sep(M1, M2, P), sep, ?, -2, 12, -0.55
# TODO: Fix this function!
Args:
M1: Mass of primary
M2: Mass of secondary
a: separation of binary
P: period of binary
amin: minimum separation of the distribution (lower bound of the range)
amax: maximum separation of the distribution (upper bound of the range)
x0: log of minimum period of the distribution (lower bound of the range)
x1: log of maximum period of the distribution (upper bound of the range)
p: slope of the distribution
Returns:
'probability' of orbital period P given the other parameters
"""
res = 0
if (M1 < 15) or (M2 / M1 < 0.1):
res = 1.0 / (math.log(amax) - math.log(amin))
else:
p1 = 1.0 + p
# For more details see the LyX document of binary_c for this distribution
# where the variables and normalisations are given
# we use the notation x=log(P), xmin=log(Pmin), x0=log(P0), ... to determine the
x = LOG_LN_CONVERTER * math.log(P)
xmin = LOG_LN_CONVERTER * math.log(calc_period_from_sep(M1, M2, amin))
xmax = LOG_LN_CONVERTER * math.log(calc_period_from_sep(M1, M2, amax))
# print("M1 M2 amin amax P x xmin xmax")
# print(M1, M2, amin, amax, P, x, xmin, xmax)
# my $x0 = 0.15;
# my $x1 = 3.5;
A1 = 1.0 / (
x0 ** p * (x0 - xmin) + (x1 ** p1 - x0 ** p1) / p1 + x1 ** p * (xmax - x1)
)
A0 = A1 * x0 ** p
A2 = A1 * x1 ** p
if x < x0:
res = 3.0 / 2.0 * LOG_LN_CONVERTER * A0
elif x > x1:
res = 3.0 / 2.0 * LOG_LN_CONVERTER * A2
else:
res = 3.0 / 2.0 * LOG_LN_CONVERTER * A1 * x ** p
return res
# print(sana12(10, 2, 10, 100, 1, 1000, math.log(10), math.log(1000), 6))
def interpolate_in_mass_izzard2012(
M: Union[int, float], high: Union[int, float], low: Union[int, float]
) -> Union[int, float]:
"""
Function to interpolate in mass
TODO: fix this function.
TODO: describe the args
high: at M=16.3
low: at 1.15
Args:
M: mass
high:
low:
Returns:
"""
log_interpolation = False
if log_interpolation:
return (high - low) / (math.log10(16.3) - math.log10(1.15)) * (
math.log10(M) - math.log10(1.15)
) + low
else:
return (high - low) / (16.3 - 1.15) * (M - 1.15) + low
def Izzard2012_period_distribution(
P: Union[int, float], M1: Union[int, float], log10Pmin: Union[int, float] = 1
) -> Union[int, float]:
"""
period distribution which interpolates between
Duquennoy and Mayor 1991 at low mass (G/K spectral type <~1.15Msun)
and Sana et al 2012 at high mass (O spectral type >~16.3Msun)
This gives dN/dlogP, i.e. DM/Raghavan's Gaussian in log10P at low mass
and Sana's power law (as a function of logP) at high mass
TODO: fix this function
Args:
P: period
M1: Primary star mass
log10Pmin: minimum period in base log10 (optional)
Returns:
'probability' of interpolated distribution function at P and M1
"""
# Check if there is input and force it to be at least 1
log10Pmin //= -1.0
log10Pmin = max(-1.0, log10Pmin)
# save mass input and limit mass used (M1 from now on) to fitted range
Mwas = M1
M1 = max(1.15, min(16.3, M1))
# print("Izzard2012 called for M={} (truncated to {}), P={}\n".format(Mwas, M1, P))
# Calculate the normalisations
# need to normalise the distribution for this mass
# (and perhaps secondary mass)
prepare_dict(distribution_constants, ["Izzard2012", M1])
if not distribution_constants["Izzard2012"][M1].get(log10Pmin):
distribution_constants["Izzard2012"][M1][
log10Pmin
] = 1 # To prevent this loop from going recursive
N = 200.0 # Resolution for normalisation. I hope 1000 is enough
dlP = (10.0 - log10Pmin) / N
C = 0 # normalisation constant.
for lP in np.arange(log10Pmin, 10, dlP):
C += dlP * Izzard2012_period_distribution(10 ** lP, M1, log10Pmin)
distribution_constants["Izzard2012"][M1][log10Pmin] = 1.0 / C
# print(
# "Normalisation constant for Izzard2012 M={} (log10Pmin={}) is\
# {}\n".format(
# M1, log10Pmin, distribution_constants["Izzard2012"][M1][log10Pmin]
# )
# )
lP = math.log10(P)
# log period
# # fits
mu = interpolate_in_mass_izzard2012(M1, -17.8, 5.03)
sigma = interpolate_in_mass_izzard2012(M1, 9.18, 2.28)
K = interpolate_in_mass_izzard2012(M1, 6.93e-2, 0.0)
nu = interpolate_in_mass_izzard2012(M1, 0.3, -1)
g = 1.0 / (1.0 + 1e-30 ** (lP - nu))
lPmu = lP - mu
# print(
# "M={} ({}) P={} : mu={} sigma={} K={} nu={} norm=%g\n".format(
# Mwas, M1, P, mu, sigma, K, nu
# )
# )
# print "FUNC $distdata{Izzard2012}{$M}{$log10Pmin} * (exp(- (x-$mu)**2/(2.0*$sigma*$sigma) ) + $K/MAX(0.1,$lP)) * $g;\n";
if (lP < log10Pmin) or (lP > 10.0):
return 0
else:
return (
distribution_constants["Izzard2012"][M1][log10Pmin]
* (math.exp(-lPmu * lPmu / (2.0 * sigma * sigma)) + K / max(0.1, lP))
* g
)
########################################################################
# Mass ratio distributions
########################################################################
def flatsections(x: float, opts: dict) -> Union[float, int]:
"""
Function to generate flat distributions, possibly in multiple sections
Args:
x: mass ratio value
opts: list containing the flat sections. Which are themselves dictionaries, with keys "max": upper bound, "min": lower bound and "height": value
Returns:
probability of that mass ratio.
"""
c = 0
y = 0
for opt in opts:
dc = (opt["max"] - opt["min"]) * opt["height"]
# print("added flatsection ({}-{})*{} = {}\n".format(
# opt['max'], opt['min'], opt['height'], dc))
c += dc
if opt["min"] <= x <= opt["max"]:
y = opt["height"]
# print("Use this\n")
c = 1.0 / c
y = y * c
# print("flatsections gives C={}: y={}\n",c,y)
return y
# print(flatsections(1, [{'min': 0, 'max': 2, 'height': 3}]))
########################################################################
# Eccentricity distributions
########################################################################
########################################################################
# Star formation histories
########################################################################
def cosmic_SFH_madau_dickinson2014(z):
"""
Cosmic star formation history distribution from Madau & Dickonson 2014 (https://arxiv.org/pdf/1403.0007.pdf)
Args:
z: redshift
Returns:
Cosmic star formation rate in Solar mass year^-1 mega parsec^-3
"""
CSFH = 0.015 * ((1 + z) ** 2.7) / (1 + (((1 + z) / 2.9) ** 5.6))
return CSFH
########################################################################
# Metallicity distributions
########################################################################
########################################################################
# Moe & DiStefano 2017 functions
#
# The code below are functions that are used to set up and interpolate
# on the Moe & DiStefano 2017 data. The interpolators take the last
# known value if we try to interpolate outside of the tables.
# There are still some open tasks and improvements that can be made:
#
# TODO: Solve the memory issues that are present.
# Are the interpolators not cleaned?
# TODO: Parallelize the setting up of the interpolators
# TODO: Generalise the code such that we can input other/newer tables
########################################################################
import py_rinterpolate
# Global dictionary to store values in
Moecache = {}
def poisson(lambda_val, n, nmax=None, verbosity=0):
"""
Function that calculates the Poisson value and normalises
TODO: improve the description
"""
cachekey = "{} {} {}".format(lambda_val, n, nmax)
if distribution_constants.get("poisson_cache", None):
if distribution_constants["poisson_cache"].get(cachekey, None):
p_val = distribution_constants["poisson_cache"][cachekey]
verbose_print(
"\tM&S: found cached value for poisson({}, {}, {}): {}".format(
lambda_val, n, nmax, p_val
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
return p_val
# Poisson distribution : note, n can be zero
#
# nmax is the truncation : if set, we normalise
# correctly.
p_val = _poisson(lambda_val, n)
if nmax:
I_poisson = 0
for i in range(nmax + 1):
I_poisson += _poisson(lambda_val, i)
p_val = p_val / I_poisson
# Add to cache
if not distribution_constants.get("poisson_cache", None):
distribution_constants["poisson_cache"] = {}
distribution_constants["poisson_cache"][cachekey] = p_val
verbose_print(
"\tM&S: Poisson({}, {}, {}): {}".format(lambda_val, n, nmax, p_val),
verbosity,
_MS_VERBOSITY_LEVEL,
)
return p_val
def _poisson(lambda_val, n):
"""
Function to return the Poisson value
"""
return (lambda_val ** n) * np.exp(-lambda_val) / (1.0 * math.factorial(n))
def get_max_multiplicity(multiplicity_array):
"""
Function to get the maximum multiplicity
"""
max_multiplicity = 0
for n in range(4):
if multiplicity_array[n] > 0:
max_multiplicity = n + 1
return max_multiplicity
def merge_multiplicities(result_array, max_multiplicity, verbosity=0):
"""
Function to fold the multiplicities higher than the max_multiplicity onto the max_multiplicity
if max_multiplicity == 1:
All the multiplicities are folded onto multiplicity == 1. This will always total to 1
if max_multiplicity == 2:
The multiplicity fractions of the triple and quadruples are folded onto that of the binary multiplicity fraction
if max_multiplicity == 3:
The multiplicity fractions of the quadruples are folded onto that of the triples
"""
if not max_multiplicity in range(1, 5):
msg = "\tM&S: merge_multiplicities: max_multiplicity has to be between 1 and 4. It is {} now".format(
max_multiplicity
)
verbose_print(
msg,
verbosity,
0,
)
raise ValueError(msg)
# Fold multiplicities:
verbose_print(
"\tM&S: merge_multiplicities: Merging multiplicities with initial array {} and max multiplicity {}".format(result_array, max_multiplicity),
verbosity,
_MS_VERBOSITY_LEVEL,
)
for i in range(max_multiplicity, len(result_array))[::-1]:
result_array[i-1] += result_array[i]
result_array[i] = 0
verbose_print(
"\tM&S: merge_multiplicities: Merging multiplicities to new array {}".format(result_array),
verbosity,
_MS_VERBOSITY_LEVEL,
)
return result_array
def normalize_dict(result_dict, verbosity=0):
"""
Function to normalise a dictionary
"""
sum_result = sum([result_dict[key] for key in result_dict.keys()])
for key in result_dict.keys():
result_dict[key] = result_dict[key] / sum_result
return result_dict
def Moe_di_Stefano_2017_multiplicity_fractions(options, verbosity=0):
"""
Function that creates a list of probability fractions and
normalises and merges them according to the users choice.
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 of mass
"""
# Use the global Moecache
global Moecache
multiplicity_modulator_array = np.array(options["multiplicity_modulator"]) # Modulator array
# Check for length
if not len(multiplicity_modulator_array)==4:
msg = "Multiplicity modulator has to have 4 elements. Now it is {}, len: {}".format(multiplicity_modulator_array, len(multiplicity_modulator_array))
verbose_print(
msg,
verbosity,
0,
)
raise ValueError(msg)
# Set up some arrays
full_fractions_array = np.zeros(4) # Meant to contain the real fractions
weighted_fractions_array = np.zeros(4) # Meant to contain the fractions multiplied by the multiplicity modulator
# Get max multiplicity
max_multiplicity = get_max_multiplicity(multiplicity_modulator_array)
# ... it's better to interpolate the multiplicity and then
# use a Poisson distribution to calculate the fractions
# (this is more accurate)
# Set up the multiplicity interpolator
if not Moecache.get("rinterpolator_multiplicity", None):
Moecache["rinterpolator_multiplicity"] = py_rinterpolate.Rinterpolate(
table=Moecache["multiplicity_table"], # Contains the table of data
nparams=1, # logM1
ndata=4, # The amount of datapoints (the parameters that we want to interpolate)
verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
)
if options["multiplicity_model"] == "Poisson":
multiplicity = Moecache["rinterpolator_multiplicity"].interpolate(
[np.log10(options["M_1"])]
)[0]
# Fill the multiplicity array
for n in range(4):
full_fractions_array[n] = poisson(multiplicity, n, 3, verbosity)
# Normalize it so it fills to one when taking all the multiplicities:
full_fractions_array = full_fractions_array/np.sum(full_fractions_array)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_multiplicity_fractions: using model {}: full_fractions_array: {}".format("Poisson", full_fractions_array),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif options["multiplicity_model"] == "data":
# use the fractions calculated from Moe's data directly
#
# note that in this case, there are no quadruples: these
# are combined with triples
# Fill with the raw values
for n in range(3):
full_fractions_array[n] = (
Moecache["rinterpolator_multiplicity"].interpolate(
[np.log10(options["M_1"])]
)[n + 1])
# Set last value
full_fractions_array[3] = 0.0 # no quadruples
verbose_print(
"\tM&S: Moe_di_Stefano_2017_multiplicity_fractions: using model {}: full_fractions_array: {}".format("data", full_fractions_array),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# Normalisation:
if options["normalize_multiplicities"] == "raw":
# Don't multiply by the multiplicity_array, but do give a fractions array
verbose_print(
"\tM&S: Moe_di_Stefano_2017_multiplicity_fractions: Not normalising (using raw results): results: {}".format(full_fractions_array),
verbosity,
_MS_VERBOSITY_LEVEL,
)
result = full_fractions_array
elif options["normalize_multiplicities"] == "norm":
# Multiply the full_multiplicity_fraction array by the multiplicity_multiplier_array, creating a weighted fractions array
weighted_fractions_array = full_fractions_array * multiplicity_modulator_array
# Normalise this so it is in total 1:
result = weighted_fractions_array/np.sum(weighted_fractions_array)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_multiplicity_fractions: Normalising with {}. result: {}".format("norm", result),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif options["normalize_multiplicities"] == "merge":
# We first take the full multiplicity array
# (i.e. not multiplied by multiplier) and do the merging
result = merge_multiplicities(full_fractions_array, max_multiplicity, verbosity=verbosity)
# Then normalise to be sure
result = result/np.sum(result)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_multiplicity_fractions: Normalising with {}, max_multiplicity={} result={}".format("merge", max_multiplicity, result),
verbosity,
_MS_VERBOSITY_LEVEL,
)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_multiplicity_fractions: {}".format(str(result)),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# return array reference
return result
def build_q_table(options, m, p, verbosity=0):
############################################################
#
# Build an interpolation table for q, given a mass and
# orbital period.
#
# $m and $p are labels which determine which system(s)
# to look up from Moe's data:
#
# $m can be M1, M2, M3, M4, or if set M1+M2 etc.
# $p can be P, P2, P3
#
# The actual values are in $opts:
#
# mass is in $opts->{$m}
# period is $opts->{$p}
#
# Since the information from the table for M&S is independent of any choice we make,
# we need to take into account that for example our choice of minimum mass leads to
# a minimum q_min that is not the same as in the table
# We should ignore those parts of the table and renormalise.
# If we are below the lowest value of qmin in the table we need to extrapolate the data
#
# Anyway, the goal of this function is to provide some extrapolated values for q when we should sample outside of the boundaries
############################################################
# We can check if we have a cached value for this already:
# TODO: fix this cache check.
incache = False
if Moecache.get("rinterpolator_q_metadata", None):
if (Moecache["rinterpolator_q_metadata"].get(m, None)) and (
Moecache["rinterpolator_q_metadata"].get(p, None)
):
if (Moecache["rinterpolator_q_metadata"][m] == options[m]) and (
Moecache["rinterpolator_q_metadata"][p] == options[p]
):
incache = True
verbose_print(
"\tM&S: build_q_table: Found cached values for m={} p={}".format(
options[m], options[p]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
else:
verbose_print(
"\tM&S: build_q_table: Cached values for different m={} p={}. Freeing current table and making new table".format(
options[m], options[p]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
#
if not incache:
# trim and/or expand the table to the range $qmin to $qmax.
# qmin is set by the minimum stellar mass : below this
# the companions are planets
# qmin = options["ranges"]["M"][
# 0
# ] # TODO: this lower range must not be lower than Mmin.
qmin = options['Mmin']/options['M_1']
verbose_print(
"\tM&S: build_q_table qmin: {}".format(
qmin,
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# qmax = maximum_mass_ratio_for_RLOF(options[m], options[p])
# TODO: change this to the above
qmax = 1
# qdata contains the table that we modify: we get
# the original data by interpolating Moe's table
qdata = {}
can_renormalize = 1
qeps = 1e-8 # small number but such that qeps+1 != 1
if qeps + 1 == 1.0:
verbose_print(
"\tM&S: build_q_table: qeps (= {}) +1 == 1. Make qeps larger".format(
qeps,
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
if qmin >= qmax:
# there may be NO binaries in this part of the parameter space:
# in which case, set up a table with lots of zero in it
qdata = {0: 0, 1: 0}
can_renormalize = 0
else:
# qmin and qmax mean we'll get something non-zero
can_renormalize = 1
# require extrapolation sets whether we need to extrapolate
# at the low and high ends
require_extrapolation = {}
if qmin >= 0.15:
# qmin is inside Moe's table : this is easy,
# we just keep points from qmin at the low
# end to qmax at the high end.
require_extrapolation["low"] = 0
require_extrapolation[
"high"
] = 1 # TODO: shouldn't the extrapolation need to happen if qmax > 0.95
qdata[qmin] = Moecache["rinterpolator_q"].interpolate(
[np.log10(options[m]), np.log10(options[p]), qmin]
)[0]
for q in np.arange(0.15, 0.950001, 0.1):
if (q >= qmin) and (q <= qmax):
qdata[q] = Moecache["rinterpolator_q"].interpolate(
[np.log10(options[m]), np.log10(options[p]), q]
)[0]
else:
require_extrapolation["low"] = 1
require_extrapolation["high"] = 1
if qmax < 0.15:
# qmax < 0.15 which is off the edge
# of the table. In this case, choose
# two points at q=0.15 and 0.16 and interpolate
# at these in case we want to extrapolate.
for q in [0.15, 0.16]:
qdata[q] = Moecache["rinterpolator_q"].interpolate(
[np.log10(options[m]), np.log10(options[p]), q]
)[0]
else:
# qmin < 0.15 and qmax > 0.15, so we
# have to generate Moe's table for
# q = 0.15 (i.e. 0.1 to 0.2) to 0.95 (0.9 to 1)
# as a function of M1 and orbital period,
# to obtain the q distribution data.
for q in np.arange(0.15, np.min([0.950001, qmax + 0.0001]), 0.1):
val = Moecache["rinterpolator_q"].interpolate(
[np.log10(options[m]), np.log10(options[p]), q]
)[0]
qdata[q] = val
# just below qmin, if qmin>qeps, we want nothing
if qmin - 0.15 > qeps:
q = qmin - qeps
qdata[q] = 0
require_extrapolation["low"] = 0
# just above qmax, if qmax<1, we want nothing
if qmax < 0.95:
q = qmax + qeps
qdata[q] = 0
require_extrapolation["high"] = 0
# sorted list of qs
qs = sorted(qdata.keys())
if len(qs) == 0:
msg = "No qs found error"
raise ValueError(msg)
elif len(qs) == 1:
# only one q value : pretend there are two
# with a flat distribution up to 1.0.
if qs[0] == 1.0:
qs[0] = 1.0 - 1e-6
qs[1] = 1
qdata[qs[0]] = 1
qdata[qs[1]] = 1
else:
qs[1] = 1.0
qdata[qs[1]] = qs[0]
# We actually should do the extrapolation now.
else:
# Loop over both the lower end and the upper end
for pre in ["low", "high"]:
if require_extrapolation[pre] == 0:
continue
else:
sign = -1 if pre == "low" else 1
end_index = 0 if pre == "low" else len(qs) - 1
indices = (
[0, 1] if pre == "low" else [len(qs) - 1, len(qs) - 2]
) # Based on whether we do the high or low end we need to use two different indices
method = options.get(
"q_{}_extrapolation_method".format(pre), None
)
qlimit = qmin if pre == "low" else qmax
verbose_print(
"\tM&S: build_q_table: Extrapolating: Q: {} method: {}, indices: {} End index: {}".format(
pre, method, indices, end_index
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# truncate the distribution
qdata[max(0.0, min(1.0, qlimit + sign * qeps))] = 0
if method == None:
# no extrapolation : just interpolate between 0.10 and 0.95
verbose_print(
"\tM&S: build_q_table: using no extrapolations".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
continue
elif method == "flat":
# use the end value and extrapolate it
# with zero slope
qdata[qlimit] = qdata[qs[end_index]]
verbose_print(
"\tM&S: build_q_table: using constant extrapolation".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif method == "linear":
qdata[qlimit] = linear_extrapolation_q(
qs=qs,
indices=indices,
qlimit=qlimit,
qdata=qdata,
end_index=end_index,
verbosity=verbosity,
)
verbose_print(
"\tM&S: build_q_table: using linear extrapolation".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
if pre=='low':
below_qlimit = qlimit-qeps
if below_qlimit > 0:
qdata[below_qlimit] = 0
qdata[0] = 0
verbose_print(
"\tM&S: build_q_table: using linear extrapolation and setting the points below the lower q bound ({}) to 0 ".format(qlimit),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif method == "plaw2":
qdata[qlimit] = powerlaw_extrapolation_q(
qs=qs, indices=indices, qdata=qdata, verbosity=verbosity
)
verbose_print(
"\tM&S: build_q_table: using powerlaw extrapolation".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif method == "nolowq":
newq = 0.05
qdata[newq] = 0
verbose_print(
"\tM&S: build_q_table: setting lowq to 0".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif method=="poly":
# TODO: consider implementing the poly method (see Perl version)
raise ValueError("M&S: build_q_table: Method 'poly' not implemented")
else:
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(msg)
# regenerate qs in new table. This is now the updated list of qs where we have some extrapolated numbers
tmp_table = []
for q in sorted(qdata.keys()):
tmp_table.append([q, qdata[q]])
# Make an interpolation table to contain our modified data
q_interpolator = py_rinterpolate.Rinterpolate(
table=tmp_table,
nparams=1,
ndata=1, # Contains the table of data # q #
verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
)
verbose_print(
"\tM&S: build_q_table: Created a new Q table",
verbosity,
_MS_VERBOSITY_LEVEL,
)
if can_renormalize:
verbose_print(
"\tM&S: build_q_table: Renormalising table",
verbosity,
_MS_VERBOSITY_LEVEL,
)
# now we integrate and renormalise (if the table is not all zero)
I = get_integration_constant_q(q_interpolator, tmp_table, qdata, verbosity=verbosity)
if I > 0:
# normalise to 1.0 by dividing the data by 1.0/$I
q_interpolator.multiply_table_column(1, 1.0 / I)
# test this
new_I = get_integration_constant_q(q_interpolator, tmp_table, qdata, verbosity=verbosity)
# fail if error in integral > 1e-6 (should be ~ machine precision)
if abs(1.0 - new_I) > 1e-6:
verbose_print(
"\tM&S: build_q_table: Error: > 1e-6 in q probability integral: {}".format(
I
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# set this new table in the cache
Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)] = q_interpolator
verbose_print(
"\tM&S: build_q_table: stored q_interpolater as {}".format(
"rinterpolator_q_given_{}_log10{}".format(m, p)
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# Store the values for which this table was set up in the dict
if not Moecache.get("rinterpolator_q_metadata", None):
Moecache["rinterpolator_q_metadata"] = {}
Moecache["rinterpolator_q_metadata"][m] = options[m]
Moecache["rinterpolator_q_metadata"][p] = options[p]
def powerlaw_extrapolation_q(qdata, qs, indices, verbosity=0):
"""
Function to do the power law extrapolation at the lower end of the q range
"""
newq = 0.05
# use a power-law extrapolation down to q=0.05, if possible
if (qdata[qs[indices[0]]] == 0.0) and (qdata[qs[indices[1]]] == 0.0):
# not possible
return 0
else:
slope = (np.log10(qdata[qs[indices[1]]]) - np.log10(qdata[qs[indices[0]]])) / (
np.log10(qs[indices[1]]) - np.log10(qs[indices[0]])
)
intercept = np.log10(qdata[qs[indices[0]]]) - slope * np.log10(qs[indices[0]])
return slope * newq + intercept
def linear_extrapolation_q(qs, indices, qlimit, qdata, end_index, verbosity=0):
"""
Function to do the linear extrapolation for q.
"""
# linear extrapolation
dq = qs[indices[1]] - qs[indices[0]]
if dq == 0:
verbose_print(
"\tM&S: build_q_table: linear dq=0".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# No change
return qs[end_index]
else:
slope = (qdata[qs[indices[1]]] - qdata[qs[indices[0]]]) / dq
intercept = qdata[qs[indices[0]]] - slope * qs[indices[0]]
qdata[qlimit] = max(0.0, slope * qlimit + intercept)
verbose_print(
"\tM&S: build_q_table: linear Slope: {} intercept: {} dn/dq({}) = {}".format(
slope, intercept, qlimit, qdata[qlimit]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
return max(0.0, slope * qlimit + intercept)
def get_integration_constant_q(q_interpolator, tmp_table, qdata, verbosity=0):
"""
Function to integrate the q interpolator and return the integration constant
"""
dq = 1e-3 # resolution of the integration/renormalisation
I = 0
# integrate: note that the value of the integral is
# meaningless to within a factor (which depends on $dq)
for q in np.arange(0, 1 + 2e-6, dq):
x = q_interpolator.interpolate([q])
if len(x) == 0:
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(msg)
else:
I += x[0] * dq
# verbose_print(
# "\tM&S: build_q_table: dn/dq ({}) = {} I -> = {}".format(q, x[0], I),
# verbosity,
# _MS_VERBOSITY_LEVEL,
# )
return I
def fill_data(sample_values, data_dict):
"""
Function that returns the normalised array of values for given logmass and logperiod
used for the e and q values
TODO: make sure we do the correct thing with the dstep
"""
data = {}
I = 0
dstep = float(sample_values[1]) - float(sample_values[0])
# Read out the data
for sample_value in sample_values:
val = data_dict[sample_value]
data[sample_value] = val
I += val
# Normalise the data
for sample_value in sample_values:
data[sample_value] = data[sample_value] / I
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 renormalise 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
):
"""
Function to calculate the P integral
We need to renormalise this because min_per > 0, and not all periods should be included
"""
global Moecache
max_logP = 10
# 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(options[mass_string], None):
I = 0
dlogP = 1e-3
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_di_Stefano_2017_pdf(options, verbosity=0):
"""
Moe & diStefano function to calculate the probability density.
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]
ecc, ecc2, ecc3 => eccentricities [numbering as for P above]
mmin => minimum allowed stellar mass (default 0.07)
mmax => maximum allowed stellar mass (default 80.0)
"""
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf with options:\n\t\t{}".format(
json.dumps(options)
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
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_di_Stefano_2017_pdf: Did not find a multiplicity value in the options dictionary"
verbose_print(
msg,
verbosity,
_MS_VERBOSITY_LEVEL,
)
raise ValueError(msg)
# multiplicity = 1
# for n in range(2, 5):
# multiplicity += 1 if options.get("M{}".format(n), None) else 0
else:
multiplicity = options["multiplicity"]
# Immediately return 0 if the multiplicity modulator is 0
if options["multiplicity_modulator"][multiplicity - 1] == 0:
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: returning 0 because of the multiplicity modulator being 0",
verbosity,
_MS_VERBOSITY_LEVEL,
)
return 0
############################################################
# multiplicity fraction
# Calculate the probability, or rather, fraction, of stars that belong to this mass
multiplicity_probability = Moe_di_Stefano_2017_multiplicity_fractions(options, verbosity)[multiplicity - 1]
prob_dict['multiplicity'] = multiplicity_probability
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob dict ({})".format(
options["M_1"], prob_dict['multiplicity'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
############################################################
# always require an IMF for the primary star
#
# NB multiply by M1 to convert dN/dM to dN/dlnM
# (dlnM = dM/M, so 1/dlnM = M/dM)
# TODO: Create an n-part-powerlaw method that can have breakpoints and slopes. I'm using a three-part power law now.
# TODO: is this actually the correct way? putting the M1 in there? Do we sample in log space?
M1_probability = Kroupa2001(options["M_1"]) * options["M_1"]
prob_dict['M_1'] = M1_probability
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob dict ({})".format(
options["M_1"], prob_dict['M_1'], prob_dict
),
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
"""
From here we go through the multiplicities.
"""
if multiplicity >= 2:
# If the multiplicity is higher than 1, we will need to construct the following tables:
# - period distribution table
# - q distribution table
# - eccentricity distribution table
# 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),
)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Created new period interpolator: {}".format(
Moecache["rinterpolator_log10P"]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# Make a table storing Moe's data for q distributions
if options.get("M_2", None) or options.get("M_3", None) or options.get("M_4", 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_di_Stefano_2017_pdf: Created new q interpolator: {}".format(
Moecache["rinterpolator_q"]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# 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_di_Stefano_2017_pdf: Created new e interpolator: {}".format(
Moecache["rinterpolator_e"]
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
###############
# Calculation for period of the binary
if options.get("M_2", None):
# Separation of the inner binary
options["sep"] = calc_sep_from_period(
options["M_1"], options["M_2"], options["P"]
)
# TODO: add check for min_P with instant RLOF?
# TODO: Actually use the value above.
# Total mass inner binary:
options["M_1+M_2"] = options["M_1"] + options["M_2"]
# Calculate P integral or use cached value
calc_P_integral(
options, 1, "P_integrals", "rinterpolator_log10P", "M_1", verbosity
)
# Set probabilty for P1
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M_1"]), np.log10(options["P"])]
)[0]
p_val = p_val / Moecache["P_integrals"][options["M_1"]]
prob_dict['P'] = p_val
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Appended period (m={}, P={}) probability ({}) to the prob list ({})".format(
options["M_1"], options["P"], prob_dict['P'], prob_dict
),
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
############################################################
# mass ratio (0 < q = M2/M1 < qmax)
#
# we need to construct the q table for the given M1
# subject to qmin = Mmin/M1
if options.get("M_2", None):
# Build the table for q
primary_mass = options["M_1"]
secondary_mass = options["M_2"]
m_label = "M_1"
p_label = "P"
# Construct the q table
build_q_table(options, m_label, p_label, verbosity=verbosity)
verbose_print(
"\tM&S: Moe_di_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
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_di_Stefano_2017_pdf: appended mass ratio (M={} P={} q={}) probability ({}) to the prob list ({}) ".format(
options["M_1"],
options["P"],
options["M_2"] / options["M_1"],
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", "M_1", "P", verbosity
)
mass_period_string = "{}_{}".format(options["M_1"], options["P"])
# Set probability for ecc
ecc_val = Moecache["rinterpolator_e"].interpolate(
[np.log10(options["M_1"]), 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_di_Stefano_2017_pdf: Appended eccentricity (m={}, P={}, ecc={}) probability ({}) to the prob list ({})".format(
options["M_1"], 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.
# For that reason, I would not advise to use things with multiplicity > 3
if multiplicity >= 3:
############################################################
# orbital period 2 =
# orbital period of star 3 (multiplicity==3) or
# the star3+star4 binary (multiplicity==4)
#
# we assume the same period distribution for star 3
# (or stars 3 and 4) but with a separation that is >10*a*(1+e)
# where 10*a*(1+e) is the maximum apastron separation of
# stars 1 and 2
# TODO: Is this a correct assumption?
max_sep = 10.0 * options["sep"] * (1.0 + options["ecc"])
min_P2 = calc_period_from_sep(options["M_1+M_2"], options["mmin"], max_sep)
if options["P2"] < min_P2:
# period is too short : system is not hierarchical
prob_dict['P2'] = 0
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: period2 is too short: {} < {}, system is not hierarchical. Added 0 to probability list".format(
options["P1"], min_P2
),
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
else:
# period is long enough that the system is hierarchical
# hence the separation between the outer star
# and inner binary
options["sep2"] = calc_sep_from_period(
options["M_3"], options["M_1+M_2"], options["P2"]
)
# Check for cached value of P integral or calculate
calc_P_integral(
options,
min_P2,
"P2_integrals",
"rinterpolator_log10P",
"M_1+M_2",
verbosity,
)
# Add the probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M_1+M_2"]), np.log10(options["P2"])]
)[0]
p_val = p_val / Moecache["P2_integrals"][options["M_1+M_2"]]
prob_dict['P2'] = p_val
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Appended period2 (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
options["M_1"], options["M_2"], 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)
# Set the variables for the masses and their names
primary_mass = options["M_1+M_2"]
secondary_mass = options["M_3"]
m_label = "M_1+M_2"
p_label = "P2"
# Build q table
build_q_table(options, m_label, p_label, verbosity=verbosity)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Called build_q_table",
verbosity,
_MS_VERBOSITY_LEVEL,
)
# Add the probability
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_di_Stefano_2017_pdf: appended mass ratio (M_1+M_2={} M_3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
options["M_1+M_2"],
options["M_3"],
options["P"],
secondary_mass / primary_mass,
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 structure of the quadruple. Is this only double binary quadruples?
############################################################
# orbital period 3
#
# we assume the same period distribution for star 4
# as for any other stars but Pmax must be such that
# sep3 < sep2 * 0.2
# TODO: fix this here
max_sep3 = 0.2 * options["sep2"] * (1.0 + options["ecc2"])
max_per3 = calc_period_from_sep(
options["M_1+M_2"], options["mmin"], max_sep3
)
# 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",
"M_1+M_2",
verbosity,
)
# Set probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M_1+M_2"]), np.log10(options["P2"])]
)[0]
p_val = p_val / Moecache["P2_integrals"][options["M_1+M_2"]]
prob_dict['P3'] = p_val
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Appended period2 (M=4) (M_1={} M_2={}, P2={}) probability ({}) to the prob list ({})".format(
options["M_1"], options["M_2"], 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
#
# 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
# Build the table for q2
primary_mass = options["M_1+M_2"]
secondary_mass = options["M_3"]
m_label = "M_1+M_2"
p_label = "P2"
# Calculate new q table
build_q_table(options, m_label, p_label, verbosity=verbosity)
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: Created q_table ".format(),
verbosity,
_MS_VERBOSITY_LEVEL,
)
# Add the probability
q3_prob = Moecache[
"rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
].interpolate([secondary_mass / primary_mass])[0]
prob_dict['q3'] = q3_prob
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: appended mass ratio (M_1+M_2={} M_3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
options["M_1+M_2"],
options["M_3"],
options["P"],
secondary_mass / primary_mass,
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
# check for input of multiplicity
elif multiplicity not in range(1, 5):
msg = "\tM&S: Moe_di_Stefano_2017_pdf: Unknown multiplicity {}".format(
multiplicity
)
verbose_print(
msg,
verbosity,
_MS_VERBOSITY_LEVEL,
)
raise ValueError(msg)
# Calculate total probdens:
prob_dict = calc_total_probdens(prob_dict)
# Some info
if multiplicity == 1:
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: M_1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
options["M_1"],
len(prob_dict),
str(prob_dict),
prob_dict['total_probdens']
),
verbosity,
_MS_VERBOSITY_LEVEL,
)
elif multiplicity == 2:
verbose_print(
"\tM&S: Moe_di_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} ({}): {} -> {}\n".format(
options["M_1"],
options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
np.log10(options["P"]),
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_di_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} M_3={} log10P2={} ecc2={} ({}): {} -> {}".format(
options["M_1"],
options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
np.log10(options["P"]),
options["ecc"] if options.get("ecc", None) else "N/A",
options["M_3"],
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_di_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} M_3={} log10P2={} ecc2={} M_4={} log10P3={} ecc3={} ({}) : {} -> {}".format(
options["M_1"],
options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
np.log10(options["P"]),
options['ecc'] if options.get("ecc", None) else "N/A",
options["M_3"],
np.log10(options["P2"]),
options["ecc2"] if options.get("ecc2", None) else "N/A",
options["M_4"],
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_dict