Source code for lenstools.simulations.limber

"""

.. module:: limber 
	:platform: Unix
	:synopsis: This module implements the tools to compute the convergence power spectrum from the 3D matter power spectrum using the Limber approximation


.. moduleauthor:: Andrea Petri <apetri@phys.columbia.edu>


"""

import numpy as np
from scipy import interpolate,integrate

from astropy import cosmology
from astropy.constants import c
from astropy.units import Mpc,def_unit

##################################################
#############Limber Integrator class##############
##################################################

[docs]class LimberIntegrator(object): """ A 3D power spectrum integrator that will compute the convergence power spectrum using the Limber approximation. :param cosmoModel: One of astropy.cosmology objects (WMAP9 cosmology is set by default) :type cosmoModel: astropy.cosmology """ def __init__(self,cosmoModel=cosmology.WMAP9): assert isinstance(cosmoModel,cosmology.FLRW),"cosmoModel should be a valid astropy cosmology instance!" self.cosmoModel = cosmoModel #Define also the Mpc/h units for convenience self.Mpc_over_h = def_unit("Mpc/h",Mpc/self.cosmoModel.h)
[docs] def load3DPowerSpectrum(self,loader,*args,**kwargs): """ Loads in the matter power spectrum from (pre-computed) external files; args and kwargs are passed to the loader :param loader: must return, in order, k, z, P(k,z) :type loader: function """ self.z,self.kappa,self.power = loader(*args,**kwargs) assert self.power.shape == self.kappa.shape + self.z.shape assert self.z[0] > 0.0,"first redshift must be >0!!" self.setUnits()
[docs] def setUnits(self,kappa_units=None,power_units=None): """ Set the physical units for wavenumber and matter power spectrum, default for length is Mpc """ if kappa_units is None: kappa_units = self.Mpc_over_h**-1 if power_units is None: power_units = self.Mpc_over_h**3 assert (power_units*(kappa_units**3)).physical_type == "dimensionless" assert hasattr(self,"kappa") assert hasattr(self,"power") self.kappa = self.kappa * kappa_units self.power = self.power * power_units
[docs] def convergencePowerSpectrum(self,l): """ Computes the convergence power spectrum with the Limber integral of the 3D matter power spectrum; this still assumes a single source redshift at z0 = max(z) :param l: multipole moments at which to compute the convergence power spectrum :type l: array :returns: array -- the convergence power spectrum at the l values specified """ z = self.z #Power spectrum normalization normalization = (9.0/4)*(self.cosmoModel.Om0)**2*(self.cosmoModel.H0/c)**4 #Compute comoving distances and integral kernel (modify in case there is an arbitrary galaxy distribution) chi = self.cosmoModel.comoving_distance(z) chi0 = chi[-1] kernel = (1.0 - chi/chi0)**2 ############################################################# #Compute the integral for lensing convergence power spectrum# ############################################################# power_interpolation = interpolate.interp1d(self.kappa.to(chi.unit**-1),self.power,axis=0,bounds_error=False,fill_value=0.0) lchi = np.outer(l,1.0/chi).reshape(len(l)*len(z)) power_integrand = power_interpolation(lchi).reshape(len(l),len(z),len(z)).diagonal(axis1=1,axis2=2) * self.power.unit full_integrand = kernel[np.newaxis,:] * (1.0 + z[np.newaxis,:])**2 * power_integrand #Finally compute the integral C = integrate.simps(full_integrand,chi,axis=1) * normalization * full_integrand.unit * chi.unit assert C.unit.physical_type == u"dimensionless" #Return the final result return C.decompose().value