diff --git a/binarycpython/utils/IMF.py b/binarycpython/utils/IMF.py
new file mode 100644
index 0000000000000000000000000000000000000000..d9b2c7de3fe45efd6afa7138846bd0a2ca28a818
--- /dev/null
+++ b/binarycpython/utils/IMF.py
@@ -0,0 +1,4 @@
+###
+# File containing functions to calculate probability of existence of star. Copied and inspired by the perl module IMF.pm or Rob Izzard
+
+# TODO: add
\ No newline at end of file
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..245da9f396098f7ae1e8974914eb4857b76750ea
--- /dev/null
+++ b/binarycpython/utils/distribution_functions.py
@@ -0,0 +1,356 @@
+import math
+import binarycpython.utils.useful_funcs
+
+###
+# File containing probability distributions
+# Mostly copied from the perl modules
+
+# TODO: make some things globally present? rob does this in his module. i guess it saves calculations but not sure if im gonna do that now
+# TODO: Add the stuff from the IMF file
+# TODO: call all of these functions to check whether they work
+# TODO: make global constants stuff
+# TODO: make description of module submodule
+
+log_ln_converter = (1.0/math.log(10.0))
+
+def flat(parameter):
+ """
+ Dummt distribution function that returns 1
+ """
+ return 1
+
+
+def number(value):
+ """
+ Dummy distribution function that returns the input
+ """
+ return value
+
+
+def powerlaw_constant(min_val, max_val, k):
+ """
+ Function that returns the constant to normalise a powerlaw
+ """
+
+ 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, max_val, k, x):
+ """
+ Single powerlaw with index k at x from min to max
+ """
+
+ # Handle faulty value
+ if k == -1:
+ print("wrong value for k")
+ raise ValueError
+
+ if (x < min_val) or (x > max_val):
+ print("input value is out of bounds!")
+ return 0
+
+ else:
+ const = powerlaw_constant(min_val, max_val, k)
+
+ # powerlaw
+ y = const * (x ** k)
+ # print(
+ # "Power law from {} to {}: const = {}, y = {}".format(
+ # min_val, max_val, const, y
+ # )
+ # )
+ return y
+
+def calculate_constants_three_part_powerlaw(m0, m1, m2, m_max, p1, p2, p3):
+ 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))
+ # local $SIG{__DIE__} = sub { Carp::confess @_ }; ?
+
+ 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)))
+
+ # print(array_constants_three_part_powerlaw)
+ 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, M0, M1, M2, M_MAX, P1, P2, P3):
+ """
+ Generalized three-part power law, usually used for mass distributions
+ """
+
+ # 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):
+ p = 0 # Below 0
+ elif (M0 < M <= M1):
+ p = three_part_powerlaw_constants[0] * (M ** P1) # Between M0 and M1
+ elif (M1 < M <= M2):
+ p = three_part_powerlaw_constants[1] * (M ** P2) # Between M1 and M2
+ elif (M2 < M <= M_MAX):
+ p = three_part_powerlaw_constants[2] * (M ** P3) # Between M2 and M_MAX
+ else:
+ p = 0 # Above M_MAX
+
+ return p
+
+def const(min_bound, max_bound, val=None):
+ """
+ a constant distribution function between min=$_[0] and max=$_[1]
+ """
+
+ if val:
+ if not (min_bound < val <= max_bound):
+ print("out of bounds")
+ p = 0
+ return p
+ p = (1.0/(min_bound - max_bound))
+ return p
+
+def set_opts(opts, newopts):
+ """
+ Function to take a default dict and override it with newer values.
+ """
+
+ # DONE: put in check to make sure that the newopts keys are contained in opts
+
+ if newopts:
+ for opt in newopts.keys():
+ if opt in opts.keys():
+ opts[key] = newopts[key]
+
+ return newopts
+
+#####
+# Mass distributions
+#####
+
+def Kroupa2001(m, newopts):
+ """
+ probability distribution function for kroupa 2001 IMF
+
+ Input: Mass, and dict of new options
+ """
+
+ # Default params and override them
+ opts = set_opts({'m0':0.1, 'm1':0.5, 'm2':1, 'mmax':100, 'p1':-1.3, 'p2':-2.3, 'p3':-2.3}, newopts)
+
+ return three_part_powerlaw(m, opts['m0'], opts['m0'], opts['m2'], opts['m0'], opts['m0'], opts['m0'], opts['m0'])
+
+def ktg93():
+ """
+ Wrapper for mass distribution of KTG93
+ """
+ # 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)
+
+ return three_part_powerlaw(m, opts['m0'], opts['m0'], opts['m2'], opts['m0'], opts['m0'], opts['m0'], opts['m0'])
+
+
+# sub ktg93_lnspace
+# {
+# # wrapper for KTG93 on a ln(m) grid
+# my $m=$_[0];
+# return ktg93(@_) * $m;
+# }
+
+# sub gaussian
+# {
+# # Gaussian distribution function e.g. for Duquennoy + Mayor 1991
+
+# # location (X value), mean and sigma, min and max range
+# my ($x,$mean,$sigma,$gmin,$gmax) = @_;
+
+# my $p;
+# if($x<$gmin || $x>$gmax)
+# {
+# $p=0.0;
+# }
+# else
+# {
+# # normalize over given range
+# my $mult= $gauss_consts{$mean}{$sigma} //
+# gaussian_normalizing_const($mean,$sigma,$gmin,$gmax);
+# $p = $mult * gaussian_func($x,$mean,$sigma);
+# }
+# return $p;
+# }
+
+# sub gaussian_normalizing_const
+# {
+# # first time: calculate multiplier for given $mean, $sigma
+# my ($mean,$sigma,$gmin,$gmax) = @_;
+# my $ptot=0.0;
+# my $d=($gmax-$gmin)/1000.0;
+# for(my $y=$gmin;$y<=$gmax;$y+=$d)
+# {
+# $ptot += $d * gaussian_func($y,$mean,$sigma);
+# }
+# $gauss_consts{$mean}{$sigma}=$ptot; # save for later
+# return $gauss_consts{$mean}{$sigma};
+# }
+
+# sub gaussian_func
+# {
+# # local function
+# # three args: x, mean, sigma
+# my $r=1.0/$_[2]; # 1/sigma
+# my $y=($_[0]-$_[1])*$r;
+# return(GAUSSIAN_PREFACTOR*$r*exp(-0.5*$y*$y));
+# }
+
+
+##
+# Binary fractions
+##
+def Arenou2010_binary_fraction(m):
+ # 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
+
+ return 0.8388 * math.tanh(0.688 * m + 0.079)
+# print(Arenou2010_binary_fraction(0.4))
+
+def raghavan2010_binary_fraction(m):
+ """
+ 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
+ """
+
+ return min(1.0,
+ max(
+ (m**0.1) * (5.12310e-01) + (-1.02070e-01),
+ (1.10450e+00) * (m**(4.93670e-01)) + (-6.95630e-01)
+ )
+ )
+##
+# Period distributions
+##
+
+def sana12(M1, M2, a, P, amin, amax, x0, x1, p): # TODO: ? wtf. vague input
+ """
+ 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
+ """
+
+ 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 normalizations 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))
+ xmin = log_ln_converter * math.log(calc_period_from_sep(m1, m2, amax))
+
+ #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
+# sana12(10, 2, 10, 100, 1, 1000, math.log(10), math.log(1000), 6)
+
+def flatsections(x, opts):
+ """
+
+ opts = list of dicts with settings for the flat sections
+ x: location to calculate the y value
+ TODO: figure out why it has to be a list of dict. why not just 1
+ """
+
+ 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}]))
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index 71590f1202f52168fcc0df66d7b56c07d9891828..779f0f7a2c9a1cff884af07fb7635814790f6f1b 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -10,7 +10,6 @@ from binarycpython.utils.custom_logging_functions import (
create_and_load_logging_function,
)
-
def parse_binary_c_version_info(version_info_string):
version_info_dict = {}
diff --git a/binarycpython/utils/probability_distributions.py b/binarycpython/utils/probability_distributions.py
deleted file mode 100644
index c546dd66b5f62bb0a6f489a91854781639a6493f..0000000000000000000000000000000000000000
--- a/binarycpython/utils/probability_distributions.py
+++ /dev/null
@@ -1,274 +0,0 @@
-# TODO: make some things globally present? rob does this in his module. i guess it saves calculations but not sure if im gonna do that now
-
-
-def flat(parameter):
- """
- Dummt distribution function that returns 1
- """
- return 1
-
-
-def number(value):
- """
- Dummy distribution function that returns the input
- """
- return value
-
-
-def powerlaw_constant(min_val, max_val, k):
- """
- Function that returns the constant to normalise a powerlaw
- """
-
- 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, max_val, k, x):
- """
- Single powerlaw with index k at x from min to max
- """
-
- # Handle faulty value
- if k == -1:
- print("wrong value for k")
- raise ValueError
-
- if (x < min_val) or (x > max_val):
- print("input value is out of bounds!")
- return 0
-
- else:
- const = powerlaw_constant(min_val, max_val, k)
-
- # powerlaw
- y = const * (x ** k)
- # print(
- # "Power law from {} to {}: const = {}, y = {}".format(
- # min_val, max_val, const, y
- # )
- # )
- return y
-
-def calculate_constants_three_part_powerlaw(m0, m1, m2, m_max, p1, p2, p3):
- 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))
- # local $SIG{__DIE__} = sub { Carp::confess @_ }; ?
-
- 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)))
-
- # print(array_constants_three_part_powerlaw)
- 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, M0, M1, M2, M_MAX, P1, P2, P3):
- """
- Generalized three-part power law, usually used for mass distributions
- """
-
- three_part_powerlaw_constants = calculate_constants_three_part_powerlaw(M0, M1, M2, M_MAX, P1, P2, P3)
-
- # Check trick with using the global value
-
- if (M < M1):
- if M < M0: # Below lower bound
- p = 0
- else:
- p = three_part_powerlaw_constants[0] * (M ** P1) # Between M0 and M1
- elif (M < M_MAX): # Below M_MAX
- if (M < M2): # Between M1 and M2
- p = three_part_powerlaw_constants[1] * (M ** P2)
- else: # Between M2 and M_MAX
- p = three_part_powerlaw_constants[2] * (M ** P3)
- else: # Above M_MAX
- p = 0
- return p
-
-
-
-
-# calculate_constants_three_part_powerlaw(0.08, 0.5, 1, 100, -1.3, -2.3, -2.3)
-# print(three_part_powerlaw(10, 0.08, 0.5, 1, 100, -1.3, -2.3, -2.3))
-
-
-# sub const
-# {
-# # a constant distribution function between min=$_[0] and max=$_[1]
-# # $_[2] is optional (if given, perform a range check)
-
-# my $C= (defined($_[2]) && (($_[2]<$_[0]) || ($_[2]>$_[1]))) ? 0.0 : (1.0/($_[1]-$_[0]));
-# #printf "CONST $_[0] to $_[1] ($_[2]) $C\n";
-# return $C;
-# }
-
-
-# sub ktg93_lnspace
-# {
-# # wrapper for KTG93 on a ln(m) grid
-# my $m=$_[0];
-# return ktg93(@_) * $m;
-# }
-
-# sub setopts
-# {
-# # wrapper to take a hash (reference) of options,
-# # override with $newopts (where appropriate) and
-# # return the hash of options
-# my $opts=$_[0];
-# my $newopts=$_[1];
-# if(defined $newopts)
-# {
-# # newopts is a hash of new options
-# foreach my $opt (keys %$newopts)
-# {
-# # overwrite opt
-# $$opts{$opt} = $$newopts{$opt};
-# }
-# }
-# return $opts;
-# }
-
-# sub Kroupa2001
-# {
-# my $m=$_[0];
-# my $newopts=$_[1];
-
-# # default parameters
-# my $opts=setopts({m0=>0.1,
-# m1=>0.5,
-# m2=>1.0,
-# mmax=>100.0,
-# p1=>-1.3,
-# p2=>-2.3,
-# p3=>-2.3},
-# $newopts);
-
-# return three_part_power_law($m,
-# $$opts{m0},
-# $$opts{m1},
-# $$opts{m2},
-# $$opts{mmax},
-# $$opts{p1},
-# $$opts{p2},
-# $$opts{p3});
-# }
-
-# sub ktg93
-# {
-# # wrapper for the mass distribution of KTG93
-# my ($m,$newopts)=@_;
-
-# 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
-# my $opts=setopts({m0=>0.1,
-# m1=>0.5,
-# m2=>1.0,
-# mmax=>80.0,
-# p1=>-1.3,
-# p2=>-2.2,
-# p3=>-2.7},
-# $newopts);
-
-# return three_part_power_law($m,
-# $$opts{m0},
-# $$opts{m1},
-# $$opts{m2},
-# $$opts{mmax},
-# $$opts{p1},
-# $$opts{p2},
-# $$opts{p3});
-# }
-
-# sub gaussian
-# {
-# # Gaussian distribution function e.g. for Duquennoy + Mayor 1991
-
-# # location (X value), mean and sigma, min and max range
-# my ($x,$mean,$sigma,$gmin,$gmax) = @_;
-
-# my $p;
-# if($x<$gmin || $x>$gmax)
-# {
-# $p=0.0;
-# }
-# else
-# {
-# # normalize over given range
-# my $mult= $gauss_consts{$mean}{$sigma} //
-# gaussian_normalizing_const($mean,$sigma,$gmin,$gmax);
-# $p = $mult * gaussian_func($x,$mean,$sigma);
-# }
-# return $p;
-# }
-
-# sub gaussian_normalizing_const
-# {
-# # first time: calculate multiplier for given $mean, $sigma
-# my ($mean,$sigma,$gmin,$gmax) = @_;
-# my $ptot=0.0;
-# my $d=($gmax-$gmin)/1000.0;
-# for(my $y=$gmin;$y<=$gmax;$y+=$d)
-# {
-# $ptot += $d * gaussian_func($y,$mean,$sigma);
-# }
-# $gauss_consts{$mean}{$sigma}=$ptot; # save for later
-# return $gauss_consts{$mean}{$sigma};
-# }
-
-# sub gaussian_func
-# {
-# # local function
-# # three args: x, mean, sigma
-# my $r=1.0/$_[2]; # 1/sigma
-# my $y=($_[0]-$_[1])*$r;
-# return(GAUSSIAN_PREFACTOR*$r*exp(-0.5*$y*$y));
-# }
diff --git a/binarycpython/utils/useful_funcs.py b/binarycpython/utils/useful_funcs.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac83dd89530d6c812b5609e16933e813b025dd09
--- /dev/null
+++ b/binarycpython/utils/useful_funcs.py
@@ -0,0 +1,93 @@
+############################################################
+# Collection of useful functions.
+
+# Part of this is copied/inspired by Robs
+# Rob's binary_stars module
+
+# Has functions to convert period to separation and vice versa.
+# calc_period_from_sep($m1,$m2,$sep) calculate the period given the separation.
+# calc_sep_from_period($m1,$m2,per) does the inverse.
+# M1,M2,separation are in solar units, period in days.
+# rzams($m,$z) gives you the ZAMS radius of a star
+# ZAMS_collision($m1,$m2,$e,$sep,$z) returns 1 if stars collide on the ZAMS
+###
+# TODO: check whether these are correct
+
+AURSUN = 2.150445198804013386961742071435e+02
+YEARDY = 3.651995478818308811241877265275e+02
+
+def calc_period_from_sep(M1, M2, sep):
+ """
+ calculate period from separation
+ args : M1, M2, separation (Rsun)
+ """
+
+ return YEARDY*(sep/AURSUN)*sqrt(sep/(AURSUN*(M1+M2)));
+
+def calc_sep_from_period(M1, M2, period):
+ """
+ inverse of the above function
+ args : M1, M2, period (days)
+ """
+
+ return AURSUN*(period*period*(M1+M2)/(YEARDY*YEARDY))**(1.0/3.0);
+
+def roche_lobe(q):
+ """
+ A function to evaluate R_L/a(q), Eggleton 1983.
+ """
+
+ p = q**(1.0/3.0);
+ return(0.49*p*p/(0.6*p*p + log(1.0+p)));
+
+def ragb(m, z):
+ """
+ Function to calculate radius of a star at first thermal pulse as a function of mass (z=0.02)
+ """
+ # Z=0.02 only, but also good for Z=0.001
+ return m*40.0+20.0; # in Rsun
+
+def zams_collission(m1, m2, sep, e, z):
+ """
+ # given m1,m2, separation and eccentricity (and metallicity)
+ # determine if two stars collide on the ZAMS
+ """
+ # calculate periastron distance
+ peri_distance=(1.0-e) * sep
+
+ # calculate star radii
+ r1=rzams(m1, z);
+ r2=rzams(m2, z);
+
+ if r1 + r2 > peri_distance:
+ return 1
+ else:
+ return 0
+
+def rzams(m, z):
+ """
+ Function to determine the radius of a ZAMS star as a function of m and z:
+ """
+ lzs = math.log10(z/0.02)
+
+ #
+ # A function to evaluate Rzams
+ # ( from Tout et al., 1996, MNRAS, 281, 257 ).
+ #
+
+ p = {}
+ xz = (666, 0.39704170, -0.32913574, 0.34776688, 0.37470851, 0.09011915, 8.52762600,-24.41225973, 56.43597107, 37.06152575, 5.45624060, 0.00025546, -0.00123461, -0.00023246, 0.00045519, 0.00016176, 5.43288900, -8.62157806, 13.44202049, 14.51584135, 3.39793084, 5.56357900,-10.32345224, 19.44322980, 18.97361347, 4.16903097, 0.78866060, -2.90870942, 6.54713531, 4.05606657, 0.53287322, 0.00586685, -0.01704237, 0.03872348, 0.02570041, 0.00383376, 1.71535900, 0.62246212, -0.92557761, -1.16996966,-0.30631491, 6.59778800, -0.42450044,-12.13339427,-10.73509484,-2.51487077, 10.08855000, -7.11727086,-31.67119479,-24.24848322,-5.33608972, 1.01249500, 0.32699690, -0.00923418, -0.03876858,-0.00412750, 0.07490166, 0.02410413, 0.07233664, 0.03040467, 0.00197741, 0.01077422, 3.08223400, 0.94472050, -2.15200882, -2.49219496,-0.63848738, 17.84778000, -7.45345690,-48.96066856,-40.05386135,-9.09331816, 0.00022582, -0.00186899, 0.00388783, 0.00142402,-0.00007671);
+
+ p['8'] = xz[36]+lzs*(xz[37]+lzs*(xz[38]+lzs*(xz[39]+lzs*xz[40])));
+ p['9'] = xz[41]+lzs*(xz[42]+lzs*(xz[43]+lzs*(xz[44]+lzs*xz[45])));
+ p['10'] = xz[46]+lzs*(xz[47]+lzs*(xz[48]+lzs*(xz[49]+lzs*xz[50])));
+ p['11'] = xz[51]+lzs*(xz[52]+lzs*(xz[53]+lzs*(xz[54]+lzs*xz[55])));
+ p['12'] = xz[56]+lzs*(xz[57]+lzs*(xz[58]+lzs*(xz[59]+lzs*xz[60])));
+ p['13'] = xz[61];
+ p['14'] = xz[62]+lzs*(xz[63]+lzs*(xz[64]+lzs*(xz[65]+lzs*xz[66])));
+ p['15'] = xz[67]+lzs*(xz[68]+lzs*(xz[69]+lzs*(xz[70]+lzs*xz[71])));
+ p['16'] = xz[72]+lzs*(xz[73]+lzs*(xz[74]+lzs*(xz[75]+lzs*xz[76])));
+
+ m195 = m**19.5
+ rad = ((p['8']*(m**2.5) + p['9']*(m**6.5) + p['10']*(m**11) + p['11']*(m**19)+p['12']*m195)/(p['13'] + p['14']*(m*m) + p['15']*(m**8.5) + (m**18.5) + p['16']*m195))
+ return rad
\ No newline at end of file
diff --git a/tests/population/plot_scaling.py b/tests/population/plot_scaling.py
index fe5b04174343082ac19815c7ebf95784b5e73eaf..f0594e227d43f032d44db91015e82b644373a474 100644
--- a/tests/population/plot_scaling.py
+++ b/tests/population/plot_scaling.py
@@ -6,12 +6,19 @@ import numpy as np
def calc_mean_and_std(arr):
return np.mean(arr), np.std(arr)
-# Readout
+
+
+# Configure
result_file = 'comparison_result.dat'
result_file = 'comparison_result_astro1.dat'
name_testcase = 'Astro1'
+
+
+
+
+#
results = []
with open(result_file, 'r') as f:
@@ -101,7 +108,9 @@ plt.ylim(0, max_speedup + 2)
plt.grid()
plt.legend()
-plt.savefig('name_testcase'+'.png')
-plt.savefig('name_testcase'+'.pdf')
+
+plt.savefig('scaling_{}.png'+'.png')
+plt.savefig('scaling_{}.png'+'.pdf')
+plt.savefig('scaling_{}.svg'+'.pdf')
plt.show()
\ No newline at end of file