diff --git a/.github/workflows/cd.yml b/.github/workflows/cd.yml index 1777a364..194e0a90 100644 --- a/.github/workflows/cd.yml +++ b/.github/workflows/cd.yml @@ -57,7 +57,7 @@ jobs: ls -ltrh ls -ltrh dist - name: Publish to Test PyPI - uses: pypa/gh-action-pypi-publish@v1.12.2 + uses: pypa/gh-action-pypi-publish@v1.12.3 with: repository-url: https://test.pypi.org/legacy/ verbose: true @@ -95,5 +95,5 @@ jobs: name: artifact path: dist - - uses: pypa/gh-action-pypi-publish@v1.12.2 + - uses: pypa/gh-action-pypi-publish@v1.12.3 if: startsWith(github.ref, 'refs/tags') diff --git a/docs/source/_toc.yml b/docs/source/_toc.yml index 3af7da8a..0d15f08e 100644 --- a/docs/source/_toc.yml +++ b/docs/source/_toc.yml @@ -7,6 +7,7 @@ chapters: - file: getting_started - file: install - file: websim + - file: gravlensingintro - file: tutorials/index sections: - file: tutorials/Introduction diff --git a/docs/source/assets/Einstein.npy b/docs/source/assets/Einstein.npy new file mode 100644 index 00000000..f239309e Binary files /dev/null and b/docs/source/assets/Einstein.npy differ diff --git a/docs/source/assets/lensconfiguration.png b/docs/source/assets/lensconfiguration.png new file mode 100644 index 00000000..f7fcd19a Binary files /dev/null and b/docs/source/assets/lensconfiguration.png differ diff --git a/docs/source/assets/multiplanelensing.png b/docs/source/assets/multiplanelensing.png new file mode 100644 index 00000000..fc39b289 Binary files /dev/null and b/docs/source/assets/multiplanelensing.png differ diff --git a/docs/source/gravlensingintro.ipynb b/docs/source/gravlensingintro.ipynb new file mode 100644 index 00000000..5c22b6ee --- /dev/null +++ b/docs/source/gravlensingintro.ipynb @@ -0,0 +1,1540 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Introduction to Gravitational Lensing\n", + "\n", + "Here we provide a primer course on gravitational lensing. It includes interactive examples using `caustics` to demo various effects.\n", + "\n", + "This gravitational lensing primer will give you a good coverage of lensing concepts and intuition. The visual examples will help drive home the mathematical concepts. However, we do not take it upon ourselves in this tutorial to derive and explore the full depths of all the lensing formalism. There are plenty of excellent resources if you would like a more formal background for these concepts. We recommend the following resources for a good start:\n", + "- Introduction to Gravitational Lensing: With Python Examples by Massimo Meneghetti\n", + "- Bartelmann and Schneider 2001\n", + "- Fleury et al. 2022\n", + "- Petkova et al. 2014\n", + "- Keeton 2001\n", + "\n", + "Quick note: we have hidden all the messy plotting commands for the sake of a clean tutorial. The code on display is the code needed to generate the data in the visualizations. If you want to see how the plots are made, just expand the bars above each figure which include the various matplotlib commands." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "import torch\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.animation as animation\n", + "from matplotlib.patches import Circle\n", + "import matplotlib\n", + "from IPython.display import HTML\n", + "import numpy as np\n", + "\n", + "import caustics\n", + "from caustics.constants import arcsec_to_rad" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# setup stuff to define a grid of points to make images\n", + "\n", + "n_pix = 256\n", + "res = 0.015\n", + "fov = n_pix * res\n", + "upsample_factor = 2\n", + "z_l = torch.tensor(0.5, dtype=torch.float32)\n", + "z_s = torch.tensor(1.0, dtype=torch.float32)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Lensing Formalism\n", + "\n", + "### Configuration of a Lens\n", + "\n", + "Before we get too deep into the weeds, lets set the stage. A photon has arrived in your telescope and you are wondering where it came from. Well if it encountered some mass on its path to you (and really what photon hasn't?) then the path it took might look something like this:\n", + "\n", + "![Lens configuration](assets/lensconfiguration.png)\n", + "\n", + "*Figure taken from Bartelmann and Schneider 2001 (figure 11)*\n", + "\n", + "The center of your image is your (0,0) coordinate and this corresponds to the *optical axis* which is the vertical dashed line. The source of the photon is said to exist in the *source plane* which is a distance of $D_s$ from the observer. The mass that deflected the photon is said to exist in the *lens plane* which is a distance $D_l$ (I use $D_l$ rather that $D_d$ for the sake of my sanity) from the observer and $D_{ls}$ from the source plane. Finally, you receive the photons in the *image plane* at the observer. Note that cosmological distances are so large that these may be very accurately treated as infinitely thin planes, even if the source and the mass are in fact 3D structures. We'll talk a little about extended mass distributions in the [Multiplane Lensing](#multiplane-lensing) section. In the source plane the photons are coming from a point $\\vec{\\eta}$ away from the optical axis, this is a physical distance you could measure in kilometers if you like. In the lens plane, the photon intersected at $\\vec{\\xi}$ away from the optical axis (again physical units). In the image plane we can say the photon was essentially at the optical axis.\n", + "\n", + "As an observer, you measure everything in angular units. The photon is the solid line which you see at some position $\\vec{\\theta}$ off from the optical axis. Had there been no mass intercepting the photon, it would have never reached your telescope, instead a photon coming from the same source would have arrived at $\\vec{\\beta}$ angular units from your optical axis. The correction at the source plane which shows how much the photon deviated from its original path is called the *deflection angle* $\\hat{\\vec{\\alpha}}$ and it is essentially a measure of the impact lensing has had on your observation. Note that in the figure the deflection angle is drawn as if the photon traveled backwards, leaving your telescope at angle $\\vec{\\theta}$ getting to the lens plane, getting deflected by $\\hat{\\vec{\\alpha}}$, and finally arriving at the source plane. This is on purpose since the geometry of lensing works in either direction, but computationally it is easier to go backwards (we'll see why in [the lens equation section](#the-lens-equation)).\n", + "\n", + "These are the names and symbols used to describe the configuration of most lensing systems. So far we haven't told you how any of this came to be, but hopefully it will give you a reference point to understand what comes next.\n", + "\n", + "### The Lens Equation\n", + "\n", + "With some simple trigonometry it is possible to relate the physical sizes and the angular sizes:\n", + "\n", + "$$\\vec{\\eta} = \\vec{\\beta} D_s$$\n", + "\n", + "$$\\vec{\\xi} = \\vec{\\theta} D_l$$\n", + "\n", + "Just note that the distances must be in angular diameter distances for these relations to hold (essentially by definition). We can also go a step further and get $\\vec{\\eta}$ using $\\vec{\\theta}$ and $\\hat{\\vec{\\alpha}}$:\n", + "\n", + "$$\\vec{\\eta} = \\vec{\\theta} D_s - \\hat{\\vec{\\alpha}} D_{ls}$$\n", + "\n", + "From here we get the lens equation using the two $\\vec{\\eta}$ formulas:\n", + "\n", + "$$\\vec{\\beta} D_s = \\vec{\\theta} D_s - \\hat{\\vec{\\alpha}} D_{ls}$$\n", + "\n", + "We are tantalizingly close to just cancelling those $D_s$ distances. How about we make a definition $\\vec{\\alpha} \\equiv \\frac{D_{ls}}{D_s}\\hat{\\vec{\\alpha}}$ then we can write the lens equation as:\n", + "\n", + "$$\\vec{\\beta} = \\vec{\\theta} - \\vec{\\alpha}$$\n", + "\n", + "Excellent! Just note that we had to change $\\hat{\\vec{\\alpha}}$ which is the real angular deflection angle that we can do geometry with, to now be $\\vec{\\alpha}$ which has a convenient refactoring such that the lens equation is easier to work with. The modified $\\vec{\\alpha}$ is called the *reduced deflection angle* and usually other quantities that have had a similar transformation to make the lens equation come out nice are called reduced quantities.\n", + "\n", + "This is a very nice and simple equation, one might be fooled into thinking it is easy to work with. Now lets explain why one typically works backwards with photons \"coming from\" the observer and going to the source plane. The deflection angle for almost all lenses is a function of the position in the sky (more deflection closer to the mass, less deflection further away) meaning that it is in fact a function $\\vec{\\alpha}(\\vec{\\theta})$. So to compute $\\vec{\\beta}$ when you know $\\vec{\\theta}$ means simply computing $\\vec{\\alpha}(\\vec{\\theta})$ and the taking the difference $\\vec{\\theta} - \\vec{\\alpha}(\\vec{\\theta})$. But if you want to compute $\\vec{\\theta}$ only knowing $\\vec{\\beta}$, then you will need to search around for a value to input into $\\vec{\\alpha}(\\vec{\\theta})$ until you can get the equation to be satisfied, which means you need to know $\\vec{\\theta}$ to compute $\\vec{\\theta}$! This is a source of great consternation in gravitational lensing researchers, but is a fact we must deal with. In the next sections we'll get into more specifics on all the interesting results that come from this simple equation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 50\n", + "\n", + "# Define a cosmology for the lensing\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "# Define a point mass for the lens plane\n", + "lens = caustics.Point(cosmology=cosmology, x0=0.0, y0=0.0, th_ein=1.0, z_l=z_l)\n", + "\n", + "# Make a bunch of theta values at which to raytrace\n", + "theta_x = torch.linspace(-fov / 2, fov / 2, F, dtype=torch.float32)\n", + "theta_y = torch.zeros(F, dtype=torch.float32)\n", + "# Evaluate the reduced deflection angles\n", + "alpha_x, alpha_y = lens.reduced_deflection_angle(theta_x, theta_y, z_s)\n", + "# Compute the lens equation at these coordinates\n", + "beta_x, beta_y = theta_x - alpha_x, theta_y - alpha_y" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This gives a nice dynamic view of what happens to light rays as they are raytraced through a gravitational lens system. Here we show two perspectives, one is from the angular coordinates where one typically works to generate images, the other is the physical raytracing of a ray in real coordinates (redshift on x-axis, kpc on y-axis). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", + "\n", + "Dl = cosmology.angular_diameter_distance(z_l) * 1000\n", + "Ds = cosmology.angular_diameter_distance(z_s) * 1000\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " ax1.clear()\n", + " ax1.scatter([z_l.item()], [0], color=\"b\", label=\"Point mass\")\n", + " ax1.plot(\n", + " [0, z_l.item()],\n", + " [theta_x[frame].item(), theta_x[frame].item()],\n", + " color=\"r\",\n", + " label=\"Light Ray\",\n", + " )\n", + " ax1.plot(\n", + " [z_l.item(), z_s.item()],\n", + " [beta_x[frame].item(), beta_x[frame].item()],\n", + " color=\"r\",\n", + " )\n", + " ax1.plot(\n", + " [z_l.item(), z_l.item()],\n", + " [theta_x[frame].item(), beta_x[frame].item()],\n", + " color=\"r\",\n", + " )\n", + " ax1.text(\n", + " z_l.item() * 0.5, theta_x[frame].item(), \"$\\\\theta$\", fontsize=16, color=\"k\"\n", + " )\n", + " ax1.text(\n", + " (z_l.item() + z_s.item()) * 0.5,\n", + " beta_x[frame].item(),\n", + " \"$\\\\beta$\",\n", + " fontsize=16,\n", + " color=\"k\",\n", + " )\n", + " ax1.text(\n", + " z_l.item(),\n", + " (theta_x[frame] + beta_x[frame]).item() * 0.5,\n", + " \"$\\\\alpha$\",\n", + " fontsize=16,\n", + " color=\"k\",\n", + " )\n", + " ax1.axvline(\n", + " x=z_l.item(), color=\"grey\", linestyle=\"--\", linewidth=2, label=\"Lens Plane\"\n", + " )\n", + " ax1.axvline(\n", + " x=z_s.item(), color=\"grey\", linestyle=\"-\", linewidth=2, label=\"Source Plane\"\n", + " )\n", + " ax1.axvline(x=0, color=\"grey\", linestyle=\":\", linewidth=2, label=\"Image Plane\")\n", + " ax1.legend(loc=\"upper left\")\n", + " ax1.set_ylim(-fov / 1.9, fov / 1.9)\n", + " ax1.set_title(\"Raytracing (Angular Coordinates)\")\n", + " ax1.set_xlabel(\"redshift z\")\n", + " ax1.set_ylabel(\"Light ray angle [arcsec]\")\n", + "\n", + " ax2.clear()\n", + " ax2.scatter([z_l.item()], [0], color=\"b\", label=\"Point mass\")\n", + " ax2.plot(\n", + " [0, z_l.item()],\n", + " [0.0, (Dl * arcsec_to_rad * theta_x[frame]).item()],\n", + " color=\"r\",\n", + " label=\"Light Ray\",\n", + " )\n", + " ax2.plot(\n", + " [z_l.item(), z_s.item()],\n", + " [\n", + " (Dl * arcsec_to_rad * theta_x[frame]).item(),\n", + " (Ds * arcsec_to_rad * beta_x[frame]).item(),\n", + " ],\n", + " color=\"r\",\n", + " )\n", + " ax2.axvline(\n", + " x=z_l.item(), color=\"grey\", linestyle=\"--\", linewidth=2, label=\"Lens Plane\"\n", + " )\n", + " ax2.axvline(\n", + " x=z_s.item(), color=\"grey\", linestyle=\"-\", linewidth=2, label=\"Source Plane\"\n", + " )\n", + " ax2.text(\n", + " z_l.item(),\n", + " (Dl * arcsec_to_rad * theta_x[frame]).item(),\n", + " \"$\\\\xi$\",\n", + " fontsize=16,\n", + " color=\"k\",\n", + " )\n", + " ax2.text(\n", + " z_s.item(),\n", + " (Ds * arcsec_to_rad * beta_x[frame]).item(),\n", + " \"$\\\\eta$\",\n", + " fontsize=16,\n", + " color=\"k\",\n", + " )\n", + " ax2.scatter([0], [0], color=\"grey\", label=\"Image Plane\")\n", + " ax2.set_ylim(\n", + " -abs(Ds * arcsec_to_rad * beta_x[0].item()) * 1.1,\n", + " abs(Ds * arcsec_to_rad * beta_x[0].item()) * 1.1,\n", + " )\n", + " ax2.set_title(\"Raytracing (Physical Coordinates)\")\n", + " ax2.set_xlabel(\"redshift z\")\n", + " ax2.set_ylabel(\"Physical distance [kpc]\")\n", + "\n", + " return\n", + "\n", + "\n", + "update(0)\n", + "\n", + "ani = animation.FuncAnimation(fig, update, frames=F, interval=80)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Lens Convergence\n", + "\n", + "Lensing is sensitive to all mass on the line of sight. Typically, we describe astronomical massive objects via a 3D density $\\rho(\\vec{\\theta},z)$, but in gravitational lensing we may describe objects as existing on an infinitely thin lens plane. The relevant quantity is then the *surface density* which we may determine by integrating over $z$ to get:\n", + "\n", + "$$\\Sigma(\\vec{\\theta}) = \\int\\rho(\\vec{\\theta}, z)dz$$\n", + "\n", + "This is often normalized by the *critical surface density* at the lens plane $\\Sigma_{cr} = \\frac{c^2D_s}{4\\pi G D_lD_{ls}}$ [^1]. This allows us to compute a dimensionless surface density or *convergence*:\n", + "\n", + "$$\\kappa(\\vec{\\theta}) = \\frac{\\Sigma(\\vec{\\theta})}{\\Sigma_{cr}}$$\n", + "\n", + "The convergence is very useful as it is a dimensionless way to describe the cumulative mass distribution along the line of sight (at least within our lens plane, we'll get more into this in the [multiplane section](#multiplane-lensing)). For some intuition on the convergence, note that a convergence at or above 1 indicates that strong lensing (multiple images) will occur. It is also the most easily translated quantity to say numerical simulations or physical models since we can simply take a projected mass histogram to get the convergence.\n", + "\n", + "Below we have plotted some convergence maps for various lensing mass distributions. In the second subfigure we included a convergence map computed from a numerical simulation (Illustris TNG) by essentially taking a 2D histogram of particle masses.\n", + "\n", + "[^1]: This is different from the cosmological critical density $\\rho_c = \\frac{3H^2}{8\\pi G}$, where $H$ is the expansion rate at a given redshift (i.e. $H_0\\approx 70 \\frac{\\rm km}{{\\rm s}\\cdot {\\rm Mpc}}$ is the expansion rate today). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "cosmology = caustics.FlatLambdaCDM()\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "\n", + "fig, axarr = plt.subplots(1, 5, figsize=(25, 5))\n", + "plt.subplots_adjust(wspace=0.01)\n", + "fig.suptitle(\"Lens Convergence Examples\")\n", + "kappa_map = np.load(\"tutorials/assets/kappa_maps.npz\")[\"kappa_maps\"][1]\n", + "lenses = [\n", + " caustics.SIE(cosmology=cosmology, x0=0, y0=0, q=0.6, phi=np.pi / 3, b=1, z_l=z_l),\n", + " caustics.PixelatedConvergence(\n", + " cosmology=cosmology,\n", + " x0=0,\n", + " y0=0,\n", + " convergence_map=kappa_map,\n", + " pixelscale=fov / kappa_map.shape[0],\n", + " z_l=z_l,\n", + " ),\n", + " caustics.NFW(cosmology=cosmology, x0=0, y0=0, m=1e12, c=5, z_l=z_l),\n", + " caustics.PseudoJaffe(\n", + " cosmology=cosmology,\n", + " x0=0,\n", + " y0=0,\n", + " mass=1e12,\n", + " core_radius=0.2,\n", + " scale_radius=1.0,\n", + " z_l=z_l,\n", + " ),\n", + " caustics.Multipole(\n", + " cosmology=cosmology, m=(3,), x0=0, y0=0, a_m=[2], phi_m=[np.pi / 3], z_l=z_l\n", + " ),\n", + "]\n", + "for i, lens in enumerate(lenses):\n", + " k = lens.convergence(theta_x, theta_y, z_s=z_s)\n", + " k = k.log() if i < 4 else k.tanh()\n", + " axarr[i].imshow(k.numpy())\n", + " axarr[i].set_title(lens.__class__.__name__)\n", + " axarr[i].axis(\"off\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Lens Potential\n", + "\n", + "There are many ways to describe a gravitational lens, with various advantages. Here we will describe the lensing potential $\\Psi$, which is much like the gravitational potential except projected into the 2D plane of our infinitely thin lens:\n", + "\n", + "$$\\hat{\\Psi}(\\vec{\\theta}) = \\frac{2D_{ls}}{D_lD_sc^2}\\int\\phi(D_l\\vec{\\theta},z)dz$$\n", + "\n", + "And the dimensionless version of this potential is:\n", + "\n", + "$$\\Psi(\\vec{\\theta}) = \\frac{D_l^2}{\\xi_0^2}\\hat{\\Psi}$$\n", + "\n", + "A convenient property of the lensing potential is that it encodes all the lens properties locally" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Relating the deflection angle, convergence, and potential\n", + "\n", + "The simplest to convert is the potential:\n", + "\n", + "$$\\vec{\\alpha}(\\vec{\\theta}) = \\vec{\\nabla}_{\\vec{\\theta}}\\Psi(\\vec{\\theta})$$\n", + "\n", + "$$\\kappa(\\vec{\\theta}) = \\frac{1}{2}\\triangle_{\\vec{\\theta}}\\Psi(\\vec{\\theta})$$\n", + "\n", + "Where $\\triangle_{\\vec{\\theta}} = \\frac{\\partial^2}{\\partial\\theta_0^2} + \\frac{\\partial^2}{\\partial\\theta_1^2}$ is the Laplacian. The potential encodes the lensing information locally, what this means is that one may take derivatives to get the other quantities. There is no need to perform any integrals over all space with the potential. This is a useful property, especially in `caustics` where all functions are automatically differentiable.\n", + "\n", + "Next we can see how the convergence transforms:\n", + "\n", + "$$\\Psi(\\vec{\\theta}) = \\frac{1}{\\pi}\\int\\kappa(\\vec{\\theta}')\\ln|\\vec{\\theta} - \\vec{\\theta}'|d\\vec{\\theta}'$$\n", + "\n", + "$$\\vec{\\alpha}(\\vec{\\theta}) = \\frac{1}{\\pi}\\int\\kappa(\\vec{\\theta}')\\frac{\\vec{\\theta}-\\vec{\\theta}'}{|\\vec{\\theta}-\\vec{\\theta}'|}d\\vec{\\theta}'$$\n", + "\n", + "Opposite to the potential, one must perform integrals over all space to convert the convergence into other quantities. One convenience is that these integrals are framed as convolutions, so one may use efficient numerical algorithms (Fast Fourier Transforms) to compute reasonable approximations of these integrals.\n", + "\n", + "It is not common to convert from the deflection angle to either the potential or convergence as this would involve inverting the above equations. This cannot be done in general and so an iterative method would likely be needed." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Deflection of a Point Mass\n", + "\n", + "In this tutorial we will not be too concerned with the derivation of lensing quantities, instead focusing on intuition and core concepts (see [Fleury et al. 2022](https://iopscience.iop.org/article/10.1088/1361-6382/abea2d), the [Meneghetti lecture notes](https://www.ita.uni-heidelberg.de/~jmerten/misc/meneghetti_lensing.pdf), and [Bartelmann and Schneider 2001](https://doi.org/10.1016/S0370-1573(00)00082-X)). So now let's see the deflection angle for a point mass:\n", + "\n", + "$$|\\hat{\\vec{\\alpha}}| = \\frac{4GM}{c^2|\\vec{\\theta}|}$$\n", + "\n", + "We use $|\\hat{\\vec{\\alpha}}|$ and $|\\vec{\\theta}|$ because angular positions on the sky are 2D vectors, but in this equation we only need to know the magnitudes. A point mass is symmetric about its center with the deflection always pointing inwards. For simplicity we have assumed the point mass is on the optical axis, you can of course translate the point mass to other parts of the field of view. Notice that the deflection drops off with distance from the point mass, it should not be too surprising that the amount of deflection decreases further from the mass.\n", + "\n", + "The above deflection is the physical deflection angle which is in regular units. If we convert to reduced deflection angle then things become simpler (as usual).\n", + "\n", + "$$|\\vec{\\alpha}| = \\frac{R_e^2}{|\\vec{\\theta}|}$$\n", + "\n", + "where $R_e$ is the Einstein radius, which is the angular radius at which the point mass has a critical curve with (technically) infinite magnification. We'll get to [magnification](#magnification-and-shear) in a later section. The Einstein radius is defined as:\n", + "\n", + "$$R_e = \\sqrt{\\frac{4GMD_{ls}}{c^2D_lD_s}}$$\n", + "\n", + "Let's see what it looks like for something to be lensed by a point mass. We will use a gaussian blob on the source plane and have it pass behind the point mass." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define a gaussian blob for the source plane\n", + "source = caustics.Sersic(q=1.0, phi=0.0, n=0.5, Re=0.3, Ie=1.0)\n", + "# Define a cosmology for the lensing\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "# Define a point mass for the lens plane\n", + "lens = caustics.Point(cosmology=cosmology, x0=0.0, y0=0.0, th_ein=1.0, z_l=z_l)\n", + "\n", + "# Make a bunch of theta values at which to raytrace\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "# Evaluate the reduced deflection angles\n", + "alpha_x, alpha_y = lens.reduced_deflection_angle(theta_x, theta_y, z_s)\n", + "# Compute the lens equation at these coordinates\n", + "beta_x, beta_y = theta_x - alpha_x, theta_y - alpha_y\n", + "\n", + "# Evaluate the brightness with the source moving along the x-axis\n", + "source_x = torch.linspace(-2, 2, 100)\n", + "source_y = torch.zeros(100)\n", + "unlensed_images = torch.vmap(lambda x: source.brightness(theta_x, theta_y, x))(\n", + " torch.stack((source_x, source_y), dim=1)\n", + ")\n", + "lensed_images = torch.vmap(lambda x: source.brightness(beta_x, beta_y, x))(\n", + " torch.stack((source_x, source_y), dim=1)\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we show the deflection angles as a function of position on the sky. The reduced deflection angle is a vector quantity, so it has two components, the x deflection and the y deflection. We plot them side by side here. You can see that the deflection is strongest near the point mass and weakest further out. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Visualize the reduced deflection angle\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4))\n", + "ax1.imshow(alpha_x, origin=\"lower\", cmap=\"seismic\")\n", + "ax1.set_title(\"Reduced deflection angle [$\\\\alpha_x$]\")\n", + "ax1.axis(\"off\")\n", + "im = ax2.imshow(alpha_y, origin=\"lower\", cmap=\"seismic\")\n", + "fig.colorbar(im, ax=ax3, label=\"deflection magnitude [arcsec]\")\n", + "ax3.axis(\"off\")\n", + "ax2.set_title(\"Reduced deflection angle [$\\\\alpha_y$]\")\n", + "ax2.axis(\"off\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we show what it looks like for a light source passing behind a point mass lens. This is just a gaussian moving from one side of the screen to the other. We have added a point where the lens mass is, and a circle to show the Einstein radius. You will see that when the lensed image reaches the full distortion (the source is exactly aligned with the point mass) the ring it creates aligns exactly with the Einstein radius." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))\n", + "\n", + "# Display the first frame of the image in the first subplot\n", + "im1 = ax1.imshow(\n", + " unlensed_images[0].numpy(), cmap=\"grey\", interpolation=\"bilinear\", vmin=0, vmax=2\n", + ")\n", + "ax1.set_title(\"unlensed\")\n", + "ax1.axis(\"off\")\n", + "im2 = ax2.imshow(\n", + " lensed_images[0].numpy(),\n", + " cmap=\"grey\",\n", + " interpolation=\"bilinear\",\n", + " vmin=0,\n", + " vmax=2,\n", + " extent=[-fov / 2, fov / 2, -fov / 2, fov / 2],\n", + ")\n", + "ax2.scatter([0], [0], color=\"r\")\n", + "c = Circle((0, 0), lens.th_ein.value.item(), fill=False, color=\"r\", linestyle=\"--\")\n", + "ax2.add_patch(c)\n", + "ax2.set_title(\"lensed\")\n", + "ax2.axis(\"off\")\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " im1.set_array(unlensed_images[frame].numpy())\n", + " im2.set_array(lensed_images[frame].numpy())\n", + "\n", + " return im1, im2\n", + "\n", + "\n", + "ani = animation.FuncAnimation(fig, update, frames=lensed_images.shape[0], interval=60)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Deflection of an Extended Mass\n", + "\n", + "Some astronomical objects like black holes and stars are often represented as point masses, but many others such as galaxies and clusters are extended mass distributions. These may be thought of as many point masses combined:\n", + "\n", + "$$\\hat{\\vec{\\alpha}}(\\vec{\\theta}) = \\sum_i\\hat{\\vec{\\alpha}}_i(\\vec{\\theta}_i - \\vec{\\theta}) = \\frac{4G}{c^2}\\sum_iM_i\\frac{\\vec{\\theta}_i - \\vec{\\theta}}{|\\vec{\\theta}_i - \\vec{\\theta}|^2}$$\n", + "\n", + "where the $\\vec{\\theta}_i$ are the positions of the point masses. In the limit that the extended mass is a continuous distribution, the $M_i$ become a density distribution $\\rho(\\vec{\\theta}, z)$ where $z$ is the dimension along the line of sight. Of course for lensing we will use the convergence $\\kappa(\\vec{\\theta})$ since a 3D density is more than we need.\n", + "\n", + "### Singular Isothermal Sphere (SIS)\n", + "\n", + "There are many mass distributions defined by different surface densities. One of the simplest comes from assuming that the mass of the lens is in the form of an ideal gas in a spherical potential. In thermal and hydrostatic equilibrium the density is given as:\n", + "\n", + "$$\\rho(r) = \\frac{\\sigma_v^2}{2\\pi Gr^2}$$\n", + "\n", + "$$\\Sigma(|\\vec{\\theta}|) = \\frac{\\sigma_v^2}{2G|\\vec{\\theta}|}$$ \n", + "\n", + "\n", + "\n", + "$$\\kappa = \\frac{1}{2|\\vec{\\theta}|}$$\n", + "\n", + "And the deflection angle is given by:\n", + "\n", + "$$\\vec{\\alpha}(\\vec{\\theta}) = \\frac{\\vec{\\theta}}{|\\vec{\\theta}|}$$\n", + "\n", + "### Elliptical lens Mass\n", + "\n", + "For any spherically symmetric mass distribution, such as the SIS it is possible to determine a change of coordinates such that the mass has been compressed along one axis, making it an elliptical mass distribution. The *axis ratio* of the ellipse $q = \\frac{b}{a}$ is the ratio of semi-minor to semi-major axis lengths and it defines the ellipse. This looks like:\n", + "\n", + "$$r\\to R = \\sqrt{\\frac{x_1^2}{q} + qx_2^2}$$\n", + "\n", + "In this transformed coordinate space the deflection angles become:\n", + "\n", + "$$\\alpha_x = \\frac{x_1}{qR}\\tilde{\\alpha}_x(R)$$\n", + "\n", + "$$\\alpha_y = \\frac{qx_2}{R}\\tilde{\\alpha}_y(R)$$\n", + "\n", + "where $\\tilde{\\alpha}$ is the unmodified deflection angle of the spherical mass distribution.\n", + "\n", + "### Singular Isothermal Ellipsoid (SIE)\n", + "\n", + "An SIE is just like an SIS except we've converted it into an elliptical mass distribution. Lets try a video like we had in the [point mass section](#deflection-of-a-point-mass). Now we will use an SIE with an axis ratio of 0.5 to see what happens to the background source." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define a gaussian blob for the source plane\n", + "source = caustics.Sersic(q=1.0, phi=0.0, n=0.5, Re=0.05, Ie=1.0)\n", + "# Define a cosmology for the lensing\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "# Define a point mass for the lens plane\n", + "lens = caustics.SIE(cosmology=cosmology, x0=0.0, y0=0.0, q=0.5, phi=0, b=1.0, z_l=z_l)\n", + "# Define a caustics simulator to handle the raytracing\n", + "sim = caustics.LensSource(lens, source, pixelscale=res, pixels_x=n_pix, z_s=z_s)\n", + "\n", + "# Evaluate the brightness with the source moving along the x-axis\n", + "source_x = torch.linspace(-2, 2, 100)\n", + "source_y = torch.zeros(100) + 0.1\n", + "unlensed_images = torch.vmap(lambda x: sim(x, lens_source=False))(\n", + " torch.stack((source_x, source_y), dim=1)\n", + ")\n", + "lensed_images = torch.vmap(sim)(torch.stack((source_x, source_y), dim=1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))\n", + "\n", + "# Display the first frame of the image in the first subplot\n", + "im1 = ax1.imshow(\n", + " unlensed_images[0].numpy(), cmap=\"grey\", interpolation=\"bilinear\", vmin=0, vmax=2\n", + ")\n", + "ax1.set_title(\"unlensed\")\n", + "ax1.axis(\"off\")\n", + "im2 = ax2.imshow(\n", + " lensed_images[0].numpy(),\n", + " cmap=\"grey\",\n", + " interpolation=\"bilinear\",\n", + " vmin=0,\n", + " vmax=2,\n", + " extent=[-2, 2, -2, 2],\n", + ")\n", + "ax2.scatter([0], [0], color=\"r\")\n", + "ax2.set_title(\"lensed\")\n", + "ax2.axis(\"off\")\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " im1.set_array(unlensed_images[frame].numpy())\n", + " im2.set_array(lensed_images[frame].numpy())\n", + "\n", + " return im1, im2\n", + "\n", + "\n", + "ani = animation.FuncAnimation(fig, update, frames=lensed_images.shape[0], interval=60)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Interestingly we see a lot more complexity in this example. The elliptical mass distribution allows for multiple images of the source to appear, rather than just a ring. If you download the tutorial for yourself, you can play around with the parameters to see all sorts of interesting configurations!\n", + "\n", + "Note we made the source very small here to make the four images appear more dramatically." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Point Mass Potential\n", + "\n", + "For a point mass the lensing potential is:\n", + "\n", + "$$\\Psi(\\vec{\\theta}) = R_e\\ln(|\\vec{\\theta}|)$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we show the potential from a series of point masses with a gaussian distribution. For comparison we show the analytic potential field of an SIE. While these aren't exactly the same distribution, it is clear that they share similar general properties." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 100\n", + "# Define a cosmology for the lensing\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "# Define a point mass for the lens plane\n", + "sie_lens = caustics.SIE(\n", + " cosmology=cosmology, x0=0.0, y0=0.0, q=0.99, phi=0, b=0.1, z_l=z_l\n", + ")\n", + "# Define a batched plane lens to combine many point masses\n", + "point_lens = caustics.BatchedPlane(\n", + " cosmology=cosmology, lens=caustics.Point(cosmology=cosmology, z_l=z_l), z_l=z_l\n", + ")\n", + "\n", + "thx, thy = caustics.utils.meshgrid(res, n_pix)\n", + "r = np.abs(np.random.randn(F)) * fov / 4\n", + "th = np.random.rand(F) * 2 * np.pi\n", + "point_locs = torch.tensor(\n", + " np.stack((r * np.cos(th), r * np.sin(th)), axis=1), dtype=torch.float32\n", + ")\n", + "images = [torch.zeros(n_pix, n_pix)]\n", + "for n in range(1, F):\n", + " images.append(\n", + " point_lens.potential(\n", + " thx,\n", + " thy,\n", + " z_s,\n", + " [point_locs[:n, 0], point_locs[:n, 1], 0.6 * torch.ones(n)],\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))\n", + "\n", + "# Display the first frame of the image in the first subplot\n", + "im1 = ax1.imshow(np.zeros((n_pix, n_pix)), cmap=\"grey\", interpolation=\"bilinear\")\n", + "ax1.set_title(\"point masses [000]\")\n", + "ax1.axis(\"off\")\n", + "im2 = ax2.imshow(\n", + " sie_lens.potential(thx, thy, z_s).numpy(),\n", + " cmap=\"grey\",\n", + " interpolation=\"bilinear\",\n", + " extent=[-2, 2, -2, 2],\n", + ")\n", + "ax2.set_title(\"SIE\")\n", + "ax2.axis(\"off\")\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " im1.set_array(images[frame].numpy())\n", + " im1.set_clim(vmin=images[frame].min(), vmax=images[frame].max())\n", + " ax1.set_title(f\"point masses [{frame:03d}]\")\n", + "\n", + " return im1\n", + "\n", + "\n", + "ani = animation.FuncAnimation(fig, update, frames=F, interval=100)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Lens Equation Jacobian\n", + "\n", + "The lens equation describes the motion of rays through a gravitational lensing system. It's jacobian encodes plenty of useful information about the system:\n", + "\n", + "$$A_{ij}(\\vec{\\theta}) = \\frac{\\partial\\beta_i(\\vec{\\theta})}{\\partial\\theta_j} = \\delta_{ij} - \\frac{\\partial\\alpha_i(\\vec{\\theta})}{\\partial\\theta_j} = \\delta_{ij} - \\frac{\\partial\\Psi(\\vec{\\theta})}{\\partial\\theta_i\\partial\\theta_j}$$\n", + "\n", + "Where we have written out multiple ways to compute this jacobian based on whether you have access to the deflection field or the potential as the more computationally efficient quantity." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmology = caustics.FlatLambdaCDM()\n", + "lens = caustics.SIE(cosmology=cosmology, x0=0.0, y0=0.0, q=0.6, phi=0, b=1.0, z_l=z_l)\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "A = lens.jacobian_lens_equation(theta_x, theta_y, z_s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "fig, axarr = plt.subplots(2, 2, figsize=(8, 8))\n", + "for i in range(2):\n", + " for j in range(2):\n", + " axarr[i, j].imshow(A[:, :, i, j].numpy(), origin=\"lower\", cmap=\"seismic\")\n", + " axarr[i, j].axis(\"off\")\n", + " axarr[i, j].set_title(f\"Jacobian component $A_{{{i}, {j}}}$\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Magnification and Shear\n", + "\n", + "The jacobian describes the local behavior of the lensing. We can break the jacobian into two parts and isotropic part and a antisymmetric part like so:\n", + "\n", + "$$ A = A_{iso} + A_{anti} = \\frac{1}{2}tr(A)\\mathbb{I} + \\left(A - \\frac{1}{2}tr(A)\\mathbb{I}\\right)$$\n", + "\n", + "Also for brevity we will use a shorthand for the components of $A$ based on the potential $\\frac{\\partial\\Psi}{\\partial\\theta_i\\partial\\theta_j} = \\Psi_{ij}$ so that now $A_{ij} = \\delta_{ij}-\\Psi_{ij}$. In this case the isotropoc component of the jacobian becomes:\n", + "\n", + "$$A_{iso} = \\left(1 - \\frac{1}{2}(\\Psi_{11} + \\Psi_{22})\\right)\\mathbb{I} = (1-\\kappa)\\mathbb{I}$$\n", + "\n", + "Which you can see is directly related to the convergence (see the [conversions](#relating-the-deflection-angle-convergence-and-potential) section). The remaining antisymmetric part of the jacobian looks like this:\n", + "\n", + "$$A_{anti} = \\left[\\begin{matrix} \n", + "-\\frac{1}{2}(\\Psi_{11} - \\Psi_{22}) & -\\Psi_{12} \\\\ \n", + "-\\Psi_{12} & \\frac{1}{2}(\\Psi_{11} - \\Psi_{22})\n", + "\\end{matrix}\\right]$$\n", + "\n", + "We define a pseudo vector $\\vec{\\gamma} = (\\gamma_1, \\gamma_2)$ where $\\gamma_1 = \\frac{1}{2}(\\Psi_{11} - \\Psi_{22})$ and $\\gamma_2 = \\Psi_{12}$ allowing us to write:\n", + "\n", + "$$A_{anti} = \\left[\\begin{matrix} -\\gamma_1 & -\\gamma_2 \\\\ -\\gamma_2 & \\gamma_1 \\end{matrix}\\right] = -|\\vec{\\gamma}|\\left[\\begin{matrix} \\cos(2\\phi) & \\sin(2\\phi) \\\\ \\sin(2\\phi) & -\\cos(2\\phi) \\end{matrix}\\right]$$\n", + "\n", + "The antisymmetric component, the shear is a $2\\times2$ tensor (hence the $2\\phi$) which represents the distortion of the source along a particular direction. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmology = caustics.FlatLambdaCDM()\n", + "lens = caustics.SIE(cosmology=cosmology, x0=0.0, y0=0.0, q=0.6, phi=0, b=1.0, z_l=z_l)\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "gamma_1, gamma_2 = lens.shear(theta_x, theta_y, z_s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "fig, axarr = plt.subplots(1, 2, figsize=(8, 4))\n", + "axarr[0].imshow(gamma_1.numpy(), origin=\"lower\", cmap=\"seismic\")\n", + "axarr[0].axis(\"off\")\n", + "axarr[0].set_title(\"Shear component $\\\\gamma_1$\")\n", + "axarr[1].imshow(gamma_2.numpy(), origin=\"lower\", cmap=\"seismic\")\n", + "axarr[1].axis(\"off\")\n", + "axarr[1].set_title(\"Shear component $\\\\gamma_2$\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let see what happens when we apply only a shear to an image. this will create distortions along a particular axis." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 100\n", + "hirez = 2\n", + "source = caustics.Pixelated(\n", + " x0=0.0, y0=0.0, image=np.load(\"assets/Einstein.npy\"), pixelscale=0.4 * res / hirez\n", + ")\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "lens = caustics.ExternalShear(cosmology=cosmology, x0=0.0, y0=0.0, z_l=z_l)\n", + "sim = caustics.LensSource(\n", + " lens, source, pixelscale=res / hirez, pixels_x=hirez * n_pix, z_s=z_s\n", + ")\n", + "\n", + "gamma_1 = torch.linspace(-0.8, 0.8, F)\n", + "gamma_2 = torch.linspace(-0.8, 0.8, F)\n", + "lensed_images1 = torch.vmap(sim)(\n", + " torch.stack((gamma_1, torch.zeros_like(gamma_1)), dim=1)\n", + ")\n", + "lensed_images2 = torch.vmap(sim)(\n", + " torch.stack((torch.zeros_like(gamma_2), gamma_2), dim=1)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))\n", + "\n", + "# Display the first frame of the image in the first subplot\n", + "im1 = ax1.imshow(lensed_images1[0].numpy(), cmap=\"grey\")\n", + "ax1.set_title(\"Shear Field $\\\\gamma_1$\")\n", + "ax1.axis(\"off\")\n", + "im2 = ax2.imshow(\n", + " lensed_images2[0].numpy(),\n", + " cmap=\"grey\",\n", + ")\n", + "ax2.set_title(\"Shear Field $\\\\gamma_2$\")\n", + "ax2.axis(\"off\")\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " im1.set_array(lensed_images1[frame].numpy())\n", + " ax1.set_title(f\"Shear Field $\\\\gamma_1 = {gamma_1[frame].item():.2f}$\")\n", + " im2.set_array(lensed_images2[frame].numpy())\n", + " ax2.set_title(f\"Shear Field $\\\\gamma_2 = {gamma_2[frame].item():.2f}$\")\n", + "\n", + " return im1, im2\n", + "\n", + "\n", + "ani = animation.FuncAnimation(fig, update, frames=F, interval=60, repeat=True)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we can look at the magnification:\n", + "\n", + "$$\\mu = \\frac{1}{det(A)} = \\frac{1}{(1-\\kappa)^2 - |\\vec{\\gamma}|^2}$$\n", + "\n", + "Where we see that the magnification is a function of both the convergence and shear, though it can be readily computed from the jacobian matrix." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmology = caustics.FlatLambdaCDM()\n", + "lens = caustics.SIE(cosmology=cosmology, x0=0.0, y0=0.0, q=0.6, phi=0, b=1.0, z_l=z_l)\n", + "mu = lens.magnification(theta_x, theta_y, z_s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "plt.imshow(mu.log().numpy(), origin=\"lower\", cmap=\"seismic\", vmin=-6, vmax=6)\n", + "plt.colorbar()\n", + "plt.axis(\"off\")\n", + "plt.title(\"Magnification $\\log(\\\\mu)$\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The magnification is formally infinite along the [critical curve](#critical-lines-and-caustics) but since we are displaying a pixelated image we simply reach large values." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Critical Lines and Caustics\n", + "\n", + "Given the formula for the magnification $\\mu = \\frac{1}{det(A)}$ it should not be too surprising that interesting things happen when $det(A)=0$. In fact, the contour of positions where this happens is called the critical curve. It is possible to plot these critical curves for a given lens. These critical curves exist in the lens plane, if we raytrace them back to the source plane then they become caustics (ba dumb tis). Lets see what this looks like for a simple SIE." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 50\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "lens = caustics.SIE(cosmology=cosmology, x0=0.0, y0=0.0, phi=0, b=0.6, z_l=z_l)\n", + "q = torch.linspace(0.1, 1.0, F)\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "A = torch.stack([lens.jacobian_lens_equation(theta_x, theta_y, z_s, [q_]) for q_ in q])\n", + "detA = torch.linalg.det(A)\n", + "images = torch.vmap(lambda q: lens.convergence(theta_x, theta_y, z_s, [q]))(q)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(7, 4))\n", + "\n", + "\n", + "def update(frame):\n", + " ax.clear()\n", + " ax.imshow(\n", + " images[frame].log(),\n", + " origin=\"lower\",\n", + " label=\"convergence\",\n", + " extent=[-fov / 2, fov / 2, -fov / 2, fov / 2],\n", + " )\n", + " ax.plot([], [], color=\"#47378D\", label=\"Convergence\", linewidth=7)\n", + " CS = ax.contour(\n", + " theta_x, theta_y, detA[frame], levels=[0.0], colors=\"orange\", linewidths=2\n", + " )\n", + " ax.plot([], [], color=\"orange\", label=\"Critical curve\", linewidth=2)\n", + " paths = CS.allsegs[0]\n", + "\n", + " for path in paths:\n", + " x1 = torch.tensor(list(float(vs[0]) for vs in path))\n", + " x2 = torch.tensor(list(float(vs[1]) for vs in path))\n", + " # raytrace the critical curve to the source plane\n", + " y1, y2 = lens.raytrace(x1, x2, z_s, [q[frame]])\n", + "\n", + " # Plot the caustic\n", + " ax.plot(y1, y2, color=\"r\", label=\"Caustic\", linewidth=2)\n", + " ax.legend()\n", + " ax.set_xlim([-2.5, 2.5])\n", + " ax.set_ylim([-1.5, 1.5])\n", + " ax.set_aspect(\"equal\")\n", + " ax.set_title(f\"Critical Curve and Caustic (q = {q[frame].item():.2f})\")\n", + "\n", + "\n", + "update(0)\n", + "ani = animation.FuncAnimation(fig, update, frames=F, interval=60, repeat=True)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Time Delay\n", + "\n", + "On the cosmic scales of lensing, the finite speed of light becomes especially apparent. This means that when lensing causes multiple images of the same object to appear, they often actually represent the object at different times. There are in fact two reasons for this. One is that the light rays will have taken different paths to reach us the observer, and those different paths are different lengths. Since light travels at a finite speed, a longer path will have taken longer to get to us. The second reason is that gravitational fields bend spacetime and so a certain time delay is induced simply by being in a strong gravitational potential. It turns out that both these effects are often of a similar magnitude and so we need to account for them both. Luckily it is fairly straightforward to do so.\n", + "\n", + "$$T = \\frac{(1+z_l)D_sD_l}{cD_{ls}}\\left(\\frac{1}{2}|\\vec{\\hat{\\alpha}}|^2 - \\Psi\\right)$$\n", + "\n", + "The first term is the geometric time delay caused by different path lengths and is proportional to the physical deflection angle. The second term comes from the strength of the gravitational field called the *Shapiro Delay* and it is proportional to the lensing potential.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 100\n", + "source = caustics.Sersic(x0=0.1, y0=0.1, q=1.0, phi=0.0, n=2, Re=0.05)\n", + "T = caustics.Param(\"T\") # time variable\n", + "source.Ie = (\n", + " lambda p: torch.sin(2 * 2 * np.pi * (p[\"T\"].value % (2 * np.pi))) ** 2\n", + ") # periodic brightness\n", + "source.Ie.link(T)\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "lens = caustics.SIE(cosmology=cosmology, x0=0.0, y0=0.0, q=0.5, phi=0, b=1.0, z_l=z_l)\n", + "t = torch.linspace(0, 1, F)\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "source_images = torch.vmap(lambda t: source.brightness(theta_x, theta_y, [t]))(t)\n", + "# Time delay fields\n", + "td_field = lens.time_delay(theta_x, theta_y, z_s) / 365\n", + "td_geometric = (\n", + " lens.time_delay(\n", + " theta_x, theta_y, z_s, shapiro_time_delay=False, geometric_time_delay=True\n", + " )\n", + " / 365\n", + ")\n", + "td_shapiro = (\n", + " lens.time_delay(\n", + " theta_x, theta_y, z_s, shapiro_time_delay=True, geometric_time_delay=False\n", + " )\n", + " / 365\n", + ")\n", + "# Raytrace including time delay\n", + "beta_x, beta_y = lens.raytrace(theta_x, theta_y, z_s)\n", + "lensed_images = torch.vmap(lambda t: source.brightness(beta_x, beta_y, [t + td_field]))(\n", + " t\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "fig, axarr = plt.subplots(1, 3, figsize=(14, 4))\n", + "plt.subplots_adjust(wspace=0.35)\n", + "im = axarr[0].imshow(td_field.numpy(), origin=\"lower\", cmap=\"grey\")\n", + "cax = fig.add_axes(\n", + " [\n", + " axarr[0].get_position().x1 + 0.01,\n", + " axarr[0].get_position().y0,\n", + " 0.02,\n", + " axarr[0].get_position().height,\n", + " ]\n", + ")\n", + "fig.colorbar(im, cax=cax, label=\"Total Time Delay [years]\")\n", + "axarr[0].axis(\"off\")\n", + "axarr[0].set_title(\"Time Delay Field\")\n", + "\n", + "im = axarr[1].imshow(td_geometric.numpy(), origin=\"lower\", cmap=\"gray\")\n", + "cax = fig.add_axes(\n", + " [\n", + " axarr[1].get_position().x1 + 0.01,\n", + " axarr[1].get_position().y0,\n", + " 0.02,\n", + " axarr[1].get_position().height,\n", + " ]\n", + ")\n", + "fig.colorbar(im, cax=cax, label=\"Geometric Time Delay [years]\")\n", + "axarr[1].axis(\"off\")\n", + "axarr[1].set_title(\"Geometric Time Delay\")\n", + "im = axarr[2].imshow(td_shapiro.numpy(), origin=\"lower\", cmap=\"gray\")\n", + "cax = fig.add_axes(\n", + " [\n", + " axarr[2].get_position().x1 + 0.01,\n", + " axarr[2].get_position().y0,\n", + " 0.02,\n", + " axarr[2].get_position().height,\n", + " ]\n", + ")\n", + "fig.colorbar(im, cax=cax, label=\"Shapiro Time Delay [years]\")\n", + "axarr[2].axis(\"off\")\n", + "axarr[2].set_title(\"Shapiro Time Delay\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))\n", + "\n", + "# Display the first frame of the image in the first subplot\n", + "im1 = ax1.imshow(\n", + " source_images[0].numpy(),\n", + " cmap=\"grey\",\n", + " interpolation=\"bilinear\",\n", + " vmin=0,\n", + " vmax=source_images.max().item(),\n", + ")\n", + "ax1.set_title(\"unlensed\")\n", + "ax1.axis(\"off\")\n", + "im2 = ax2.imshow(\n", + " lensed_images[0].numpy(),\n", + " cmap=\"grey\",\n", + " interpolation=\"bilinear\",\n", + " vmin=0,\n", + " vmax=2,\n", + " extent=[-fov / 2, fov / 2, -fov / 2, fov / 2],\n", + ")\n", + "ax2.set_title(\"lensed\")\n", + "ax2.axis(\"off\")\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " im1.set_array(source_images[frame].numpy())\n", + " im2.set_array(lensed_images[frame].numpy())\n", + "\n", + " return im1, im2\n", + "\n", + "\n", + "ani = animation.FuncAnimation(fig, update, frames=lensed_images.shape[0], interval=60)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Multiplane Lensing\n", + "\n", + "The Universe has more than one mass in it, in fact the Universe has many masses in it. A typical gravitational lensing system may be dominated by a single large mass, but it is very common for there so be many *line of sight halos* along the path which perturb a light ray. In some cases there may even be multiple objects of similar mass along the line of sight which each induce their own strong lensing effect. In all of these situations, the simplifying assumptions we made for the above analysis become invalid and we need to be a bit more careful about how we treat the lensing, this is called *multiplane lensing*. Even including many lenses, it is very accurate to consider each mass as individually existing on a thin lensing plane, but now we will have many such planes. Here is a schematic of a general light ray path through a multiplane lensing system.\n", + "\n", + "![Multiplane Lensing](assets/multiplanelensing.png)\n", + "\n", + "*Credit: Petkova et al. 2013, GLAMER Part II: Multiple Plane Gravitational Lensing*\n", + "\n", + "In this figure, the observer is on the left, and the $D_{i}$ are transverse comoving distances to each plane, $D_{i,j}$ is the transverse comovong distance between planes $i$ and $j$, The $\\vec{X}_{i}$ are the comoving coordinates of where the ray strikes each plane, the $\\vec{\\alpha}_{i}$ are the physical deflection angles at each plane. $\\vec{\\theta}_i$ will be the current angle of the ray at each plane.\n", + "\n", + "To perform raytracing in this scenario we will again start with the ray spreading out from a single point at the observer. Thus the ray will all have comoving coordinates $\\vec{X}_0 = 0$ (starting from a point), and a starting angle of $\\vec{\\theta}_0$. Since the ray may take a circuitous path through all the planes, we need to keep track of both the $\\vec{X}_i$ at each plane and the $\\vec{\\theta}_i$ as it updates at each plane. Note that we use comoving coordinates so we don't have to account for the expansion of the Universe as a function of redshift. This is unlike the single plane case where we could use angular diameter distances since only one lensing and one source plane were involved and the terms could cancel out.\n", + "\n", + "The equivalent of the lens equation for multiplane lensing is very similar:\n", + "\n", + "$$\\vec{\\theta}_s = \\vec{\\theta}_0 - \\sum_i\\frac{\\mathcal{D}_{is}}{\\mathcal{D}_s}\\vec{\\hat{\\alpha}}_i(\\vec{\\theta}_i)$$\n", + "\n", + "where a subscript $s$ indicates the source plane quantity, and $\\mathcal{D}$ is an angular diameter distance. The challenge with this equation is that $\\vec{\\theta}_i$ is an argument of the deflection angle, making the expression cumbersome to evaluate. We now make the transition to comoving distances $D_i = (1+z_i)\\mathcal{D}_i$ and work in comoving coordinates. We can write an expression for the comoving coordinates of the next plane using only the planes before it like so:\n", + "\n", + "$$\\vec{X}_{i+1} = \\vec{X}_i + D_{i+1,i}\\left(\\vec{\\theta}_0 - \\sum_j^{i-1}\\vec{\\hat{\\alpha}}_j(\\vec{X}_j)\\right)$$\n", + "\n", + "Which we can further break down into two recursion functions to be applied at each plane.\n", + "\n", + "$$\\vec{\\theta}_{i+1} = \\vec{\\theta}_i - \\vec{\\hat{\\alpha}}_i(\\vec{X}_i / D_i)$$\n", + "$$\\vec{X}_{i+1} = \\vec{X}_i + D_{i,i+1}\\vec{\\theta}_{i}$$\n", + "\n", + "where we see the first equation essentially breaks down the sum of the angles at each plane (only computing them as we need them, as we figure out the path one plane at a time), the second equation gives the recursion relation but with the angles already accumulated.\n", + "\n", + "Note that the indexing for $\\vec{X}_i$ and $\\vec{\\theta}_i$ is slightly different. It is best for $\\vec{X}_i$ to think of each plane (including the image plane and source plane) as being the index for $i$, so $i=0$ is the image plane, $i=1$ is the first lens plane, until $i=s$ is the source plane. For $\\vec{\\theta}_i$ it is best to think of is as indexing between planes. So for $\\vec{\\theta}_i$ the $i=0$ index is for angles between the image plane and the first lens plane, $i=1$ would be between the first lens plane and the second, and so on. This means that the indexing for $\\vec{X}_i$ ultimately goes one higher than for $\\vec{\\theta}_i$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 200\n", + "\n", + "z_ls = torch.tensor([0.2, 0.4, 0.6, 0.8], dtype=torch.float32)\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "lenses = [\n", + " caustics.Point(\n", + " cosmology=cosmology,\n", + " x0=0.0,\n", + " y0=0.0,\n", + " th_ein=1.0 / np.sqrt(len(z_ls)),\n", + " z_l=z,\n", + " s=1e-3,\n", + " )\n", + " for z in z_ls\n", + "]\n", + "\n", + "theta_x = [torch.linspace(-fov / 2, fov / 2, F, dtype=torch.float32)]\n", + "theta_y = [torch.zeros(F, dtype=torch.float32)]\n", + "D01 = cosmology.transverse_comoving_distance(z_ls[0]) * arcsec_to_rad * 1000\n", + "X = [torch.zeros_like(theta_x[0]), theta_x[0] * D01]\n", + "Y = [torch.zeros_like(theta_y[0]), theta_y[0] * D01]\n", + "\n", + "for i, lens in enumerate(lenses):\n", + " # Next plane\n", + " z_next = z_ls[i + 1] if i + 1 < len(z_ls) else z_s\n", + " # Relevant distances\n", + " Di = cosmology.transverse_comoving_distance(z_ls[i]) * arcsec_to_rad * 1000\n", + " Dij = (\n", + " cosmology.transverse_comoving_distance_z1z2(z_ls[i], z_next)\n", + " * arcsec_to_rad\n", + " * 1000\n", + " )\n", + " # Deflection\n", + " alpha_x, alpha_y = lens.physical_deflection_angle(X[-1] / Di, Y[-1] / Di, z_s)\n", + " # Update angles\n", + " theta_x.append(theta_x[i] - alpha_x)\n", + " theta_y.append(theta_y[i] - alpha_y)\n", + " # Update positions\n", + " X.append(X[-1] + theta_x[i + 1] * Dij)\n", + " Y.append(Y[-1] + theta_y[i + 1] * Dij)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 4))\n", + "Ds = cosmology.transverse_comoving_distance(z_s) * arcsec_to_rad * 1000\n", + "\n", + "\n", + "def update(frame):\n", + " ax.clear()\n", + " ax.scatter(z_ls.numpy(), [0] * len(z_ls), color=\"b\", label=\"Point mass\")\n", + " ax.plot(\n", + " [0, z_ls[0].item()], [X[0][frame], X[1][frame]], color=\"r\", label=\"Light Ray\"\n", + " )\n", + " for i in range(len(z_ls) - 1):\n", + " ax.plot(\n", + " [z_ls[i].item(), z_ls[i + 1].item()],\n", + " [X[i + 1][frame], X[i + 2][frame]],\n", + " color=\"r\",\n", + " )\n", + " ax.plot([z_ls[-1].item(), z_s.item()], [X[-2][frame], X[-1][frame]], color=\"r\")\n", + " for i, z in enumerate(z_ls):\n", + " ax.axvline(\n", + " x=z.item(),\n", + " color=\"grey\",\n", + " linestyle=\"--\",\n", + " linewidth=2,\n", + " label=\"Lens Plane\" if i == 0 else None,\n", + " )\n", + " ax.axvline(\n", + " x=z_s.item(), color=\"grey\", linestyle=\"-\", linewidth=2, label=\"Source Plane\"\n", + " )\n", + " ax.scatter([0], [0], color=\"grey\", label=\"Observer\")\n", + " ax.set_ylim(-Ds * fov / 1.9, Ds * fov / 1.9)\n", + " ax.legend(loc=\"upper left\")\n", + " ax.set_xlabel(\"redshift z\")\n", + " ax.set_ylabel(\"Light ray position [kpc/(1+z)]\")\n", + " ax.set_title(\"Multiplane Raytracing (Physical Coordinates)\")\n", + "\n", + "\n", + "update(0)\n", + "ani = animation.FuncAnimation(fig, update, frames=F, interval=60, repeat=True)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It's neat to watch the path of a ray through space, but now lets see what an image looks like!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = 100\n", + "hirez = 2\n", + "source = caustics.Pixelated(\n", + " x0=0.0,\n", + " y0=0.0,\n", + " image=np.flip(np.load(\"assets/Einstein.npy\"), 0).copy(),\n", + " pixelscale=0.4 * res / hirez,\n", + ")\n", + "z_ls = torch.tensor([0.2, 0.4, 0.6, 0.8], dtype=torch.float32)\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "lenses = [\n", + " caustics.Point(cosmology=cosmology, th_ein=1.0 / np.sqrt(len(z_ls)), z_l=z, s=1e-3)\n", + " for z in z_ls\n", + "]\n", + "lens = caustics.Multiplane(cosmology=cosmology, lenses=lenses)\n", + "sim = caustics.LensSource(\n", + " lens, source, pixelscale=res / hirez, pixels_x=hirez * n_pix, z_s=z_s\n", + ")\n", + "\n", + "# Set the point mases on random rotations\n", + "T = caustics.Param(\"T\") # time variable\n", + "lens_pos_funcs = [\n", + " lambda p: p.meta[\"r\"]\n", + " * torch.cos(p.meta[\"sign\"] * p.meta[\"w\"] * p[\"T\"].value + p.meta[\"phi\"]),\n", + " lambda p: p.meta[\"r\"]\n", + " * torch.sin(p.meta[\"sign\"] * p.meta[\"w\"] * p[\"T\"].value + p.meta[\"phi\"]),\n", + "]\n", + "for sublens in lenses:\n", + " r = np.random.rand() * fov / 3\n", + " sign = np.random.choice([-1, 1])\n", + " w = 2 * np.pi * (np.random.rand() * 1.5 + 0.5)\n", + " phi = np.random.rand() * 2 * np.pi\n", + " sublens.x0.meta = {\"r\": r, \"sign\": sign, \"w\": w, \"phi\": phi}\n", + " sublens.x0 = lens_pos_funcs[0]\n", + " sublens.y0.meta = {\"r\": r, \"sign\": sign, \"w\": w, \"phi\": phi}\n", + " sublens.y0 = lens_pos_funcs[1]\n", + " sublens.x0.link(T)\n", + " sublens.y0.link(T)\n", + "\n", + "t = torch.linspace(0, 1, F)\n", + "images = torch.vmap(sim)(t.unsqueeze(1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# Create animation\n", + "fig, ax = plt.subplots(1, 1, figsize=(5, 5))\n", + "\n", + "\n", + "def update(frame):\n", + " \"\"\"Update function for the animation.\"\"\"\n", + " # Update the image in the first subplot\n", + " T.value = t[frame].item()\n", + " ax.clear()\n", + " ax.imshow(\n", + " images[frame].numpy(),\n", + " cmap=\"grey\",\n", + " extent=[-fov / 2, fov / 2, -fov / 2, fov / 2],\n", + " origin=\"lower\",\n", + " )\n", + " ax.axis(\"off\")\n", + " ax.set_title(f\"Multiplane Lensing (t={t[frame].item():.2f})\")\n", + " for sublens, c in zip(lenses, [\"b\", \"g\", \"orange\", \"r\"]):\n", + " ax.scatter(\n", + " [lens_pos_funcs[0](sublens.x0)], [lens_pos_funcs[1](sublens.y0)], color=c\n", + " )\n", + " return im1, im2\n", + "\n", + "\n", + "update(0)\n", + "ani = animation.FuncAnimation(fig, update, frames=F, interval=60, repeat=True)\n", + "\n", + "plt.close()\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we said earlier, multiplane lenses break a lot of the conventions set out for single plane lenses. However, it is still helpful for some of these quantities to compute effective versions which have similar qualitative behavior. For example, the lens equation is a very convenient construct for computing various quantities. With $\\vec{\\beta} = \\vec{\\theta} - \\vec{\\alpha}(\\vec{\\theta})$ the $\\vec{\\beta}$ term describes the angular position of where the ray strikes the source plane and $\\vec{\\theta}$ describes the angular position in the image plane, while $\\vec{\\alpha}(\\vec{\\theta})$ can be considered as the accumulation of the deflection after a ray has been traced through the multiplane lensing system. This is an effective reduced deflection angle since it actually just hides a more complex recursive deflection angle calculation. Now that we have a lens equation, we can also compute a jacobian of the lens equation and compute quantities from that. For example, the convergence in a single plane lens is the surface mass density divided by the critical surface density, the catch for multiplane is that the critical surface density is defined relative to a single lensing plane. In multiplane lensing there is no single plane that we can use to normalize the surface mass density, and even the surface mass density itself is spread over cosmological distances! There is, however, a way to compute an effective convergence for a multiplane lensing system. Recall that we can relate the convergence to the isotropic component of $A$ thus we can say:\n", + "\n", + "$$\\kappa_{eff, div} = \\frac{1}{2}\\left(\\frac{\\partial\\alpha_1}{\\partial\\theta_1} + \\frac{\\partial\\alpha_2}{\\partial\\theta_2}\\right)$$\n", + "\n", + "Now in single plane lensing recall that $A_{12} = A_{21}$, which is connected to the fact that the deflection angles are curl free, and only have divergence. In multiplane lensing, this is no longer the case. This means that we can't just replace all the complicated multiplane lensing with a single plane using $\\kappa_{eff,div}$ (though that is sometimes a nice approximation). There is now a curl based effective convergence as well: \n", + "\n", + "$$\\kappa_{eff,curl} = \\frac{1}{2}\\left(\\frac{\\partial\\alpha_2}{\\partial\\theta_1} - \\frac{\\partial\\alpha_1}{\\partial\\theta_2}\\right)$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_ls = torch.tensor([0.2, 0.4, 0.6, 0.8], dtype=torch.float32)\n", + "cosmology = caustics.FlatLambdaCDM()\n", + "lenses = [\n", + " caustics.SIE(\n", + " cosmology=cosmology,\n", + " x0=(np.random.rand() - 0.5) * fov * 0.6,\n", + " y0=(np.random.rand() - 0.5) * fov * 0.6,\n", + " q=np.random.rand() * 0.4 + 0.3,\n", + " phi=np.random.rand() * np.pi,\n", + " b=0.5 / np.sqrt(len(z_ls)),\n", + " z_l=z,\n", + " s=1e-3,\n", + " )\n", + " for z in z_ls\n", + "]\n", + "lens = caustics.Multiplane(cosmology=cosmology, lenses=lenses)\n", + "theta_x, theta_y = caustics.utils.meshgrid(res, n_pix)\n", + "kappa_div = lens.effective_convergence_div(theta_x, theta_y, z_s)\n", + "kappa_curl = lens.effective_convergence_curl(theta_x, theta_y, z_s)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the divergence component of the effective convergence is much larger than the curl component. This is typically the case except in very extreme scenarios." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "fig, axarr = plt.subplots(1, 2, figsize=(8, 4))\n", + "im = axarr[0].imshow(\n", + " kappa_div.numpy(), origin=\"lower\", cmap=\"grey\", norm=matplotlib.colors.LogNorm()\n", + ")\n", + "cax = fig.add_axes(\n", + " [\n", + " axarr[0].get_position().x1 + 0.01,\n", + " axarr[0].get_position().y0,\n", + " 0.02,\n", + " axarr[0].get_position().height,\n", + " ]\n", + ")\n", + "fig.colorbar(im, cax=cax, label=\"convergence (div)\")\n", + "axarr[0].axis(\"off\")\n", + "axarr[0].set_title(\"Effective Convergence (divergence)\")\n", + "im = axarr[1].imshow(\n", + " kappa_curl.numpy(),\n", + " origin=\"lower\",\n", + " cmap=\"seismic\",\n", + " vmin=-kappa_curl.abs().max(),\n", + " vmax=kappa_curl.abs().max(),\n", + ")\n", + "cax = fig.add_axes(\n", + " [\n", + " axarr[1].get_position().x1 + 0.01,\n", + " axarr[1].get_position().y0,\n", + " 0.02,\n", + " axarr[1].get_position().height,\n", + " ]\n", + ")\n", + "fig.colorbar(im, cax=cax, label=\"convergence (curl)\")\n", + "axarr[1].axis(\"off\")\n", + "axarr[1].set_title(\"Effective Convergence (curl)\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PY39", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/src/caustics/lenses/base.py b/src/caustics/lenses/base.py index c1420481..46a16747 100644 --- a/src/caustics/lenses/base.py +++ b/src/caustics/lenses/base.py @@ -1076,7 +1076,7 @@ def time_delay( potential = self.potential(x, y, z_s) TD = TD - potential if geometric_time_delay: - ax, ay = self.physical_deflection_angle(x, y, z_s) + ax, ay = self.reduced_deflection_angle(x, y, z_s) fp = 0.5 * (ax**2 + ay**2) TD = TD + fp diff --git a/tests/test_sie.py b/tests/test_sie.py index 16f7a879..3579dea1 100644 --- a/tests/test_sie.py +++ b/tests/test_sie.py @@ -1,8 +1,8 @@ from math import pi from io import StringIO -import lenstronomy.Util.param_util as param_util import torch +import lenstronomy.Util.param_util as param_util from lenstronomy.LensModel.lens_model import LensModel from utils import lens_test_helper diff --git a/tests/test_time_delay.py b/tests/test_time_delay.py new file mode 100644 index 00000000..8396586d --- /dev/null +++ b/tests/test_time_delay.py @@ -0,0 +1,50 @@ +import torch +import numpy as np +import lenstronomy.Util.param_util as param_util +from lenstronomy.LensModel.lens_model import LensModel + +import caustics + +import pytest + + +@pytest.mark.parametrize("q", [0.5, 0.7, 0.9]) +@pytest.mark.parametrize("phi", [0.0, np.pi / 3, np.pi / 2]) +@pytest.mark.parametrize("bx,by", [(0.1, -0.05), (0.2, 0.1), (0.0, 0.0)]) +def test_time_delay_pointsource(q, phi, bx, by): + + # configuration parameters + bx = torch.tensor(bx) + by = torch.tensor(by) + z_l = torch.tensor(0.5) + z_s = torch.tensor(1.0) + + # Define caustics lens + cosmo = caustics.FlatLambdaCDM(name="cosmo") + lens = caustics.SIE(cosmology=cosmo, z_l=z_l, x0=0.0, y0=0.0, q=q, phi=phi, b=1.0) + x, y = lens.forward_raytrace(bx, by, z_s) + + # Define lenstronomy lens + lens_model_list = ["SIE"] + lens_ls = LensModel( + lens_model_list=lens_model_list, z_lens=z_l.item(), z_source=z_s.item() + ) + e1, e2 = param_util.phi_q2_ellipticity(phi=phi, q=q) + kwargs_ls = [{"theta_E": 1.0, "e1": e1, "e2": e2, "center_x": 0.0, "center_y": 0.0}] + + # Compute time delay caustics + tdc = lens.time_delay(x, y, z_s).detach().cpu().numpy() + tdc = tdc - np.min(tdc) + np.sort(tdc) + + # Compute time delay lenstronomy + time_delays = lens_ls.arrival_time( + x.detach().cpu().numpy(), + y.detach().cpu().numpy(), + kwargs_ls, + ) + time_delays = time_delays - np.min(time_delays) + np.sort(time_delays) + + # Compare time delays + assert np.allclose(tdc, time_delays, atol=1e-3)