notebooks/WFC3/point_spread_function/hst_point_spread_function.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "32375168-9e76-49dc-a5c0-f52ef009e1dc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "<a id=\"top\"></a>\n",
    "# HST WFC3 Point Spread Function Modeling   <img style=\"float: right;\" src=\"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"Space Telescope Logo\" width=\"200px\"/>\n",
    "\n",
    "***\n",
    "<a id=\"learning\"></a>\n",
    "## Learning Goals\n",
    "This notebook demonstrates how to create Point Spread Function (PSF) models for HST [Wide Field Camera 3 (WFC3)](https://www.stsci.edu/hst/instrumentation/wfc3) observations.\n",
    "\n",
    "The tutorial will allow users to complete the following tasks:\n",
    "\n",
    " - Perform simple aperture photometry on science images for comparison with PSF photometry.\n",
    " - Retrieve or create a PSF model for WFC3 single exposures (FLCs) and drizzled (DRZs) images.\n",
    " - Complete basic science workflows including PSF photometry, subtraction, and decomposition.\n",
    "\n",
    "<a id=\"toc\"></a>\n",
    "## Table of Contents\n",
    "\n",
    "**&ensp;[Learning Goals](#learning)** <br>\n",
    "**&ensp;[1. Introduction](#introduction)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[1.1 Environment Setup](#environment)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[1.2 Import Packages](#import)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[1.3 Download Data](#data)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[1.4 Aperture Photometry](#aperture)** <br>\n",
    "**&ensp;[2. Exposure PSF Models](#mainc)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[2.1 Analytical PSF Models](#analytical)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[2.2 Empirical PSF Models](#empirical)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[2.2.1 Default Models](#empirical-1)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[2.2.2 Focus Perturbation](#empirical-2)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[2.2.3 Spatial Perturbation](#empirical-3)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[2.3 Stacked PSF Models](#stacked)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[2.3.1 Stellar Stacks](#stellar)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[2.3.2 MAST Stacks](#mast)** <br>\n",
    "**&ensp;[3. Drizzle PSF Models](#drizzles)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[3.1 Create Drizzle](#create)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[3.2 Stack Stars](#stack)** <br>\n",
    "**&ensp;[Conclusions](#onclusions)** <br>\n",
    "**&ensp;[About this Notebook](#about)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[Additional Resources](#additional)** <br>\n",
    "**&ensp;&nbsp;&nbsp;&nbsp;&nbsp;[Software Citations](#cite)** <br>\n",
    "\n",
    "**Acronyms:**\n",
    "- Point Spread Function (PSF)\n",
    "- Hubble Space Telescope (HST)\n",
    "- Wide Field Camera 3 (WFC3)\n",
    "- WFC3 InfraRed detector (IR)\n",
    "- WFC3 Ultraviolet and VISable detector (UVIS)\n",
    "- Advanced Camera for Surveys (ACS)\n",
    "- Wide Field and Planetary Camera 2 (WFPC2)\n",
    "- Mikulski Archive for Space Telescopes (MAST)\n",
    "- Flexible Image Transport System (FITS)\n",
    "- Instrument Science Report (ISR)\n",
    "- Full Width at Half Maximum (FWHM)\n",
    "- Signal-to-Noise (S/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d1d8205-0a56-4aca-bf73-6955a1c7078b",
   "metadata": {},
   "source": [
    "<a id=\"introduction\"></a>\n",
    "## 1. Introduction\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "**What is the Point Spread Function and why do we want to model it?**\n",
    "\n",
    "The purpose of a telescope is to collect, magnify, and record the intensity of light from astronomical sources. A telescope's ability to resolve the size of objects is limited by its spatial resolution, which is determined by its aperture and the wavelength of light being observed, such that every telescope is fundamentally limited by the diffraction nature of light. The minimum size of an image that the light can be focused onto is known as the Airy disk. In addition, the shape of the image will be altered by any optical aberrations in the telescope's system of mirrors, lenses, filters, and camera detectors. When observing unresolved point sources such as stars, these factors determine the size and shape of the focused image. This distribution of light as a function of position on the telescopes detectors is known as the Point Spread Function (PSF). More simply, the PSF describes how the light from an unresolved source is distributed across multiple pixels when it is recorded by the telescope's camera.\n",
    "\n",
    "While a well-focused telescope produces a very small image of point sources, the shape will be distorted and the light is effectively scattered to an infinite radius. Therefore, a model of the PSF is required to accurately measure astrometry (positions) and photometry (brightnesses) of objects in astronomical images. In the case of the HST WFC3 this is especially important because the images are undersampled, which means the telescope is capable of delivering higher resolution than the detector pixels can record. The PSF will vary depending on the telescopes focus, and also significantly with location on the WFC3 detector due to optical distortions in the light path that are designed to minimize light loss by reducing the number of optical reflections. For these reasons, the WFC3 team has produced a library of PSF models and tools for users on the [WFC3 PSF webpage](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf).\n",
    "\n",
    "Finally, astrometric and photometric measurements can generally be made either from individual calibrated exposures known as [FLTs and FLCs](https://hst-docs.stsci.edu/wfc3dhb/chapter-2-wfc3-data-structure/2-1-types-of-wfc3-files) for WFC3, or exposures that are combined into a mosaic through drizzling. In theory, modeling individual exposures is always preferable because these images are on the telescope's intrinsic pixel grid without any modification. However, this requires bright unresolved sources that have been observed with a relatively short exposure time, so that images are not strongly affected by cosmic rays. In addition, many astronomical sources are very faint, requiring a combination of multiple dithered exposures to obtain sufficient signal-to-noise (S/N) for analysis. For these reasons, we will demonstrate techniques for constructing PSF models for both individual calibrated exposures and drizzled mosaics.\n",
    "\n",
    "Users looking for a fundamental introduction to [aperture photometry](https://boyce-astro.org/aperture-photometry/) and [PSF photometry](https://boyce-astro.org/point-spread-function-psf/) may find this excellent series of [video tutorials](https://boyce-astro.org/videos/photometry/) by Boyce-Astro helpful, as well as this description of the [Airy Disk](https://www.edmundoptics.com/knowledge-center/application-notes/imaging/limitations-on-resolution-and-contrast-the-airy-disk/) by Edmund Optics. In addition, detailed information on WFC3 can be found on the [WFC3 Website](https://www.stsci.edu/hst/instrumentation/wfc3), [WFC3 PSF Webpage](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf), [WFC3 Instrument Handbook](https://hst-docs.stsci.edu/wfc3ihb), and [WFC3 Data Handbook](https://hst-docs.stsci.edu/wfc3dhb). In addition, the publication by [Anderson & King 2000](https://ui.adsabs.harvard.edu/abs/2000PASP..112.1360A/abstract) provides information about the effects of undersampling on the HST PSF. Finally, we note that the flux measured for each source will depend on the radial extent of the PSF model. Specifically, models truncated at a given radius will encompass a fraction of the total electrons from each star, which can be accounted for by using an encircled energy ratio correction. See [Section 9.1.8](https://hst-docs.stsci.edu/wfc3dhb/chapter-9-wfc3-data-analysis/9-1-photometry#id-9.1Photometry-9.1.89.1.8ApertureCorrections) of the WFC3 Data Handbook, [WFC3 ISR 2022-02](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2022/WFC3-ISR-2022-02.pdf), and the WFC3 [UVIS](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/photometric-calibration/uvis-encircled-energy) and [IR](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/photometric-calibration/ir-encircled-energy) encircled energy webpages for more information.\n",
    "\n",
    "We hope this notebook provides valuable information and examples for constructing HST WFC3 PSF models. While the focus in this notebook is for the WFC3 instrument, the majority of the code also applies to HST Advanced Camera for Surveys (ACS) observations, as well as archival Wide Field and Planetary Camera 2 (WFPC2) observations. See [ACS ISR 2006-01](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/acs/documentation/instrument-science-reports-isrs/_documents/isr0601.pdf) for more information on the ACS PSF models. Finally, users may find several helpful resources on the [WFC3 Software webpage](https://www.stsci.edu/hst/instrumentation/wfc3/software-tools)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88701fd4-f3c7-46ff-9824-f382c4a0b037",
   "metadata": {},
   "source": [
    "<a id=\"environment\"></a>\n",
    "### 1.1 Environment Setup\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "This Jupyter Notebook requires users to install the ```stenv``` Python environment maintained by STScI.\n",
    "\n",
    "Detailed installation instructions can be found on the stenv [readthedocs](https://stenv.readthedocs.io/en/latest/getting_started.html). For completeness, the terminal commands for MacOS with ARM chips are included below:\n",
    "\n",
    "<code>conda update conda\n",
    "conda env create --name stenv2024 --file https://github.com/spacetelescope/stenv/releases/download/2024.02.05/stenv-macOS-ARM64-py3.11-2024.02.05.yaml\n",
    "conda activate stenv2024\n",
    "conda install -c conda-forge gfortran\n",
    "conda clean -t\n",
    "jupyter lab</code>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><b>NOTE:</b> The following non-Python tools are needed to run some advanced features of the Jupyter Notebook. Specifically, a FORTRAN compiler and the hst1pass code are only required for PSF perturbation with individual exposures (FLCs). These tools are not required for stacking stars or creating PSF models for drizzled images (DRZs). <b>The code will automatically compile hst1pass. If this step encounters errors, then manual instructions for installing FORTRAN and hst1pass are included below:</b> </div>\n",
    "\n",
    " - [**hst1pass**](https://www.stsci.edu/~jayander/HST1PASS/) - FORTRAN code for performing photometry and PSF perturbation as documented in [WFC3 ISR 2022-05](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2022/WFC3-ISR-2022-05.pdf).\n",
    "\n",
    "FORTRAN can be installed by typing the command ```conda install -c conda-forge gfortran``` into a terminal. If the installation is successful, then typing ```gfortran --version``` should return a message similar to the following: ```GNU Fortran (GCC) 13.2.0 Copyright (C) 2023 Free Software Foundation, Inc.```.\n",
    "\n",
    "Next, to install hst1pass proceed to [https://www.stsci.edu/~jayander/HST1PASS/CODE/hst1pass/](https://www.stsci.edu/~jayander/HST1PASS/CODE/hst1pass/) and download the latest version of ```hst1pass.F```, e.g. ```hst1pass.2023.11.07_v1f.F```, or use the version in this notebook's repository. Save this to your computer and compile the code using the command: ```gfortran hst1pass.2023.11.07_v1f.F -o hst1pass.e```."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "922995ca-d521-4b8e-af77-09ea5e054648",
   "metadata": {},
   "source": [
    "<a id=\"import\"></a>\n",
    "### 1.2 Import Packages\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "The following Python packages are required to run the Jupyter Notebook:\n",
    " - [**os**](https://docs.python.org/3/library/os.html) - change and make directories\n",
    " - [**glob**](https://docs.python.org/3/library/glob.html) - gather lists of filenames\n",
    " - [**shutil**](https://docs.python.org/3/library/shutil.html) - copy files between folders\n",
    " - [**tarfile**](https://docs.python.org/3/library/tarfile.html) - read and write compressed tar files\n",
    " - [**urllib**](https://docs.python.org/3/library/urllib.html) - download files hosted online\n",
    " - [**ipython**](https://ipython.org) - cell formatting and interactives\n",
    "     - [display.clear_output](http://ipython.org/ipython-doc/dev/api/generated/IPython.display.html#IPython.display.clear_output) - clear text from cells\n",
    " - [**numpy**](https://numpy.org) - math and array functions\n",
    "     - [percentile](https://numpy.org/doc/stable/reference/generated/numpy.percentile.html) - lognorm percentile clipping\n",
    " - [**matplotlib**](https://matplotlib.org) - create graphics\n",
    "     - [pyplot](https://matplotlib.org/stable/tutorials/pyplot.html) - make figures and graphics\n",
    "     - [colors.LogNorm](https://matplotlib.org/stable/api/_as_gen/matplotlib.colors.LogNorm.html) - display image normalization\n",
    "     - [patches.Rectangle](https://matplotlib.org/stable/api/_as_gen/matplotlib.patches.Rectangle.html) - overlay rectangular regions\n",
    " - [**astroquery**](https://astroquery.readthedocs.io/en/latest/) - download astronomical data\n",
    "     - [mast.Observations](https://astroquery.readthedocs.io/en/latest/mast/mast_obsquery.html) - query MAST database\n",
    " - [**astropy**](https://www.astropy.org) - model fitting and file handling\n",
    "     - [io.fits](https://docs.astropy.org/en/stable/io/fits/) - import FITS files\n",
    "     - [table.QTable](https://docs.astropy.org/en/stable/api/astropy.table.QTable.html) - tables with physical units\n",
    "     - [modeling.models](https://docs.astropy.org/en/stable/modeling/models.html) - Gaussian and Moffat models\n",
    "     - [modeling.fitting](https://docs.astropy.org/en/stable/modeling/fitting.html) - fitting methods for models\n",
    "     - [visualization.simple_norm](https://docs.astropy.org/en/stable/api/astropy.visualization.mpl_normalize.simple_norm.html) - display image normalization\n",
    "     - [stats.sigma_clipped_stats](https://docs.astropy.org/en/stable/api/astropy.stats.sigma_clipped_stats.html) - mean, median, and standard deviation\n",
    "     - [stats.SigmaClip](https://docs.astropy.org/en/stable/api/astropy.stats.SigmaClip.html) - sigma-clip provided data\n",
    " - [**photutils**](https://photutils.readthedocs.io/en/stable/index.html) - aperture and PSF photometry tools\n",
    "     - [detection.DAOStarFinder](https://photutils.readthedocs.io/en/stable/api/photutils.detection.DAOStarFinder.html#photutils.detection.DAOStarFinder) - finding stars in images\n",
    "     - [detection.find_peaks](https://photutils.readthedocs.io/en/stable/api/photutils.detection.find_peaks.html#photutils.detection.find_peaks) - finding peaks in images\n",
    "     - [photutils.centroids](https://photutils.readthedocs.io/en/stable/centroids.html) - finding stellar centroids\n",
    "     - [photutils.aperture](https://photutils.readthedocs.io/en/stable/aperture.html) - aperture photometry tools\n",
    "       - [aperture_photometry](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.aperture_photometry.html#photutils.aperture.aperture_photometry) - measure aperture photometry\n",
    "       - [CircularAperture](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAperture.html) - measure photometry in circle\n",
    "       - [CircularAnnulus](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAnnulus.html) - measure photometry in annulus\n",
    "       - [ApertureStats](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.ApertureStats.html) - calculate aperture statistics\n",
    "     - [photutils.psf](https://photutils.readthedocs.io/en/stable/psf.html) - PSF photometry tools\n",
    "         - [PSFPhotometry](https://photutils.readthedocs.io/en/stable/api/photutils.psf.PSFPhotometry.html#photutils.psf.PSFPhotometry) - perform PSF photometry\n",
    "         - [SourceGrouper](https://photutils.readthedocs.io/en/stable/api/photutils.psf.SourceGrouper.html#photutils.psf.SourceGrouper) - group sources for fitting\n",
    "         - [IntegratedGaussianPRF](https://photutils.readthedocs.io/en/stable/api/photutils.psf.IntegratedGaussianPRF.html) - Gaussian model for PSF fitter\n",
    "         - [GriddedPSFModel](https://photutils.readthedocs.io/en/stable/api/photutils.psf.GriddedPSFModel.html#photutils.psf.GriddedPSFModel) - enables spatially-dependent PSF grids\n",
    "         - [STDPSFGrid](https://photutils.readthedocs.io/en/stable/api/photutils.psf.STDPSFGrid.html#photutils.psf.STDPSFGrid) - reads STScI standard PSF models into a grid\n",
    "         - [stdpsf_reader](https://photutils.readthedocs.io/en/stable/api/photutils.psf.stdpsf_reader.html#photutils.psf.stdpsf_reader) - reads STScI standard PSF format models\n",
    "         - [FittableImageModel](https://photutils.readthedocs.io/en/stable/api/photutils.psf.FittableImageModel.html#photutils.psf.FittableImageModel) - convert stacked stars to fitable model\n",
    " - [**scipy**](https://scipy.org) - mathematical interpolation and shift functions\n",
    "     - [ndimage.shift](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.shift.html) - fractional array shifts\n",
    " - [**drizzlepac**](https://www.stsci.edu/scientific-community/software/drizzlepac) - combine HST images into mosaics\n",
    "     - [astrodrizzle](https://drizzlepac.readthedocs.io/en/deployment/astrodrizzle.html) - combine exposures by drizzling\n",
    "     - [tweakreg](https://drizzlepac.readthedocs.io/en/deployment/tweakreg.html) - align exposures to a common WCS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acfcc04d-596f-41a6-a2fa-bce20ee2cfbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import shutil\n",
    "import tarfile\n",
    "from IPython.display import clear_output\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "from matplotlib.patches import Rectangle\n",
    "from astroquery.mast import Observations\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table, QTable\n",
    "from astropy.visualization import simple_norm\n",
    "from astropy.stats import sigma_clipped_stats, SigmaClip\n",
    "from photutils.detection import DAOStarFinder\n",
    "from photutils.aperture import (ApertureStats, CircularAnnulus,\n",
    "                                CircularAperture, aperture_photometry)\n",
    "from photutils.psf import (FittableImageModel, GriddedPSFModel, IntegratedGaussianPRF, \n",
    "                           PSFPhotometry, SourceGrouper, STDPSFGrid)\n",
    "from drizzlepac import tweakreg, astrodrizzle\n",
    "\n",
    "# Custom functions written in psf_utilities.py.\n",
    "# Contains functions for making various figures and maniputlating science images.\n",
    "from psf_utilities import save_figure, setup_matplotlib, create_mask, plot_apertures, \\\n",
    "    plot_psf_results, download_psf_model, make_cutouts, stack_cutouts, plot_cutout_grid\n",
    "\n",
    "# Custom functions written in mast_api_psf.py.\n",
    "# Contains functions for downloading image cutouts from the MAST PSF library server.\n",
    "import mast_api_psf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "622da46b-bf0f-473f-88f3-3669dba40783",
   "metadata": {},
   "source": [
    "Note that we can access the helpful documentation string for any of the custom-defined functions using help(). For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7999490-ece9-4728-af45-664194ebf04c",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(save_figure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "316b41c6-fd8d-4218-a548-0cdbb8b84500",
   "metadata": {},
   "source": [
    "Set the main code, data, PSF, plot, AstroDrizzle, and TweakReg directories. Set the backend configuration to 'retina', which improves the resolution of figures in notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69b74d05-5a9d-4101-88ed-0cd0ee5c635a",
   "metadata": {},
   "outputs": [],
   "source": [
    "code_dir = os.getcwd()+'/'\n",
    "data_dir = code_dir+'data/'\n",
    "psfs_dir = code_dir+'psfs/'\n",
    "plot_dir = code_dir+'plot/'\n",
    "driz_dir = code_dir+'driz/'\n",
    "tweak_dir = code_dir+'tweak/'\n",
    "\n",
    "os.chdir(code_dir)\n",
    "\n",
    "print('\\nThe code_dir is:', code_dir)\n",
    "\n",
    "for folder in ['data', 'psfs', 'plot', 'driz', 'tweak']:\n",
    "    if not os.path.exists(folder):\n",
    "        os.mkdir(folder)\n",
    "    else:\n",
    "        print(f'{folder} directory exists.')\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c186b65-9c1e-4f08-ac94-51f43d467c35",
   "metadata": {},
   "source": [
    "<a id=\"data\"></a>\n",
    "### 1.3 Download Data\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "In this example, we will be using HST WFC3 imaging of the globular cluster Omega Centauri (NGC 5139). Specifically, utilizing exposures that are not centered on the cluster core so there is a moderate density of bright stars. This source was chosen because it serves as a frequent target for both science and calibration exposures, and has been observed with every UVIS and IR filter onboard WFC3. A complete list of WFC3 filters is available in the Data Handbook for the [UVIS](https://hst-docs.stsci.edu/wfc3ihb/chapter-6-uvis-imaging-with-wfc3/6-5-uvis-spectral-elements) and [IR](https://hst-docs.stsci.edu/wfc3ihb/chapter-7-ir-imaging-with-wfc3/7-5-ir-spectral-elements) channels. First, we use astroquery to retrieve the calibrated exposure from MAST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "808a45b7-d743-4f10-b027-cb22472de4d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "hst_wfc3_f606w_file = 'id8048dyq_flc.fits'\n",
    "\n",
    "Observations.download_file('mast:HST/product/'+hst_wfc3_f606w_file, local_path=data_dir+hst_wfc3_f606w_file, cache=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0860c860-a5f3-4559-86dc-7da555e920f3",
   "metadata": {},
   "source": [
    "Next, let's display the image of Omega Centauri and add a rectangle showing the region of the image that we will focus on in this analysis. We select a small region of the image that contains isolated and blended stars to highlight key concepts in aperture and PSF photometry. The WFC3/UVIS pixel scale is 0.03962 arcseconds per pixel so selecting a region of 75 pixels corresponds to approximately 3 arcseconds. The UVIS detector has two chips that can be accessed in the FITS files from extensions ```['SCI', 1]``` for UVIS2 and ```['SCI', 2]``` for UVIS1. The lower-left corner of the detector is the location (0, 0) pixels, so we use the ```origin = 'lower'``` option. Finally, we display the image on a logarithmic scale to view both very faint and bright stars simultaneously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bd92a08-2877-40ac-b3ea-a5224a01e05e",
   "metadata": {},
   "outputs": [],
   "source": [
    "with fits.open(data_dir+hst_wfc3_f606w_file) as hdul:\n",
    "    sci_data = hdul['SCI', 1].data\n",
    "print('The image dimensions are:', np.shape(sci_data))\n",
    "\n",
    "my_figure_size, my_fontsize = setup_matplotlib('notebook', 1.2)\n",
    "figure, axes = plt.subplots(1, 1, figsize=my_figure_size)\n",
    "cutout_size = 75 # 75 WFC3 UVIS pixels = 3 arcseconds.\n",
    "\n",
    "axes.add_patch(Rectangle((1818-cutout_size/2, 986-cutout_size/2), width=cutout_size, height=cutout_size, fill=None, color='w'))\n",
    "norm = LogNorm(vmin=np.percentile(sci_data, 10.0), vmax=np.percentile(sci_data, 99.9))\n",
    "axes.imshow(sci_data, cmap='gray', origin='lower', aspect='equal', norm=norm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e46e40d0-cfdf-4426-9737-2ce89bc8ab47",
   "metadata": {},
   "source": [
    "<a id=\"aperture\"></a>\n",
    "### 1.4 Aperture Photometry\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "First, we can perform simple aperture photometry on the image. This consists of detecting stars in the image, assigning them an aperture within which to measure the enclosed flux, and creating a background annulus that is used to subtract the background level. In ground-based astronomy this is commonly referred to as the \"sky\" level. In this example, we create a mask so that we only find and measure the properties of stars that are within our selected box region.\n",
    "\n",
    "The cell below first calculates the mean, median, and standard deviation for the background, initializes the star finder function with a threshold of 500 electrons and a FWHM of three pixels, and then creates a mask so that the star finding algorithm will only search within our zoomed-in box region. Next, the `DAOStarFinder` function searches for peaks and returns a list of sources. We then assign to the position of each source a circular annulus designed to measure the flux of the star, and a background annulus that measures the local background flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9b4311-e2ab-4da0-945b-0cfe07a04c99",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean, median, std = sigma_clipped_stats(sci_data, sigma=3.0)\n",
    "daofind = DAOStarFinder(threshold=500, fwhm=3.0, roundhi=0.5)\n",
    "mask = create_mask(sci_data, cutout_size, 1818, 986)\n",
    "sources = daofind(sci_data - median, mask=mask)\n",
    "positions = np.transpose((sources['xcentroid'], sources['ycentroid']))\n",
    "apertures_stellar = CircularAperture(positions, r=5.0)\n",
    "apertures_annulus = CircularAnnulus(positions, r_in=9, r_out=12)\n",
    "phot_table = aperture_photometry(sci_data, apertures_stellar)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78576b80-1a6c-456d-a37b-39ff4e727c75",
   "metadata": {},
   "source": [
    "Next, we use the aperture and annulus measurement for each star to calculate the sigma-clipped background flux, and subtract it to obtain only the flux of each star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5195b262-92a2-4cc1-b41a-e19d2398c63d",
   "metadata": {},
   "outputs": [],
   "source": [
    "aperture_stats = ApertureStats(sci_data, apertures_stellar, sigma_clip=None)\n",
    "background_stats = ApertureStats(sci_data, apertures_annulus, sigma_clip=SigmaClip(sigma=3.0, maxiters=10))\n",
    "total_background = background_stats.median * aperture_stats.sum_aper_area.value # area of annulus.\n",
    "flux_background_sub = aperture_stats.sum - total_background\n",
    "phot_table['total_background'] = total_background\n",
    "phot_table['flux_background_sub'] = flux_background_sub\n",
    "\n",
    "for col in phot_table.colnames:  \n",
    "    if col not in ('id', 'npix'):\n",
    "        phot_table[col].info.format = '%.2f'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e637c0d-33de-4860-b891-d27f6deb9ac0",
   "metadata": {},
   "source": [
    "Finally, we display the flux measurements in a table and create a figure showing the science data, stellar apertures (green circles) and background annuli (white circles)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8cea058-fb2a-4a85-90a8-90c563d34fed",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('\\n', phot_table, '\\n')\n",
    "\n",
    "plot_apertures(data=sci_data, \n",
    "               xcenter=1818, \n",
    "               ycenter=986, \n",
    "               cutout_size=cutout_size, \n",
    "               apertures_stellar=apertures_stellar, \n",
    "               apertures_annulus=apertures_annulus)\n",
    "\n",
    "save_figure(figure, plot_dir+'fig1_aperture_photometry.pdf', True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "325396c3-178f-4b7d-b63c-adf6e22ec4c1",
   "metadata": {},
   "source": [
    "<a id=\"mainc\"></a>\n",
    "## 2. Exposure PSF Models\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "As described in the [Introduction](#introduction), we may perform PSF modeling, photometry, and astrometry on either individual exposures or mosaics that combine exposures through drizzling. The sections below contain several different workflows with various options depending on your science goals and the available data. The flowchart below is meant as a general guide for deciding which sections are most relevant for your project. We highlight that several different paths may be applicable to a single dataset and users are encouraged to experiment with the options below.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><b>NOTE:</b> Users that are creating PSF models for drizzled images and not single exposures can proceed to the section on Drizzle PSF Models. </div>\n",
    "\n",
    "![This flowchart provides users with a general guide for deciding which PSF modeling technique is the most relevant based on their dataset and science goals. The top of the chart asks users if they are modeling at the exposure or drizzle level. The drizzle workflow forks to the right, and asks users how many bright, isolated stars are in their image. If there are fewer than 10 stars, they can fit an analytical profile or download custom tools as described below. If they have 10 or more stars, they can extract and stack stars as described in Section 3. If users are working at the exposure level, then the diagram forks to the left. If the user's goal is focused on high-accuracy photometry and astrometry of stars, then they follow the left fork. If their goals are to model the extended wings and diffraction spikes of the PSF, they follow the fork right. In following the left fork for precision photometry and astrometry, users are asked if an empirical model from the WFC3 team exists. If it does, then they can run hst1pass and examine the residuals. If there is not an empirical model, they can use the model for the closest filter and perturb it with hst1pass if there are at least 10 stars. In following the right fork for modeling extended wing emission and diffraction spikes, users are asked if they have more or less than 10 stars in their image. If they have fewer than 10, they can use the MAST cutout service to retrieve and stack stars from other observations. If they have at least 10 stars, then they can stack stars directly from their observations. This flowchart is a general guide, and we note that several different paths may be applicable to a single dataset. Users are encouraged to experiment with the options below.](workflow.png)\n",
    "\n",
    "Note that the left fork under the drizzle level branch is an advanced workflow that involves drizzling the empirical models provided by the WFC3 team into FLC files. This requires installing the [wfc3_photometry](https://github.com/spacetelescope/wfc3_photometry) package, running the `make_model_star_image()` function within `PSFUtils.py`, and then stacking the drizzled models. This functionality is under development and not yet officially supported."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e16eb70c-277e-4a6a-a285-a2e987760873",
   "metadata": {},
   "source": [
    "<a id=\"analytical\"></a>\n",
    "### 2.1 Analytical PSF Models\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "The example above demonstrates how to perform basic aperture photometry of sources. However, we can see that two of the stars overlap in the image and their fluxes are blended. Therefore, aperture photometry cannot be used reliably for crowded fields, and requires corrections when using multiple filters as the apertures do not account for changing resolution with wavelength. Without any knowledge of the intrinsic shape of the PSF, we can first fit purely analytical functions such as a Gaussian, Moffat, Lorentzian, or Voigt profile to each star. While the majority of our functions are specified in the ```psf_utilities.py``` module, we explicitly define two common photometry functions here so users can easily experiment and change common options depending on their data and goals.\n",
    "\n",
    "Critically, we note that the size of the `fit_shape` variable (in pixels) determines the number of pixels used to constrain the fit. If the value is small, the core pixels of the source will dominate the fit, while larger values will weight more wing emission when determining the best fit model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a79caaa4-196a-46bd-934d-78eab9c66406",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_photometry(sci_data, sources, psf_model):\n",
    "    \n",
    "    \"\"\"\n",
    "    A function that performs PSF photometry given a science image, \n",
    "    a list of sources, and a PSF model in the PSFPhotometry format. \n",
    "    See the photutils documentation below for additional examples:\n",
    "    \n",
    "    https://photutils.readthedocs.io/en/stable/psf.html\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    sci_data : np.ndarray\n",
    "        The np.ndarray containing the array of science data.\n",
    "    sources : astropy.table.table.QTable\n",
    "        The astropy table of sources from the daofind function.\n",
    "    psf_model : photutils.psf.GriddedPSFModel\n",
    "        An astropy compatible model PSF or grid of PSFs.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    phot : astropy.table.table.QTable\n",
    "        An astropy quantity table with measured source photometry.\n",
    "    psfphot : photutils.psf.PSFPhotometry\n",
    "        Contains results of PSF fitting for use with other functions.\n",
    "    \"\"\"\n",
    "\n",
    "    # Perform PSF fitting using an 9x9 pixel region.    \n",
    "    fit_shape = (9, 9)\n",
    "    grouper = SourceGrouper(min_separation=20)\n",
    "    finder = DAOStarFinder(threshold=500, fwhm=3.0, roundhi=0.5)\n",
    "    psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder, aperture_radius=5, grouper=grouper)\n",
    "    \n",
    "    init_params = QTable()\n",
    "    init_params['x'] = sources['xcentroid']\n",
    "    init_params['y'] = sources['ycentroid']\n",
    "    phot = psfphot(sci_data, init_params=init_params)\n",
    "    phot['x_fit'].info.format = '.2f'\n",
    "    phot['y_fit'].info.format = '.2f'\n",
    "    phot['flux_fit'].info.format = '.2f'\n",
    "    phot['qfit'].info.format = '.3f'\n",
    "    print('\\033[1m'+'\\nPSF Photometry (before subtraction):\\n'+'\\033[0m')\n",
    "    print(phot[('id', 'x_fit', 'y_fit', 'flux_fit', 'qfit')])\n",
    "    print('')\n",
    "    \n",
    "    return phot, psfphot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adba3a96-ff9b-4b32-8c72-600514e3ea70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_residuals(resid, sources):\n",
    "    \n",
    "    \"\"\"\n",
    "    A function that calculates the residual flux at the \n",
    "    locations of stars in a given PSF-subtracted image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    resid : np.ndarray\n",
    "        A residual image created by psfphot.make_residual_image().\n",
    "    sources : astropy.table.table.QTable\n",
    "        The astropy table of sources from the daofind function.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    phot_table : astropy.table.table.QTable\n",
    "        Prints an astropy table with measured source photometry.\n",
    "    \"\"\"\n",
    "\n",
    "    # Use the PSF fitting residuals to calculate remaining fluxes.\n",
    "    phot_table = aperture_photometry(resid, apertures_stellar)\n",
    "    aperture_stats = ApertureStats(resid, apertures_stellar, sigma_clip=None)\n",
    "    background_stats = ApertureStats(resid, apertures_annulus, sigma_clip=SigmaClip(sigma=3.0, maxiters=10))\n",
    "    total_background = background_stats.median * aperture_stats.sum_aper_area.value # area of annulus.\n",
    "    flux_background_sub = aperture_stats.sum - total_background\n",
    "    phot_table['total_background'] = total_background\n",
    "    phot_table['flux_background_sub'] = flux_background_sub\n",
    "    \n",
    "    for col in phot_table.colnames:  \n",
    "        if col not in ('id', 'npix'):\n",
    "            phot_table[col].info.format = '%.2f'\n",
    "\n",
    "    print('\\033[1m'+'\\nAperture Photometry (after subtraction):\\n'+'\\033[0m')\n",
    "    \n",
    "    return print(phot_table, '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c34b63-871b-4444-98e7-010dcb59574c",
   "metadata": {},
   "source": [
    "In this first analytical example, we fit Gaussian profiles to the four stars in our subimage and display the data, model, and residuals. We provide an initial guess for the width (sigma) using an approximate full width at half maximum (FWHM) of 2 pixels, where $\\sigma$ = FWHM/2.355. Users are encouraged to experiment with this value to see how the resulting fit and residuals change. First, the `calculate_photometry` function uses the provided PSF model with the photutils `PSFPhotometry` function to fit the model to the stars. In essence, the flux of the normalized PSF model is scaled up to match the flux of each star. In this case, neighboring stars are fit simultaneously. The `make_residual_image` function then subtracts the best fit model so the flux residuals can be examined. Finally, we call the `calculate_residuals` function with the exact same stellar and background photometry apertures used in [Section 1.4](#aperture) in order to measure the residual fluxes. A smaller residual flux generally implies a better fit, but the residuals should be inspected as over and undersubtractions from different portions of the star can also yield a small net residual flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb319d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf_model = IntegratedGaussianPRF(flux=1000, sigma=2.0/2.355)\n",
    "photometry, psfphot = calculate_photometry(sci_data, sources, psf_model)\n",
    "resid = psfphot.make_residual_image(sci_data, (19, 19))\n",
    "figure = plot_psf_results(sci_data, resid, 1818, 986, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig2_gaussian_model.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04e4cef7-8f74-4008-9d44-bd04d70bab82",
   "metadata": {},
   "source": [
    "In this case, the Gaussian fit leaves significant residual flux in the stellar wings, and has oversubtracted the emission in the core. The three panels above are shown with the same minimum and maximum flux values, so the residuals can be directly compared with the data and model. In the next section, we improve this using models of HST data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "692251c0-7004-414c-a82b-d6937ac44f57",
   "metadata": {},
   "source": [
    "<a id=\"empirical\"></a>\n",
    "### 2.2 Empirical PSF Models\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "The example above demonstrates how a Gaussian profile can determine reasonable centroids for stellar sources, but they fail to accurately model the radial light distribution. While it is possible to use more advanced analytical forms such as [Moffat profiles](https://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Moffat2D.html), the overall shape can be more accurately represented by models constructed from bright stars in the field of view. The examples below demonstrate how to use pre-generated PSF models provided by the WFC3 team, as well as how to extract and stack stars from your observations. As described below, the pre-generated models are excellent for precision astrometry and photometry of unresolved sources, while stacking stars can be more beneficial for accurately modeling the PSF wings and diffraction spikes.\n",
    "\n",
    "First, we display the pre-generated PSF models provided by the WFC3 team, which are stored in a standard format known as \"STDPSF_WFC3\". The FITS file for each filter contains a grid of models that represents the PSF at regularly spaced intervals of a few hundred pixels across the detector, as shown below. The models are oversampled by a factor of 4x to account for the fact that the morphology of the PSF changes depending on where a star falls within a given pixel. These models are tapered to zero flux at large radii to prevent discontinuities in image subtraction, at the expense of PSF wing accuracy. The specific details of these models are provided on the [WFC3 PSF webpage](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf), as well as in [WFC3 ISR 2016-12](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2016/WFC3-2016-12.pdf). The models are retrieved from [https://www.stsci.edu/~jayander/HST1PASS/LIB/PSFs/STDPSFs/](https://www.stsci.edu/~jayander/HST1PASS/LIB/PSFs/STDPSFs/), and we save them to the main directory for access by hst1pass.\n",
    "\n",
    "There are also \"focus-diverse\" models for commonly used filters with an additional dimension containing the same grids at different focus intervals. The best fit model can be determined experimentally, or by running the hst1pass code to measure the focus. However, if there are a sufficient number of stars to measure the focus then hst1pass can \"perturb\" the pre-generated PSF models to best match the data. As this approach is generalized for any filter we do not use the focus-diverse models in this notebook. Users can find information in [WFC3 ISR 2018-14](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2018/WFC3-2018-14.pdf). In addition, similar focus-diverse models exist for [ACS](http://acspsf.stsci.edu) as described in [ACS ISR 2023-06](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/acs/documentation/instrument-science-reports-isrs/_documents/isr2306.pdf).\n",
    "\n",
    "Finally, in the examples below we complete a single iteration in the PSF fitting, but it is possible to recover faint blended stars iteratively using [IterativePSFPhotometry](https://photutils.readthedocs.io/en/latest/api/photutils.psf.IterativePSFPhotometry.html#photutils.psf.IterativePSFPhotometry)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a74aa76a-1241-4d19-9aa8-a138f412c040",
   "metadata": {},
   "outputs": [],
   "source": [
    "hst_wfc3_f606w_psf_model = download_psf_model(psfs_dir, 'WFC3UV', 'F606W')\n",
    "hdul = fits.open(psfs_dir+hst_wfc3_f606w_psf_model, ignore_missing_end=True)\n",
    "setup_matplotlib('notebook', 1.4)\n",
    "psfgrid = STDPSFGrid(psfs_dir+hst_wfc3_f606w_psf_model).plot_grid()\n",
    "save_figure(psfgrid, plot_dir+'fig3_empirical_grid.pdf', True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "563bde21-c750-4c89-8981-e5a5816d480f",
   "metadata": {},
   "source": [
    "<a id=\"empirical-1\"></a>\n",
    "#### 2.2.1 Default Models\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "First, we can simply use the pre-generated PSF models provided by the WFC3 team. In many cases, these models will be sufficient for measuring very accurate astrometry and photometry of stellar sources in individual exposures (FLTs/FLCs) without further modifications. We utilize functions within the Astropy-affiliated [photutils](https://photutils.readthedocs.io/en/stable/) package to read the models into a recognized format for the PSF fitter. As noted in the packages [documentation](https://photutils.readthedocs.io/en/stable/api/photutils.psf.stdpsf_reader.html), for ACS/WFC and WFC3/UVIS the detector value should be \"1\" for WFC2, UVIS2 (sci, 1), and \"2\" for WFC1, UVIS1 (sci, 2).\n",
    "\n",
    "In the examples below we only utilize the `['SCI', 1]` extension of the science image, so the hst1pass and array image coordinates will match. When analyzing `['SCI', 2]` (`k = 1` in hst1pass catalogs) users must subtract 2048 from the hst1pass `y` coordinates so they match the array coordinates. This is because hst1pass detects over the full 4096x2048 dimensions of the detector, but each chip is read in as a 2048x2048 array.\n",
    "\n",
    "We now read in the model and use it with the PSF fitter to calculate the best fit model and residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee854918-cb7c-4668-bbc1-abc925485392",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_uvis1 = GriddedPSFModel.read(psfs_dir+hst_wfc3_f606w_psf_model, format='stdpsf', detector_id=2)\n",
    "# model_uvis2 = GriddedPSFModel.read(psfs_dir+hst_wfc3_f606w_psf_model, format='stdpsf', detector_id=1) # Would be used for UVIS2.\n",
    "photometry, psfphot = calculate_photometry(sci_data, sources, psf_model=model_uvis1)\n",
    "resid = psfphot.make_residual_image(sci_data, (9, 9))\n",
    "figure = plot_psf_results(sci_data, resid, 1818, 986, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig4_empirical_default.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b09ed551-fc60-40a7-8038-109e8938b480",
   "metadata": {},
   "source": [
    "<a id=\"empirical-2\"></a>\n",
    "#### 2.2.2 Focus Perturbation\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "In the case above, the pre-generated PSF model subtraction has substantial flux remaining in the residual image. We can use the perturbation feature within hst1pass to modify the shape of the PSF to account for the focus of these specific observations, as well as other effects such as minor modulations due to guiding accuracy. As described in Section 3.3 of [WFC3 ISR 2022-05](https://www.stsci.edu/~jayander/JWST1PASS/DOCs/ISR_HST1PASS.pdf), the routine has a `PERT` option that can modify the library PSF if there are at least 10 relatively bright, isolated stars in the field. If `PERT=1`, then hst1pass constructs a single PSF that applies across the entire image. This often provides the best results, but `PERT=2` can be specified to generate an NxN array of PSFs spaced evenly across the detector. This requires ample bright, isolated stars across the detector. In general, it is recommended that users begin with the unaltered PSFs before utilizing the focus or spatially-dependent perturbed PSFs. The perturbed models have a small effect on astrometry and moderate on photometry.\n",
    "\n",
    "The commands below assume that you have followed the instructions for installing hst1pass given in the [Environment Setup](#environment) section. The executable can be in the same directory as assumed below, or you can provide the full path to the location of the executable. We now attempt to compile and call hst1pass with `PERT=1` to construct a single perturbed PSF model that applies at all locations on the detector. We make a copy of the model because it will be overwritten by hst1pass, but call the original file for compatibility with `GriddedPSFModel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62acd32-175d-4d0b-8946-d02edfe7ab91",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(code_dir+'hst1pass.e'):\n",
    "    os.chdir(code_dir)\n",
    "    !gfortran hst1pass.2023.11.07_v1f.F -o hst1pass.e\n",
    "    print('The hst1pass.e has been compiled.')\n",
    "else:\n",
    "    print('The executable already exists.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e6768d0-9f21-486d-b36b-db6d9e027bd5",
   "metadata": {},
   "source": [
    "We now call hst1pass and specify the minimum isolation distance in pixels from the next brightest pixel (`HMIN`), the minimum source flux in counts in the central pixel (`FMIN`), and the maximum source flux in the central pixel (`PMAX`). These arguments are followed by the `OUT` output columns of the resulting catalog, followed by the PSF model to use, and the FLC file to be analyzed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d6b9c86-1a77-4f50-b4c2-a8626e90543e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(data_dir)\n",
    "!./../hst1pass.e HMIN=5 FMIN=1000.0 PMAX=66000.0 PERT=1 OUT=xympqks PSF=./../psfs/STDPSF_WFC3UV_F606W.fits id8048dyq_flc.fits\n",
    "shutil.copy('id8048dyq_psf.fits', psfs_dir+'id8048dyq_psf_pert1.fits')\n",
    "shutil.copy('id8048dyq_flc.xympqks', 'id8048dyq_flc_pert1.xympqks')\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68559dfb-05d4-4744-9fda-d433d92b1fe1",
   "metadata": {},
   "source": [
    "Read in the perturbed model and examine the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e3a176a-ec16-4ee6-8ea4-27ea99736dc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf_model_perturbed_1 = data_dir+'id8048dyq_psf.fits'\n",
    "model_uvis1 = GriddedPSFModel.read(psf_model_perturbed_1, format='stdpsf', detector_id=2)\n",
    "# model_uvis2 = GriddedPSFModel.read(psf_model_perturbed_1, format='stdpsf', detector_id=1) # Would be used for UVIS2.\n",
    "photometry, psfphot = calculate_photometry(sci_data, sources, psf_model=model_uvis1)\n",
    "resid = psfphot.make_residual_image(sci_data, (9, 9))\n",
    "figure = plot_psf_results(sci_data, resid, 1818, 986, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig5_empirical_focus.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbfd9a1a-fde7-4458-9210-da9075665ab0",
   "metadata": {},
   "source": [
    "We see that the residual flux is smaller than when using the default models, and visually the subtraction is improved by accounting for the focus of these observations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69c8d145-f94a-49e0-a315-0ceb636e6c28",
   "metadata": {},
   "source": [
    "<a id=\"empirical-3\"></a>\n",
    "#### 2.2.3 Spatial Perturbation\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "In general, using hst1pass with `PERT=1` to construct a single, focus-dependent PSF model that applies across the detector will provide the best results. However, for specialized use cases with many bright, isolated stars that are located near one location on the detector, or dispersed across the entire detector, hst1pass can create a spatially-dependent PSF model in the same NxN array format as the STDPSF_WFC3 standard models. We now use hst1pass with `PERT=2` to generate a grid of spatially-dependent PSF models optimized for the focus of this exposure at each location on the detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30ec74bf-ee16-4829-82e8-eb56a3a4776b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(data_dir)\n",
    "!./../hst1pass.e HMIN=5 FMIN=1000.0 PMAX=66000.0 PERT=2 OUT=xympqks PSF=./../psfs/STDPSF_WFC3UV_F606W.fits id8048dyq_flc.fits\n",
    "shutil.copy('id8048dyq_psf.fits', psfs_dir+'id8048dyq_psf_pert2.fits')\n",
    "shutil.copy('id8048dyq_flc.xympqks', 'id8048dyq_flc_pert2.xympqks')\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b97646c-2d6a-4f61-ad12-8c3c915a368b",
   "metadata": {},
   "source": [
    "Read in the perturbed model and examine the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3cd48ed-a0ad-4e58-ad8b-4f297d65bbd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf_model_perturbed_2 = data_dir+'id8048dyq_psf.fits'\n",
    "model_uvis1 = GriddedPSFModel.read(psf_model_perturbed_2, format='stdpsf', detector_id=2)\n",
    "# model_uvis2 = GriddedPSFModel.read(psf_model_perturbed_2, format='stdpsf', detector_id=1) # Would be used for UVIS 2.\n",
    "photometry, psfphot = calculate_photometry(sci_data, sources, psf_model=model_uvis1)\n",
    "resid = psfphot.make_residual_image(sci_data, (9, 9))\n",
    "figure = plot_psf_results(sci_data, resid, 1818, 986, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig6_empirical_spatial.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14ba694b-ebe4-432b-9856-ff4ea9c9e060",
   "metadata": {},
   "source": [
    "Interestingly, the residual flux is minutely larger when using the spatially-dependent model with this dataset, which can be seen by comparing the `flux_background_sub` values for each star in the above tables. This is primarily due to the fact that there are an insufficient number of bright, isolated stars to accurately constrain the PSF model at each location on the detector. Thus, the best solution is produced by using `PERT=1` with hst1pass.\n",
    "\n",
    "In summary, the default models will work well for many science cases. If there are substantial residuals, then hst1pass can apply a perturbation that slightly changes the PSF model shape to compensate for different focus positions (`PERT=1`). It is also possible to account for both the focus and spatial variation of the PSF across the detector (`PERT=2`). However, changes in the PSF due to focus are generally more significant than the spatial variations unless sources are close to the edges of the detector, and accounting for both effects requires a very large number of bright, isolated stars.\n",
    "\n",
    "Finally, we copy the hst1pass catalog generated with `PERT=1` so the best catalog is used later in the stacking analysis. We also copy the catalog to the main directory so users who are unable to run hst1pass can complete the tutorial. Users can change this cell to experiment with catalogs made by running hst1pass with different settings.\n",
    "\n",
    "This concludes the section on using empirical PSF models provided by the WFC3 team."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95ac8546-a490-40d9-a3db-6b1fe3d2d3dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.exists(data_dir+'id8048dyq_flc_pert1.xympqks'):\n",
    "    shutil.copy(data_dir+'id8048dyq_flc_pert1.xympqks', data_dir+'id8048dyq_flc.xympqks') # Used below for source selection.\n",
    "else:\n",
    "    shutil.copy(code_dir+'id8048dyq_flc.xympqks', data_dir+'id8048dyq_flc.xympqks') # Included in the Github repository."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb77493-1a3f-41a0-b7ca-362be7b5499e",
   "metadata": {},
   "source": [
    "<a id=\"stacked\"></a>\n",
    "### 2.3 Stacked PSF Models\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "The empirical models used above that are provided by the WFC3 team can provide excellent fits at the exposure level for very accurate astrometry and photometry. However, these models have three limitations for specific science cases. First, the models are intentionally tapered to zero flux at large radii to prevent oversubtraction and fitting of noise at large radii, which means they cannot accurately capture very extended emission of bright sources. Second, while these models are available for the majority of commonly used filters on WFC3, models are not available for more specialized narrowband and quadrant filters. Third, these models are generated in the exposure frame, and so do not accurately represent the PSF in drizzled image mosaics, which will differ depending on how the drizzle was created.\n",
    "\n",
    "With these considerations, we now explore extracting and stacking stars from both individual exposures and drizzled mosaics. Stacking a large number of sources allows for efficient outlier rejection and can produce a high quality PSF that extends to large radii, but can inherently smear or broaden the PSF at a miniscule level depending on how the sources are aligned. In this example, we will first extract stars from our image of Omega Cen and stack them to produce a PSF model. Next, we will examine the case where there are insufficient stars in an image, and retrieve stellar cutouts from the MAST database to construct a model. **In this case the images have a very low background flux, but images with a high background would require users to background subtract their data and/or PSF model for proper fitting and photometric accuracy.**\n",
    "\n",
    "In summary, users may wish to stack stars if 1) there is no empirical model available for their filter, 2) they are interested in extended wing emission and diffraction spikes, or 3) they need to generate a PSF model for a drizzled mosaic, including those with diffraction spikes at multiple angles. Users can stack stars extracted from their images, or use the MAST database cutout service to retrieve observations of a similar focus and quality for their science. We will use the hst1pass catalogs produced above for examples on the exposure level, which are the same as the centroids provided by the [MAST](https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html) cutout service, but users can provide coordinates from an external catalog (e.g. SourceExtractor, photutils) if desired."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a74124b-8cc6-4170-b62b-95506cbc37d2",
   "metadata": {},
   "source": [
    "<a id=\"stellar\"></a>\n",
    "#### 2.3.1 Stellar Stacks\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "In this example, we use the hst1pass catalog to extract and stack stars from the FLC exposure of Omega Centauri. The hst1pass catalog columns are described in Section 2.3 of [WFC3 ISR 2022-05](https://www.stsci.edu/~jayander/JWST1PASS/DOCs/ISR_HST1PASS.pdf) and include `x` and `y` detector coordinates, `m`:magnitude, `p`:peak flux, `q`:quality of fit, `k`:chip number, and `s`:sky value. We note that running hst1pass is not required, and users can use an external catalog generated with photutils or other packages provided the catalog contains accurate `x` and `y` coordinates. An example of this is shown in the section on [Drizzle PSF Models](#drizzles).\n",
    "\n",
    "We first read the hst1pass catalog and select bright stars with excellent qfit values on chip two. We also require a low sky value to remove some cutouts with many neighboring stars, and create a buffer around the edges of the detector so all cutouts contain a full array of data. The `rpix` parameter sets the radius of the cutout in pixels. For example, `rpix = 19` produces a 39x39 pixel PSF model, corresponding to a size of 1.54$\\times$1.54 arcseconds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d52d0ad4-83c7-4da0-9ea7-ccf7b1d23a7a",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\"><b>NOTE:</b> Users who are unable to run hst1pass can continue the tutorial using the included catalog file id8048dyq_flc.xympqks that is automatically copied to the data_dir. </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4c06875-e535-4b3a-a9c7-3fd43eb91564",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read the hst1pass catalog and select stars.\n",
    "hst1pass_cat = Table.read(data_dir+'id8048dyq_flc.xympqks', \n",
    "                          format='ascii', \n",
    "                          guess=False, \n",
    "                          fast_reader=False, \n",
    "                          data_start=1, \n",
    "                          names=['x', 'y', 'm', 'p', 'q', 'k', 's'])\n",
    "\n",
    "print('The original catalog has', len(hst1pass_cat), 'entries.')\n",
    "\n",
    "rpix = 19 # Set the radius in pixels for the cutouts and resulting PSF.\n",
    "\n",
    "# Select sources with specific brightness, qfit, and location criteria.\n",
    "mask = (\n",
    "    (hst1pass_cat['x'] > rpix) &\n",
    "    (hst1pass_cat['x'] < 4096.0-rpix) &\n",
    "    (hst1pass_cat['y'] > rpix) &\n",
    "    (hst1pass_cat['y'] < 2048.0-rpix) &\n",
    "    (hst1pass_cat['p'] >= 10000.0) &\n",
    "    (hst1pass_cat['p'] <= 63300.0) &\n",
    "    (hst1pass_cat['q'] <= 0.035) &\n",
    "    (hst1pass_cat['k'] == 2) &\n",
    "    (hst1pass_cat['s'] <= 100.0)\n",
    ")\n",
    "\n",
    "hst1pass_cat = hst1pass_cat[mask]\n",
    "\n",
    "star_ids = list(range(len(hst1pass_cat)))\n",
    "\n",
    "hst1pass_cat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40365740-0ec4-4264-84b4-37d14a34995d",
   "metadata": {},
   "source": [
    "The `make_cutouts` function uses the science image and catalog coordinates to extract cutouts of the stars meeting the criteria set by the mask in the cell above. The `scale_stars` option will normalize the flux of each star, which is generally recommended for proper averaging. The `sub_pixel` option uses the user-supplied `x` and `y` coordinates to align the cutouts at the subpixel level and should always be used except for testing. When `verbose` is True, additional diagnostic messages will be printed. Specifically, after performing the subpixel alignment, the function will calculate the center of mass using the inner five pixels. This value should generally be very close to the user-suppled `rpix` value, but will deviate minutely depending on how the `x` and `y` centroids were calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63f590e3-cf20-49ad-8aa1-7ae7dfe61c71",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "star_cutouts = make_cutouts(image=sci_data, \n",
    "                            star_ids=star_ids, \n",
    "                            xis=hst1pass_cat['x'], \n",
    "                            yis=hst1pass_cat['y'], \n",
    "                            rpix=rpix, \n",
    "                            scale_stars=True, \n",
    "                            sub_pixel=True, \n",
    "                            show_figs=True, \n",
    "                            verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bceefc1a-693a-4b9b-b151-cbe119ae9cf0",
   "metadata": {},
   "source": [
    "The `stack_cutouts` function then combines the array of cutouts into a single mean or median stack, with the option to normalize the flux. In general, using a `stack_type` of median is recommended to reject contaminating sources. Similarly, setting `scale_flux` to True is recommended so that all sources have normalized fluxes, which prevents biased weighting by S/N and produces a normalized PSF, which is typically required by many photometric tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d99c6b3e-be88-4b92-b2f6-aa7fcbb17c51",
   "metadata": {},
   "outputs": [],
   "source": [
    "star_stack = stack_cutouts(star_cutouts, \n",
    "                           rpix=rpix, \n",
    "                           stack_type='median', \n",
    "                           scale_flux=True, \n",
    "                           export_file=psfs_dir+'hst_wfc3_uvis_f606w_star_psf_model.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11cf0eb3-a232-4223-814c-4e8f8e0be8e2",
   "metadata": {},
   "source": [
    "Finally, we provide our stacked model to the photutils [FittableImageModel](https://photutils.readthedocs.io/en/stable/api/photutils.psf.FittableImageModel.html) function, which allows it to be used with the PSF fitter. While the core emission was better fit using the empirical models in the previous section, using the stacked model that is not tapered to zero flux allows us to fit out to large radii, which is set to 19 pixels below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9597fecd-f991-4b57-baf8-1808d215c453",
   "metadata": {},
   "outputs": [],
   "source": [
    "star_stack_psf = FittableImageModel(star_stack[0])\n",
    "photometry, psfphot = calculate_photometry(sci_data, sources, psf_model=star_stack_psf)\n",
    "resid = psfphot.make_residual_image(sci_data, (19, 19))\n",
    "figure = plot_psf_results(sci_data, resid, 1818, 986, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig7_star_stack.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdb48d6f-07bb-4896-9325-3753f45f312c",
   "metadata": {},
   "source": [
    "<a id=\"mast\"></a>\n",
    "#### 2.3.2 MAST Stacks\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "The example above demonstrates how to extract and stack stars from a provided exposure. However, in many cases there may be an insufficient number of stars to create a high-quality median stack. In those cases, users can utilize the MAST cutout database described in [WFC3 ISR 2021-12](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2021/ISR_2021_12.pdf), with access instructions provided on the [PSF Search webpage](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf/psf-search). The MAST API commands are listed on [this webpage](https://mast.stsci.edu/api/v0/pyex.html#MastCatalogsFilteredWfc3PsfUvisPy), with parameter options on [this webpage](https://mast.stsci.edu/api/v0/_w_f_c3__p_s_ffields.html). The examples below utilize sections of code from the `download_psf_cutouts.ipynb` developed by Dauphin et al. (*in preparation*). In this case, we set limits on the x and y detector locations, quality of fit, exposure time, isolation index, integrated flux (in electrons for WFC3/UVIS and WFPC2, and in electrons per second for WFC3/IR), the central pixel flux, sky flux, and exclude subarrays. The stellar centroids provided in the MAST database were calculated using hst1pass and the empirical PSF models described in [Section 2.2](#empirical). See [WFC3 ISR 2021-12](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2021/ISR_2021_12.pdf) for additional information.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><b>NOTE:</b> We explicitly assume that users have familiarized themselves with the contents of the \"download_psf_cutouts.ipynb\" notebook for this section. </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a906a387-96ef-4465-ad77-3881cac8d570",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(code_dir)\n",
    "\n",
    "# Searching a larger radius around the cutout for a sufficient number of stars.\n",
    "x_min_max = mast_api_psf.set_min_max(1818-cutout_size*5, 1818+cutout_size*5)\n",
    "y_min_max = mast_api_psf.set_min_max(986-cutout_size*5, 986+cutout_size*5)\n",
    "qfit_min_max = mast_api_psf.set_min_max(0.0, 0.05)\n",
    "exp_time = mast_api_psf.set_min_max(10, 30000)\n",
    "iso_index = mast_api_psf.set_min_max(20, 4096)\n",
    "psf_flux = mast_api_psf.set_min_max(2000, 1000000000) # This is an integrated value so set an arbitarily large upper bound.\n",
    "pixc = mast_api_psf.set_min_max(59000, 63300)\n",
    "sky_flux = mast_api_psf.set_min_max(0, 50)\n",
    "focus = mast_api_psf.set_min_max(-10, 10) # Change these values to see the effects of focus.\n",
    "subarray = mast_api_psf.set_min_max(0, 0)\n",
    "\n",
    "detector = 'UVIS'\n",
    "\n",
    "# Set the parameters, noting that the last four are not set above.\n",
    "parameters = {\n",
    "    'psf_x_center': x_min_max, \n",
    "    'psf_y_center': y_min_max, \n",
    "    'qfit': qfit_min_max, \n",
    "    'exptime': exp_time, \n",
    "    'iso_index': iso_index, \n",
    "    'psf_flux': psf_flux, \n",
    "    'pixc': pixc, \n",
    "    'sky': sky_flux, \n",
    "    'focus': focus, \n",
    "    'subarray': ['0'],\n",
    "    'n_sat_pixels': ['0'], \n",
    "    'filter': ['F606W'], \n",
    "    'mtflag': ['0'],\n",
    "}\n",
    "\n",
    "filts = mast_api_psf.set_filters(parameters)\n",
    "columns = ['id', 'rootname', 'filter', 'chip', 'exptime', 'psf_x_center', 'psf_y_center', 'pixc', 'sky', 'qfit', 'iso_index', 'subarray']\n",
    "obs = mast_api_psf.mast_query_psf_database(detector=detector, filts=filts, columns=columns)\n",
    "obs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55d3e120-393d-45e4-937a-0d9f6f7613ef",
   "metadata": {},
   "source": [
    "As described and detailed in the `download_psf_cutouts.ipynb`, the below cell constructs the filepaths for the cutouts, requests to download them from the MAST cutout service, and then extracts the files from a compressed tar folder. Finally, the filepaths for each cutout are saved to a list in `path_data` and passed to an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c70c2789-8257-4efc-b8e6-d565a5262af8",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(data_dir)\n",
    "file_suffix = ['flc']\n",
    "dataURIs = mast_api_psf.make_dataURIs(obs, detector=detector, file_suffix=file_suffix)\n",
    "filename = mast_api_psf.download_request(dataURIs, filename='mastDownload.tar.gz', download_type='bundle.tar.gz')\n",
    "with tarfile.open(filename, 'r:gz') as tar:\n",
    "    path_mast = tar.getnames()[0]\n",
    "    tar.extractall(filter='data')\n",
    "path_data = sorted(glob.glob(f'{path_mast}/*/*'))\n",
    "mast_cutouts = np.array([fits.getdata(path) for path in path_data])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c4d1aa5-ff40-44f9-8bbf-6ffae0420c14",
   "metadata": {},
   "source": [
    "Next we visually inspect all of the downloaded cutouts, where each panel shows the cutout, and any additional sources identified at > 10$\\sigma$ above the background in green circles. Users may wish to manually or automatically filter the results to create a sample with fewer contaminants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8b7af6d-ebb5-4e5c-9841-b4e3085ceab8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(data_dir)\n",
    "plot_cutout_grid(cutouts=mast_cutouts, path_data=path_data, max_cutouts=50, verbose=False)\n",
    "save_figure(figure, plot_dir+'fig8_mast_grid.pdf', True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54c012be-5451-434a-927e-8e6a47554e7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mast_stack = stack_cutouts(mast_cutouts, \n",
    "                           rpix=25, \n",
    "                           stack_type='median', \n",
    "                           scale_flux=True, \n",
    "                           export_file=psfs_dir+'hst_wfc3_uvis_f606w_mast_psf_model.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e59d071-5586-44de-ab60-a589e524abc9",
   "metadata": {},
   "source": [
    "We now use this model to perform PSF fitting and subtraction on the science image. We see that the cutouts retrieved and stacked from MAST only marginally fit the shapes of the stars in this particular exposure, similar to the results that we obtained with the empirical model before perturbation. The model could be revised by querying the MAST cutout service with different parameters and/or carefully selecting cutouts based on their FWHM. The optimal workflow will depend on the data and science goals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5581e94a-da1c-42e5-819b-e7476c31e04d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mast_stack_psf = FittableImageModel(mast_stack[0])\n",
    "photometry, psfphot = calculate_photometry(sci_data, sources, psf_model=mast_stack_psf)\n",
    "resid = psfphot.make_residual_image(sci_data, (9, 9))\n",
    "figure = plot_psf_results(sci_data, resid, 1818, 986, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig9_mast_stack.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa93c28f-887f-446b-b521-49ae16fd65e4",
   "metadata": {},
   "source": [
    "<a id=\"drizzles\"></a>\n",
    "## 3. Drizzle PSF Models\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "The examples in Section 2 apply to PSF modeling at the exposure level. However, in many cases with faint astronomical sources, multiple exposures are combined through drizzling to yield an image mosaic with increased signal-to-noise that is free from artifacts such as cosmic rays and hot pixels. This process changes the shape of the PSF due to resampling, so dedicated models are required. The process of drizzling is beyond the scope of this notebook, and users can find more information on the [DrizzlePac Website](https://www.stsci.edu/scientific-community/software/drizzlepac), including a series of [tutorial notebooks](https://spacetelescope.github.io/hst_notebooks/notebooks/DrizzlePac/README.html) available on [Github](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/DrizzlePac) with useful examples, and a brief explanation [here](https://www.stsci.edu/~fruchter/dither/drizzle.html?fbclid=IwAR3mrcUR7kxAZuDk7l6GUQtuJKrs-ecvv5YHSmWYfM64EgXP3unqp02Rdpk). In this case, it is sufficient to note that drizzling allows us to account for undersampling of the WFC3 detectors to recover the full spatial resolution of HST, and account for distortion across the field of view to provide a flux-conserved, geometrically rectified mosaic generated from the individual exposures. In this case, we combine three F606W exposures covering the same area of Omega Cen, and focused on the same four stars shown in [Section 1.4](#aperture)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91fdd13a-9da5-41b2-8ed9-46b9a2ce5a5a",
   "metadata": {},
   "source": [
    "<a id=\"create\"></a>\n",
    "### 3.1 Create Drizzle\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "In the cell below, we first download the three F606W exposures and copy them to a new directory because their headers will be modified in the alignment and drizzling process. Next, we use `TweakReg` to find sources in common between the images and align them to the first exposure. See the `TweakReg` [readthedocs](https://drizzlepac.readthedocs.io/en/deployment/tweakreg.html) for more information on each parameter. Note also that users can utilize hst1pass catalogs of their FLC files to perform high-accuracy alignments if there are a sufficient number of stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53daf06d-3c42-4298-bdea-614695edb78c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(data_dir)\n",
    "\n",
    "tweak_files = ['id8048e0q_flc.fits', 'id8048dyq_flc.fits', 'id8048e3q_flc.fits']\n",
    "\n",
    "for tweak_file in tweak_files:\n",
    "    print('Copying', tweak_file)\n",
    "    Observations.download_file('mast:HST/product/'+tweak_file, local_path=data_dir+tweak_file, cache=True)\n",
    "    shutil.copy2(tweak_file, tweak_dir)\n",
    "    \n",
    "os.chdir(tweak_dir)\n",
    "\n",
    "# Call TweakReg.\n",
    "tweakreg.TweakReg(tweak_files,\n",
    "                  refimage=tweak_files[0],\n",
    "                  expand_refcat=True,\n",
    "                  clean=True,\n",
    "                  interactive=False,\n",
    "                  updatehdr=True,\n",
    "                  reusename=False,\n",
    "                  wcsname='TWEAK', \n",
    "                  descrip='HST WFC3 PSF Notebook Example',\n",
    "                  catalog=tweak_files[0],\n",
    "                  shiftfile=True,\n",
    "                  outshifts='tweakreg_shifts.dat',\n",
    "                  minobj=100,\n",
    "                  searchrad=0.25,\n",
    "                  searchunits='arcseconds',\n",
    "                  fitgeometry='rscale',\n",
    "                  nclip=5,\n",
    "                  sigma=3.0,\n",
    "                  imagefindcfg={'threshold': 50, 'peakmax': 62000, 'use_sharp_round': True},\n",
    "                  refimagefindcfg={'threshold': 50, 'peakmax': 62000, 'use_sharp_round': True})\n",
    "\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "188a4348-67db-44de-ba71-139dc65c883b",
   "metadata": {},
   "source": [
    "Now, we examine the `TweakReg` residuals to confirm they are small. The columns are filename, shift in x, shift in y, rotation in degrees, plate scale, and the x and y RMS errors on the shift. In general, RMS values below $\\approx$ 0.1 pixels are desired for a quality alignment solution. A value of NaN indicates a failure and the search parameters must be adjusted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4d1a1b0-65f3-4877-9db1-c7812115f303",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(tweak_dir)\n",
    "shift_table = Table.read('tweakreg_shifts.dat', \n",
    "                         format='ascii.no_header',\n",
    "                         names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])\n",
    "\n",
    "formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']\n",
    "for i, col in enumerate(shift_table.colnames[1:]):\n",
    "    shift_table[col].format = formats[i]\n",
    "print(shift_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0e43e17-fb42-46e0-8168-2898da04bb3d",
   "metadata": {},
   "source": [
    "In this case, the exposures are already well aligned as the shifts are less than $\\approx$ 0.1 pixels. Because the errors on the shifts are generally smaller than the measured shifts, we adopt the shifted values. Next, we use `AstroDrizzle` to combine the aligned exposures into a single mosaic. See the `AstroDrizzle` [readthedocs](https://drizzlepac.readthedocs.io/en/deployment/astrodrizzle.html) for more information on each parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "393a94bc-4a13-4ff0-b895-1713b64ff0bc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(tweak_dir)\n",
    "\n",
    "for driz_file in tweak_files:\n",
    "    print('Copying', driz_file, 'to the drizzle directory.')\n",
    "    shutil.copy2(driz_file, driz_dir)\n",
    "    \n",
    "os.chdir(driz_dir)\n",
    "\n",
    "# Run AstroDrizzle.\n",
    "astrodrizzle.AstroDrizzle(tweak_files,\n",
    "                          output='omega_cen_f606w_drz.fits', # The output filename.\n",
    "                          wcskey='TWEAK',              # Use tweakreg WCS solution.\n",
    "                          build=True,                  # Save combined SCI and WHT.\n",
    "                          updatewcs=False,             # Don't run updatewcs again.\n",
    "                          stepsize=1,                  # Don't interpolate the WCS.\n",
    "                          num_cores=1,                 # Single core saves memory.\n",
    "                          preserve=False,              # We save the files separately.\n",
    "                          clean=True,                  # We don't need the mask files.\n",
    "                          skysub=True,                 # Definitely subtract the sky.\n",
    "                          skymethod='globalmin+match', # Recommended for diffuse data.\n",
    "                          combine_type='iminmed',      # 'iminmed' for < 10 images.\n",
    "                          driz_cr_corr=False,          # Save the cosmic-ray masks.\n",
    "                          driz_combine=True,           # Use all images in drizzle.\n",
    "                          final_wht_type='EXP',        # EXP weight for bright objects.\n",
    "                          final_pixfrac=0.8,           # Fraction to shrink drops by.\n",
    "                          final_units='cps')           # Set units to counts per second.\n",
    "\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd42691d-fb68-4e6c-b206-fabb9e89343d",
   "metadata": {},
   "source": [
    "Finally, we can display the drizzled image and highlight the same boxed region shown in Sections 1-2. The mosaic displays a significantly smoother background than the individual exposures, without significant cosmic rays or hot pixels. Stars that were saturated in the original FLCs contain artifacts and therefore should not be stacked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a2b747e-2aad-4705-a6b5-cdb4eb19af3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "with fits.open(driz_dir+'omega_cen_f606w_drz.fits') as hdul:\n",
    "    sci_driz = np.nan_to_num(hdul['SCI'].data)\n",
    "print('The image dimensions are:', np.shape(sci_driz))\n",
    "\n",
    "my_figure_size, my_fontsize = setup_matplotlib('notebook', 1.2)\n",
    "figure, axes = plt.subplots(1, 1, figsize=my_figure_size)\n",
    "cutout_size = 75 # 75 WFC3 UVIS pixels = 3 arcseconds.\n",
    "\n",
    "# The same box is at: 1823.5781, 1096.5125 in the drizzled mosaic.\n",
    "axes.add_patch(Rectangle((1823.5781-cutout_size/2, 1096.5125-cutout_size/2), width=cutout_size, height=cutout_size, fill=None, color='w'))\n",
    "norm = simple_norm(sci_driz, 'sqrt', percent=99.8-0.2)\n",
    "axes.imshow(sci_driz, cmap='gray', origin='lower', aspect='equal', interpolation='nearest', norm=norm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fec5bebd-c522-4e5c-938a-0050ac1756b9",
   "metadata": {},
   "source": [
    "<a id=\"stack\"></a>\n",
    "### 3.2 Stack Stars\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "Next, we want to detect stars that surround the location of interest, and stack them to construct a PSF model. In this case we have a sufficient number of stars (> 10) that are near the location on the detector for which we want to construct a PSF model. In this case we select a 525$\\times$525 pixel region near the four stars shown above. The functionality of the code below is described in more detail in [Section 1.4](#aperture), and the resulting stellar flux measurements can be shown by adding `print(phot_table)`. **In this case the drizzles have been background-subtracted, but images with a high background would require users to background subtract their data and/or PSF model for proper fitting and photometric accuracy.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e73741c0-aaa3-4692-be2a-dec575e6f8f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate statistics for the background, initialize the star finder, locate the stars, and assign flux apertures.\n",
    "mean, median, std = sigma_clipped_stats(sci_driz, sigma=3.0)\n",
    "daofind = DAOStarFinder(threshold=2.5, fwhm=3.0, roundhi=0.5)\n",
    "mask = create_mask(sci_driz, cutout_size=cutout_size*7.0, xcenter=1823, ycenter=1096)\n",
    "sources = daofind(sci_driz - median, mask=mask)\n",
    "positions = np.transpose((sources['xcentroid'], sources['ycentroid']))\n",
    "apertures_stellar = CircularAperture(positions, r=5.0)\n",
    "apertures_annulus = CircularAnnulus(positions, r_in=9, r_out=12)\n",
    "phot_table = aperture_photometry(sci_driz, apertures_stellar)\n",
    "\n",
    "# Calculate the sigma-clipped background levels and subtract that from the flux measured within the stellar aperture.\n",
    "aperture_stats = ApertureStats(sci_driz, apertures_stellar, sigma_clip=None)\n",
    "background_stats = ApertureStats(sci_driz, apertures_annulus, sigma_clip=SigmaClip(sigma=3.0, maxiters=10))\n",
    "total_background = background_stats.median * aperture_stats.sum_aper_area.value # area of annulus.\n",
    "flux_background_sub = aperture_stats.sum - total_background\n",
    "phot_table['total_background'] = total_background\n",
    "phot_table['flux_background_sub'] = flux_background_sub\n",
    "\n",
    "# Display the science image and photometric apertures.\n",
    "plot_apertures(data=sci_driz, \n",
    "               xcenter=1823, \n",
    "               ycenter=1096, \n",
    "               cutout_size=cutout_size*7.0, \n",
    "               apertures_stellar=apertures_stellar, \n",
    "               apertures_annulus=apertures_annulus)\n",
    "\n",
    "save_figure(figure, plot_dir+'fig10_drizzle_sources.pdf', True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eca9bbc-c924-41b0-8bdd-f4c659c04065",
   "metadata": {},
   "source": [
    "Now we use the photometric catalog to select stars from the drizzled image for stacking. In this case we shift the catalog coordinates by one pixel so that the image and catalog both have the same origin of one, rather than being zero based. Here we use the centroids obtained directly from the photometric catalog, which yields coordinates that agree with the center of mass within $\\lesssim$ 0.05 pixels. This is sufficient for many science cases, but it may be possible to improve the alignment by using the stacked model to fit the stars with the [PSFPhotometry](https://photutils.readthedocs.io/en/stable/api/photutils.psf.PSFPhotometry.html#photutils.psf.PSFPhotometry) package and iterating the fits until the centroids converge to stable values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4d695cd-8894-4d1a-aeec-f2c4c50e7b5b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "star_ids = list(phot_table['id'].value)\n",
    "xis = [x+1 for x in list(phot_table['xcenter'].value)]\n",
    "yis = [y+1 for y in list(phot_table['ycenter'].value)]\n",
    "\n",
    "driz_cutouts = make_cutouts(image=sci_driz, \n",
    "                            star_ids=star_ids, \n",
    "                            xis=xis, \n",
    "                            yis=yis, \n",
    "                            rpix=rpix, \n",
    "                            scale_stars=True, \n",
    "                            show_figs=True, \n",
    "                            sub_pixel=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32e05757-18e5-4430-be1e-30da71207aca",
   "metadata": {},
   "outputs": [],
   "source": [
    "driz_stack = stack_cutouts(driz_cutouts, \n",
    "                           rpix=rpix, \n",
    "                           stack_type='median', \n",
    "                           scale_flux=True, \n",
    "                           export_file=psfs_dir+'hst_wfc3_uvis_f606w_driz_psf_model.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d537deb-ba6a-4bf1-ab86-0ed80df1b7dd",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "We now isolate the same four stars show at the exposure level in Sections 1-2. The drizzled image has a substantially smoother background with cosmic ray and hot pixel artifacts removed by the drizzling process. In addition, a faint fifth star that was barely detectable at the exposure level is now readily visible and can also be analyzed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b0d7442-8574-4d27-a874-c2674ca04122",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate statistics for the background, initialize the star finder, locate the stars, and assign flux apertures.\n",
    "mean, median, std = sigma_clipped_stats(sci_driz, sigma=3.0)\n",
    "daofind = DAOStarFinder(threshold=0.1, fwhm=3.0, roundhi=0.5)\n",
    "mask = create_mask(sci_driz, cutout_size, 1823, 1096)\n",
    "sources = daofind(sci_driz - median, mask=mask)\n",
    "positions = np.transpose((sources['xcentroid'], sources['ycentroid']))\n",
    "apertures_stellar = CircularAperture(positions, r=5.0)\n",
    "apertures_annulus = CircularAnnulus(positions, r_in=9, r_out=12)\n",
    "phot_table = aperture_photometry(sci_driz, apertures_stellar)\n",
    "\n",
    "# Calculate the sigma-clipped background levels and subtract that from the flux measured within the stellar aperture.\n",
    "aperture_stats = ApertureStats(sci_driz, apertures_stellar, sigma_clip=None)\n",
    "background_stats = ApertureStats(sci_driz, apertures_annulus, sigma_clip=SigmaClip(sigma=3.0, maxiters=10))\n",
    "total_background = background_stats.median * aperture_stats.sum_aper_area.value # area of annulus.\n",
    "flux_background_sub = aperture_stats.sum - total_background\n",
    "phot_table['total_background'] = total_background\n",
    "phot_table['flux_background_sub'] = flux_background_sub\n",
    "\n",
    "for col in phot_table.colnames:  \n",
    "    if col not in ('id', 'npix'):\n",
    "        phot_table[col].info.format = '%.2f'\n",
    "print('\\n', phot_table, '\\n')\n",
    "\n",
    "plot_apertures(data=sci_driz, \n",
    "               xcenter=1823, \n",
    "               ycenter=1096, \n",
    "               cutout_size=cutout_size, \n",
    "               apertures_stellar=apertures_stellar, \n",
    "               apertures_annulus=apertures_annulus)\n",
    "\n",
    "save_figure(figure, plot_dir+'fig11_aperture_photometry.pdf', True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4caca377-2d38-44e1-8f1b-30975c973f8c",
   "metadata": {},
   "source": [
    "Finally, we use the stacked model from the drizzled image to perform a PSF fit and subtraction of the same stars examined earlier, with the addition of a faint fifth star on the left of the image. Overall, the model accurately encapsulates the total flux of each star, but some artifacts remain because drizzled images resample the observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6a19b3e-3584-4651-9b2c-cc9fd26db6dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "driz_stack_psf = FittableImageModel(driz_stack[0])\n",
    "photometry, psfphot = calculate_photometry(sci_driz, sources, psf_model=driz_stack_psf)\n",
    "resid = psfphot.make_residual_image(sci_driz, (17, 17))\n",
    "figure = plot_psf_results(sci_driz, resid, 1823, 1096, cutout_size)\n",
    "save_figure(figure, plot_dir+'fig12_driz_stack.pdf', True)\n",
    "calculate_residuals(resid, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba2cf78e-01bb-4465-858e-4548aceeeba1",
   "metadata": {},
   "source": [
    "If users have fewer than approximately 5-10 stars in their drizzled images and cannot create a stacked model, then it is possible to drizzle the empirical models provided by the WFC3 team into FLC files. This requires installing the [wfc3_photometry](https://github.com/spacetelescope/wfc3_photometry) package, running the included `make_model_star_image()` function within `PSFUtils.py` on the FLC files, and then extracting and stacking the drizzled models. This functionality is under development and not yet officially supported."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dc383b4-a3fa-45bb-b116-91049b317bd6",
   "metadata": {},
   "source": [
    "<a id=\"onclusions\"></a>\n",
    "## Conclusions\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "We hope this tutorial provides you with valuable tools and techniques for PSF modeling of HST observations. In summary, we have shown several different workflows for generating PSF models depending on the available observations. First, we demonstrated in [Section 2.1](#analytical) that analytical models such as Gaussian profiles fail to characterize the complex shape of the HST PSF. In [Section 2.2](#empirical), we used empirical models provided by the WFC3 team to accurately characterize the shape of the PSF. We explored using the [default models](#empirical-1), as well as those modified to match the [focus](#empirical-2) of specific observations, and characterize [spatial variations](#empirical-3) of the PSF across the detector. These models provide high-accuracy astrometry and photometry of stellar sources, but do not account for extended wing emission or diffraction spikes. We tackled this issue in [Section 2.3](#stacked) by stacking stars [extracted from](#stellar) the observations, or by [retrieving cutouts](#mast) from MAST. Stacking a large number of sources allows for efficient outlier rejection and can produce a high quality PSF that extends to large radii, but can also broaden the PSF core depending on how the sources are aligned. Finally, we extended this process to drizzled images in [Section 3](#drizzles). In some cases, more than one workflow is possible for a given set of observations, and users are encouraged to experiment with these tools to develop the best PSF models for their science goals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4115421-039e-4232-ad06-fd25c38adda7",
   "metadata": {},
   "source": [
    "<a id=\"about\"></a>\n",
    "## About this Notebook\n",
    "[Table of Contents](#toc)\n",
    "\n",
    "***\n",
    "\n",
    "**Author:** Mitchell Revalski<br>\n",
    "**Created:** 15 Apr 2024<br>\n",
    "**Updated:** 05 Jun 2024<br>\n",
    "**Source:** [https://github.com/spacetelescope/hst_notebooks](https://github.com/spacetelescope/hst_notebooks)\n",
    "\n",
    "<a id=\"additional\"></a>\n",
    "### Additional Resources\n",
    "\n",
    "Below are some additional resources that may be helpful. Please send any questions to the [HST Help Desk](https://stsci.service-now.com/hst).\n",
    "\n",
    "- [WFC3 Website](https://www.stsci.edu/hst/instrumentation/wfc3)\n",
    "- [WFC3 PSF Webpage](https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf)\n",
    "- [WFC3 Data Handbook](https://hst-docs.stsci.edu/wfc3dhb)\n",
    "- [WFC3 Instrument Handbook](https://hst-docs.stsci.edu/wfc3ihb)\n",
    "- [WFC3 Software Tools Webpage](https://www.stsci.edu/hst/instrumentation/wfc3/software-tools)\n",
    "- [Introduction to PSF photometry](https://boyce-astro.org/point-spread-function-psf/)\n",
    "- HST PSF: [Anderson & King (2000)](https://ui.adsabs.harvard.edu/abs/2000PASP..112.1360A/abstract)\n",
    "\n",
    "<a id=\"cite\"></a>\n",
    "### Software Citations\n",
    "If you use Python packages such as `astropy`, `astroquery`, `drizzlepac`, `matplotlib`, or `numpy` for published research, please cite the authors. Follow the links below for more information about citing various packages:\n",
    "\n",
    "* [Citing `astropy`](https://www.astropy.org/acknowledging.html)\n",
    "* [Citing `astroquery`](https://github.com/astropy/astroquery/blob/main/astroquery/CITATION)\n",
    "* [Citing `drizzlepac`](https://drizzlepac.readthedocs.io/en/latest/getting_started/citing_drizzlepac.html)\n",
    "* [Citing `matplotlib`](https://matplotlib.org/stable/users/project/citing.html)\n",
    "* [Citing `numpy`](https://numpy.org/citing-numpy/)\n",
    "* [Citing `photutils`](https://photutils.readthedocs.io/en/stable/getting_started/citation.html)\n",
    "***\n",
    "\n",
    "[Top of Page](#top)\n",
    "<img style=\"float: right;\"\n",
    "src = \"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"Space Telescope Logo\" width=\"200px\"/>"
   ]
  }
 ],
 "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}