KSZ work

python_sample/build_dipole_ksz_from_galaxies.py

@@ -100,13 +100,35 @@ def build_unit_vectors(N):
 
   return ux,uy,uz
 
-def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_hubble=False):
+def compute_vcmb(l, b):
+    # Motion is obtained from Tully (2007): sun_vs_LS + LS_vs_CMB
+    motion = [-25.,-246.,277.];
+
+    x = np.cos(l*np.pi/180) * np.cos(b*np.pi/180)
+    y = np.sin(l*np.pi/180) * np.cos(b*np.pi/180)
+    z = np.sin(b*np.pi/180)
+
+    return x*motion[0] + y*motion[1] + z*motion[2]
+
+
+def compute_vlg(l,b):
+
+    motion = [-79,296,-36]; # [-86, 305, -33];
+
+    x = np.cos(l*np.pi/180) * np.cos(b*np.pi/180)
+    y = np.sin(l*np.pi/180) * np.cos(b*np.pi/180)
+    z = np.sin(b*np.pi/180)
+
+    return x*motion[0] + y*motion[1] + z*motion[2]
+
+
+def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_hubble=False,x=2.37,bright=-np.inf,bright_list=[],use_vlg=True,sculpt=-1):
   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))]
+#  cat = cat[np.where((cat['distance']>100*dmin)*(cat['distance']<dmax*100))]
 
   deg2rad = np.pi/180
   Npix = 12*Nside**2
@@ -119,14 +141,19 @@ def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_h
     ksz_hubble_template = np.zeros(ksz_template.size, dtype=np.float64)
 
   for i in pbar.ProgressBar(maxval = cat.size, widgets=[pbar.Bar(), pbar.ETA()])(cat):
+    # Skip too point sources
+    if i['name'] in bright_list:
+      print("Object %s is in bright list" % i['name'])
+      continue
+
     if do_random:
       l = np.random.rand()*360
       b = np.arcsin(2*np.random.rand()-1)*180/np.pi
     else:
-      l,b=i['gal_long'],i['gal_lat']
+      l0,b0=i['gal_long'],i['gal_lat']
 
-    l=ne.evaluate('l*deg2rad')
+    l=ne.evaluate('l0*deg2rad')
-    b=ne.evaluate('b*deg2rad')
+    b=ne.evaluate('b0*deg2rad')
 
     dtheta,dphi = np.random.randn(2)*perturb
     theta,l=move_direction(dtheta,dphi,0.5*np.pi - b, l)
