Initialize the python environment with the scientific toolkits

In [1]:

```
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
#astropy package for astronomy
from astropy.units import deg,Mpc,m,lyr
from astropy.constants import G,c
from astropy.cosmology import WMAP9
```

In [2]:

```
G
```

Out[2]:

In [3]:

```
c
```

Out[3]:

Let's start by importing the LensTools basic functionalities

In [4]:

```
from lenstools.simulations import PotentialPlane,RayTracer
```

In [5]:

```
tracer = RayTracer(lens_mesh_size=512)
#Add the lenses to the system
for i in range(11,57):
tracer.addLens(PotentialPlane.load("LensTools/Test/Data/lensing/planes/snap{0}_potentialPlane0_normal0.fits".format(i)))
tracer.lens[-1].data *= 20
```

**Solely** for demonstration purposes, we boosted the lensing potential by a factor of 20 to emphatize the lensing effect; in a real world example, the ray deflections are **much** smaller, this is why this kind of gravitational lensing is called *weak*. Next we make sure that the lensed are ordered according to the distance to the observer

In [6]:

```
#Rearrange the lenses according to redshift and roll them randomly along the axes
tracer.reorderLenses()
tracer.randomRoll()
```

In [7]:

```
WMAP9.comoving_distance(z=2.0)
```

Out[7]:

In [8]:

```
#These are the initial ray positions as they hit the first lens
sample_ray_bucket = np.array([[0.0,-0.1,0.1,-0.2,0.2],[0.0,0.0,0.0,0.0,0.0]]) * deg
#This method computes the light ray deflections through the lenses and displays them schematically
tracer.displayRays(sample_ray_bucket,z=2.0)
```

The length measure units we are using, Mpc, translates into more familiar ones as follows

In [9]:

```
(1.0 * Mpc).to(m)
```

Out[9]:

In [10]:

```
(1.0*Mpc).to(lyr)
```

Out[10]:

Einstein General Relativity tells us that the light ray path is the solution of the the following differential equation...$$\frac{d^2 \mathbf{x}(w)}{dw^2} = -\frac{2}{c^2}\nabla_{\mathbf{x}}\Phi(\mathbf{x}(w),w)$$

where $\mathbf{x}=(x,y)$ are the transverse coordinates and $d$ is the longitudinal coordinate. $\Phi(\mathbf{x},w)$ is the three dimensional gravitational potential that describes the spacetime distortions; in our framework the full profile of $\Phi$ is approximated by a discrete set of finely spaced two dimensional lens planes, for which we can define a two dimensional projected lensing potential $$\Psi(\pmb{\beta},w)\approx \frac{3H_0\Omega_m w\Delta}{2c^2}\Phi(\pmb{\beta}w,w)$$

In [11]:

```
(WMAP9.H0,WMAP9.Om0)
```

Out[11]:

In [12]:

```
fig,ax = plt.subplots(1,2,figsize=(16,8))
#Compute and visualize the density and corresponding lensing potential
tracer.lens[30].density().visualize(fig,ax[0],colorbar=True)
tracer.lens[30].visualize(fig,ax[1],colorbar=True)
ax[0].set_title(r"Projected matter density $\sigma$")
ax[1].set_title(r"Lensing potential $\Psi$")
```

Out[12]:

The equation for the light trajectories in angular coordinates looks something like $$\frac{d^2\pmb{\beta}(w)}{dw^2} = -\frac{2}{w}\frac{d\pmb{\beta}(w)}{dw} - \frac{1}{w\Delta}\nabla_{\beta}\Psi(\pmb{\beta}(w),w)$$

Which has to be integrated with the following initial conditions

$$\pmb{\beta}(0) = \pmb{\theta} \,\,\, \frac{d\pmb{\beta}}{dw}(0)=0$$What we are interested in is the final angle, after the final lens has been crossed $\pmb{\beta}_f = \pmb{\beta}(w_f;\pmb{\theta})$

**tracer.shoot()** method will do precisely that: it takes in an initial position $\pmb{\theta}$ and integrates the lensing ODE returning the final angle $\pmb{\beta}_f$

In [13]:

```
theta = np.array([[1.0],[1.5]]) * deg
tracer.shoot(theta,z=2.0)
```

Out[13]:

In [14]:

```
#load unlensed image from png file and discretize it in points
image_unlensed = plt.imread("imageNY.png")[:,:,0]
pos_original = (np.array(np.where(image_unlensed>0)) * 2.0/image_unlensed.shape[0]) * deg
pos_original = np.roll(pos_original,1,axis=0)
pos_original[1] *= -1
pos_original[1] += 2.0*deg
#Visualize the unlensed distribution of sources
fig,ax = plt.subplots()
ax.scatter(pos_original[0],pos_original[1])
ax.set_xlabel(r"$\beta_x$({0})".format(pos_original.unit.to_string()))
ax.set_ylabel(r"$\beta_y$({0})".format(pos_original.unit.to_string()))
```

