notebooks/overview/overview_1_lensing.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Overview: Lensing\n",
        "-----------------\n",
        "\n",
        "When two galaxies are aligned perfectly down the line-of-sight to Earth, the background galaxy's light is bent by the\n",
        "intervening mass of the foreground galaxy. Its light can be fully bent around the foreground galaxy, traversing multiple\n",
        "paths to the Earth, meaning that the background galaxy is observed multiple times. This by-chance alignment of two\n",
        "galaxies is called a strong gravitational lens and a two-dimensional scheme of such a system is pictured below.\n",
        "\n",
        "PyAutoLens is software designed for modeling these strong lensing systems!\n",
        "\n",
        "To begin, lets import autolens and the plot module."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib inline\n",
        "from pyprojroot import here\n",
        "workspace_path = str(here())\n",
        "%cd $workspace_path\n",
        "print(f\"Working Directory has been set to `{workspace_path}`\")\n",
        "\n",
        "import autolens as al\n",
        "import autolens.plot as aplt"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Grids__\n",
        "\n",
        "To describe the deflection of light, **PyAutoLens** uses `Grid2D` data structures, which are two-dimensional\n",
        "Cartesian grids of (y,x) coordinates. \n",
        "\n",
        "Below, we make and plot a uniform Cartesian grid in units of arcseconds. \n",
        "\n",
        "All quantities which are distance units (e.g. coordinate centre's radii) are in units of arc-seconds, as this is the\n",
        "most convenient unit to represent lensing quantities."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "grid = al.Grid2D.uniform(\n",
        "    shape_native=(100, 100),\n",
        "    pixel_scales=0.05,  # The pixel-scale describes the conversion from pixel units to arc-seconds.\n",
        ")\n",
        "\n",
        "grid_plotter = aplt.Grid2DPlotter(grid=grid)\n",
        "grid_plotter.set_title(label=\"Cartesian (y,x) Grid (arcsec)\")\n",
        "grid_plotter.figure_2d()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Light Profiles__\n",
        "\n",
        "We will ray-trace this `Grid2D`'s coordinates to calculate how the lens galaxy's mass deflects the source \n",
        "galaxy's light. We therefore need analytic functions representing a galaxy's light and mass distributions. \n",
        "\n",
        "This requires analytic functions representing the light and mass distributions of galaxies, for example the \n",
        "elliptical `Sersic` `LightProfile`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "sersic_light_profile = al.lp.Sersic(\n",
        "    centre=(0.0, 0.0),\n",
        "    ell_comps=(0.2, 0.1),\n",
        "    intensity=0.005,\n",
        "    effective_radius=2.0,\n",
        "    sersic_index=2.5,\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By passing this profile a `Grid2D`, we can evaluate the light at every (y,x) coordinate on the `Grid2D` and create an \n",
        "image of the Sersic.\n",
        "\n",
        "All images in **PyAutoLens** are in units of electrons per second."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "image = sersic_light_profile.image_2d_from(grid=grid)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The **PyAutoLens** plot module provides methods for plotting objects and their properties, like light profile's image."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "light_profile_plotter = aplt.LightProfilePlotter(\n",
        "    light_profile=sersic_light_profile, grid=grid\n",
        ")\n",
        "light_profile_plotter.set_title(label=\"Image of Sersic Light Profile\")\n",
        "light_profile_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Mass Profiles__\n",
        "\n",
        "**PyAutoLens** uses `MassProfile` objects to represent a galaxy's mass distribution and perform ray-tracing\n",
        "calculations. \n",
        "\n",
        "Below we create an `Isothermal` mass profile and compute its deflection angles on our Cartesian grid, which describe\n",
        "how the source galaxy's light rays are deflected as they pass this mass distribution."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "isothermal_mass_profile = al.mp.Isothermal(\n",
        "    centre=(0.0, 0.0), ell_comps=(0.1, 0.0), einstein_radius=1.6\n",
        ")\n",
        "deflections = isothermal_mass_profile.deflections_yx_2d_from(grid=grid)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets plot the isothermal mass profile's deflection angle map.\n",
        "\n",
        "The black curve on the figure is the tangential critical curve of the mass profile, if you do not know what this is\n",
        "don't worry about it for now!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "mass_profile_plotter = aplt.MassProfilePlotter(\n",
        "    mass_profile=isothermal_mass_profile, grid=grid\n",
        ")\n",
        "mass_profile_plotter.set_title(label=\"Isothermal Deflection Angles (y)\")\n",
        "mass_profile_plotter.figures_2d(\n",
        "    deflections_y=True,\n",
        ")\n",
        "mass_profile_plotter.set_title(label=\"Isothermal Deflection Angles (x)\")\n",
        "mass_profile_plotter.figures_2d(\n",
        "    deflections_x=True,\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "There are many other lensing quantities which can be plotted, for example the convergence and gravitational\n",
        "potential.\n",
        "\n",
        "If you are not familiar with gravitational lensing and therefore are unclear on what the convergence and potential \n",
        "are, don't worry for now!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "mass_profile_plotter.set_title(label=\"Isothermal Mass Convergence\")\n",
        "mass_profile_plotter.figures_2d(\n",
        "    convergence=True,\n",
        ")\n",
        "mass_profile_plotter.set_title(label=\"Isothermal Mass Potential\")\n",
        "mass_profile_plotter.figures_2d(\n",
        "    potential=True,\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Galaxies__\n",
        "\n",
        "A `Galaxy` object is a collection of `LightProfile` and `MassProfile` objects at a given redshift. \n",
        "\n",
        "The code below creates two galaxies representing the lens and source galaxies shown in the strong lensing diagram above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lens_galaxy = al.Galaxy(\n",
        "    redshift=0.5, bulge=sersic_light_profile, mass=isothermal_mass_profile\n",
        ")\n",
        "\n",
        "source_light_profile = al.lp.Exponential(\n",
        "    centre=(0.3, 0.2), ell_comps=(0.1, 0.0), intensity=0.1, effective_radius=0.5\n",
        ")\n",
        "\n",
        "source_galaxy = al.Galaxy(redshift=1.0, bulge=source_light_profile)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The geometry of the strong lens system depends on the cosmological distances between the Earth, the lens galaxy and \n",
        "the source galaxy. It there depends on the redshifts of the `Galaxy` objects. \n",
        "\n",
        "By passing these `Galaxy` objects to a `Tracer` with a `Cosmology` object, **PyAutoLens** uses these galaxy redshifts \n",
        "and a cosmological model to create the appropriate strong lens system."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "tracer = al.Tracer.from_galaxies(\n",
        "    galaxies=[lens_galaxy, source_galaxy], cosmology=al.cosmo.Planck15()\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Ray Tracing__\n",
        "\n",
        "We can now create the image of the strong lens system! \n",
        "\n",
        "When calculating this image, the `Tracer` performs all ray-tracing for the strong lens system. This includes using the \n",
        "lens galaxy's total mass distribution to deflect the light-rays that are traced to the source galaxy. As a result, \n",
        "the source's appears as a multiply imaged and strongly lensed Einstein ring."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "image = tracer.image_2d_from(grid=grid)\n",
        "\n",
        "tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid)\n",
        "tracer_plotter.set_title(label=\"Image of Strong Lens System\")\n",
        "tracer_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `TracerPlotter` includes the mass quantities we plotted previously, which can be plotted as a subplot \n",
        "that plots all these quantities simultaneously.\n",
        "\n",
        "The black and white lines in the source-plane image are the tangential and radial caustics of the mass, which again\n",
        "you do not need to worry about for now if you don't know what that is!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "tracer_plotter.set_title(label=None)\n",
        "tracer_plotter.subplot_tracer()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The tracer is composed of planes. The system above has two planes, an image-plane (at redshift=0.5) and a \n",
        "source-plane (at redshift=1.0). \n",
        "\n",
        "When creating an image via a Tracer, the mass profiles are used to ray-trace the image-plane grid (plotted above) \n",
        "to a source-plane grid, via the mass profile's deflection angles.\n",
        "\n",
        "We can use the tracer`s `traced_grid_2d_list_from` method to calculate and plot the image-plane and source-plane grids."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "traced_grid_list = tracer.traced_grid_2d_list_from(grid=grid)\n",
        "\n",
        "grid_plotter = aplt.Grid2DPlotter(grid=traced_grid_list[0])\n",
        "grid_plotter.set_title(label=\"Image Plane Grid\")\n",
        "grid_plotter.figure_2d()\n",
        "\n",
        "grid_plotter = aplt.Grid2DPlotter(grid=traced_grid_list[1])\n",
        "grid_plotter.set_title(label=\"Source Plane Grid\")\n",
        "grid_plotter.figure_2d()  # Source-plane grid."
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Extending Objects__\n",
        "\n",
        "The **PyAutoLens** API has been designed such that all of the objects introduced above are extensible. `Galaxy` \n",
        "objects can take many `LightProfile`'s and `MassProfile`'s. `Tracer`' objects can take many `Galaxy`'s. \n",
        "\n",
        "If the galaxies are at different redshifts a strong lensing system with multiple lens planes will be created, \n",
        "performing complex multi-plane ray-tracing calculations.\n",
        "\n",
        "To finish, lets create a `Tracer` with 3 galaxies at 3 different redshifts, forming a system with two distinct Einstein\n",
        "rings! The mass distribution of the first galaxy also has separate components for its stellar mass and dark matter."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lens_galaxy_0 = al.Galaxy(\n",
        "    redshift=0.5,\n",
        "    bulge=al.lmp.Sersic(\n",
        "        centre=(0.0, 0.0),\n",
        "        ell_comps=(0.0, 0.05),\n",
        "        intensity=0.5,\n",
        "        effective_radius=0.3,\n",
        "        sersic_index=3.5,\n",
        "        mass_to_light_ratio=0.6,\n",
        "    ),\n",
        "    disk=al.lmp.Exponential(\n",
        "        centre=(0.0, 0.0),\n",
        "        ell_comps=(0.0, 0.1),\n",
        "        intensity=1.0,\n",
        "        effective_radius=2.0,\n",
        "        mass_to_light_ratio=0.2,\n",
        "    ),\n",
        "    dark=al.mp.NFWSph(centre=(0.0, 0.0), kappa_s=0.08, scale_radius=30.0),\n",
        ")\n",
        "\n",
        "lens_galaxy_1 = al.Galaxy(\n",
        "    redshift=1.0,\n",
        "    bulge=al.lp.Exponential(\n",
        "        centre=(0.00, 0.00),\n",
        "        ell_comps=(0.05, 0.05),\n",
        "        intensity=1.2,\n",
        "        effective_radius=0.1,\n",
        "    ),\n",
        "    mass=al.mp.Isothermal(\n",
        "        centre=(0.0, 0.0), ell_comps=(0.05, 0.05), einstein_radius=0.3\n",
        "    ),\n",
        ")\n",
        "\n",
        "source_galaxy = al.Galaxy(\n",
        "    redshift=2.0,\n",
        "    bulge=al.lp.Sersic(\n",
        "        centre=(0.0, 0.0),\n",
        "        ell_comps=(0.0, 0.111111),\n",
        "        intensity=1.4,\n",
        "        effective_radius=0.1,\n",
        "        sersic_index=1.5,\n",
        "    ),\n",
        ")\n",
        "\n",
        "tracer = al.Tracer.from_galaxies(galaxies=[lens_galaxy_0, lens_galaxy_1, source_galaxy])"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is what the lens looks like. \n",
        "\n",
        "Note how crazy the critical curves are!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid)\n",
        "tracer_plotter.set_title(label=\"Image of Complex Strong Lens\")\n",
        "tracer_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Wrap Up__\n",
        "\n",
        "If you are unfamiliar with strong lensing and not clear what some of the above quantities or plots mean, fear not, \n",
        "in chapter 1 of the **HowToLens** lecture series we'll take you through strong lensing theory in detail, whilst teaching\n",
        "you how to use **PyAutoLens** at the same time! \n",
        "\n",
        "Checkout the\n",
        "`tutorials <https://pyautolens.readthedocs.io/en/latest/tutorials/howtolens.html>`_ section of the readthedocs!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "kernelspec": {
      "display_name": "Python 3",
      "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}