@@ -137,10 +164,23 @@ def generate_from_catalog(dmin,dmax,Nside,perturb=0.0,y=0.0,do_random=False,do_h
     y0 = np.sin(l)*np.cos(b)
     z0 = np.sin(b)
 
-    DA =i['distance']/100
+    if use_vlg:
+      vlg = i['best_velcmb'] - compute_vcmb(l0, b0) + compute_vlg(l0, b0)
+      DA = vlg/100
+    else:
+      DA = i['best_velcmb'] / 100
+
+    if DA < dmin or DA > dmax:
+      continue
+
     Lgal = DA**2*10**(0.4*(tmpp_cat['Msun']-i['K2MRS']+25))
 
-    profiler = ksz.KSZ_Isothermal(Lgal, 2.37, y=y)
+    M_K=i['K2MRS']-5*np.log10(DA)-25
+    # Skip too bright galaxies
+    if M_K < bright:
+      continue
+
+    profiler = ksz.KSZ_Isothermal(Lgal, x, y=y, sculpt=sculpt)
 
     idx0 = hp.query_disc(Nside, (x0,y0,z0), 3*profiler.rGalaxy/DA)
 
@@ -177,12 +217,17 @@ def get_args():
   parser.add_argument('--depth_max', type=float, default=60)
   parser.add_argument('--ksz_map', type=str, required=True)
   parser.add_argument('--base_fig', type=str, default="kszfig.png")
-  parser.add_argument('--build_dipole', type=bool, default=False)
+  parser.add_argument('--build_dipole', action='store_true')
   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('--x',type=float,default=2.37)
+  parser.add_argument('--random', action='store_true')
   parser.add_argument('--perturb', type=float, default=0)
-  parser.add_argument('--hubble_monopole', type=bool, default=False)
+  parser.add_argument('--hubble_monopole', action='store_true')
+  parser.add_argument('--remove_bright', type=float, default=-np.inf)
+  parser.add_argument('--bright_file', type=str)
+  parser.add_argument('--lg', action='store_true')
+  parser.add_argument('--sculpt_beam', type=float, default=-1)
   return parser.parse_args()
 
 def main():
@@ -195,7 +240,16 @@ def main():
 
   print("Generating map...")
 
-  r = generate_from_catalog(args.depth_min,args.depth_max,args.Nside,perturb=args.perturb,y=args.y,do_random=args.random,do_hubble=args.hubble_monopole)
+  with open("crap.txt", mode="r") as f:
+    bright_list = [l.split('#')[0].strip(" \t\n\r") for l in f]
+
+  if args.bright_file:
+    with open(args.bright_file, mode="r") as f:
+      idx_name = f.readline().split(',').index('name_2')
+      bright_list = bright_list + [l.split(',')[idx_name] for l in f]
+
+  print("Built bright point source list: " + repr(bright_list))
+  r = generate_from_catalog(args.depth_min,args.depth_max,args.Nside,perturb=args.perturb,y=args.y,do_random=args.random,do_hubble=args.hubble_monopole,x=args.x,bright=args.remove_bright,use_vlg=args.lg,bright_list=bright_list,sculpt=args.sculpt_beam)
   hubble_map = None
   if args.hubble_monopole:
     proj,mask,hubble_map = r

python_sample/icgen/borgicgen.py

@@ -4,10 +4,11 @@ import cosmolopy as cpy
 from cosmogrowth import *
 import borgadaptor as ba
 
-def gen_posgrid(N, L):
+@ct.timeit
+def gen_posgrid(N, L, delta=1, dtype=np.float32):
     """ Generate an ordered lagrangian grid"""
 
-    ix = (np.arange(N)*L/N).astype(np.float32)
+    ix = (np.arange(N)*(L/N*delta)).astype(dtype)
 
     x = ix[:,None,None].repeat(N, axis=1).repeat(N, axis=2)
     y = ix[None,:,None].repeat(N, axis=0).repeat(N, axis=2)
@@ -147,6 +148,8 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, supergenerate
     
     D2 = -3./7 * D1_0**2
 
+    if do_lpt2:
+      psi2 = lpt.lpt2('all')
     for j in xrange(3):
         # Generate psi_j (displacement along j)
         print("LPT1 axis=%d" % j)
@@ -154,8 +157,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, supergenerate
         psi = psi.reshape((psi.size,))
         if do_lpt2:
           print("LPT2 axis=%d" % j)
-          psi2 = lpt.lpt2(j)
+          psi += D2 * psi2[j].reshape((psi2[j].size,))
-          psi += D2 * psi2.reshape((psi2.size,))
         # Generate posx
         posx.append(((posq[j] + psi)%L).astype(np.float32))
         # Generate vel

python_sample/icgen/cosmogrowth.py

@@ -4,6 +4,7 @@ import pyfftw
 import weakref
 import numpy as np
 import cosmolopy as cpy
+import cosmotool as ct
 
 class CubeFT(object):
   def __init__(self, L, N, max_cpu=-1):
@@ -98,6 +99,7 @@ class LagrangianPerturbation(object):
         self._kz = self.ik[None,None,:(self.N/2+1)]
         self.cache = {}#weakref.WeakValueDictionary()
 
+    @ct.timeit_quiet
     def upgrade_sampling(self, supersample):
         N2 = self.N * supersample
         N = self.N
@@ -113,14 +115,26 @@ class LagrangianPerturbation(object):
         self.N = N2
         self.cube = CubeFT(self.L, self.N, max_cpu=self.max_cpu)
         
+    @ct.timeit_quiet
     def _gradient(self, phi, direction):
-        ne.evaluate('phi_hat * i * kv / (kx**2 + ky**2 + kz**2)', out=self.cube.dhat,
+        if direction == 'all':
-              local_dict={'i':-1j, 'phi_hat':phi, 'kv':self._kdir(direction),
+          dirs = [0,1,2]
+          copy = True
+        else:
+          dirs = [direction]
+          copy = False
+        ret=[]
+        for dir in dirs:
+          ne.evaluate('phi_hat * i * kv / (kx**2 + ky**2 + kz**2)', out=self.cube.dhat,
+              local_dict={'i':-1j, 'phi_hat':phi, 'kv':self._kdir(dir),
                           'kx':self._kx, 'ky':self._ky, 'kz':self._kz},casting='unsafe')
 #        self.cube.dhat = self._kdir(direction)*1j*phi 
-        self.cube.dhat[0,0,0] = 0
+          self.cube.dhat[0,0,0] = 0
-        return self.cube.irfft()
+          x = self.cube.irfft()
+          ret.append(x.copy() if copy else x)
+        return ret[0] if len(ret)==1 else ret
             
+    @ct.timeit_quiet
     def lpt1(self, direction=0):
         return self._gradient(self.dhat, direction)
         
@@ -155,6 +169,7 @@ class LagrangianPerturbation(object):
           self.cube.density = array
         return self.cube.rfft()
 
+    @ct.timeit_quiet
     def lpt2(self, direction=0):
 #        k2 = self._get_k2()
 #        k2[0,0,0] = 1
@@ -170,12 +185,13 @@ class LagrangianPerturbation(object):
             for j in xrange(3):
                 q = self._do_irfft( potgen0(j) ).copy()
                 for i in xrange(j+1, 3):
-                    ne.evaluate('div + q * pot', out=div_phi2, 
+                    with ct.time_block("LPT2 elemental (%d,%d)" %(i,j)):
-                          local_dict={'div':div_phi2, 'q':q,'pot':self._do_irfft( potgen0(i), copy=False ) } 
+                      ne.evaluate('div + q * pot', out=div_phi2, 
-                          )
+                            local_dict={'div':div_phi2, 'q':q,'pot':self._do_irfft( potgen0(i), copy=False ) } 
-                    ne.evaluate('div - pot**2',out=div_phi2,
+                            )
-                          local_dict={'div':div_phi2,'pot':self._do_irfft(potgen(i,j), copy=False) }
+                      ne.evaluate('div - pot**2',out=div_phi2,
-                          )
+                            local_dict={'div':div_phi2,'pot':self._do_irfft(potgen(i,j), copy=False) }
+                            )
 
             phi2_hat = self._do_rfft(div_phi2)
             #self.cache['lpt2_potential'] = phi2_hat

python_sample/ksz/constants.py

@@ -27,6 +27,7 @@ tmpp_cat={'Msun':3.29,
 baryon_fraction = Omega_baryon / Omega_matter
 
 ksz_normalization = -T_cmb*sigmaT*v_unit/(lightspeed*mu*mp) * baryon_fraction
+assert ksz_normalization < 0
 rho_mean_matter = Omega_matter * (3*(100e3/Mpc)**2/(8*np.pi*G)) 
 Lbar = tmpp_cat['lbar'] / Mpc**3
 M_over_L_galaxy = rho_mean_matter / Lbar

python_sample/ksz/gal_prof.py

@@ -10,7 +10,14 @@ class KSZ_Profile(object):
   R_star= 0.0 # 15 kpc/h 
   L_gal0 = 10**(0.4*(tmpp_cat['Msun']-tmpp_cat['Mstar']))
 
-  def __init__(self):
+  def __init__(self,sculpt):
+    """Base class for KSZ profiles
+
+    Arguments:
+      sculpt (float): If negative, do not sculpt. If positive, there will be a 2d 
+                      suppression of the profile with a radius given by sculpt (in arcmins).
+    """
+    self.sculpt = sculpt * np.pi/180/60.
     self.rGalaxy = 1.0
 
   def evaluate_profile(self, r):
@@ -37,6 +44,11 @@ class KSZ_Profile(object):
     if tan_theta_2.size > 0:
       idx_mask = np.append(idx_mask,idx[tan_theta_2.argmin()])
 
+    if self.sculpt > 0:
+      theta = np.arctan(tan_theta)
+      cond = theta < self.sculpt
+      m[cond] *= (theta[cond]/self.sculpt)**2
+
     return idx,idx_mask,m
      
 
@@ -48,10 +60,20 @@ class KSZ_Isothermal(KSZ_Profile):
   sigma_FP=160e3  #m/s
   R_innergal = 0.030
 
-  def __init__(self, Lgal, x, y=0.0):
+  def __init__(self, Lgal, x, y=0.0, sculpt=-1):
-    "Support for Isothermal profile"
+    """Support for Isothermal profile
+
+    Arguments:
+      Lgal (float): Galaxy luminosity in solar units
+      x (float): extent of halo in virial radius units
+
+    Keyword arguments:
+      y (float): Inner part where there is no halo
+      sculpt (float): If negative, do not sculpt. If positive, there will be a 2d 
+                      suppression of the profile with a radius given by sculpt (in arcmins).
+    """
     
-    super(KSZ_Isothermal,self).__init__()
+    super(KSZ_Isothermal,self).__init__(sculpt)
     
     self.R_gal = 0.226 * x
     self.R_innergal *= y
@@ -71,20 +93,20 @@ class KSZ_Isothermal(KSZ_Profile):
     
     Q = np.zeros(r.size)
 
-    cond = (r<=0)
+    cond = (r<=1e-10)
     Q[cond] = rho0*2/Mpc * (rGalaxy-rInner)/(rGalaxy*rInner)
 
     cond = (r>0)*(r <= rInner)
     D['r'] = r[cond]
     Q[cond] = ne.evaluate('rho0*2/(Mpc*r) * (arctan(sqrt( (rGalaxy/r)**2 -1 ))  - arctan(sqrt( (rInner/r)**2 - 1 )))',
                 local_dict=D)
-    
+ 
     cond = (r > rInner)*(r <= rGalaxy)
     D['r'] = r[cond]
     Q[cond] = ne.evaluate('rho0*2/(Mpc*r) * arctan(sqrt( (rGalaxy/r)**2 -1 ))',
                 local_dict=D)
-    
+
-    return Q,np.where(r<rInner)[0]
+    return Q,[] #np.where(r<rInner)[0]
 
 
 # -----------------------------------------------------------------------------