Population class extension: distribution_functions module

Module containing the predefined distribution functions

The user can use any of these distribution functions to generate probability distributions for sampling populations

To add custom functions you can take any function and add it to the class instance before running the code. See https://stackoverflow.com/a/28060251 for some tips on how to do that

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, Mathieu eccentricity
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

class binarycpython.utils.population_extensions.distribution_functions.distribution_functions(**kwargs)[source]

Bases: object

Extension for the Population grid object that contains the distribution functions

Arenou2010_binary_fraction(m)[source]

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

Parameters: m (Union[int, float]) – mass to evaluate the distribution at
Return type: Union[int, float]
Returns: binary fraction at m

Izzard2012_period_distribution(P, M1, log10Pmin=- 1.0)[source]

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

Parameters

P (Union[int, float]) – period
M1 (Union[int, float]) – Primary star mass
log10Pmin (Union[int, float]) – minimum period in base log10 (optional)

Return type

Union[int, float]

Returns

‘probability’ of interpolated distribution function at P and M1

Kroupa2001(m, newopts=None)[source]

Probability distribution function for Kroupa 2001 IMF.

The (default) values for this is:

default = {
   "m0": 0.1,
   "m1": 0.5,
   "m2": 1,
   "mmax": 100,
   "p1": -1.3,
   "p2": -2.3,
   "p3": -2.3
}

Parameters

m (Union[int, float]) – mass to evaluate the distribution at
newopts (Optional[dict]) – optional dict to override the default values.

Return type

Union[int, float]

Returns

‘probability’ of distribution function evaluated at m

Moe_di_Stefano_2017_multiplicity_fractions(options, verbosity=0)[source]

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 quite low.

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

Moe_di_Stefano_2017_pdf(options, verbosity=0)[source]

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)

build_q_table(options, m, p, verbosity=0)[source]

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 Moe and di Stefano 2017 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 TODO: fix description to be correct for python

calc_P_integral(options, min_logP, max_logP, integrals_string, interpolator_name, mass_string, verbosity=0)[source]

Function to calculate the P integral

We need to renormalise this because min_per > 0, and not all periods should be included

calc_e_integral(options, integrals_string, interpolator_name, mass_string, period_string, verbosity=0)[source]

Function to calculate the e integral

We need to renormalise this because min_per > 0, and not all periods should be included

calc_total_probdens(prob_dict)[source]: Function to calculate the total probability density

calculate_constants_three_part_powerlaw(m0, m1, m2, m_max, p1, p2, p3)[source]

Function to calculate the constants for a three-part power law

TODO: use the power law_constant function to calculate all these values

Parameters

m0 (Union[int, float]) – lower bound mass
m1 (Union[int, float]) – second boundary, between the first slope and the second slope
m2 (Union[int, float]) – third boundary, between the second slope and the third slope
m_max (Union[int, float]) – upper bound mass
p1 (Union[int, float]) – first slope
p2 (Union[int, float]) – second slope
p3 (Union[int, float]) – third slope

Return type

Union[int, float]

Returns

array of normalisation constants

const_distribution(min_bound, max_bound, val=None)[source]

a constant distribution function between min=min_bound and max=max_bound.

Parameters

min_bound (Union[int, float]) – lower bound of the range
max_bound (Union[int, float]) – upper bound of the range

Returns

returns 0

Return type

returns the value of 1/(max_bound-min_bound). If val is provided, it will check whether min_bound < val <= max_bound. if not

cosmic_SFH_madau_dickinson2014(z)[source]

