blob: 8cc46ee01ba20118b7f52694662988ef61418929 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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)
|
