docs/photometry_and_petrosian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fa820f90",
   "metadata": {},
   "source": [
    "# Photometry\n",
    "\n",
    "To measure the Petrosian properties of galaxies, we construct a photometric curve of growth using a set of concentric apertures with varying radii to measure the flux. From this curve of growth, we will measure different properties of the galaxy including the total flux, characteristic radii such as the Petrosian radius and half-light radius, and the concentration index. \n",
    "\n",
    "To achieve an accurate measurement, we must remove the signature of any background flux and minimize the effects of neighboring sources by masking or removing them.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2f32b9c",
   "metadata": {},
   "source": [
    "## Loading Example Data\n",
    "\n",
    "The following data is a cutout of a group of bright galaxies in Abell 2744. The original data was acquired by the [Hubble Frontier Fields](https://frontierfields.org) team via the WFC3 instrument in the `F105W` filter and can be directly downloaded from the [Mikulski Archive for Space Telescopes](https://archive.stsci.edu/pub/hlsp/frontier/abell2744/images/hst/v1.0/hlsp_frontier_hst_wfc3-60mas_abell2744_f105w_v1.0_drz.fits). \n",
    "\n",
    "We first use `astropy`'s ``CCDData`` to load the example data and visualize it through `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cca798ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.nddata import CCDData\n",
    "\n",
    "image = CCDData.read('data/abell_2744_dwarf_galaxy_f105w.fits.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0df0ebd",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "# Hidden cell\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d936492",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [10, 10]\n",
    "plt.rcParams['image.origin'] = 'lower'\n",
    "\n",
    "vmax = image.data.std() # Use the image std as max and min of all plots \n",
    "vmin = - vmax \n",
    "\n",
    "plt.imshow(image.data, vmin=0, vmax=image.data.std())\n",
    "plt.title(\"Galaxy in Abell 2744\")\n",
    "plt.xlabel(\"Pixels\")\n",
    "plt.ylabel(\"Pixels\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20444fd6",
   "metadata": {},
   "source": [
    "### Estimate data noise at dark area"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34743c6c",
   "metadata": {},
   "source": [
    "In this section, we estimate the noise levels in the image. This estimate is used as a threshold for the detection and segmentation steps. To estimate the background noise, we find a dark region in the image (in this case the corner of the image) and make a cutout using `Cutout2D`. We then compute the statistics for that cutout region.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "101d412e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.nddata import Cutout2D\n",
    "\n",
    "# Estimate data noise at dark area\n",
    "# --------------------------------\n",
    "noise_cutout_pos = (50, 50)\n",
    "noise_cutout_size = 50 \n",
    "noise_cutout = Cutout2D(image, noise_cutout_pos, noise_cutout_size)\n",
    "\n",
    "noise_mean = noise_cutout.data.mean()\n",
    "noise_sigma = noise_cutout.data.std()\n",
    "noise_3_sigma = noise_sigma * 3.\n",
    "noise_8_sigma = noise_sigma * 8.\n",
    "\n",
    "print(noise_mean, noise_3_sigma, noise_8_sigma)\n",
    "\n",
    "# Plot image and noise distribution\n",
    "# ---------------------------------\n",
    "plt.imshow(noise_cutout.data, vmax=noise_mean+noise_3_sigma, vmin=noise_mean-noise_3_sigma)\n",
    "plt.title(\"Dark Patch\")\n",
    "plt.xlabel(\"Pixels\")\n",
    "plt.ylabel(\"Pixels\")\n",
    "plt.show()\n",
    "\n",
    "n, bins, patches = plt.hist(noise_cutout.data.flatten(), bins=35, align='left', color='black')\n",
    "plt.plot(bins[:-1], n, c='r', linewidth=3)\n",
    "plt.axvline(noise_mean, label=\"noise_mean\", linestyle=\"--\")\n",
    "\n",
    "plt.xlabel('Flux Bins [{}]'.format(str(image.unit)))\n",
    "plt.ylabel('Count')\n",
    "plt.title('Noise Histogram')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9088b8c7",
   "metadata": {},
   "source": [
    "We can also use the standard deviation from the `astropy` `sigma_clipped_stats` function to estimate the detection threshold. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec9baf08",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.stats import sigma_clipped_stats\n",
    "image_mean, image_median, image_stddev = sigma_clipped_stats(image.data, sigma=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36d0015a",
   "metadata": {},
   "source": [
    "## Catalogs\n",
    "\n",
    "Before we can perform photometry, we need to first identify the sources in the image. To do this, we use `photutils` and the wrappers in `petrofit`. \n",
    "\n",
    "### Make Catalog \n",
    "\n",
    "We first start by defining the detection threshold and we select this value to be the data standard dev or `noise_8_sigma` computed in the last section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dc2e8af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define detect threshold\n",
    "threshold = noise_8_sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36f2709b",
   "metadata": {},
   "source": [
    "Next, we define the parameters of a Gaussian kernel that is used to smooth the image before segmentation. The `kernel_size` defines the dimensions of the smoothing kernel and `fwhm` which is used to compute the `sigma` value of the Gaussian distribution. The `fwhm` is used to compute the `sigma` value as follows: \n",
    "\n",
    "```\n",
    "sigma = fwhm * gaussian_fwhm_to_sigma\n",
    "kernel = Gaussian2DKernel(sigma, x_size=kernel_size, y_size=kernel_size)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebb05648",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define smoothing kernel\n",
    "kernel_size = 3\n",
    "fwhm = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95102409",
   "metadata": {},
   "source": [
    "We also have to define the number of pixels that make up the smallest object. `npixels` is the number of connected pixels, each greater than the threshold value, that an object must have to be detected. `npixels` must be a positive integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b042b37",
   "metadata": {},
   "outputs": [],
   "source": [
    "npixels = 4**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0d7d08c",
   "metadata": {},
   "source": [
    "\n",
    "To make a catalog of sources and segmentation maps, we use ` petrofit`'s `make_catalog` function. The function returns a source catalog (`photutils.SourceCatalog`), a segmentation image, and, if `deblend=True`, a deblended segmentation image. \n",
    "\n",
    "The `make_catalog` function wraps three steps into one function:\n",
    "\n",
    "1. **Segmentation:**\n",
    "\n",
    "    To identify sources in the image, we first segment the image. The image is smoothed with a gaussian kernel if ` kernel_size ` is provided and clipped at the threshold specified (after smoothing if applicable). The image is then segmented using the ` petrofit.segmentation.make_segments` function, which is a wrapper for `photutil`’s ` detect_sources` functionality. \n",
    "\n",
    "\n",
    "2. **Deblending:**\n",
    "\n",
    "    To further distinguish the sources, we use `photutils.deblend_sources` to deblend the sources into individual galaxies. The `contrast` parameter is the fraction of the total (blended) source flux that a local peak must have (at any one of the multi-thresholds) to be considered as a separate object and `nlevels` is the number of multi-thresholding levels to use. If `deblend` is set to ``False``, ``None`` is returend for `segm_deblend`.\n",
    "\n",
    "3. **Source Catalog:**\n",
    "\n",
    "    Now that we have deblended the sourecs into individual sources, the next step is to create a source catalog (`photutils.SourceCatalog`) that contains properties like `(xcentroid, ycentroid)`, `eccentricity` and `area`. Note that the deblended map is used to make the source catalog but if `deblend` is set to ``False``, the segmentation map is used instead.\n",
    "    \n",
    "By setting the `plot` flag to `True`, we see plots of the segmentation and deblended segmentation maps. The image is plotted along with a color-coded overplot of the segmentation map (each source is a different color). `vmax` and `vmin` can be used the same way as in `plt.imshow`. \n",
    "\n",
    "\n",
    "After the computation, the following objects are returned:\n",
    "\n",
    "- `cat` : A catalog of sources.\n",
    "\n",
    "- `segm` : Segmentation map.\n",
    "\n",
    "- `segm_deblend` : Deblended segmentation map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "574d4356",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.segmentation import make_catalog, plot_segments\n",
    "\n",
    "cat, segm, segm_deblend = make_catalog(    \n",
    "    image.data, \n",
    "    threshold, \n",
    "    deblend=True,                 \n",
    "    kernel_size=kernel_size,                  \n",
    "    fwhm=fwhm, \n",
    "    npixels=npixels,\n",
    "    plot=True, vmax=vmax, vmin=vmin\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d85df3a",
   "metadata": {},
   "source": [
    "To demonstrate the useful information in the catalog, we convert the `SourceCatalog` to an `astropy.table.Table` and display the first 10 objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eabaace",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display source properties\n",
    "print(\"Num of Targets:\", len(cat))\n",
    "\n",
    "# Convert to table\n",
    "cat_table = cat.to_table()\n",
    "\n",
    "cat_table[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0712b1dc",
   "metadata": {},
   "source": [
    "### Plotting Segmentation Maps\n",
    "\n",
    "To plot segmentations, you can use the `plot_segments` function included in PetroFit as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf7e6e2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.segmentation import plot_segments\n",
    "\n",
    "plot_segments(segm, image=image.data, vmax=vmax, vmin=vmin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f1bfef8",
   "metadata": {},
   "source": [
    "As you can see, the segmentation resulted in the identification sources but the sources at the center were classified as a single object because they have interconnecting pixels that are above the threshold.\n",
    "\n",
    "Next we use the same functiton to plot the deblended segmentation map, notice how the central sources are now deblended into individual sources:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e928bfa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_segments(segm_deblend, image=image.data, vmax=vmax, vmin=vmin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab074429",
   "metadata": {},
   "source": [
    "We can also plot the background pixels that are not a part of a source’s segmentation footprint using the `plot_segment_residual` function. We significantly lower the `vmax` and `vmin` values so the background pixels become more apparent. This plot can be used to see if the threshold used to segment the image was too high. If the threshold is high, we would notice bright pixels that are part of the source in this plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8e230ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.segmentation import plot_segment_residual\n",
    "plot_segment_residual(segm, image.data, vmax=vmax/5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62b1115d",
   "metadata": {},
   "source": [
    "## Photometry on Single Source"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cc2bd90",
   "metadata": {},
   "source": [
    "The purpose of this step is to perform aperture photometry to construct a curve of growth that we can use for the Petrosian measurements. \n",
    "\n",
    "### Source Selection\n",
    "\n",
    "For this example, we will focus on a single source. We have included a helper function `order_cat` that will produce a list of indices sorted by a key (default is 'area'). We use the `order_cat` function to identify the source of interest and perform photometry on its cutout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfa11b07",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.photometry import order_cat\n",
    "\n",
    "# Sort and get the largest object in the catalog\n",
    "sorted_idx_list = order_cat(cat, key='area', reverse=True)\n",
    "idx = sorted_idx_list[1] # index 0 is largest \n",
    "source = cat[idx]  # get source from the catalog "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcc8ce3c",
   "metadata": {},
   "source": [
    "### Aperture Radii \n",
    "\n",
    "To construct the curve of growth, we measure the photometry using circular and/or elliptical apertures of varying, concentric radii that are centered on the source. Before we perform the photometry, we need to provide a list of radii that will be used to construct the circular and elliptical apertures. To achieve this, we have provided a helper function that takes in the max radius in pixels (`max_pix`) and the number radii (`n`). The function will return a list of radii by dividing the range `(max_pix/n, max_pix)` into n  equally spaced integers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e59cc07",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.photometry import make_radius_list\n",
    "\n",
    "r_list = make_radius_list(\n",
    "    max_pix=50, # Max pixel to go up to\n",
    "    n=50 # the number of radii to produce \n",
    ")\n",
    "\n",
    "print(repr(r_list))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aae24fef",
   "metadata": {},
   "source": [
    "### Photometry Calculation\n",
    "\n",
    "The `source_photometry` function is used to perform photometry on a single source (`photutils.segmentation.catalog.SourceCatalog`). In this example, we provide the function with the source object, the raw image (note that this is the 2D array and not `CCDData`), segmentation from the last step (`SegmentationImage`) and the list of radii we made using `make_radius_list`. Given these parameters and the options below the function performs photometry as follows:\n",
    "\n",
    "- The positions (max value of source in pixels), elongation and position-angle of the source are determined from the `SourceCatalog` object. The position will be used to center the apertures and the elongation and position angles will be used as parameters of the elliptical apertures. Each of the radii will be assigned an elliptical aperture with these parameters. \n",
    "\n",
    "\n",
    "- If `cutout_size` is defined, the code will use it to make a cutout of that size with the source centered. If the `cutout_size` is larger than the image or contains pixels outside the image, those pixels outside of the image are replaced by `np.nan`.\n",
    "\n",
    "\n",
    "- Once the `cutout_size` is determined, cutouts of the error map (if provided) and image are produced. Before the raw image is cutout, sources that are not the source of interest are masked using the segmentation map. The `mask_background` option gives us the ability to also mask pixels that are considered to be background pixels because they do not belong to any source’s segmentation map. All masked pixels are replaced by `np.nan` and are not counted in the returned area array.  \n",
    "\n",
    "\n",
    "- If `bkg_sub` is set to true, a 2D plane is used to fit pixels that are below a specified sigma from the mean using the `petrofit.fitting.fit_plane` function. The sigma `sigma` value is used to determine noise pixels. Once the pixels above this value are masked, a 2D plane is fit to determine the background. The 2D plane model is then converted into an image and subtracted from the cutout of the target source. `sigma_type` is used to set how this `sigma` value will be used. The `sigma_type` options are `'clip'` and `'bound'`:\n",
    "    - ``'clip'`` (default): Uses `astropy.stats.sigma_clipping.sigma_clip` to clip at the provided `sigma` std value. Note that `sigma` in this case is the number of stds above the mean.\n",
    "    - ``'bound'`` : After computing the mean of the image, clip at `mean - sigma` and `mean + sigma`. Note that `sigma` in this case is a value and not the number of stds above the mean.\n",
    "\n",
    "\n",
    "- The resulting image (after being noise subtracted if `bkg_sub` is set to true) is passed to the `petrofit.photometry.photometry_step` which constructs the apertures and performs photometry. \n",
    "\n",
    "After calculating the photometry at each radius, three arrays are returned:\n",
    "\n",
    "* `flux_arr`: Photometric sum in aperture.\n",
    "\n",
    "* `area_arr`: Exact area of the aperture.\n",
    "\n",
    "* `error_arr`: if error map is provided, error of measurements.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a6edf4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.photometry import source_photometry\n",
    "\n",
    "# Photomerty \n",
    "flux_arr, area_arr, error_arr = source_photometry(\n",
    "    \n",
    "    # Inputs \n",
    "    source, # Source (`photutils.segmentation.catalog.SourceCatalog`)\n",
    "    image.data, # Image as 2D array \n",
    "    segm_deblend, # Deblended segmentation map of image\n",
    "    r_list, # list of aperture radii  \n",
    "    \n",
    "    # Options \n",
    "    cutout_size=max(r_list)*2, # Cutout out size, set to double the max radius  \n",
    "    bkg_sub=True, # Subtract background  \n",
    "    sigma=3, sigma_type='clip', # Fit a 2D plane to pixels within 3 sigma of the mean\n",
    "    plot=True, vmax=vmax, vmin=vmin, # Show plot with max and min defined above\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "531fde92",
   "metadata": {},
   "source": [
    "If `source_photometry`’s `plot` option is set to True, four plots are displayed: \n",
    "- The top left plot shows the cutout with the 2D plane background subtraction and the surrounding sources masked (replaced by `np.nan`).\n",
    "- The top right plot shows the curve of growth with increasing radius. The red lines represent the aperture radius.\n",
    "- The bottom two plots show the source profile sliced at the center of the image in the y and x direction respectively.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd59c51",
   "metadata": {},
   "source": [
    "### Save Photometry Arrays\n",
    "\n",
    "There are many ways of saving the photometry results. For example adding the photometry to the source catalog with the radii as columns. One simple way to save the results for a single source is to save it as a `csv` file using AstroPy `Table`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f675dd4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "\n",
    "t = Table(\n",
    "    data=[r_list, flux_arr, area_arr, error_arr],\n",
    "    names=['r_list', 'flux_arr', 'area_arr', 'error_arr'],   \n",
    ")\n",
    "\n",
    "t.write('temp/abell_2744_galaxy_f105w_photometry.csv', overwrite=True)\n",
    "\n",
    "t[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bf8ff3d",
   "metadata": {},
   "source": [
    "# Petrosian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e375cfc3",
   "metadata": {},
   "source": [
    "## Construct Petrosian from Photometry\n",
    "\n",
    "In this section, we use photometric values stored in flux, aperture area, and radii arrays to construct a `Petrosian` object. The following inputs are needed as inputs:\n",
    "\n",
    "* `r_list`: Array of radii in pixels.\n",
    "\n",
    "* `area_list`: Array of aperture areas.\n",
    "\n",
    "* `flux_list` : Array of photometric flux values.\n",
    "\n",
    "These values should represent  the curve of growth and can be computed by using the [PetroFit photometry tools](#Photometry).\n",
    "\n",
    "The user can also specify the `eta` and `epsilon` values.\n",
    "\n",
    "* `eta` (default=0.2) is the Petrosian value that defines the Petrosian radius.\n",
    "\n",
    "\n",
    "* `epsilon` (default=2) is used to determine the radius of total flux.\n",
    "\n",
    "    * `r_total_flux = r_petrosian * epsilon`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30054dac",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.petrosian import Petrosian\n",
    "\n",
    "p = Petrosian(r_list, area_arr, flux_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa5649d5",
   "metadata": {},
   "source": [
    "## Petrosian Radii\n",
    "\n",
    "PetroFit uses the curve of growth of a galaxy’s flux to compute its Petrosian properties such as Petrosian radius and concentration index."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74bf2a89",
   "metadata": {},
   "source": [
    "### Petrosian Radius\n",
    "\n",
    "The Petrosian radius is defined as the radius at which the Petrosian profile reaches the Eta (`eta`, default=0.2) value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b795a20",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.r_petrosian # in pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a33476cc",
   "metadata": {},
   "source": [
    "### Petrosian Total Flux Radius\n",
    "\n",
    "The Petrosian flux or total flux radius is the radius that ideally encloses all the flux of the galaxy. The Petrosian total flux radius is estimated by multiplying `r_petrosian` with `epsilon` (default=2). \n",
    "\n",
    "`r_total_flux = r_petrosian * epsilon`\n",
    "\n",
    "We can use the `r_total_flux_arcsec` function, by passing it a WCS object, to compute the total flux radius in arcsec."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c81e7bfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.r_total_flux # pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ee5e5da",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.r_total_flux_arcsec(image.wcs) # arcsec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae23476b",
   "metadata": {},
   "source": [
    "### Petrosian Half-Light Radius\n",
    "\n",
    "The half-light radius contains half of the galaxy's total flux. To compute the half-light radius, we find the total flux (flux at `r_total_flux`) and divide it by half to find the “half flux” or “half-light”. We then find the pixel closest to the half-light value and define it as the half-light radius. Please note that interpolation is used between the inputted flux radii to find the radius that best matches the half-light flux.   \n",
    "\n",
    "We can use the `r_half_light_arcsec` function, by passing it a WCS object, to compute the half-light radius in arcsec.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87b5963c",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.r_half_light # pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb618407",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.r_half_light_arcsec(image.wcs) # arcsec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dbab0bc",
   "metadata": {},
   "source": [
    "### Fraction of Flux Radius\n",
    "\n",
    "We can compute a radius that contains a specific fraction of the total flux using the `fraction_flux_to_r` function. For example we can compute the radius that contains 60% of the total flux as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1229844c",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.fraction_flux_to_r(fraction=0.6) # pixels "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6004f5f1",
   "metadata": {},
   "source": [
    "### Concentration Index\n",
    "\n",
    "The concentration index is the ratio of two aperture radii that contain a fraction (percent) of the total flux. It is computed as follows  \n",
    "\n",
    "`concentration_index = 5 * np.log10( r(fraction_2) / r(fraction_1) )`\n",
    "\n",
    "The default is set to `fraction_1 = 0.2` and `fraction_2 = 0.8`. The `concentration_index` function returns the `r_fraction_1`, `r_fraction_2` and `concentration_index`.\n",
    "\n",
    "In these examples, we comput the default `C2080` and `C5090` concentration indices for the input galaxy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "450a0274",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_20, r_80, c2080 = p.concentration_index()  # defualt c2080\n",
    "\n",
    "r_20, r_80, c2080 # Radii in pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cde8f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_50, r_90, c5090 = p.concentration_index(\n",
    "    fraction_1=0.5, \n",
    "    fraction_2=0.9\n",
    ")\n",
    "\n",
    "r_50, r_90, c5090 # Radii in pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bda683b9",
   "metadata": {},
   "source": [
    "### Total Petrosian Flux \n",
    "\n",
    "We can also use `Petrosian` to compute the total Petrosian flux, which is defined as the flux at `r_total_flux`. If the `r_total_flux` is outside the photometric aperture radii, ``np.nan`` is returned.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "769cd8cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "p.total_flux"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3737b4e6",
   "metadata": {},
   "source": [
    "For Hubble data, we can use the `flux_to_abmag` function to convert flux values into `mags` by providing a header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5e01dc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.photometry import flux_to_abmag\n",
    "\n",
    "flux_to_abmag(p.total_flux, header=image.header)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a74a785",
   "metadata": {},
   "source": [
    "## Petrosian Plots\n",
    "\n",
    "### Profile Plot \n",
    "The  Petrosian plot shows the Petrosian profile, the `eta` valued to define the Petrosian radius and the  Petrosian radius. The blue points are the data points and the orange curve is the interpolated values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61f0afdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the Petrosian profile\n",
    "p.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26626023",
   "metadata": {},
   "source": [
    "We can overplot the half-light radius and total flux radius by setting `plot_r=True`. We can also overplot a normalized flux curve of growth by setting `plot_normalized_flux=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed484e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot with radii and overplot normalized flux curve of growth\n",
    "p.plot(plot_r=True, plot_normalized_flux=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbda2822",
   "metadata": {},
   "source": [
    "### Image Overplot\n",
    "\n",
    "Another way to visualize the radii is to overplot them over an image. To do this we first plot the image as usual and use the ` Petrosian.imshow` function to overplot the `r_half_light`, `r_total_flux`, `r_20` and `r_80`. The ` Petrosian.imshow` requires the center of the apertures and plots the radii in pixels. Since elliptical apertures were used, we also provide the `elongation` and orientation (`theta`) of the apertures. We get these values from the source object and use utility functions (`get_source_position`,  `get_source_elong`, `get_source_theta`) to extract them: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d51ca62c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.segmentation import get_source_position, get_source_elong, get_source_theta\n",
    "\n",
    "position = get_source_position(source)\n",
    "elong = get_source_elong(source)\n",
    "theta = get_source_theta(source)\n",
    "\n",
    "p.imshow(position=position, elong=elong, theta=theta, lw=1.25)\n",
    "\n",
    "plt.imshow(image.data, vmax=vmax, vmin=vmin)\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a02e00f3",
   "metadata": {},
   "source": [
    "### Change eta and epsilon\n",
    "\n",
    "We can change the Petrosian `eta` and `epsilon` values after the `Petrosian` object has been initiated by setting their respective attributes. After setting the attributes, all other calculations and plots will use the new values.\n",
    "\n",
    "In this example we copy the `Petrosian` object and change the `eta` and `epsilon` values to see how the radii change. Note how `r_half_light` and `r_total_flux` changed. Since `epsilon = 1`, `r_total_flux = r_petrosian` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9e5c432",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from copy import copy\n",
    "\n",
    "p_copy = copy(p)\n",
    "p_copy.eta = 0.3\n",
    "p_copy.epsilon = 1\n",
    "\n",
    "print('eta =', p_copy.eta)\n",
    "print('epsilon =', p_copy.epsilon)\n",
    "print('r_half_light (old vs new) = {:0.2f} vs {:0.2f}'.format(p.r_half_light, p_copy.r_half_light))\n",
    "print('r_total_flux (old vs new) = {:0.2f} vs {:0.2f}'.format(p.r_total_flux, p_copy.r_total_flux))\n",
    "p_copy.plot(plot_r=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e1fe6dd",
   "metadata": {},
   "source": [
    "## Petrosian Corrections\n",
    "\n",
    "The most accurate value for the total flux is given when the `epsilon` value is a function of the Sersic index. As sersic index has been shown to be correlated with the concentration index,  we can use the measured concentration index to determine the correct value of `epsilon` to return the the radii at which to measure the total flux. \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffce04ba",
   "metadata": {},
   "source": [
    "### Load Correction Grid \n",
    "\n",
    "As a convenience, a grid of values for `epsilon` for different concentration indices and half-light radii is provided such that the radius given by `r_total_flux = r_petrosian * epsilon` will theoretically return 99% of the flux for the source. \n",
    "\n",
    "To load a grid, simply pass its path to the `PetrosianCorrection` class as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98ef1852",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.petrosian import PetrosianCorrection\n",
    "\n",
    "pc = PetrosianCorrection(\"data/concentration_index_grid_f105w_60mas.yaml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b204bb75",
   "metadata": {},
   "source": [
    "### Estimate Sersic Index\n",
    "\n",
    "We can estimate the Sersic index `n` by passing the `estimate_n` function the uncorrected `r_half_light` and `C2080`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5837b99",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimated_n = pc.estimate_n(\n",
    "    p.r_half_light,\n",
    "    p.concentration_index()[-1]\n",
    ")\n",
    "\n",
    "estimated_n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecf650fc",
   "metadata": {},
   "source": [
    "### Estimate Corrected Epsilon\n",
    "\n",
    "`epsilon` is the factor we multiply the Petrosian radius with to find the total flux radius and in turn the half-light radius. Here we estimate the corrected value that will result in a better estimate of the total flux radius and half-light radius by passing the `estimate_epsilon` function the uncorrected `r_half_light` and `C2080`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8caf22d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimated_epsilon = pc.estimate_epsilon(\n",
    "    p.r_half_light,\n",
    "    p.concentration_index()[-1]\n",
    ")\n",
    "\n",
    "estimated_epsilon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ae7114",
   "metadata": {},
   "source": [
    "### Corrected Petrosian\n",
    "\n",
    "Now that we have a better estimate of `epsilon`, we override the default `epsilon` value. We recompute the half-light radius and expect it to match the r_eff parameter above (r_eff is the half-light radius). To save the original `Petrosian` object, we make a new instance with the corrected `epsilon`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83994012",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_corrected = Petrosian(\n",
    "    p.r_list,\n",
    "    p.area_list,\n",
    "    p.flux_list,\n",
    "    epsilon=estimated_epsilon,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29fad565",
   "metadata": {},
   "source": [
    "Plot the corrected Petrosian radii:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d04a9d9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_corrected.plot(plot_r=True, plot_normalized_flux=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66d5484c",
   "metadata": {},
   "source": [
    "We use the same `position`, `elong`, and `theta` values from the [Image Overplot](#Image-Overplot) section to overplot corrected and uncorrected (in black) radii over the image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f19de61",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncorrected (Black)\n",
    "p.imshow(position=position, elong=elong, theta=theta, color='black')\n",
    "\n",
    "# Corrected \n",
    "p_corrected.imshow(position=position, elong=elong, theta=theta)\n",
    "\n",
    "plt.imshow(image.data, vmax=vmax, vmin=vmin)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6cb52e6",
   "metadata": {},
   "source": [
    "### Corrected Petrosian Total Flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db00c715",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Uncorrected Flux = {}\".format(p.total_flux * image.unit))\n",
    "print(\"Corrected Flux = {}\".format(p_corrected.total_flux * image.unit))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}