cosmotool/python_sample/build_nbody_ksz_from_galaxies.py

168 lines
4.3 KiB
Python
Raw Normal View History

2014-06-18 16:13:16 +02:00
import healpy as hp
import numpy as np
import cosmotool as ct
import argparse
import h5py as h5
2014-07-03 15:37:37 +02:00
import matplotlib
matplotlib.use('Agg')
2014-06-18 16:13:16 +02:00
from matplotlib import pyplot as plt
import ksz
from ksz.constants import *
from cosmotool import interp3d
def wrapper_impulsion(f):
class _Wrapper(object):
def __init__(self):
pass
def __getitem__(self,direction):
if 'velocity' in f:
return f['velocity'][:,:,:,direction]
n = "p%d" % direction
return f[n]
return _Wrapper()
parser=argparse.ArgumentParser(description="Generate Skymaps from CIC maps")
parser.add_argument('--boxsize', type=float, required=True)
parser.add_argument('--Nside', type=int, default=128)
parser.add_argument('--base_h5', type=str, required=True)
parser.add_argument('--base_fig', type=str, required=True)
parser.add_argument('--start', type=int, required=True)
parser.add_argument('--end', type=int, required=True)
parser.add_argument('--step', type=int, required=True)
parser.add_argument('--minval', type=float, default=-0.5)
parser.add_argument('--maxval', type=float, default=0.5)
parser.add_argument('--depth_min', type=float, default=10)
parser.add_argument('--depth_max', type=float, default=60)
parser.add_argument('--iid', type=int, default=0)
parser.add_argument('--ksz_map', type=str, required=True)
args = parser.parse_args()
L = args.boxsize
Nside = args.Nside
def build_unit_vectors(N):
ii = np.arange(N)*L/N - 0.5*L
d = np.sqrt(ii[:,None,None]**2 + ii[None,:,None]**2 + ii[None,None,:]**2)
d[N/2,N/2,N/2] = 1
ux = ii[:,None,None] / d
uy = ii[None,:,None] / d
uz = ii[None,None,:] / d
return ux,uy,uz
def build_radial_v(v):
2014-07-03 15:37:37 +02:00
N = v[0].shape[0]
2014-06-18 16:13:16 +02:00
u = build_unit_vectors(N)
2014-07-03 15:37:37 +02:00
vr = v[0] * u[2]
2014-06-18 16:13:16 +02:00
vr += v[1] * u[1]
2014-07-03 15:37:37 +02:00
vr += v[2] * u[0]
2014-06-18 16:13:16 +02:00
2014-07-03 15:37:37 +02:00
return vr.transpose()
def generate_from_catalog(vfield,Boxsize,dmin,dmax):
2014-06-18 16:13:16 +02:00
import progressbar as pbar
cat = np.load("2m++.npy")
cat['distance'] = cat['best_velcmb']
cat = cat[np.where((cat['distance']>100*dmin)*(cat['distance']<dmax*100))]
deg2rad = np.pi/180
2014-07-03 15:37:37 +02:00
Npix = 12*Nside**2
2014-06-18 16:13:16 +02:00
xp,yp,zp = hp.pix2vec(Nside, np.arange(Npix))
N2 = np.sqrt(xp**2+yp**2+zp**2)
2014-07-03 15:37:37 +02:00
ksz_template = np.zeros(Npix, dtype=np.float64)
ksz_mask = np.zeros(Npix, dtype=np.uint8)
2014-06-18 16:13:16 +02:00
2014-07-03 15:37:37 +02:00
pb = pbar.ProgressBar(maxval = cat.size, widgets=[pbar.Bar(), pbar.ETA()]).start()
2014-06-18 16:13:16 +02:00
for k,i in np.ndenumerate(cat):
2014-07-03 15:37:37 +02:00
pb.update(k[0])
2014-06-18 16:13:16 +02:00
l,b=i['gal_long'],i['gal_lat']
2014-07-03 15:37:37 +02:00
ra,dec=i['ra'],i['dec']
2014-06-18 16:13:16 +02:00
l *= deg2rad
b *= deg2rad
2014-07-03 15:37:37 +02:00
ra *= deg2rad
dec *= deg2rad
# x0 = np.cos(l)*np.cos(b)
# y0 = np.sin(l)*np.cos(b)
# z0 = np.sin(b)
x0 = xra = np.cos(ra)*np.cos(dec)
y0 = yra = np.sin(ra)*np.cos(dec)
z0 = zra = np.sin(dec)
2014-06-18 16:13:16 +02:00
DA =i['distance']/100
Lgal = DA**2*10**(0.4*(tmpp_cat['Msun']-i['K2MRS']+25))
2014-07-03 15:37:37 +02:00
profiler = ksz.KSZ_Isothermal(Lgal, 2.37)
2014-06-18 16:13:16 +02:00
2014-07-03 15:37:37 +02:00
idx0 = hp.query_disc(Nside, (x0,y0,z0), 3*profiler.rGalaxy/DA)
vr = interp3d(DA * xra, DA * yra, DA * zra, vfield, Boxsize)
2014-06-18 16:13:16 +02:00
xp1 = xp[idx0]
yp1 = yp[idx0]
zp1 = zp[idx0]
N2_1 = N2[idx0]
cos_theta = x0*xp1+y0*yp1+z0*zp1
cos_theta /= np.sqrt(x0**2+y0**2+z0**2)*(N2_1)
idx,idx_masked,m = profiler.projected_profile(cos_theta, DA)
idx = idx0[idx]
idx_masked = idx0[idx_masked]
ksz_template[idx] += vr * m
ksz_mask[idx_masked] = 0
pb.finish()
return ksz_template, ksz_mask
for i in xrange(args.start,args.end,args.step):
ff=plt.figure(1)
plt.clf()
v=[]
2014-07-03 15:37:37 +02:00
fname = args.base_h5 % i
if False:
print("Opening %s..." % fname)
with h5.File(fname, mode="r") as f:
p = wrapper_impulsion(f)
for j in xrange(3):
v.append(p[j] / f['density'][:])
print("Building radial velocities...")
2014-06-18 16:13:16 +02:00
# Rescale by Omega_b / Omega_m
2014-07-03 15:37:37 +02:00
vr = build_radial_v(v)
# _,_,vr = build_unit_vectors(128)
# vr *= 1000 * 500
vr *= ksz_normalization*1e6
del v
# np.save("vr.npy", vr)
else:
print("Loading vrs...")
vr = np.load("vr.npy")
print("Generating map...")
2014-06-18 16:13:16 +02:00
2014-07-03 15:37:37 +02:00
proj,mask = generate_from_catalog(vr,args.boxsize,args.depth_min,args.depth_max)
2014-06-18 16:13:16 +02:00
hp.write_map(args.ksz_map % i, proj)
2014-07-03 15:37:37 +02:00
hp.write_map((args.ksz_map % i) + "_mask", mask)
2014-06-18 16:13:16 +02:00
hp.mollview(proj, fig=1, coord='CG', cmap=plt.cm.coolwarm, title='Sample %d' % i, min=args.minval,
max=args.maxval)
ff.savefig(args.base_fig % i)