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 G c from lenstools.simulations import PotentialPlane,RayTracer 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 #Rearrange the lenses according to redshift and roll them randomly along the axes tracer.reorderLenses() tracer.randomRoll() WMAP9.comoving_distance(z=2.0) #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) (1.0 * Mpc).to(m) (1.0*Mpc).to(lyr) (WMAP9.H0,WMAP9.Om0) 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$") theta = np.array([[1.0],[1.5]]) * deg tracer.shoot(theta,z=2.0) #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())) #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())) #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) 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) 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())) #Perform forward ray tracing with grid interpolation to compute the image distortion pos_apparent = tracer.shootForward(pos_original,z=2.0,save_intermediate=True) #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])) #Arrange the plot panels on the figure fig = plt.figure(figsize=(12,6)) gs = plt.GridSpec(2,2,height_ratios=[2,1]) ax0 = plt.subplot(gs[0,0],axisbg="black") ax1 = plt.subplot(gs[0,1],axisbg="black") ax2 = plt.subplot(gs[1,:],axisbg="black") #Plot the rays tracer.reset() tracer.displayRays(sample_ray_bucket,z=tracer.redshift[2],fig=fig,ax=[ax0,ax1]) #Plot the sources image image = ax2.imshow(densityEstimator(pos_apparent[0],64),origin="lower",extent=(0.4,1.6,0.8,1.2)) ax2.set_xlabel(r"$x$({0})".format(pos_apparent.unit.to_string())) ax2.set_ylabel(r"$y$({0})".format(pos_apparent.unit.to_string())) ax2.set_title("z={0:.2f}".format(tracer.redshift[2])) fig.tight_layout() import matplotlib.animation as animation #Define the function that updates the animation and makes the movie def make_lensing_movie(i,tracer,image): tracer.displayRays(sample_ray_bucket,z=tracer.redshift[2+i],fig=fig,ax=[ax0,ax1]) image.set_array(densityEstimator(pos_apparent[i],64)) image.axes.set_title("z={0:.2f}".format(tracer.redshift[2+i])) fig.savefig("frame{0}.png".format(i)) #Instantiate the animation ani = animation.FuncAnimation(fig,make_lensing_movie,frames=range(pos_apparent.shape[0]-2),fargs=(tracer,image)) ani.save("lens.m4v",codec="libx264") from IPython.display import HTML from base64 import b64encode with open("lens.m4v","rb") as video: video_encoded = b64encode(video.read()) video_tag = '