Cosmic star formation history distribution from Madau & Dickonson 2014 (https://arxiv.org/pdf/1403.0007.pdf)

Parameters: z – redshift
Returns: Cosmic star formation rate in Solar mass year^-1 mega parsec^-3

duquennoy1991(logper)[source]

Period distribution from Duquennoy + Mayor 1991. Evaluated the function self.gaussian(logper, 4.8, 2.3, -2, 12)

Parameters: logper (Union[int, float]) – logarithm of period to evaluate the distribution at
Return type: Union[int, float]
Returns: ‘probability’ at self.gaussian(logper, 4.8, 2.3, -2, 12)

fill_data(sample_values, data_dict)[source]

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

flat()[source]

Dummy distribution function that returns 1

Returns: 1
Return type: a flat uniform distribution

flatsections(x, opts)[source]

Function to generate flat distributions, possibly in multiple sections

Parameters

x (float) – mass ratio value
opts (dict) – list containing the flat sections. Which are themselves dictionaries, with keys “max”: upper bound, “min”: lower bound and “height”: value

Return type

Union[float, int]

Returns

probability of that mass ratio.

gaussian(x, mean, sigma, gmin, gmax)[source]

Gaussian distribution function. used for e.g. Duquennoy + Mayor 1991

Parameters

x (Union[int, float]) – location at which to evaluate the distribution
mean (Union[int, float]) – mean of the Gaussian
sigma (Union[int, float]) – standard deviation of the Gaussian
gmin (Union[int, float]) – lower bound of the range to calculate the probabilities in
gmax (Union[int, float]) – upper bound of the range to calculate the probabilities in

Return type

Union[int, float]

Returns

‘probability’ of the Gaussian distribution between the boundaries, evaluated at x

gaussian_func(x, mean, sigma)[source]

Function to evaluate a Gaussian at a given point, but this time without any boundaries.

Parameters

x (Union[int, float]) – location at which to evaluate the distribution
mean (Union[int, float]) – mean of the Gaussian
sigma (Union[int, float]) – standard deviation of the Gaussian

Return type

Union[int, float]

Returns

value of the Gaussian at x

gaussian_normalizing_const(mean, sigma, gmin, gmax)[source]

Function to calculate the normalisation constant for the Gaussian

Parameters

mean (Union[int, float]) – mean of the Gaussian
sigma (Union[int, float]) – standard deviation of the Gaussian
gmin (Union[int, float]) – lower bound of the range to calculate the probabilities in
gmax (Union[int, float]) – upper bound of the range to calculate the probabilities in

Return type

Union[int, float]

Returns

normalisation constant for the Gaussian distribution(mean, sigma) between gmin and gmax

get_integration_constant_q(q_interpolator, tmp_table, qdata, verbosity=0)[source]: Function to integrate the q interpolator and return the integration constant

get_max_multiplicity(multiplicity_array)[source]: Function to get the maximum multiplicity

imf_chabrier2003(m)[source]

Probability distribution function for IMF of Chabrier 2003 PASP 115:763-795

Parameters: m (Union[int, float]) – mass to evaluate the distribution at
Return type: Union[int, float]
Returns: ‘probability’ of distribution function evaluated at m

imf_scalo1986(m)[source]

Probability distribution function for Scalo 1986 IMF (defined up until 80Msol): self.three_part_powerlaw(m, 0.1, 1.0, 2.0, 80.0, -2.35, -2.35, -2.70)

Parameters: m (Union[int, float]) – mass to evaluate the distribution at
Return type: Union[int, float]
Returns: ‘probability’ of distribution function evaluated at m

imf_scalo1998(m)[source]

From Scalo 1998

Probability distribution function for Scalo 1998 IMF (defined up until 80Msol): self.three_part_powerlaw(m, 0.1, 1.0, 10.0, 80.0, -1.2, -2.7, -2.3)

Parameters: m (Union[int, float]) – mass to evaluate the distribution at
Return type: Union[int, float]
Returns: ‘probability’ of distribution function evaluated at m

imf_tinsley1980(m)[source]

Probability distribution function for Tinsley 1980 IMF (defined up until 80Msol): self.three_part_powerlaw(m, 0.1, 2.0, 10.0, 80.0, -2.0, -2.3, -3.3)

Parameters: m (Union[int, float]) – mass to evaluate the distribution at
Return type: Union[int, float]
Returns: ‘probability’ of distribution function evaluated at m

interpolate_in_mass_izzard2012(M, high, low)[source]

Function to interpolate in mass

TODO: fix this function. TODO: describe the args high: at M=16.3 low: at 1.15

Parameters

M (Union[int, float]) – mass
high (Union[int, float]) –
low (Union[int, float]) –

Returns:

Return type: Union[int, float]

ktg93(m, newopts=None)[source]

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}

