Added CAMBTransferFunction class

docs/source/code.rst

@@ -109,6 +109,12 @@ Nicaea bindings
 .. autoclass:: lenstools.simulations.Nicaea
 	:members: fromCosmology,convergencePowerSpectrum,shearTwoPoint
 
+CAMB
+====
+
+.. autoclass:: lenstools.simulations.camb.CAMBTransferFunction
+	:inherited-members:
+
 N-body simulation snapshot handling
 ===================================
 
@@ -155,7 +161,7 @@ Directory tree handling
 	:members: path,mkdir,newCollection,getCollection,collections,newTelescopicMapSet,getTelescopicMapSet,telescopicmapsets
 
 .. autoclass:: lenstools.pipeline.simulation.SimulationCollection
-	:members: newRealization,getRealization,realizations,newMapSet,getMapSet,writeCAMB,camb2ngenic
+	:members: newRealization,getRealization,realizations,newMapSet,getMapSet,writeCAMB,loadTransferFunction,camb2ngenic
 
 .. autoclass:: lenstools.pipeline.simulation.SimulationIC
 	:members: newPlaneSet,getPlaneSet,writeNGenIC,writeGadget2

lenstools/pipeline/simulation.py

@@ -18,6 +18,8 @@
 from .settings import *
 
 from .deploy import JobHandler
+
+from ..simulations.camb import CAMBTransferFunction
 from ..simulations import Gadget2SnapshotDE
 from ..simulations.raytracing import PotentialPlane
 
@@ -2305,7 +2307,9 @@ def catalogs(self):
 			return catalogs
 
 
-	#################################################################################################################################
+	###################################################################
+	##################Useful methods for CAMB I/O######################
+	###################################################################
 
 	def writeCAMB(self,z,settings):
 
@@ -2336,8 +2340,32 @@ def writeCAMB(self,z,settings):
 		with self.syshandler.open(os.path.join(self.home_subdir,"camb.p"),"w") as settingsfile:
 			self.syshandler.pickledump(settings,settingsfile)
 
+	def loadTransferFunction(self):
 
-	################################################################################################################################
+		"""
+		Loads in the CDM transfer function calculated with CAMB (k,delta(k)/k^2) for a unit super-horizon perturbation
+
+		:rtype: :py:class:`~lenstools.simulations.camb.CAMBTransferFunction`
+
+		"""
+
+		#Look in the home subdirectory for files camb_transferfunc_zxxxxxxxx.dat
+		tfr_files = self.ls("camb_transferfunc_z*.dat","home")
+		if not len(tfr_files):
+			raise ValueError("Transfer functions have not been computed yet!")
+
+		#Look for the wavenumbers in the first file
+		k = np.loadtxt(tfr_files[0],usecols=(0,))*self.cosmology.h*(self.Mpc_over_h**-1)
+		tfr = CAMBTransferFunction(k)
+
+		#Add transfer function information from each redshift
+		for f in tfr_files:
+			z, = re.search(r"z([0-9.]+)",f).groups()
+			transfer_values = np.loadtxt(f,usecols=(1,))
+			tfr.add(float(z.rstrip(".")),transfer_values)
+
+		#Return to user
+		return tfr
 
 	def camb2ngenic(self,z):
 
@@ -2361,6 +2389,7 @@ def camb2ngenic(self,z):
 
 		print("[+] CAMB matter power spectrum at {0} converted into N-GenIC readable format at {1}".format(camb_ps_file,ngenic_ps_file))
 
+################################################################################################################################
 
 
 ##########################################################

lenstools/simulations/camb.py

@@ -2,6 +2,7 @@
 import StringIO
 
 import numpy as np
+from scipy import interpolate
 import astropy.units as u
 from astropy.cosmology import FLRW
 
@@ -315,3 +316,66 @@ def write(self,output_root,cosmology,redshifts):
 		s.seek(0)
 		return s.read()
 
+#########################################################################################################################################
+
+#CAMB transfer function
+class CAMBTransferFunction(object):
+
+	def __init__(self,k):
+
+		"""
+		:param k: wavenumbers at which the transfer function is computed at
+		:type k: quantity
+
+		"""
+
+		assert k.unit.physical_type=="wavenumber"
+		self._k = k.to((u.Mpc)**-1)
+		self._transfer = dict()
+		self._interpolated = dict()
+
+	def add(self,z,T):
+
+		"""
+		Add transfer function information at redshift z
+
+		:param z: redshift
+		:type z: float.
+
+		:param T: CDM transfer function from CAMB output
+		:type T: array 
+
+		"""
+
+		assert T.shape==k.shape,"There should be exactly one transfer function value for each wavenumber! len(T)={0} len(k)={1}".format(len(T),len(self.k))
+		self._transfer[z] = T 
+
+	def compute(self,z,k):
+
+
+		"""
+		Compute the transfer function at redshift z by linear interpolation
+
+		:param z: redshift
+		:type z: float.
+
+		:param k: wavenumbers at which to compute the transfer function (linearly interpolated with scipy.interp1d)
+		:type k: quantity
+
+		:returns: transfer function at k
+		:rtype: array 
+
+		"""
+
+		assert k.unit.physical_type=="wavenumber"
+
+		#If interpolator has not been built yet for the current redshift, build it
+		if z not in self._interpolated:
+
+			if z not in self._transfer:
+				raise ValueError("There is no information at redshift {0:.3f}".format(z))
+
+			self._interpolated[z] = interpolate.interp1d(self._k.value,self._transfer[z])
+
+		#Use interpolator to compute the transfer function
+		return self._interpolated[z](k.to((u.Mpc)**-1).value)