diff --git a/binarycpython/utils/probability_distributions.py b/binarycpython/utils/probability_distributions.py
index 5506cf4eba7dd6e11a0aaa3b6bfcbe98ff2661c9..c546dd66b5f62bb0a6f489a91854781639a6493f 100644
--- a/binarycpython/utils/probability_distributions.py
+++ b/binarycpython/utils/probability_distributions.py
@@ -21,11 +21,11 @@ def powerlaw_constant(min_val, max_val, k):
"""
k1 = k + 1.0
- print(
- "Powerlaw consts from {} to {}, k={} where k1={}".format(
- min_val, max_val, k, k1
- )
- )
+ # 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
@@ -42,6 +42,7 @@ def powerlaw(min_val, max_val, k, x):
raise ValueError
if (x < min_val) or (x > max_val):
+ print("input value is out of bounds!")
return 0
else:
@@ -49,13 +50,72 @@ def powerlaw(min_val, max_val, k, x):
# powerlaw
y = const * (x ** k)
- print(
- "Power law from {} to {}: const = {}, y = {}".format(
- min_val, max_val, const, y
- )
- )
+ # 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
# {
@@ -67,29 +127,6 @@ def powerlaw(min_val, max_val, k, x):
# return $C;
# }
-# sub initialize_three_part_power_law_consts
-# {
-# local $SIG{__DIE__} = sub { Carp::confess @_ };
-# # calculate constants for this power law
-# my ($m0,$m1,$m2,$mmax,$p1,$p2,$p3)=@_;
-# #print "INIT three-part power law m0=$m0 m1=$m1 m2=$m2 mmax=$mmax p1=$p1 p2=$p2 p3=$p3\n";
-# my $array=[];
-
-# $$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];
-# }
# sub ktg93_lnspace
# {
@@ -191,36 +228,6 @@ def powerlaw(min_val, max_val, k, x):
# $$opts{p3});
# }
-# sub three_part_power_law
-# {
-# # generalized three part power law, usually used for mass distributions
-# # args are M, M0, M1, M2, MMAX, P1, P2, P3 (powers)
-# my $m=shift;
-
-# # use pre-calculated consts if possible, otherwise calculate them
-# my $consts=$threepart_powerlaw_consts{"@_"} //
-# initialize_three_part_power_law_consts(@_);
-
-# my $p;
-# if($m<$_[1])
-# {
-# $p = $m<$_[0]
-# ? 0.0
-# : $$consts[0]*($m**$_[4]);
-# }
-# elsif($m<$_[3])
-# {
-# $p = $m<$_[2]
-# ? $$consts[1]*($m**$_[5])
-# : $$consts[2]*($m**$_[6]);
-# }
-# else
-# {
-# $p=0.0;
-# }
-# return $p;
-# }
-
# sub gaussian
# {
# # Gaussian distribution function e.g. for Duquennoy + Mayor 1991