diff options
Diffstat (limited to 'rand')
| -rw-r--r-- | rand/luminosity_func.py | 96 | ||||
| -rw-r--r-- | rand/pointsrc_coord.py | 98 | ||||
| -rw-r--r-- | rand/sphere.py | 57 |
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 + + |