fixes to kappa profiles

hmvec/cosmology.py

@@ -17,7 +17,7 @@
 
 class Cosmology(object):
 
-    def __init__(self,params,halofit=None,engine='camb'):
+    def __init__(self,params=None,halofit=None,engine='camb'):
         assert engine in ['camb','class']
         if engine=='class': raise NotImplementedError
 
@@ -29,6 +29,15 @@ def __init__(self,params,halofit=None,engine='camb'):
         self._init_cosmology(self.p,halofit)
 
 
+    def sigma_crit(self,zlens,zsource):
+        Gval = 4.517e-48 # Newton G in Mpc,seconds,Msun units
+        cval = 9.716e-15 # speed of light in Mpc,second units
+        Dd = self.angular_diameter_distance(zlens)
+        Ds = self.angular_diameter_distance(zsource)
+        Dds = np.asarray([self.results.angular_diameter_distance2(zl,zsource) for zl in zlens])
+        return cval**2 * Ds / 4 / np.pi / Gval / Dd / Dds
+
+
     def P_mm_linear(self,zs,ks):
         pass
 
@@ -38,6 +47,9 @@ def P_mm_nonlinear(self,ks,zs,halofit_version='mead'):
     def comoving_radial_distance(self,z):
         return self.results.comoving_radial_distance(z)
 
+    def angular_diameter_distance(self,z):
+        return self.results.angular_diameter_distance(z)
+
     def hubble_parameter(self,z):
         # H(z) in km/s/Mpc
         return self.results.hubble_parameter(z)

hmvec/hmvec.py

@@ -579,51 +579,52 @@ def _get_term(iname):
             print("Two-halo consistency2: " , consistency2,integral2)
         return self.Pzk * (integral+b1-consistency1)*(integral2+b2-consistency2)
 
-    def sigma_1h_profiles(self,thetas,Ms,concs,sig_theta=None,delta=200,rho='critical',rho_at_z=True):
-        import cluster-lensing as cl
+    def sigma_1h_profiles(self,thetas,Ms,concs,sig_theta=None,delta=200,rho='mean',rho_at_z=True):
+        import clusterlensing as cl
         zs = self.zs
-        chis = self.comoving_radial_distance(zs)
+        Ms = np.asarray(Ms)
+        concs = np.asarray(concs)
+        chis = self.angular_diameter_distance(zs)
         rbins = chis * thetas
         offsets = chis * sig_theta if sig_theta is not None else None
         if rho=='critical': rhofunc = self.rho_critical_z 
         elif rho=='mean': rhofunc = self.rho_matter_z
         rhoz = zs if rho_at_z else zs * 0
-        Rdeltas = R_from_M(Ms,rhofunc(rhoz)[:,None],delta=delta)
+        Rdeltas = R_from_M(Ms,rhofunc(rhoz),delta=delta)
         rs = Rdeltas / concs
         rhocrits = self.rho_critical_z(zs)
         delta_c =  Ms / 4 / np.pi / rs**3 / rhocrits / Fcon(concs)
-        smd = cl.nfw.SurfaceMassDensity(rs, delta_c, rho_crit,rbins=rbins,offsets=offsets)
+        smd = cl.nfw.SurfaceMassDensity(rs, delta_c, rhocrits,rbins=rbins,offsets=offsets)
         sigma = smd.sigma_nfw()
         return sigma
 
-    def sigma_crit(self,zsource):
-        zlens = self.zs
-        Gval = 4.517e-48 # Newton G in Mpc,seconds,Msun units
-        cval = 9.716e-15 # speed of light in Mpc,second units
-        Dd = self.comoving_radial_distance(zlens)
-        Ds = self.comoving_radial_distance(zsource)
-        Dds = Ds - Dd
-        return cval**2 * Ds / 4 / np.pi / Gval / Dd / Dds
-
-    def kappa_1h_profiles(self,thetas,Ms,concs,zsource,sig_theta=None,delta=200,rho='critical',rho_at_z=True):
+    def kappa_1h_profiles(self,thetas,Ms,concs,zsource,sig_theta=None,delta=200,rho='mean',rho_at_z=True):
         sigma = self.sigma_1h_profiles(thetas,Ms,concs,sig_theta=sig_theta,delta=delta,rho=rho,rho_at_z=rho_at_z)
-        sigmac = self.sigma_crit(zsource)
+        sigmac = self.sigma_crit(self.zs,zsource)
         return sigma / sigmac
 
-    def kappa_2h_profiles(self,thetas,Ms,concs,zsource,delta=200,rho='critical',rho_at_z=True,lmin=2,lmax=10000):
+    def kappa_2h_profiles(self,thetas,Ms,concs,zsource,delta=200,rho='mean',rho_at_z=True,lmin=100,lmax=10000):
         from scipy.special import j0
         zlens = self.zs
-        sigmac = self.sigma_crit(zsource)
+        sigmac = self.sigma_crit(zlens,zsource)
         rhomz = self.rho_matter_z(zlens)
         chis = self.comoving_radial_distance(zlens)
+        DAz = self.results.angular_diameter_distance(zlens)
         ells = self.ks*chis
         sel = np.logical_and(ells>lmin,ells<lmax)
         ells = ells[sel]
         #Array indexing is as follows:
         #[z,M,k/r]
         Ps = self.Pzk[:,sel]
-        integrand = rhomz * self.bh * Ps / (1+zlens)**3. / sigmac / DAz**2 * j0(ells*thetas) * ells / 2./ np.pi
-        return np.trapz(integrand,ells)
+        bhs = []
+        for i in range(zlens.shape[0]): # vectorize this
+            bhs.append( interp1d(self.ms,self.bh[i])(Ms))
+        bhs = np.asarray(bhs)
+        ints = []
+        for theta in thetas: # vectorize
+            integrand = rhomz * bhs * Ps / (1+zlens)**3. / sigmac / DAz**2 * j0(ells*theta) * ells / 2./ np.pi
+            ints.append( np.trapz(integrand,ells) )
+        return np.asarray(ints)
 
 """
 Mass function

hmvec/tinker.py

@@ -1,5 +1,6 @@
 import numpy as np
 from scipy.interpolate import interp1d
+import os
 
 """
 Implements Tinker et al 2010 and Tinker et al 2008
@@ -59,7 +60,8 @@ def f_nu(nu,zs,delta=200.,norm_consistency=True,
     gamma = gamma0 * (1+zs)**(-0.01)
     unnormalized = (1. + (beta*nu)**(-2.*phi))*(nu**(2*eta))*np.exp(-gamma*nu**2./2.)
     if norm_consistency:
-        izs,ialphas = np.loadtxt("alpha_consistency.txt",unpack=True) # FIXME: hardcoded
+        aroot = os.path.dirname(__file__)+"/../data/alpha_consistency.txt"
+        izs,ialphas = np.loadtxt(aroot,unpack=True) # FIXME: hardcoded
         alpha = interp1d(izs,ialphas,bounds_error=True)(zs)
     return alpha * unnormalized