{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4bf8ff3d",
   "metadata": {},
   "source": [
    "# Petrosian\n",
    "\n",
    "In this section, we use the photometric measurements (curve of growth) that were made in the [Photometry](./photometry.ipynb#Photometry) section to construct a Petrosian profile. We use the Petrosian profile to measure various radii and concentrations. For a quick guide on how to construct curves of growth and Petrosian profiles, please see the [Making a Photutils Source Catalog](./quick_start.ipynb#Making-a-Photutils-Source-Catalog) and [Curve of Growth and Petrosian Radii](./quick_start.ipynb#Curve-of-Growth-and-Petrosian-Radii) sections in the [Quick Start](./quick_start.ipynb#Quick-Start) guide.\n",
    "\n",
    "To start with `PetroFit`, simply import it as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad41b77d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import petrofit as pf"
   ]
  },
  {
   "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 (located at `(3.596248, -30.388517)`). 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). The cutout image used in this documentation can be found in the git repository at the following path `petrofit/docs/data/abell_2744_dwarf_galaxy_f105w.fits.gz`.\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",
    "from astropy.io import fits \n",
    "\n",
    "image = CCDData.read('data/abell_2744_dwarf_galaxy_f105w.fits.gz')\n",
    "rms = fits.getdata('data/abell_2744_dwarf_galaxy_f105w_rms.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'] = [6, 6]\n",
    "plt.rcParams['image.origin'] = 'lower'\n",
    "plt.rcParams['font.size'] = 12\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=vmin, vmax=vmax)\n",
    "plt.title(\"Galaxy in Abell 2744\")\n",
    "plt.xlabel(\"Pixels\")\n",
    "plt.ylabel(\"Pixels\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b38c1691",
   "metadata": {},
   "source": [
    "## Load Photometry\n",
    "\n",
    "In the [Photometry Chapter](./photometry.ipynb#Photometry) we constructed a curve of growth for the football shaped galaxy displayed there:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff790ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.plot_target(image.data, (138.5, 89.3), size=100,\n",
    "               vmin=vmin, vmax=vmax, lw=2)\n",
    "plt.xlabel(\"Pixels\")\n",
    "plt.ylabel(\"Pixels\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3922f03e",
   "metadata": {},
   "source": [
    "We load the photometry as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d25db5a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "phot_table = Table.read('data/abell_2744_galaxy_f105w_photometry.ecsv') # Read table \n",
    "\n",
    "# Load data\n",
    "r_list = np.array(phot_table['r_list'])\n",
    "flux_arr = np.array(phot_table['flux_arr'])\n",
    "area_arr = np.array(phot_table['area_arr'])\n",
    "error_arr = np.array(phot_table['error_arr'])"
   ]
  },
  {
   "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 within apertures.\n",
    "\n",
    "* `error_arr`: Array of flux errors.\n",
    "\n",
    "These values should represent  the curve of growth and can be computed by using the [PetroFit photometry tools](./photometry.ipynb#Photometry).\n",
    "\n",
    "We 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": [
    "p = pf.Petrosian(r_list, area_arr, flux_arr, flux_err=error_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": [
    "print(\"{:0.4f} ± {:0.4f} pix\".format(p.r_petrosian, p.r_petrosian_err))"
   ]
  },
  {
   "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": [
    "print(\"{:0.4f} ± {:0.4f} pix\".format(p.r_total_flux, p.r_total_flux_err))"
   ]
  },
  {
   "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 input 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": [
    "print(\"{:0.4f} ± {:0.4f} pix\".format(p.r_half_light, p.r_half_light_err))"
   ]
  },
  {
   "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": [
    "r_60 = p.fraction_flux_to_r(fraction=0.6) # pixels \n",
    "r_60_err = p.fraction_flux_to_r_err(fraction=0.6) # pixels \n",
    "print(\"{:0.4f} ± {:0.4f} pix\".format(r_60, r_60_err))"
   ]
  },
  {
   "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": [
    "print(\"{:0.4f} ± {:0.4f} pix\".format(p.total_flux, p.total_flux_err))"
   ]
  },
  {
   "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": [
    "pf.hst_flux_to_abmag(p.total_flux, header=image.header) # Mag"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a74a785",
   "metadata": {},
   "source": [
    "## Petrosian and COG 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."
   ]
  },
  {
   "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": "4c7a1de7",
   "metadata": {},
   "source": [
    "Much in the same way we can plot the curve of growth as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed484e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot curve of growth \n",
    "p.plot_cog(flux_unit=image.unit)\n",
    "\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. \n",
    "To do this we first plot the image as usual and use the ` Petrosian.imshow` \n",
    "function to overplot the `r_half_light`, `r_total_flux`, `r_20` and `r_80`. \n",
    "The ` Petrosian.imshow` requires the center of the apertures and plots the \n",
    "radii in pixels. Since elliptical apertures were used, we also provide the `elongation` \n",
    "and orientation (`theta`) of the apertures. We get these values from the source object \n",
    "and use utility functions (`get_source_position`,  `get_source_elong`, `get_source_theta`) as showen in\n",
    "the [Photometry](./photometry.ipynb#Photometry) section. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d51ca62c",
   "metadata": {},
   "outputs": [],
   "source": [
    "position = phot_table.meta['position']\n",
    "elong = phot_table.meta['elong']\n",
    "theta = phot_table.meta['theta']\n",
    "\n",
    "p.imshow(position=position, elong=elong, theta=theta, lw=1.25)\n",
    "\n",
    "plt.imshow(image.data, vmax=vmax, vmin=vmin)\n",
    "plt.xlabel(\"Pixels\")\n",
    "plt.ylabel(\"Pixels\")\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. To review:\n",
    "\n",
    "- `eta`: The Petrosian index, which dictates the value at which the Petrosian radius is determined. \n",
    "  Typically, the default value is 0.2.\n",
    "\n",
    "- `epsilon`: A multiplication factor used to scale the Petrosian radius to compute the \"epsilon radius\" (`r_epsilon`). \n",
    "  If epsilon is `0.275`, for instance, then `r_epsilon` would be `0.275 * r_p`.\n",
    "\n",
    "- `r_epsilon`: Refers to the \"epsilon radius\". It's the radius that encompasses an epsilon_fraction of the galaxy's total flux. \n",
    "  This radius is derived from the Petrosian radius by scaling it with the epsilon factor: `r_epsilon = epsilon * r_p`.\n",
    "\n",
    "- `epsilon_fraction`: Represents the fraction of the total flux contained within the `r_epsilon`. For instance, a value of 0.5 suggests that 50% of the galaxy's total flux is inside this radius.\n",
    "\n",
    "By default `epsilon_fraction = total_flux_fraction = 0.99` or `0.98`, therefore the defualt is `r_epsilon = r_total_flux`.\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.13\n",
    "p_copy.epsilon = 0.275\n",
    "p_copy.epsilon_fraction = 0.5\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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2bfedfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig, axs = plt.subplots(1,2, figsize=(12, 6))\n",
    "plt.sca(axs[0])\n",
    "p_copy.plot(plot_r=True)\n",
    "\n",
    "plt.sca(axs[1])\n",
    "p_copy.plot_cog(plot_r=True, flux_unit=image.unit)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8e64814",
   "metadata": {},
   "source": [
    "## Corrections and Approximations\n",
    "\n",
    "### Important Note\n",
    "\n",
    "Before going into corrections and approximations of Petrosian profiles, \n",
    "it's important to understand the broader context in which they're utilized. \n",
    "A significant portion of astronomical literature, including foundational analyses from the Sloan Digital Sky Survey (SDSS), \n",
    "predominantly relies on standard parameters: $\\eta=0.2$ and $\\epsilon=2$ for deriving the total flux. \n",
    "The rationale behind this prevalent practice is two-fold:\n",
    "\n",
    "1. **Noise Considerations**: In real-world observations, the faint extremities of galaxies often lie beneath the instrument's noise threshold. \n",
    "   This can render intricate corrections superfluous since the adjustments might not significantly alter the results.\n",
    "  \n",
    "2. **Practicality Over Precision**: While the discussed corrections provide higher precision, \n",
    "   their application demands optimal background subtraction. Hence, they're often reserved for \n",
    "   scenarios that have ideal subtraction for flux measurements or for measuring a galaxy's size at higher precision \n",
    "   given low signal-to-noise.\n",
    "\n",
    "In other words, the default settings have historically been good enough. \n",
    "Given this backdrop, the forthcoming sub-section goes into specific techniques \n",
    "for those seeking to venture beyond the default parameters, either out of \n",
    "necessity or research specificity.\n",
    "\n",
    "\n",
    "### Techniques\n",
    "\n",
    "For accurately determining the total flux, the epsilon value should ideally be a function of the Sérsic index. As the Sérsic index is correlated with the concentration index, we can leverage the measured concentration index to derive the appropriate value of epsilon, which, in turn, helps pinpoint the radii for total flux measurements.\n",
    "\n",
    "To estimate epsilon, a few methods are at our disposal:\n",
    "\n",
    "**Utilizing Approximations**:\n",
    "\n",
    "- Using the uncorrected concentration index (`U2080`), we can approximate both the Sérsic index and the corrected epsilon. This method relies on approximations derived from fitting these relationships against a standard Sérsic profile.\n",
    "  \n",
    "  **Pros**:\n",
    "  \n",
    "  - A swift and straightforward means to approximate the Sérsic index and the corrected epsilon.\n",
    "  \n",
    "  - Only requires a moderate level of background subtraction.\n",
    "  \n",
    "  **Cons**:\n",
    "  \n",
    "  - Does not account for PSF, which notably affects high Sérsic index profiles with a small $r_{50}$.\n",
    "  \n",
    "  - Assumes an ideal singular Sérsic component profile (no noise limit). This could lead to a total flux radius that's significantly larger than what is practically measurable.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5da316b",
   "metadata": {},
   "outputs": [],
   "source": [
    "U2080 = p.c2080 # Uncorrected concentration index\n",
    "C2080 = pf.PetroApprox.u2080_to_c2080(U2080) # Corrected concentration index\n",
    "sersic_index = pf.PetroApprox.c2080_to_n(C2080) # Sersic index\n",
    "corrected_epsilon = pf.PetroApprox.n_to_epsilon(sersic_index) # Corrected epsilon\n",
    "\n",
    "# Correct:\n",
    "p_approx = copy(p) # Copy \n",
    "p_approx.epsilon = corrected_epsilon\n",
    "\n",
    "print(corrected_epsilon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72d5eecc",
   "metadata": {},
   "source": [
    "**Adopting `epsilon=0.5` & `epsilon_fraction=0.5`**:\n",
    "\n",
    "- This method leverages the observation that for a wide range of Sérsic indices, \n",
    "  the half-light radius ($r_{50}$) is approximately half of the Petrosian radius. \n",
    "  This approximation is particularly helpful if our interest is primarily in determining $r_{50}$ \n",
    "  and if the curve of growth does not extend to the total flux radius of the uncorrected petrosian profile.\n",
    "\n",
    "  **Pros**:\n",
    "  \n",
    "  - Simplifies the process by using a consistent approximation across different Sérsic indices.\n",
    "  \n",
    "  - Does not require the curve of growth to extend to the total flux radius.\n",
    "  \n",
    "  - Can be faster and more direct when only $r_{50}$ is of interest.\n",
    "  \n",
    "  **Cons**:\n",
    "  \n",
    "  - Might not be the most accurate method for determining total flux or for profiles with \n",
    "    significant deviations from the average. This is because the slope of the curve of growth is steep near $r_{50}$.\n",
    "  \n",
    "  - Assumes the relationship between the half-light radius and petrosian radius is consistent across different galaxy profiles.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf68df37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correct:\n",
    "p_r_50 = copy(p) # Copy \n",
    "p_r_50.epsilon = 0.5  # epsilon\n",
    "p_r_50.epsilon_fraction = 0.5 # epsilon corrsponds to 50% of total flux\n",
    "\n",
    "print(p.r_50, p_approx.r_50, p_r_50.r_50)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1057bce3",
   "metadata": {},
   "source": [
    "This results in a $r_{50}$ that is about a pixel off from the approximated and uncorrected profiles in this example. \n",
    "But the total flux radii will be poorly approximated by this method (~30% error):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72fd23ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p.r_total_flux, p_approx.r_total_flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ebad44",
   "metadata": {},
   "source": [
    "**Utilizing Correction Grids**:\n",
    "\n",
    "- This approach involves the simulation of a grid of Sérsic profiles, each varying in terms of Sérsic index and effective radius. This simulation results in a lookup table that associates the petrosian radius, the uncorrected half light radius, and the `C2080` value with specific Sérsic index and corrected epsilon values.\n",
    "  \n",
    "  **Pros**:\n",
    "  \n",
    "  - Amongst available methods, this stands as the most precise, especially when focusing on single Sérsic component profiles.\n",
    "  \n",
    "  - Takes into account the PSF, which is important for high Sérsic index profiles due to it \"smearing out\" radii.\n",
    "  \n",
    "  **Cons**:\n",
    "  \n",
    "  - Conducting simulations can be time-consuming.\n",
    "  \n",
    "  - Vastly different PSFs across multiple bands necessitate individual grids for each.\n",
    "  \n",
    "  - Effective corrections demand accurate image background subtraction.\n",
    "\n",
    "Plese see the [Correction Grids](./correction_grids.ipynb#Correction-Grids) section for more details on how to generate correction grids. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f45a2830",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read F105W grid:\n",
    "pc = pf.PetrosianCorrection.read('./data/f105w_psf_corr.ecsv')\n",
    "\n",
    "# Pass uncorrected p to the correct function\n",
    "p_corrected = pc.correct(p)\n",
    "\n",
    "print(p.r_50, p_approx.r_50, p_r_50.r_50, p_corrected.r_50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5e78d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig, axs = plt.subplots(1,2, figsize=(12, 6))\n",
    "plt.sca(axs[0])\n",
    "ax = p.plot_cog(plot_r=True, flux_unit=image.unit)\n",
    "ax.set_title('Before Correction')\n",
    "\n",
    "plt.sca(axs[1])\n",
    "ax = p_corrected.plot_cog(plot_r=True, flux_unit=image.unit)\n",
    "ax.set_title('After Correction')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5df49f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(image.data, vmax=vmax, vmin=vmin)\n",
    "p_corrected.imshow(position=position, elong=elong, theta=theta)\n",
    "plt.xlabel(\"Pixels\")\n",
    "plt.ylabel(\"Pixels\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}