adding distribution functions

binarycpython/utils/distribution_functions.py

@@ -149,6 +149,53 @@ def set_opts(opts, newopts):
 
     return newopts
 
+def gaussian(x, mean, sigma, gmin, gmax):
+    """
+    Gaussian distribution function. used for e..g Duquennoy + Mayor 1991
+    
+    Input: location, mean, sigma, min and max:
+    """
+    # # location (X value), mean and sigma, min and max range
+    # my ($x,$mean,$sigma,$gmin,$gmax) = @_;
+
+    if ((x < gmin) or (x > gmax)):
+        p = 0
+    else:
+        # normalize over given range
+        # TODO: add loading into global var
+        normalisation = gaussian_normalizing_const(mean, sigma, gmin, gmax)
+        p = normalisation * gaussian_func(x, mean, sigma)
+
+    return p
+
+def gaussian_normalizing_const(mean, sigma, gmin, gmax):
+    """
+    Function to calculate the normalisation constant for the gaussian
+    """
+
+    # First time; calculate multipllier 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, mean, sigma):
+    """
+    Function to evaluate a gaussian at a given point
+    """
+    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)
+
 #####
 # Mass distributions
 ##### 
@@ -208,52 +255,6 @@ def ktg93():
 #     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
 ## 
@@ -285,6 +286,8 @@ def raghavan2010_binary_fraction(m):
             (1.10450e+00) * (m**(4.93670e-01)) + (-6.95630e-01)
             )
         )
+# print(raghavan2010_binary_fraction(2))
+
 ## 
 # Period distributions
 ## 
@@ -329,7 +332,7 @@ def sana12(M1, M2, a, P, amin, amax, x0, x1, p): # TODO: ? wtf. vague input
             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)
+# print(sana12(10, 2, 10, 100, 1, 1000, math.log(10), math.log(1000), 6))
 
 def flatsections(x, opts):
     """
@@ -353,4 +356,4 @@ def flatsections(x, opts):
     y = y * c
     # print("flatsections gives C={}: y={}\n",c,y)
     return y
-print(flatsections(1, [{'min': 0, 'max': 2, 'height': 3}]))
+# print(flatsections(1, [{'min': 0, 'max': 2, 'height': 3}]))