add all scripts/tools

3 files changed, 251 insertions, 0 deletions

diff --git a/rand/luminosity_func.py b/rand/luminosity_func.py
new file mode 100644
index 0000000..8cc46ee
--- /dev/null
+++ b/rand/luminosity_func.py
@@ -0,0 +1,96 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+#
+# Aaron LI
+# 2015/07/01
+#
+
+"""
+Generate random numbers (i.e., fluxes) with respect to the
+provided luminosity function.
+"""
+
+import numpy as np
+import random
+
+def luminosity_func(Lx, N0=1.0):
+    """
+    The *cumulative* luminosity function: N(>=L)
+    The number of objects with luminosities >= L(x) for each L(x).
+    """
+    # broken power-law model (Xu et al. 2005)
+    # Nx = (1) N0 * (Lx/L_b)^(-alpha_l);  for Lx <= L_b
+    #      (2) N0 * (Lx/L_b)^(-alpha_h);  for Lx > L_b
+    L_b = 4.4e38 # break point (erg/s) (+2.0/-1.4)
+    alpha_h = 2.28 # (+1.72/-0.53)
+    alpha_l = 1.08 # (+0.15/-0.33)
+    if isinstance(Lx, np.ndarray):
+        Nx = np.zeros(Lx.shape)
+        Nx[Lx <= 0] = 0.0
+        Nx[Lx <= L_b] = N0 * (Lx[Lx <= L_b] / L_b)**(-alpha_l)
+        Nx[Lx > L_b] = N0 * (Lx[Lx > L_b] / L_b)**(-alpha_h)
+    else:
+        # Lx is a single number
+        if Lx <= 0.0:
+            Nx = 0.0
+        elif Lx <= L_b:
+            Nx = N0 * (Lx/L_b)**(-alpha_l)
+        else:
+            Nx = N0 * (Lx/L_b)**(-alpha_h)
+    return Nx
+
+
+def luminosity_density(Lx, N0=1.0):
+    """
+    Function of number density at luminosity at Lx. => PDF
+
+    PDF(Lx) = - d(luminosity_func(Lx) / d(Lx)
+    """
+    L_b = 4.4e38 # break point (erg/s) (+2.0/-1.4)
+    alpha_h = 2.28 # (+1.72/-0.53)
+    alpha_l = 1.08 # (+0.15/-0.33)
+    if isinstance(Lx, np.ndarray):
+        Px = np.zeros(Lx.shape)
+        Px[Lx<=0] = 0.0
+        Px[Lx<=L_b] = N0 * (alpha_l/L_b) * (Lx[Lx<=L_b] / L_b)**(-alpha_l-1)
+        Px[Lx>L_b] = N0 * (alpha_h/L_b) * (Lx[Lx>L_b] / L_b)**(-alpha_h-1)
+    else:
+        # Lx is a single number
+        if Lx <= 0.0:
+            Px = 0.0
+        elif Lx <= L_b:
+            Px = N0 * (alpha_l/L_b) * (Lx/L_b)**(-alpha_l-1)
+        else:
+            Px = N0 * (alpha_h/L_b) * (Lx/L_b)**(-alpha_h-1)
+    return Px
+
+
+def luminosity_pdf(Lx):
+    """
+    Probability density function
+    """
+    h = 1e-5 * Lx # step size for numerical deviation
+    p = - (luminosity_func(Lx+0.5*h) - luminosity_func(Lx-0.5*h)) / h
+    return p
+
+
+def sampler(min, max, number=1):
+    """
+    Generate a sample of luminosity values within [min, max] from
+    the above luminosity distribution.
+    """
+    # Get the maximum value of the density function
+    M = luminosity_density(min)
+    results = []
+    for i in range(number):
+        while True:
+            u = random.random() * M
+            y = random.random() * (max-min) + min
+            if u <= luminosity_density(y):
+                results.append(y)
+                break
+    if len(results) == 1:
+        return results[0]
+    else:
+        return np.array(results)
+
diff --git a/rand/pointsrc_coord.py b/rand/pointsrc_coord.py
new file mode 100644
index 0000000..1da9be2
--- /dev/null
+++ b/rand/pointsrc_coord.py
@@ -0,0 +1,98 @@
+# -*- coding: utf-8 -*-
+#
+# Aaron LI
+# 2015/07/01
+#
+
+"""
+Generate random coordinates for point sources with respect to the r^{1/4}
+distribution.
+"""
+
+import numpy as np
+import random
+
+
+def cdf(r, N0=1.0):
+    """
+    Cumulative distribution function of the number of point sources.
+
+    r^{1/4} distribution law: de Vaucouleurs 1948
+    """
+    return N0 * r**(1.0/4.0)
+
+
+def pdf(r, N0=1.0):
+    """
+    Density function of the number of point sources.
+
+    pdf = d(pdf) / d(r)
+    """
+    if isinstance(r, np.ndarray):
+        p = np.zeros(r.shape)
+        p[r<=0.0] = 0.0
+        p[r>0.0] = 0.25 * N0 * r[r>0.0]**(-3.0/4.0)
+    else:
+        if r <= 0.0:
+            p = 0.0
+        else:
+            p = 0.25 * N0 * r**(-3.0/4.0)
+    return p
+
+
+def sampler(min, max, number=1):
+    """
+    Generate a sample of coordinates (only r) within [min, max] from
+    the above density distribution.
+
+    min, max: the minimum and maximum r values (in degree)
+    """
+    # Get the maximum value of the density function
+    M = pdf(min)
+    results = []
+    for i in range(number):
+        while True:
+            u = random.random() * M
+            y = random.random() * (max-min) + min
+            if u <= pdf(y):
+                results.append(y)
+                break
+    if len(results) == 1:
+        return results[0]
+    else:
+        return np.array(results)
+
+
+def add_angle(r):
+    """
+    Add angle for each r value to make up a coordinate of a polar coordinate.
+    """
+    coords = []
+    for ri in r:
+        theta = random.random() * 360
+        coords.append((ri, theta))
+    if len(coords) == 1:
+        return coords[0]
+    else:
+        return coords
+
+
+def to_radec(coords, xc=0, yc=0):
+    """
+    Convert the generated coordinates to (ra, dec) (unit: degree).
+
+    xc, yc: the center coordinate (ra, dec)
+    """
+    results = []
+    for r, theta in coords:
+        # FIXME: spherical algebra should be used!!!
+        dx = r * np.cos(theta*np.pi/180)
+        dy = r * np.sin(theta*np.pi/180)
+        x = xc + dx
+        y = yc + dy
+        results.append((x, y))
+    if len(results) == 1:
+        return results[0]
+    else:
+        return results
+
diff --git a/rand/sphere.py b/rand/sphere.py
new file mode 100644
index 0000000..4220766
--- /dev/null
+++ b/rand/sphere.py
@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Randomly pick point on the sphere surface.
+#
+# References:
+# [1] Shpere Poin Picking -- from Wolfram MathWorld
+#     http://mathworld.wolfram.com/SpherePointPicking.html
+# [2] Random Points on a Sphere
+#     https://www.jasondavies.com/maps/random-points/
+#
+# Aaron LI
+# 2015/06/18
+
+__version__ = "0.1.0"
+__date__ = "2015/06/16"
+
+import math
+import random
+
+def sphere_point(n=1, unit="rad"):
+    """
+    Randomly uniformly pick a point on the sphere surface.
+    Using the method "Sphere Point Picking" from Wolfram MathWorld.
+
+    Arguments:
+        n: number of points to be generated
+        unit: unit of output values: rad/deg
+
+    Return:
+        (theta, phi): spherical coordinate (unit: rad).
+                      theta: [0, 2\pi); phi: [0 - \pi]
+        If n > 1, then return a list of (theta, phi)
+    """
+    points = []
+    for i in range(n):
+        u = random.random()
+        v = random.random()
+        theta = 2.0 * math.pi * u
+        phi = math.acos(2.0*v - 1.0)
+        if unit == "deg":
+            theta = rad2deg(theta)
+            phi = rad2deg(phi)
+        points.append((theta, phi))
+    if n == 1:
+        return points[0]
+    else:
+        return points
+
+
+def rad2deg(x):
+    return x * 180.0 / math.pi
+
+def deg2rad(x):
+    return x * math.pi / 180.0
+
+