allow for convergence_ks at the end of raytracing

lenstools/data/lens_default.ini

@@ -43,5 +43,6 @@ first_realization = 1
 
 #Which lensing quantities do we need?
 convergence = True
+convergence_ks = True
 shear = True
 omega = True

lenstools/image/shear.py

@@ -619,7 +619,7 @@ class ShearMap(Spin2):
 	def convergence(self):
 
 		"""
-		Reconstructs the convergence from the E component of the shear
+		Reconstructs the convergence from the E component of the shear (KS method)
 
 		:returns: new ConvergenceMap instance 
 
@@ -632,6 +632,7 @@ def convergence(self):
 		conv = fftengine.irfft2(ft_E)
 
 		#Return the ConvergenceMap instance
-		return ConvergenceMap(conv,self.side_angle)
+		kwargs = dict((k,getattr(self,k)) for k in self._extra_args)
+		return ConvergenceMap(conv,self.side_angle,**kwargs)
 
 

lenstools/pipeline/settings.py

@@ -401,6 +401,7 @@ def _init_commmon(self):
 		#Which lensing quantities do we need?
 		self.tomographic_convergence = False
 		self.convergence = True
+		self.convergence_ks = False
 		self.shear = False
 		self.omega = False
 
@@ -466,18 +467,35 @@ def _read_common(self,options,section):
 
 		self.seed = options.getint(section,"seed")
 
+		###########################################################################################
+
 		try:
 			self.tomographic_convergence = options.getboolean(section,"tomographic_convergence")
 		except NoOptionError:
 			pass
 
 		try:
 			self.convergence = options.getboolean(section,"convergence")
+		except NoOptionError:
+			pass
+
+		try:
+			self.convergence_ks = options.getboolean(section,"convergence_ks")
+		except NoOptionError:
+			pass
+
+		try:
 			self.shear = options.getboolean(section,"shear")
+		except NoOptionError:
+			pass
+
+		try:
 			self.omega = options.getboolean(section,"omega")
 		except NoOptionError:
 			pass
 
+		###########################################################################################
+
 		try:
 			self.integration_type = options.get(section,"integration_type")
 		except NoOptionError:

lenstools/scripts/raytracing.py

@@ -319,12 +319,20 @@ def singleRedshift(pool,batch,settings,batch_id):
 
 			##############################################################################################################################
 
-			if settings.shear:
+			if settings.shear or settings.convergence_ks:
 
-				shearMap = ShearMap(data=np.array([0.5*(jacobian[3]-jacobian[0]),-0.5*(jacobian[1]+jacobian[2])]),angle=map_angle)
-				savename = batch.syshandler.map(os.path.join(save_path,"WLshear_z{0:.2f}_{1:04d}r.{2}".format(source_redshift,r+1,settings.format)))
-				logdriver.info("Saving shear map to {0}".format(savename))
-				shearMap.save(savename) 
+				shearMap = ShearMap(data=np.array([0.5*(jacobian[3]-jacobian[0]),-0.5*(jacobian[1]+jacobian[2])]),angle=map_angle,cosmology=map_batch.cosmology,redshift=source_redshift)
+
+				if settings.shear:
+					savename = batch.syshandler.map(os.path.join(save_path,"WLshear_z{0:.2f}_{1:04d}r.{2}".format(source_redshift,r+1,settings.format)))
+					logdriver.info("Saving shear map to {0}".format(savename))
+					shearMap.save(savename)
+
+				if settings.convergence_ks:
+					convMap = shearMap.convergence() 
+					savename = batch.syshandler.map(os.path.join(save_path,"WLconv-ks_z{0:.2f}_{1:04d}r.{2}".format(source_redshift,r+1,settings.format)))
+					logdriver.info("Saving convergence (KS) map to {0}".format(savename))
+					convMap.save(savename)
 
 			##############################################################################################################################
 

lenstools/simulations/raytracing.py

@@ -897,7 +897,7 @@ def convergence(self):
 
 		#Compute the jacobian and the convergence by tracing it
 		jacobian = self.jacobian()
-		convergence = 0.5*(jacobian[0]+jacobian[3])
+		convergence = -0.5*(jacobian[0]+jacobian[3])
 
 		#Instantiate the convergence map
 		return ConvergenceMap(convergence,angle=self.side_angle)

lenstools/tests/test_ray.py

@@ -63,6 +63,7 @@ def test_ray_simple():
 	conv = dfl.convergence()
 	shear = dfl.shear()
 	omega = dfl.omega()
+	conv_ks = shear.convergence()
 
 	now = time.time()
 	logging.info("Weak lensing calculations completed in {0:.3f}s".format(now-last_timestamp))
@@ -71,6 +72,8 @@ def test_ray_simple():
 	#Finally visualize the result
 	conv.visualize(colorbar=True)
 	conv.savefig("raytraced_convergence.png")
+	conv_ks.visualize(colorbar=True)
+	conv_ks.savefig("raytraced_convergence_ks.png")
 	omega.visualize(colorbar=True)
 	omega.savefig("raytraced_omega.png")
 	shear.visualize(colorbar=True)