Skip to content
Snippets Groups Projects
Commit ec888d03 authored by David Hendriks's avatar David Hendriks
Browse files

added three-part-powerlaw distribution function

parent 74349adf
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment