-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathprofile_1D.py
308 lines (251 loc) · 10.7 KB
/
profile_1D.py
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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
#-------------------------------------------------------------------------------
#Author: Karen Yin-Yee Ng <[email protected]>
#Purpose for plotting 1D reduced shear profile of a NFW halo
#Date: 07/04/2013
#License: BSD
#-------------------------------------------------------------------------------
from __future__ import division
import cosmo
import profiles
import numpy as np
def azimuthal_avg_shear_in_bins(k,x,startbin=0,endbin=10,b_width=1):
'''
Purpose:
average the value of shear within a certain radius
Input:
k = shear that you want to bin
startbin = the value at the starting bin
endbin = the value at the ending bin
b_width = the width of the bin
'''
bins = int((endbin-startbin)/b_width)
#print '# of bins is ',bins
#initialize bin edge
bin_edge=np.arange(startbin, endbin+b_width,b_width)
#print 'bin_edge is ',bin_edge
#initialize the average value of each bin
bin_array = np.arange(startbin+b_width/2.0 , endbin+b_width/2, b_width)
#print 'bin array is', bin_array
#exclude values outside the desired range
mask = np.logical_and(x > startbin, x<= endbin)
red_x = x[mask]
red_k = k[mask]
azim_kappa = np.ones(bin_array.size)
for i in range(bins):
#only look at values within each bin within each iteration
temp_mask = np.logical_and(red_x>= bin_edge[i], red_x<bin_edge[i+1])
if red_k[temp_mask].size==0:
temp_k = 0
else:
temp_k = red_k[temp_mask]
#average the e1_pix value within this bin
if np.isnan(np.mean(temp_k)):
raise ValueError('NaN encountered')
azim_kappa[i] = np.mean(temp_k)
#print '# of data point in this bin = ',len(temp_k), ',bin = ',(bin_edge[i]+bin_edge[i+1] )/2., ' bin_array=',bin_array[i]
return azim_kappa, bin_array
def ext_1D_reduced_shear(theta, conc, r_s):
'''
function for calculating the tangential shear given by a NFW profile
for fitting the concentration parameter and the scale radius r_s
This is written by adopting the expression in Umetsu
theta = range of radius of the halo that we 're considering,
unit in arcsec
conc = concentration parameter at the given radius
r_s = scale radius of the NFW profile (Mpc)
z_halo = halo redshift
### WARNING
### cosmological parameters are written within the function
h_scale = hubble scale H = h*100 km/s / Mpc
Om = matter energy density
Ol = dark energy density
Or radiation energy density
halo coord 1
halo coord 2
z_halo = redshift of the lens
z_source = redshift of the source galaxies
beta = ratio of D_LS and D_S see James Jee 's paper for exact def
'''
##latest parameters for El Gordo
z_halo = 0.87
z_source = 1.318
beta = 0.258
#old parameters
#z_halo =0.89
#z_source = 1.22
#beta = 0.216
#Cosmological parameters
#print "profile_1D.tan_1D_reduced_shear: this func uses its own "
#print "cosmological parameters"
Om = 0.3
Ol = 0.7
Or = 0.0
c = 3e5 #speed of light units km/s
G = 6.673*10**(-11) # m^3/kg/s^2
h_scale=0.7
kminMpc = 3.08568025*10**19 # km in a Megaparsec
minMpc = 3.08568025*10**22 #m in a Megaparsec
kginMsun = 1.988e30 #kg
#angular diameter distance to the lens
dl = cosmo.Da(z_halo,h_scale,Om,Ol) #in Mpc
del_c = 200/3.*conc**3/(np.log(1+conc)-conc/(1+conc)) #unitless
nfwSigmacr = c**2/(4*np.pi*G*dl*beta)*1000/kminMpc #in units of kg/m^2
#print 'nfwSigmacr in units of M_sun /pc^2 is ', nfwSigmacr/kginMsun*(3.086*10**16)**2
# in units of kg/m^3
rho_cr = cosmo.rhoCrit(z_halo,h_scale,Om,Ol,Or)#/kginMsun*minMpc**3
#if scale radius is not given then infer from Duffy et al.
#assuming full halo
#if r_s == np.nan:
#r_s = scale_radius(conc,z_halo,Om,Ol,Or)
kappa_s = 2*del_c*rho_cr*(r_s*minMpc)/nfwSigmacr
#just do a simple conversion from arcsec to physical distance using
#angular diameter distance
#make sure that the units is in arcsec
theta_s = r_s /dl*180.0/np.pi*60.*60.
x = theta / theta_s
#write an array with finer bins than x
fx = x #np.arange(np.min(x),np.max(x),(x[1]-x[0])/10.)
#code the expression for kappa as a function of x
#print kappa_s
func_kappa0 = lambda fx: kappa_s/(1-fx**2.)*\
(-1.+2./np.sqrt(1.-fx**2.)*\
np.arctanh(np.sqrt(1.-fx)/np.sqrt(1.+fx)))
func_kappa1 = lambda fx: kappa_s*1./3.
func_kappa2 = lambda fx: kappa_s/(fx**2.-1)*\
(1.-2./np.sqrt(fx**2.-1.)*\
np.arctan(np.sqrt(fx-1.)/np.sqrt(fx+1.)))
kappa = np.piecewise(fx,[fx<1.0, fx==1.0, fx>1.0],
[func_kappa0,func_kappa1,func_kappa2])
g_theta = np.piecewise(x, [x<1.0, x==1.0, x>1.0],
[lambda x: np.log(x/2.)+2./np.sqrt(1-x**2)*\
np.arctanh(np.sqrt(1-x)/np.sqrt(1+x)),
lambda x: np.log(x/2.)+1.,
lambda x: np.log(x/2.)+2./np.sqrt(x**2.-1.)*\
np.arctan(np.sqrt(x-1.)/np.sqrt(x+1.))])
kappa_bar = 2.*kappa_s/x**2.*g_theta
azim_kappa = kappa #no need to azimuthally smooth function since
#the function did n't contain azimuthal angular info to begin with
#, azim_bins = azimuthal_avg_shear_in_bins(kappa, fx, np.min(x), np.max(x)+x[1]-x[0], x[1]-x[0])
#print 'bins are ', x
#print 'azim_bins are ', azim_bins
#print 'size of kappa_bar is ', kappa_bar.size
#print 'size of azim_kappa is ', azim_kappa.size
#print 'kappa_bar is ', kappa_bar
#print 'azim_kappa is ', azim_kappa
red_shear= (kappa_bar-azim_kappa)/(1-azim_kappa)
return red_shear
def tan_1D_reduced_shear(theta, conc):
'''
function for calculating the tangential shear given by a NFW profile
for fitting the concentration parameter and the scale radius r_s
This is written by adopting the expression in Umetsu
theta = range of radius of the halo that we 're considering, unit in arcsec
conc = concentration parameter at the given radius
z_halo = halo redshift
### WARNING
### cosmological parameters are written within the function
h_scale = hubble scale H = h*100 km/s / Mpc
Om = matter energy density
Ol = dark energy density
Or radiation energy density
halo coord 1
halo coord 2
z_halo = redshift of the lens
z_source = redshift of the source galaxies
beta = ratio of D_LS and D_S see James Jee 's paper for exact def
'''
##latest parameters for El Gordo
z_halo = 0.87
z_source = 1.318
beta = 0.258
#old parameters
#z_halo =0.89
#z_source = 1.22
#beta = 0.216
#Cosmological parameters
#print "profile_1D.tan_1D_reduced_shear: this func uses its own "
#print "cosmological parameters"
Om = 0.3
Ol = 0.7
Or = 0.0
c = 3e5 #speed of light units km/s
G = 6.673*10**(-11) # m^3/kg/s^2
h_scale=0.7
kminMpc = 3.08568025*10**19 # km in a Megaparsec
minMpc = 3.08568025*10**22 #m in a Megaparsec
kginMsun = 1.988e30 #kg
#angular diameter distance to the lens
dl = cosmo.Da(z_halo,h_scale,Om,Ol) #in Mpc
del_c = 200/3.*conc**3/(np.log(1+conc)-conc/(1+conc)) #unitless
nfwSigmacr = c**2/(4*np.pi*G*dl*beta)*1000/kminMpc #in units of kg/m^2
#print 'nfwSigmacr in units of M_sun /pc^2 is ', nfwSigmacr/kginMsun*(3.086*10**16)**2
# in units of kg/m^3
rho_cr = cosmo.rhoCrit(z_halo,h_scale,Om,Ol,Or)#/kginMsun*minMpc**3
#if scale radius is not given then infer from Duffy et al.
#assuming full halo
#if r_s == np.nan:
r_s = scale_radius(conc,z_halo,Om,Ol,Or)
kappa_s = 2*del_c*rho_cr*(r_s*minMpc)/nfwSigmacr
#just do a simple conversion from arcsec to physical distance using
#angular diameter distance
#make sure that the units is in arcsec
theta_s = r_s /dl*180.0/np.pi*60.*60.
x = theta / theta_s
#write an array with finer bins than x
fx = x #np.arange(np.min(x),np.max(x),(x[1]-x[0])/10.)
#code the expression for kappa as a function of x
#print kappa_s
func_kappa0 = lambda fx: kappa_s/(1-fx**2.)*\
(-1.+2./np.sqrt(1.-fx**2.)*\
np.arctanh(np.sqrt(1.-fx)/np.sqrt(1.+fx)))
func_kappa1 = lambda fx: kappa_s*1./3.
func_kappa2 = lambda fx: kappa_s/(fx**2.-1)*\
(1.-2./np.sqrt(fx**2.-1.)*\
np.arctan(np.sqrt(fx-1.)/np.sqrt(fx+1.)))
kappa = np.piecewise(fx,[fx<1.0, fx==1.0, fx>1.0],
[func_kappa0,func_kappa1,func_kappa2])
g_theta = np.piecewise(x, [x<1.0, x==1.0, x>1.0],
[lambda x: np.log(x/2.)+2./np.sqrt(1-x**2)*\
np.arctanh(np.sqrt(1-x)/np.sqrt(1+x)),
lambda x: np.log(x/2.)+1.,
lambda x: np.log(x/2.)+2./np.sqrt(x**2.-1.)*\
np.arctan(np.sqrt(x-1.)/np.sqrt(x+1.))])
kappa_bar = 2.*kappa_s/x**2.*g_theta
azim_kappa = kappa #no need to azimuthally smooth function since
#the function did n't contain azimuthal angular info to begin with
#, azim_bins = azimuthal_avg_shear_in_bins(kappa, fx, np.min(x), np.max(x)+x[1]-x[0], x[1]-x[0])
#print 'bins are ', x
#print 'azim_bins are ', azim_bins
#print 'size of kappa_bar is ', kappa_bar.size
#print 'size of azim_kappa is ', azim_kappa.size
#print 'kappa_bar is ', kappa_bar
#print 'azim_kappa is ', azim_kappa
red_shear= (kappa_bar-azim_kappa)/(1-azim_kappa)
return red_shear
def scale_radius(conc, z_halo,Om=0.3, Ol=0.7, Or=0.0,
A_200=5.71, B_200=-0.084, C_200=-0.47, h_scale=0.7):
'''
purpose: compute the scale radius given the concentration parameter
default parameters based on full profiles of Duffy et al. 2008
input:
conc = concentration parameter at r200
z_halo = redshift of that halo
Om, Ol, Or = cosmological paramters
A200 = duffy et al parameter for concentration radius relationship
B200 = "
C200 = "
h_scale = reduced hubble parameter
output:
scale radius in Mpc
'''
#unit conversion values
minMpc = 3.08568025*10**22 #m in a Megaparsec
kginMsun = 1.988e30 #kg
rho_cr = cosmo.rhoCrit(z_halo,h_scale,Om,Ol,Or) #in kg/m**3
#the h_scale is multiplied because the pivotal mass
#m_pivotal = 2e12 is in unit of M_sun h_scale^(-1)
m_200 = profiles.nfwM200(conc,z_halo, A_200,B_200,C_200,h_scale)*\
kginMsun
r_200 = (m_200/(4*np.pi/3*200*rho_cr))**(1/3.) #in m
r_s = r_200 / conc / minMpc
return r_s