Out[14]:

**tracer.shoot()** method calculates $\pmb{\beta}_f(\pmb{\theta})$, now we are interested in the inverse problem); this could be quite hard to solve exactly, this is why we adopt an approximate approach: we start from a grid of points at the observer position, $\pmb{\theta}_G$ and we use the ODE solver to compute $\pmb{\beta}_G(\pmb{\theta}_G)$.

In [15]:

```
#Regular grid at the observer location
s = np.linspace(0.0,2.0,16)
theta_G = np.array(np.meshgrid(s,s)) * deg
beta_G = tracer.shoot(theta_G,z=2.0)
#This is how the grids look like
fig,ax = plt.subplots(1,2,figsize=(16,8))
ax[0].scatter(theta_G[0],theta_G[1],color="black")
ax[0].set_title(r"$\theta_G$",fontsize=18)
ax[1].scatter(beta_G[0],beta_G[1],color="red")
ax[1].set_title(r"$\beta_G(z=2.0)$",fontsize=18)
for i in range(2):
ax[i].set_xlabel(r"$x$({0})".format(theta_G.unit.to_string()))
ax[i].set_ylabel(r"$y$({0})".format(theta_G.unit.to_string()))
```

**tracer.shootForward()** which will take care of the ODE solver and the KD tree operations

In [16]:

```
#Perform forward ray tracing with grid interpolation to compute the image distortion
pos_apparent = tracer.shootForward(pos_original,z=2.0)
fig,ax = plt.subplots(1,2,figsize=(16,8))
#Unlensed positions
ax[0].scatter(pos_original[0],pos_original[1])
ax[0].set_xlabel(r"$\beta_x$({0})".format(pos_original.unit.to_string()),fontsize=16)
ax[0].set_ylabel(r"$\beta_y$({0})".format(pos_original.unit.to_string()),fontsize=16)
#Apparent positions
ax[1].scatter(pos_apparent[0],pos_apparent[1])
ax[1].set_xlabel(r"$\theta_x$({0})".format(pos_apparent.unit.to_string()),fontsize=16)
ax[1].set_ylabel(r"$\theta_y$({0})".format(pos_apparent.unit.to_string()),fontsize=16)
```

Out[16]:

This looks a lot nicer if we use a density estimator for the sources instead

In [17]:

```
from lenstools.extern import _gadget
from scipy.ndimage import filters
def densityEstimator(points,grid_size=64):
#Use the LensToold gridding functionality written in C
pt = np.vstack((points,np.zeros(points.shape[1]))).transpose().astype(np.float32)
grid_binning = (np.linspace(0.0,2.0,grid_size+1),np.linspace(0.5,1.5,grid_size+1),np.array([-0.1,0.1]))
density = _gadget.grid3d(pt,grid_binning)
#Smooth the image
return filters.gaussian_filter(density.sum(-1).transpose(),1)
```

In [18]:

```
fig,ax = plt.subplots(1,2,figsize=(16,8))
ax[0].imshow(densityEstimator(pos_original,64),origin="lower",extent=(0.0,2.0,0.5,1.5))
ax[1].imshow(densityEstimator(pos_apparent,64),origin="lower",extent=(0.0,2.0,0.5,1.5))
for i in range(2):
ax[i].set_xlabel(r"$x$({0})".format(pos_apparent.unit.to_string()))
ax[i].set_ylabel(r"$y$({0})".format(pos_apparent.unit.to_string()))
```

**tracer.shootForward()** method allows also to save the intermediate image distortions, adding the lenses from the observer to the source one by one, by specifying the save_intermediate keyword (this takes a little bit more time to run though)

In [19]:

```
#Perform forward ray tracing with grid interpolation to compute the image distortion
pos_apparent = tracer.shootForward(pos_original,z=2.0,save_intermediate=True)
```

In [20]:

```
#Plot some of the distorted images
fig,ax = plt.subplots(2,5,figsize=(40,8))
for i in range(2):
for j in range(5):
ax[i,j].imshow(densityEstimator(pos_apparent[(5*i + j)*4],64),origin="lower",extent=(0.4,1.6,0.8,1.2))
ax[i,j].set_xlabel(r"$x$({0})".format(pos_apparent.unit.to_string()))
ax[i,j].set_ylabel(r"$y$({0})".format(pos_apparent.unit.to_string()))
ax[i,j].set_title("z={0:.2f}".format(tracer.redshift[(5*i + j)*4]))
```