import numexpr as ne import healpy as hp import numpy as np import cosmotool as ct import argparse import h5py as h5 import matplotlib matplotlib.use('Agg') from matplotlib import pyplot as plt import ksz from ksz.constants import * from cosmotool import interp3d import numpy as np from scipy.special import sinc def move_direction_new(d_theta, d_phi, theta, phi): cos=np.cos sin=np.sin sqrt=np.sqrt amplitude = sqrt(d_theta*d_theta + d_phi*d_phi); cos_alpha = d_theta/amplitude; sin_alpha = d_phi/amplitude; if (amplitude == 0): return theta,phi cos_d = cos(amplitude); sin_d = sin(amplitude); cos_theta = cos(theta); sin_theta = sin(theta); cos_phi = cos(phi); sin_phi = sin(phi); basis = [ [ cos_phi * sin_theta, sin_phi * sin_theta, cos_theta ], [ cos_phi * cos_theta, sin_phi * cos_theta, -sin_theta ], [ -sin_phi, cos_phi, 0 ] ] np0 = [ cos_d, cos_alpha*sin_d, sin_alpha*sin_d ] np1 = [ sum([basis[j][i] * np0[j] for j in xrange(3)]) for i in xrange(3) ] dnp = sqrt(sum([np1[i]**2 for i in xrange(3)])) theta = np.arccos(np1[2]/dnp); phi = np.arctan2(np1[1], np1[0]) % (2*np.pi); return theta,phi def move_direction_new2(delta_theta, delta_phi, theta, phi): cos,sin,sqrt=np.cos,np.sin,np.sqrt grad_len = sqrt(delta_theta**2 + delta_phi**2) if grad_len==0: return theta,phi cth0 = cos(theta) sth0 = sin(theta) topbottom = 1 if (theta < 0.5*np.pi) else -1 sinc_grad_len = sinc(grad_len) cth = topbottom*cos(grad_len) * cth0 - sinc_grad_len*sth0*delta_theta sth = max(1e-10, sqrt((1.0-cth)*(1.0+cth)) ) phi = phi + np.arcsin(delta_phi * sinc_grad_len / sth) theta = np.arccos(cth) return theta,phi move_direction = move_direction_new 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() def build_unit_vectors(N): ii = np.arange(N,dtype=np.float64)/N - 0.5 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 generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False): import progressbar as pbar cat = np.load("2m++.npy") cat['distance'] = cat['best_velcmb'] cat = cat[np.where((cat['distance']>100*dmin)*(cat['distance'] 0: proj *= mask proj = hp.ud_grade(proj, nside_out=args.degrade) mask = hp.ud_grade(mask, nside_out=args.degrade) Nside = args.degrade else: Nside = args.Nside hp.write_map(args.ksz_map + ".fits", proj) hp.write_map(args.ksz_map + "_mask.fits", mask) if args.build_dipole: x,y,z=hp.pix2vec(Nside, np.arange(hp.nside2npix(Nside))) hp.write_map(args.ksz_map + "_x.fits", proj*x) hp.write_map(args.ksz_map + "_y.fits", proj*y) hp.write_map(args.ksz_map + "_z.fits", proj*z) hp.mollview(proj*100*1e6, fig=1, coord='GG', cmap=plt.cm.coolwarm, title='', min=args.minval, max=args.maxval) ff.savefig(args.base_fig) main()