docs/correction_grids.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "09eb9c87",
   "metadata": {},
   "source": [
    "# Correction Girds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb595677",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "# Hidden cell\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c20bce2",
   "metadata": {},
   "source": [
    "## Define Sampling Points \n",
    "\n",
    "The correction grid is generated along with a set of half-light radii and Sersic indices. The generator loops through the Sersic indices for each half-light radius in the radius list. For this documentation, we define a small set of values, with the Sersic indices of a Gaussian and de Vaucouleurs’ profile.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcec33cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "r_eff_list = np.array([15, 20, 25])\n",
    "n_list = np.array([1, 2, 4])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f5cf7ae",
   "metadata": {},
   "source": [
    "## Correction Grid\n",
    "\n",
    "We do a simple API call to generate the correction grid. We provide a PSF and oversampling rule as well. Oversampling becomes important for small half-light radii since the model needs to be sampled well in the center. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef5216be",
   "metadata": {},
   "source": [
    "### Load PSF \n",
    "\n",
    "We load an HST `F105W` PSF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "potential-princeton",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Load PSF image (2D array)\n",
    "PSF = fits.getdata('data/f105w_psf.fits.gz')\n",
    "\n",
    "# Normalize PSF \n",
    "PSF = PSF / PSF.sum()\n",
    "\n",
    "# Note that the PSF shape is odd on all sides\n",
    "print(\"PSF Shape = {}\".format(PSF.shape))\n",
    "\n",
    "# Plot PSF and use vmax and vmin to show difraction spikes\n",
    "plt.imshow(PSF, vmin=0, vmax=5e-4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "562ba847",
   "metadata": {},
   "source": [
    "### Generate Correction Grid\n",
    "\n",
    "The `generate_petrosian_sersic_correction` grid follows the following steps to generate the correction grid:\n",
    "\n",
    "- Computes the total flux (`total_flux`) of an ideal Sersic profile with the sampling points using the ` petrofit.modeling.models.sersic_enclosed` function and setting the radius to `np.inf`.\n",
    "- Computes the radius equal to `total_flux * 0.99`.\n",
    "- Makes a PSF convolved Sersic Model image and measures the photometry. \n",
    "- Measures the uncorrected Petrosian radius (). \n",
    "- Computes the correct `epsilon` value as ` corrected_epsilon = r_total_flux / r_petrosian`.\n",
    "- It also saves other values, such as the uncorrected `C2080` which can be used to mapto the Sersic index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dedicated-helen",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.correction import generate_petrosian_sersic_correction\n",
    "\n",
    "petrosian_grid = generate_petrosian_sersic_correction(\n",
    "    output_yaml_name='temp/example_correction_gid.yaml',\n",
    "    psf=PSF,\n",
    "    r_eff_list=r_eff_list,\n",
    "    n_list=n_list,\n",
    "    oversample=('x_0', 'y_0', 10, 5),\n",
    "    plot=False,\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85317016",
   "metadata": {},
   "source": [
    "## Test Grid\n",
    "\n",
    "We test the correction grid by loading it into a ` PetrosianCorrection` and using it estimate values.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "222dde66",
   "metadata": {},
   "outputs": [],
   "source": [
    "from petrofit.petrosian import PetrosianCorrection\n",
    "\n",
    "pc = PetrosianCorrection('temp/example_correction_gid.yaml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aecb98d8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "corrected_epsilon = pc.estimate_epsilon(\n",
    "    r_hl_pet=15, \n",
    "    c2080pet=3,  \n",
    "    verbose=True\n",
    ")\n",
    "\n",
    "print(corrected_epsilon)"
   ]
  }
 ],
 "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
}