In [1]:
%matplotlib inline
import numpy as np
from math import *
import matplotlib.pyplot as pyplt
from matplotlib.colors import LogNorm
import scipy.stats as stats
In [2]:
import triNFW
In [3]:
# Set up toy models
halo1 = triNFW.pNFW(mass = 3E14,
concen = 3.0,
zd = 0.3,
overden = 200,
wrt = "crit",
qa = 1.0,
qb = 1.0,
theta = 0.0,
phi = 0.0,)
halo2 = triNFW.pNFW(mass = 3E14,
concen = 3.0,
zd = 0.3,
overden = 200,
wrt = "crit",
qa = 0.4,
qb = 0.7,
theta = pi / 3.2,
phi = -pi/ 1.2,)
In [4]:
# Set up boundary
xmpc_edges = np.linspace(-2,2,1000)
ympc_edges = np.linspace(-2,2,1000)
xmpc_bins = 0.5 * (xmpc_edges[1:] + xmpc_edges[:-1])
ympc_bins = 0.5 * (ympc_edges[1:] + ympc_edges[:-1])
xmpc_mesh, ympc_mesh = np.meshgrid(xmpc_bins, ympc_bins)
In [5]:
# calculate maps
map1 = halo1.Sigma_XY(XX = xmpc_mesh, YY = ympc_mesh)
map2 = halo2.Sigma_XY(XX = xmpc_mesh, YY = ympc_mesh)
In [6]:
pyplt.imshow(map1[::-1],
interpolation = "nearest",
norm = LogNorm(vmin=1E12, vmax=1E16),
extent = [xmpc_edges[0], xmpc_edges[-1], ympc_edges[0], ympc_edges[-1]])
Out[6]:
In [7]:
pyplt.imshow(map2[::-1],
interpolation = "nearest",
norm = LogNorm(vmin=1E12, vmax=1E16),
extent = [xmpc_edges[0], xmpc_edges[-1], ympc_edges[0], ympc_edges[-1]])
Out[7]:
In [8]:
# calculate the psi angle (the third Euler angle) on the plan of sky
print halo2.psi * 180 / pi
In [9]:
# construct P(qa) - eq~17 from Jing & Suto 2002
def Pqa(qa):
# sanitize
qa = np.array(qa, ndmin=1)
return stats.norm.pdf(x = qa, loc = 0.54, scale = 0.113)
# construct P(a/b|qa) - eq~19 from Jing & Suto 2002
def P_aoverb_given_qa(a_over_b, qa):
# sanitize
qa = np.array(qa , ndmin=1)
a_over_b = np.array(a_over_b, ndmin=1)
# mesh it - in the shape of (len(a_over_b), len(qa))
qa_mesh, a_over_b_mesh = np.meshgrid(qa, a_over_b)
# for the case of qa >= 0.5
case1 = 3.0 / (2.0 * (1.0 - qa_mesh)) * ( 1.0 - ( (2 * a_over_b_mesh - 1.0 - qa_mesh)/(1.0 - qa_mesh) )**2 )
case1[ (a_over_b_mesh < qa_mesh) ] = 0.0
# for the case of qa < 0.5
case2 = 3.0 / (2.0 * (1.0 - 0.5 )) * ( 1.0 - ( (2 * a_over_b_mesh - 1.0 - 0.5 )/(1.0 - 0.5 ) )**2 )
case2[ (a_over_b_mesh < 0.5 ) ] = 0.0
# return_me
return_me = np.nan * np.ones(np.shape(qa_mesh))
return_me[ (qa_mesh >= 0.5) ] = np.copy( case1[ (qa_mesh >= 0.5) ] )
return_me[ (qa_mesh < 0.5) ] = np.copy( case2[ (qa_mesh < 0.5) ] )
# return
return return_me
In [10]:
# check whether the pdf seems right
qa_edges = np.linspace(1E-3, 1.0, 1000)
qa_bins = 0.5 * (qa_edges[1:] + qa_edges[:-1])
qa_steps = (qa_edges[1:] - qa_edges[:-1])
pyplt.plot(qa_bins, Pqa(qa_bins) * qa_steps, "k--")
pyplt.show()
a_over_b_edges = np.linspace(1E-3, 1.0, 1000)
a_over_b_bins = 0.5 * ( a_over_b_edges[1:] + a_over_b_edges[:-1] )
a_over_b_steps = ( a_over_b_edges[1:] - a_over_b_edges[:-1] )
pyplt.plot(a_over_b_bins, P_aoverb_given_qa(a_over_b = a_over_b_bins, qa = 0.34) * a_over_b_steps, "r--")
pyplt.show()
In [11]:
# create a pool of clusters
nsample = 10000
# draw qa from a distribution
qa_sample = np.random.normal(loc = 0.54, scale = 0.113, size = nsample)
# draw qb for a given qa - tweaking the shape of array here
qb_sample = np.array([
np.random.choice(a = a_over_b_bins,
p = P_aoverb_given_qa(a_over_b = a_over_b_bins, qa = aa).T[0] / \
np.sum(P_aoverb_given_qa(a_over_b = a_over_b_bins, qa = aa).T[0]),
size = 1)
for aa in qa_sample ]).T[0]
# random angles
theta_sample = np.random.uniform(-pi/2, pi/2, size = nsample)
phi_sample = np.random.uniform(0, pi, size = nsample)
In [12]:
# create halos - let's see how fast it can be
import time
# JS02 halos
t1 = time.time()
halos = [ triNFW.pNFW(mass = 3E14,
concen = 3.0,
zd = 0.3,
overden = 200,
wrt = "crit",
qa = qa_sample[nn],
qb = qb_sample[nn],
theta = theta_sample[nn],
phi = phi_sample[nn],
) for nn in range(nsample) ]
t2 = time.time()
print "#", "time:", t2 - t1, "s"
# Prolate halos
t1 = time.time()
prolate_halos = [ triNFW.pNFW(mass = 3E14,
concen = 3.0,
zd = 0.3,
overden = 200,
wrt = "crit",
qa = qa_sample[nn],
qb = qa_sample[nn],
theta = theta_sample[nn],
phi = phi_sample[nn],
) for nn in range(nsample) ]
t2 = time.time()
print "#", "time:", t2 - t1, "s"
In [13]:
# qproj_sample
qproj_sample = np.array([ the_halo.q_proj for the_halo in halos ])
prolate_qproj_sample = np.array([ the_halo.q_proj for the_halo in prolate_halos ])
# fgoe_sample
fgeo_sample = np.array([the_halo.fgeo for the_halo in halos])
prolate_fgeo_sample = np.array([the_halo.fgeo for the_halo in prolate_halos])
In [15]:
# quantify the distribution
import Qpdf
print "#" * 10
print "#", "JS02 qproj"
Qpdf.quan_sample(qproj_sample, verbose = True)
print "#" * 10
print "#" * 10
print "#", "prolate qproj"
Qpdf.quan_sample(prolate_qproj_sample, verbose = True)
print "#" * 10
print "#" * 10
print "#", "JS02 fgeo"
Qpdf.quan_sample(fgeo_sample, verbose = True)
print "#" * 10
print "#" * 10
print "#", "prolate fgeo"
Qpdf.quan_sample(prolate_fgeo_sample, verbose = True)
print "#" * 10
In [36]:
axis_ratio_edges = np.linspace(0.0, 1.0, 100)
axis_ratio_bins = 0.5 * (axis_ratio_edges[1:] + axis_ratio_edges[:-1])
axis_ratio_steps = (axis_ratio_edges[1:] - axis_ratio_edges[:-1])
hist_qa = np.histogram(a = qa_sample, bins = axis_ratio_edges)[0]
hist_qb = np.histogram(a = qb_sample, bins = axis_ratio_edges)[0]
hist_2d = np.histogram(a = qproj_sample, bins = axis_ratio_edges)[0]
pyplt.plot(axis_ratio_bins, hist_qa * 1.0 / nsample, "g-", label = "a/c")
pyplt.plot(axis_ratio_bins, hist_qb * 1.0 / nsample, "r-", label = "b/c")
pyplt.plot(axis_ratio_bins, hist_2d * 1.0 / nsample, "k-", label = "axis ratio in 2d")
pyplt.xlabel("axis ratios")
pyplt.ylabel("Prob.")
pyplt.legend(loc=2, numpoints = 1)
fgeo_edges = np.linspace(0.0, 5.0, 500)
fgeo_bins = 0.5 * (fgeo_edges[1:] + fgeo_edges[:-1])
fgeo_steps = (fgeo_edges[1:] - fgeo_edges[:-1])
hist_fgeo = np.histogram(fgeo_sample, bins = fgeo_edges)[0]
hist_fgeo_prolate= np.histogram(prolate_fgeo_sample, bins = fgeo_edges)[0]
pyplt.figure("fgeo")
pyplt.plot(fgeo_bins, hist_fgeo * 1.0 / nsample, "k-", label = "JS02 halos")
pyplt.plot(fgeo_bins, hist_fgeo_prolate * 1.0 / nsample, "r-", label = "prolate")
pyplt.legend(loc = 1, numpoints = 1)
pyplt.xlabel("fgeo")
pyplt.ylabel("Prob.")
Out[36]:
In [15]:
In [15]:
In [ ]: