Added support for perturbed position

python_sample/build_dipole_ksz_from_galaxies.py

@@ -10,6 +10,67 @@ 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):
 
@@ -38,7 +99,7 @@ def build_unit_vectors(N):
 
   return ux,uy,uz
 
-def generate_from_catalog(dmin,dmax,Nside,y=0.0,do_random=False):
+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")
@@ -57,13 +118,18 @@ def generate_from_catalog(dmin,dmax,Nside,y=0.0,do_random=False):
   for i in pbar.ProgressBar(maxval = cat.size, widgets=[pbar.Bar(), pbar.ETA()])(cat):
     if do_random:
       l = np.random.rand()*360
-      b = np.arcsin(2*np.random.rand()-1)*90/np.pi
+      b = np.arcsin(2*np.random.rand()-1)*180/np.pi
     else:
       l,b=i['gal_long'],i['gal_lat']
 
     l=ne.evaluate('l*deg2rad')
     b=ne.evaluate('b*deg2rad')
 
+    dtheta,dphi = np.random.randn(2)*perturb
+    theta,l=move_direction(dtheta,dphi,0.5*np.pi - b, l)
+	
+    b = 0.5*np.pi-theta
+
     x0 = np.cos(l)*np.cos(b)
     y0 = np.sin(l)*np.cos(b)
     z0 = np.sin(b)
@@ -105,6 +171,7 @@ def get_args():
   parser.add_argument('--degrade', type=int, default=-1)
   parser.add_argument('--y',type=float,default=0.0)
   parser.add_argument('--random', type=bool, default=False)
+  parser.add_argument('--perturb', type=float, default=0)
   return parser.parse_args()
 
 def main():
@@ -117,7 +184,7 @@ def main():
 
   print("Generating map...")
 
-  proj,mask = generate_from_catalog(args.depth_min,args.depth_max,args.Nside,args.y,do_random=args.random)
+  proj,mask = generate_from_catalog(args.depth_min,args.depth_max,args.Nside,perturb=args.perturb,y=args.y,do_random=args.random)
 
   if args.degrade > 0:
     proj *= mask