Parameters

m (Union[int, float]) – mass to evaluate the distribution at
newopts (Optional[dict]) – optional dict to override the default values.

Return type

Union[int, float]

Returns

‘probability’ of distribution function evaluated at m

linear_extrapolation_q(qs, indices, qlimit, qdata, end_index, verbosity=0)[source]: Function to do the linear extrapolation for q.

merge_multiplicities(result_array, max_multiplicity, verbosity=0)[source]

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

number(value)[source]

Dummy distribution function that returns the input

Parameters: value (Union[int, float]) – the value that will be returned by this function.
Return type: Union[int, float]
Returns: the value that was provided

poisson(lambda_val, n, nmax=None)[source]: Function that calculates the Poisson value and normalises TODO: improve the description

powerlaw(min_val, max_val, k, x)[source]

Single power law with index k at x from min to max

Parameters

min_val (Union[int, float]) – lower bound of the power law
max_val (Union[int, float]) – upper bound of the power law
k (Union[int, float]) – slope of the power law
x (Union[int, float]) – position at which we want to evaluate

Return type

Union[int, float]

Returns

probability at the given position(x)

powerlaw_constant(min_val, max_val, k)[source]

Function that returns the constant to normalise a power law

TODO: what if k is -1?

Parameters

min_val (Union[int, float]) – lower bound of the range
max_val (Union[int, float]) – upper bound of the range
k (Union[int, float]) – power law slope

Return type

Union[int, float]

Returns

constant to normalise the given power law between the min_val and max_val range

powerlaw_constant_nocache(min_val, max_val, k)[source]

Function that returns the constant to normalise a power law

TODO: what if k is -1?

Parameters

min_val (Union[int, float]) – lower bound of the range
max_val (Union[int, float]) – upper bound of the range
k (Union[int, float]) – power law slope

Return type

Union[int, float]

Returns

constant to normalise the given power law between the min_val and max_val range

powerlaw_extrapolation_q(qdata, qs, indices)[source]: Function to do the power-law extrapolation at the lower end of the q range

raghavan2010_binary_fraction(m)[source]

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

Parameters: m (Union[int, float]) – mass to evaluate the distribution at
Return type: Union[int, float]
Returns: binary fraction at m

sana12(M1, M2, a, P, amin, amax, x0, x1, p)[source]

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! Half of the input here can be taken out and calculated within the function itself.

Parameters

M1 (Union[int, float]) – Mass of primary
M2 (Union[int, float]) – Mass of secondary
a (Union[int, float]) – separation of binary
P (Union[int, float]) – period of binary
amin (Union[int, float]) – minimum separation of the distribution (lower bound of the range)
amax (Union[int, float]) – maximum separation of the distribution (upper bound of the range)
x0 (Union[int, float]) – log of minimum period of the distribution (lower bound of the range)
x1 (Union[int, float]) – log of maximum period of the distribution (upper bound of the range)
p (Union[int, float]) – slope of the distribution

Return type

Union[int, float]

Returns

‘probability’ of orbital period P given the other parameters

three_part_powerlaw(m, m0, m1, m2, m_max, p1, p2, p3)[source]

Generalised three-part power law, usually used for mass distributions

Parameters

m (Union[int, float]) – mass at which we want to evaluate the distribution.
m0 (Union[int, float]) – lower bound mass
m1 (Union[int, float]) – second boundary, between the first slope and the second slope
m2 (Union[int, float]) – third boundary, between the second slope and the third slope
m_max (Union[int, float]) – upper bound mass
p1 (Union[int, float]) – first slope
p2 (Union[int, float]) – second slope
p3 (Union[int, float]) – third slope

Return type

Union[int, float]

Returns

‘probability’ at given mass m