In [1]:
%matplotlib inline
import evillens as evil
import matplotlib.pyplot as plt
Build analytic SIE lens, and compute its kappa distribution.
In [2]:
SIElensA = evil.AnalyticSIELens(0.4,1.5)
SIElensA.setup_grid(NX = 800, NY = 800, pixscale = 0.025, n = 8)
SIElensA.build_kappa_map(centroid = [0.01,0.01])
In [3]:
SIElensA.plot("kappa")
print('Dimensions of Convergence map:', SIElensA.kappa.shape)
Compute deflection angles, we'll need those later. Then write kappa map out to file so that we can read it as basic evil gravitational lens object.
In [4]:
SIElensA.deflect()
SIElensA.write_kappa_to("analytic_SIE_kappa.fits")
Now, start a gravitational lens object, and read in the kappa from the fits file we just created, and compute deflection angles numerically.
In [5]:
SIElensB = evil.GravitationalLens(0.4,1.5)
SIElensB.read_kappa_from("analytic_SIE_kappa.fits")
SIElensB.plot("kappa")
In [6]:
SIElensB.deflect()
Now let's plot maps of the deflection angles for each lens, analytic and numerical, for comparison. Recall, SIElensA is the analytic SIE lens, and SIElensB was read in as a kappa map from the SIElensA output kappa map FITS file."
In [7]:
SIElensA.plot("alpha")
SIElensB.plot("alpha")
"These look similar, but let's check the difference maps to make sure. The rms fractional difference is a useful statistic for summarizing the accuracy."
In [8]:
import numpy as np
SIElensC = SIElensB-SIElensA
SIElensC.plot('alpha')
print("RMS accuracy (%):",100.0*np.std(SIElensC.alpha_x/SIElensA.alpha_x), 100.0*np.std(SIElensC.alpha_y/SIElensA.alpha_y))
The peak accuracy is lower in the very center of the maps (the lens center), but unless we are studying central images this may not be important. The next test could be to assess deflection angle accuracy when a small substructure is present - but this will be best done by using the recovered substructure mass as a metric.
Just for completeness, the difference maps are plotted below on a different color scale. We have also plotted the total error magnitude of the form
$\vec{\epsilon} = \hat{\epsilon}\sqrt{\epsilon_{x}^2+\epsilon_{y}^{2}}$
where $\hat{\epsilon}$ denotes the direction of the error. The vectors plotted above the error magnitudes show the direction of the uncertainty.
In [20]:
plt.imshow(SIElensC.alpha_y,origin='lower',vmin=-0.03,vmax=0.03,extent = (np.min(SIElensC.image_x),np.max(SIElensC.image_x),np.min(SIElensC.image_y),np.max(SIElensC.image_y)))
plt.xlabel('x / arcsec')
plt.ylabel('y / arcsec')
plt.colorbar()
Out[20]:
In [22]:
plt.imshow(SIElensC.alpha_x,vmin=-0.03,vmax=0.03,extent = (np.min(SIElensC.image_x),np.max(SIElensC.image_x),np.min(SIElensC.image_y),np.max(SIElensC.image_y)))
plt.xlabel('x / arcsec')
plt.ylabel('y / arcsec')
plt.colorbar()
Out[22]:
In [21]:
plt.imshow(np.sqrt(SIElensC.alpha_x**2+SIElensC.alpha_y**2),extent = (np.min(SIElensC.image_x),np.max(SIElensC.image_x),np.min(SIElensC.image_y),np.max(SIElensC.image_y)),vmin=0,vmax=0.02)
plt.colorbar()
plt.xlabel('x / arcsec')
plt.ylabel('y / arcsec')
plt.quiver(SIElensC.image_x[::6,::6], SIElensC.image_y[::6,::6], SIElensC.alpha_x[::6,::6],SIElensC.alpha_y[::6,::6])
Out[21]:
In [12]:
In [ ]: