examples/NIRCam_TSO_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create NIRCam TSO Simulated Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to run [Mirage](https://mirage-data-simulator.readthedocs.io/en/latest/) to create TSO data.\n",
    "\n",
    "It uses an xml and a pointing file exported from an [APT](https://jwst-docs.stsci.edu/jwst-astronomers-proposal-tool-overview) file for a proposal containing both Grism and Imaging mode time series observations. From these files, along with several other user-inputs, Mirage will create a series of input yaml files. One yaml file is produced for each detector in each exposure. Mirage then creates the simulated data corresponding to each yaml file.\n",
    "\n",
    "Note that in a TSO exposure, users will typically get data from 3 [detectors](https://jwst-docs.stsci.edu/near-infrared-camera/nircam-instrumentation/nircam-detector-overview) (A1 and A3 for shortwave data, and A5 for longwave data) in each exposure.\n",
    "\n",
    "In this example, we will use an example APT program that looks at WASP-79. The proposal contains 3 observations:\n",
    "\n",
    "* Grism TSO - F444W, with accompanying Weak Lens +8 and F182M shortwave imaging observation\n",
    "* Grism TSO - F322W2, with accompanying Weak Lens +8 and F210M shortwave imaging observation\n",
    "* Imaging TSO - F470N in the longwave channel, and Weak Lens +8 and F210M in the shortwave channel\n",
    "\n",
    "We will create yaml files for all of these observations (including the accompanying Target Acquisition exposures). We will then use Mirage to create simulated data for a longwave Grism Time Series observation, one of the Imaging observations that accompany the Grism data, as well as a standalone Imaging Time Series Observation.\n",
    "\n",
    "As written, the notebook can be executed from any directory and will find the example xml and pointing files in the Mirage repository. Outputs are all saved in the working directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [Imports](#imports)\n",
    "* [Inputs](#inputs)\n",
    "  * [Stellar Spectrum](#stellar_spectrum)\n",
    "  * [Batman parameters](#batman_parameters)\n",
    "  * [Transmission Spectrum - needed for GrTSO](#transmission_spectrum)\n",
    "  * [Grism TSO catalog](#grism_tso_catalog)\n",
    "  * [Lightcurve File - needed for Imaging TSO](#lightcurve_file)\n",
    "  * [Imaging TSO catalog](#imaging_tso_catalog)\n",
    "  * [Catalog of Background Sources](#background_catalog)\n",
    "  * [Create Input Yaml Files](#yaml_files)\n",
    "* [Create Simulated Data](#create_simulated_data)\n",
    "  * [Grism TSO](#grism_data)\n",
    "  * [Accompanying Imaging TSO](#accompanying_imaging_data)\n",
    "  * [Imaging TSO](#imaging_data)\n",
    "* [Run Calibration Pipeline](#calibration_pipeline)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set Environment Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have not yet set your CRDS-related environment variables, do that here. This must be done prior to importing the CRDS software package. Your MIRAGE_DATA environment variable should also be pointing to the location of Mirage's collection of reference files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "#os.environ[\"MIRAGE_DATA\"] = \"/my_files/jwst/simulations/mirage_data\"\n",
    "os.environ[\"CRDS_PATH\"] = os.path.expandvars(\"$HOME/crds_cache\")\n",
    "os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds.stsci.edu\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"imports\"></a>\n",
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits, ascii\n",
    "from astropy.table import Table\n",
    "from astropy.visualization import simple_norm, imshow_norm\n",
    "from astropy import units as u\n",
    "import batman\n",
    "import h5py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "import matplotlib.cm as cmx\n",
    "import stsynphot as stsyn\n",
    "from synphot import SourceSpectrum, SpectralElement\n",
    "from synphot import units\n",
    "import yaml\n",
    "\n",
    "from mirage.catalogs.hdf5_catalog import save_tso\n",
    "from mirage.catalogs.catalog_generator import GrismTSOCatalog, ImagingTSOCatalog, PointSourceCatalog\n",
    "from mirage.catalogs.catalog_generator import TSO_GRISM_INDEX\n",
    "from mirage.grism_tso_simulator import GrismTSO\n",
    "from mirage.imaging_simulator import ImgSim\n",
    "from mirage.seed_image.catalog_seed_image import Catalog_seed\n",
    "from mirage.utils.utils import ensure_dir_exists\n",
    "from mirage.yaml import yaml_generator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define paths to help organize inputs and outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = os.path.dirname(yaml_generator.__file__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_data_path = os.path.abspath(os.path.join(path, '../../examples/tso_example_data'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_data_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dir = './'\n",
    "output_yaml_dir = os.path.abspath('./yaml_files')\n",
    "ensure_dir_exists(output_yaml_dir)\n",
    "output_data_dir = os.path.abspath('./sim_data')\n",
    "ensure_dir_exists(output_data_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The yaml files that will serve as the inputs to Mirage will be saved in this directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_yaml_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulated data produced by Mirage will be saved in this directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_data_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for viewing simulated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show(array, title, min=0, max=1000):\n",
    "    \"\"\"Quick view of an array.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    array : numpy.ndimage\n",
    "        2D array\n",
    "        \n",
    "    title : str\n",
    "        Title to place on image\n",
    "        \n",
    "    min : int\n",
    "        Signal level corresponding to bottom of color scale\n",
    "        \n",
    "    max : int\n",
    "        Signal corresponding to top of color scale\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots(figsize=(10, 10))\n",
    "    norm = simple_norm(array, stretch='log', min_cut=min, max_cut=max)\n",
    "    cax = ax.imshow(array, norm=norm, origin='lower')\n",
    "    cbar = fig.colorbar(cax)\n",
    "    plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"prepare_inputs\"></a>\n",
    "## Prepare Inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prior to simulated data creation, there are a number of input parameters to set and files to be created. To start with, define the xml and pointing files associated with the APT file. These files can be saved from within APT by selecting them from table that appears when you select **Export** from the **File** menu. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xml_file = os.path.join(input_data_path, 'wasp-79_example_TSO.xml')\n",
    "pointing_file = xml_file.replace('.xml', '.pointing')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"stellar_spectrum\"></a>\n",
    "### Stellar spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the spectrum of the unocculted star associated with the TSO object. For this example, we'll use the Castelli & Kurucz models in stsynphot to generate a spectrum that is similar to WASP-79. You can generate your spectrum using any tool you like, as it will eventually be saved in an hdf5 file below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_eff = 6750  # surface temperature\n",
    "metallicity = 0.03\n",
    "log_g = 4.26  # surface gravity = 182 m/s^2\n",
    "sp = stsyn.grid_to_spec('ck04models', t_eff, metallicity, log_g) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalize the spectrum to be k = 9.06. There are two ways you can scale your spectrum. This first is shown here, where you manually scale the spectrum before saving to the hdf5 file. The second way is to leave the scaling to Mirage. In that case, you save the spectrum as-is, and set the flux units in the hdf5 file to 'normalized'. With that option, Mirage will automatically scale the spectrum to the magnitude indicated in the grism TSO source catalog. If you choose to scale the spectrum manually (and use flux units of 'flam' in the hdf5 file), then Mirage will ignore the source magnitude listed in the grism TSO catalog, and use the saved spectrum with no changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize the spectrum\n",
    "bp = SpectralElement.from_filter('johnson_k')\n",
    "vega = SourceSpectrum.from_vega()\n",
    "sp_norm = sp.normalize(9.06 * units.VEGAMAG, bp, vegaspec=vega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get wavelengths and flux densities of the spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelengths = sp_norm.waveset.to(u.micron)\n",
    "fluxes = sp_norm(wavelengths, flux_unit='flam')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the spectrum in the wavelengths of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, a = plt.subplots()\n",
    "a.plot(wavelengths, fluxes)\n",
    "a.set_xlim(1, 5)\n",
    "a.set_ylim(0, 1e-13)\n",
    "a.set_xlabel('Wavelength (microns)')\n",
    "a.set_ylabel('Flux density (FLAM)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the units for the wavelength and flux density arrays. It's generally recommended to use flux denisty units of FLAM (erg / s / cm^2 / 𝐴˚). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelength_units = 'microns'\n",
    "flux_units = 'flam'\n",
    "#flux_units = 'normalized'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Name of the file in which to save the spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sed_file = os.path.join(output_dir, 'test_grism_tso_sed_file_wasp79.hdf5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HDF5 files can contain multiple datasets, so make the flux and wavelength variables into lists, even though for TSO observations there will be only one source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes = [fluxes]\n",
    "wavelengths = [wavelengths]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the spectrum in the hdf5 file. Give the dataset a special index number, defined by TSO_GRISM_INDEX, to help Mirage keep this source separate from any provided by other catalogs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(sed_file, \"w\") as file_obj:\n",
    "    for i in range(len(fluxes)):\n",
    "        dset = file_obj.create_dataset(str(i+TSO_GRISM_INDEX), data=[wavelengths[i].value, fluxes[i].value],\n",
    "                                       dtype='f', compression=\"gzip\", compression_opts=9)\n",
    "        dset.attrs[u'wavelength_units'] = wavelength_units\n",
    "        dset.attrs[u'flux_units'] = flux_units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"batman_parameters\"></a>\n",
    "### Batman Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model lightcurves are generated using the [Batman](https://github.com/lkreidberg/batman) package. Here, set the parameters that Batman needs. See the [Batman documentation](https://www.cfa.harvard.edu/~lkreidberg/batman/) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this case, we'll use parameters based on WASP-79b, but modified for a shorter exposure time\n",
    "# and a deeper eclipse, in order to save processing time and make the resulting data easier to\n",
    "# examine.\n",
    "params = batman.TransitParams()       # object to store transit parameters \n",
    "params.t0 = 280.                        # time of inferior conjunction \n",
    "params.per = 3162.24                       # orbital period \n",
    "params.rp = 0.723                       # planet radius (in units of stellar radii) \n",
    "params.a = 9.37                        # semi-major axis (in units of stellar radii) \n",
    "params.inc = 83.3                      # orbital inclination (in degrees) \n",
    "params.ecc = 0.                       # eccentricity \n",
    "params.w = 90.                        # longitude of periastron (in degrees) \n",
    "params.limb_dark = \"nonlinear\"        # limb darkening model \n",
    "params.u = [0.5, 0.1, 0.1, -0.1]      # limb darkening coefficients [u1, u2, u3, u4] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the array of times at which the lightcurve will be calculated. In this case we know that the total exposure time is about 570 seconds, so we extend the array of times just beyond that range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = np.linspace(0, 580, 1000)  # times at which to calculate light curve "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"transmission_spectrum\"></a>\n",
    "### Create a Transmission Spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transmission spectrum is the wavelength-dependent effective radius of the planet, in units of the stellar radius. It must be saved in an ascii file, and will be used in the creation of the Grism TSO data. For this example, we use the simple case of a flat spectrum. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "waves = np.linspace(0.9, 5.5, 1000)  # microns\n",
    "trans = np.repeat(params.rp, 1000)  # R_planet / R_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tran_spec_file = os.path.join(output_dir,'transmission_spectrum.txt')\n",
    "tab = Table()\n",
    "tab['Wavelength'] = waves\n",
    "tab['Transmission'] = trans\n",
    "tab.write(tran_spec_file, format='ascii', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot transmission spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, a = plt.subplots()\n",
    "a.scatter(waves, trans, color='red', marker='o')\n",
    "a.set_xlabel('Wavelength (microns)')\n",
    "a.set_ylabel('Transmission')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"grism_tso_catalog\"></a>\n",
    "### Create Grism TSO catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the stellar spectrum and transmission spectrum saved, we can now create the Grism TSO source catalog that Mirage will use when creating the simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Name of catalog file to hold information on the TSO source\n",
    "grism_tso_catalog = os.path.join(output_dir,'tso_grism_source.cat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic information on the source. Note that the object magnitude will be ignored if the saved stellar spectrum is in units of FLAM. Conversely, if the stellar spectrum has units of 'normalized', then Mirage will scale the spectrum to the magnitude indicated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "object_ra = 66.37090333\n",
    "object_dec = -30.60044722\n",
    "object_f444w_mag = 9.0\n",
    "object_f322w2_mag = 9.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the Grism TSO catalog object and populate RA, Dec, Batman parameters, times, and the name of the transmission spectrum file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grism_cat = GrismTSOCatalog(ra=[object_ra], dec=[object_dec], semimajor_axis=[params.a],\n",
    "                            orbital_inclination=[params.inc], eccentricity=[params.ecc],\n",
    "                            orbital_period=[params.per], longitude_of_periastron=[params.w],\n",
    "                            limb_dark_model=[params.limb_dark], limb_dark_coeffs=[params.u],\n",
    "                            time_units=['second'], start_time=[np.min(times)],\n",
    "                            end_time=[np.max(times)], inferior_conj=[params.t0],\n",
    "                            transmission_spectrum=[tran_spec_file])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the source magnitudes to the catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add source magnitude\n",
    "grism_cat.add_magnitude_column([object_f444w_mag], magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f444w')\n",
    "grism_cat.add_magnitude_column([object_f322w2_mag], magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f322w2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grism_cat.save(grism_tso_catalog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the contents of the catalog. There should be only a single source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grism_cat.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"lightcurve_file\"></a>\n",
    "### Create Lightcurve File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the imaging time series observations, we need to provide a file that contains the lightcurve to use when simulating the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the batman model using the parameters specified in the [Batman Parameters](batman_parameters) section and generate a lightcurve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = batman.TransitModel(params, times)\n",
    "flux = m.light_curve(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the lightcurve to be used to generate the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, a = plt.subplots()\n",
    "a.scatter(times, flux, color='red', marker='v')\n",
    "a.set_xlabel('Time (sec)')\n",
    "a.set_ylabel('Normalized Signal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightcurve_file = os.path.join(output_dir, 'example_lightcurve.hdf5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Place the lightcurve into a dictionary to prepare for saving. The keys are object indexes corresponding to objects in the Mirage input catalogs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contents = {}\n",
    "contents['1'] = {'times': times,\n",
    "                 'fluxes': flux}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save\n",
    "save_tso(contents, lightcurve_file, time_unit='second')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"imaging_tso_catalog\"></a>\n",
    "### Create Imaging TSO catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, rather than the `tso_grism_catalog` in the yaml file, the user must supply a `tso_imaging_catalog`. This catalog will contain the list of sources whose flux will be varying with time. As in the grism TSO case, objects listed in the other types of catalogs will be added to the simulation in order to create a more realistic scene.\n",
    "\n",
    "See the [notebook on catalog creation](https://github.com/spacetelescope/mirage/blob/master/examples/Catalog_Generation_Tools.ipynb) in the examples directory of the Mirage repository for examples of how to create source catalogs. Catalogs may also be created manually.\n",
    "\n",
    "Mirage does not yet support the use of input spectra in order to determine filter-based magnitudes for imaging simulations. Until that ability is added, users must estimate the source magnitude in each filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imaging_tso_catalog = os.path.join(output_dir, 'tso_imaging_source.cat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsimg_cat = ImagingTSOCatalog(ra=[object_ra], dec=[object_dec], lightcurve_file=[lightcurve_file])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "object_f182m_mag = 10.0\n",
    "object_f210m_mag = 9.5\n",
    "object_f470n_mag = 9.0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add source magnitudes\n",
    "tsimg_cat.add_magnitude_column([object_f182m_mag], magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f182m/wlp8')\n",
    "tsimg_cat.add_magnitude_column([object_f210m_mag], magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f210m')\n",
    "tsimg_cat.add_magnitude_column([object_f470n_mag], magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f470n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsimg_cat.save(imaging_tso_catalog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsimg_cat.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"background_catalog\"></a>\n",
    "### Catalog of Background Sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a catalog containing nearby sources. Mirage will add these to the data. This can help inform desired roll angles in order to avoid contamination.\n",
    "\n",
    "Mirage accepts a series of ascii source catalogs. In this case, the options that can be used are:\n",
    "\n",
    "1. Point source catalog\n",
    "2. Galaxy catalog\n",
    "3. Extended source catalog\n",
    "\n",
    "The point source, galaxy, and extended source catalogs contain \"background\" sources. That is, sources other than the primary TSO target. All catalogs are input to Mirage as entries within the input yaml file.\n",
    "\n",
    "See the notebook on [catalog creation](https://github.com/spacetelescope/mirage/blob/master/examples/Catalog_Generation_Tools.ipynb) in the examples directory of the Mirage repository for examples of how to create source catalogs. Catalogs may also be created manually. For this example, we will create a point source catalog with a single star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_sources_ra = [66.378244]\n",
    "bkgd_sources_dec = [-30.598312]\n",
    "bkgd_sources_f444w_mag = [12.50]\n",
    "bkgd_sources_f322w2_mag = [12.75]\n",
    "bkgd_sources_f182m_mag = [13.75]\n",
    "bkgd_sources_f210m_mag = [13.39]\n",
    "bkgd_sources_f470n_mag = [12.40]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_cat_file = os.path.join(output_dir, 'ptsrcs.cat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_cat = PointSourceCatalog(ra=bkgd_sources_ra, dec=bkgd_sources_dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add source magnitudes\n",
    "bkgd_cat.add_magnitude_column(bkgd_sources_f182m_mag, magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f182m/wlp8')\n",
    "bkgd_cat.add_magnitude_column(bkgd_sources_f444w_mag, magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f444w')\n",
    "bkgd_cat.add_magnitude_column(bkgd_sources_f322w2_mag, magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f322w2')\n",
    "bkgd_cat.add_magnitude_column(bkgd_sources_f210m_mag, magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f210m')\n",
    "bkgd_cat.add_magnitude_column(bkgd_sources_f470n_mag, magnitude_system='vegamag',\n",
    "                               instrument='nircam', filter_name='f470n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_cat.save(bkgd_cat_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_cat.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"yaml_files\"></a>\n",
    "### Create input yaml files for Mirage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mirage requires one input yaml file for each detector and exposure. Using the files exported from ATP combined with the various files produced above, Mirage can automatically generate all of the yaml files necessary to produce all of the data in the APT file. See the Mirage documentation on [creating yaml files from APT inputs](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html) for details on the input parameters used below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Populate the input catalog dictionary with the names of the source catalogs created above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogs = {'WASP-79': {'point_source': bkgd_cat_file,\n",
    "                        'tso_imaging_catalog': imaging_tso_catalog,\n",
    "                        'tso_grism_catalog': grism_tso_catalog,\n",
    "                        }\n",
    "           }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the desired background level for the observations. Most common will to specify 'low', 'medium', or 'high'. The definitions of these three levels match those used in the ETC. Note that for the grism observations, if the `dateobs_for_background` is set to True, the background will be created based on the given observation date rather than the level given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "background = 'medium'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the desired telescope roll angle and date for the observation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pav3 = 0.\n",
    "dates = '2021-10-25'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "yam = yaml_generator.SimInput(xml_file, pointing_file, catalogs=catalogs, verbose=True,\n",
    "                              output_dir=output_yaml_dir, simdata_output_dir=output_data_dir,\n",
    "                              background=background, roll_angle=pav3,\n",
    "                              dates=dates, datatype='linear, raw', dateobs_for_background=True,\n",
    "                              reffile_defaults='crds')\n",
    "\n",
    "yam.use_linearized_darks = True\n",
    "yam.create_inputs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List the yaml files created by the generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yam.yaml_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a table showing details of what exposure each yaml file describes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basenames = []\n",
    "modes = []\n",
    "apertures = []\n",
    "filters = []\n",
    "pupils = []\n",
    "for yfile in yam.yaml_files:\n",
    "    with open(yfile) as file_obj:\n",
    "        info = yaml.safe_load(file_obj)\n",
    "    basenames.append(os.path.basename(yfile))\n",
    "    modes.append(info['Inst']['mode'])\n",
    "    apertures.append(info['Readout']['array_name'])\n",
    "    filters.append(info['Readout']['filter'])\n",
    "    pupils.append(info['Readout']['pupil'])\n",
    "info_tab = Table()\n",
    "info_tab['Filename'] = basenames\n",
    "info_tab['Mode'] = modes\n",
    "info_tab['Aperture'] = apertures\n",
    "info_tab['Filter'] = filters\n",
    "info_tab['Pupil'] = pupils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exposures using `imaging` mode are for the target acquisition exposures. Those with `ts_grism` mode describe the Grism time series observations. Those with `ts_imaging` mode are for Imaging time series observations. This includes the data collected with the shortwave channel detectors while the longwave detector is producing the grism time series observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"create_simulated_data\"></a>\n",
    "## Create Simulated Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate data from some of the files above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"grism_data\"></a>\n",
    "### Grism TSO Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's simulate a grism time series exposure using the F444W filter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr_tso_yaml_file = os.path.join(output_yaml_dir, 'jw88888001001_01101_00002_nrca5.yaml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr_tso_yaml_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# NOTE: This cell will take a while (~30 min) to run \n",
    "gr_f444w = GrismTSO(gr_tso_yaml_file, SED_file=sed_file, SED_normalizing_catalog_column=None,\n",
    "                    final_SED_file=None, save_dispersed_seed=True, source_stamps_file=None,\n",
    "                    extrapolate_SED=True, override_dark=None, disp_seed_filename=None,\n",
    "                    orders=[\"+1\", \"+2\"])\n",
    "gr_f444w.create()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examine the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mirage will split the exposure into mulitple files if it is too large, in the same way that the JWST calibration pipeline will do. Exposures will be split around integrations. That is, an integration will not be split between two files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, look at the seed image, which is the noiseless scene that contains only signal from the astronomical sources, including background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gr_f444w.seed_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the first seed file. Checking the array dimensions, we see that it holds the first 60 integrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_444 = fits.getdata(gr_f444w.seed_files[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_444.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the seed image. The TSO source is the lower source of the two. The upper source is the background star contained in the point source catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(data_444[0, 4, :, :], 'F444W: Final group of integration 0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the seed image, plot the signal along the trace of the TSO object for each integration, as a sanity check to see that the signal does change with time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multi_spec_plot(array, yval, xstart):\n",
    "    \"\"\"Plot the signal along the `yval` row from `xstart` to the end,\n",
    "    in the final group of each integration. This will be a quick\n",
    "    sanity check to show how the signal is changing with wavelength\n",
    "    across the integrations. Step along the color table and use a\n",
    "    different color for each integration.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    array : numpy.ndarray\n",
    "        4D array containing the spectrum in each group and integration\n",
    "        \n",
    "    yval : int\n",
    "        Row number to plot\n",
    "        \n",
    "    xstart : int\n",
    "        First column number to plot. Plot from here to the end of the row\n",
    "    \"\"\"\n",
    "    num_integrations = array.shape[0]\n",
    "    num_grps = array.shape[1]\n",
    "    clrs = cmx.rainbow(np.linspace(0, 1, num_integrations))\n",
    "    f, a = plt.subplots()\n",
    "    f.set_size_inches(11.5,9)\n",
    "    for integration, color in zip(range(num_integrations), clrs): \n",
    "        slopes = array[integration, num_grps - 1, yval, xstart:] \n",
    "        a.plot(np.arange(len(slopes)), slopes, color=color)\n",
    "        a.set_xlabel('Column Number')\n",
    "        a.set_ylabel('Signal (DN)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plotting function above uses the `rainbow` color table, which starts with blue and moves through the spectrum to red. Below we plot the signal in the 51st row of the aperture in the final group of each of the 60 integrations. We see that the initial integrations (blue, but covered and not visible) and the final integrations (red) overlay each other, as expected since they occur before and after the transit, respectively. In the middle of the exposure, the signal in the spectrum decreases with time (light blue -> green) down to a minimum and then increases again (orange), as expected. The F444W filter covers a wavelength range of approximately 3.8 - 5.1 microns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multi_spec_plot(data_444, 51, 700)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the raw, or uncalibrated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uncal_file = os.path.join(output_data_dir, 'jw88888001001_01101_00002-seg001_nrca5_uncal.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uncal = fits.getdata(uncal_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at a difference image between the last and first groups of an integration. By using the difference, we remove bias offsets and make the sources more visible. Raw files have signals in integers, so multiply by 1.0 in order to make them floats before subtracting. In this case, the signal from the background star is faint enough that it is not visible here above the noise.\n",
    "\n",
    "The elevated background signal in the 512-column wide area is the result of each 512-column block being read out by a different amplifier. These differences will be minimized once the data have been run through the calibration pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = 1. * uncal[0, 4, :, :] - 1. * uncal[0, 0, :, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(diff, \"Uncalibrated Data\", min=0, max=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create simulated data for the same observation, but in the F322W2 filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr_tso_f322w2_yaml_file = os.path.join(output_yaml_dir, 'jw88888002001_01101_00002_nrca5.yaml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gr_f322w2 = GrismTSO(gr_tso_f322w2_yaml_file, SED_file=sed_file, SED_normalizing_catalog_column=None,\n",
    "                     final_SED_file=None, save_dispersed_seed=True, source_stamps_file=None,\n",
    "                     extrapolate_SED=True, override_dark=None, disp_seed_filename=None,\n",
    "                     orders=[\"+1\", \"+2\"])\n",
    "gr_f322w2.create()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in one of the seed image files and examine the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_322 = fits.getdata(gr_f322w2.seed_files[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how with F322W2 as the crossing filter, the location of the dispersed objects have shifted to the left relative to what they were with the F444W filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "show(data_322[0, 4, :, :], 'F322W2: Final group of integration 0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the signal along the trace, as with the F444W data. Again we see that the initial (blue) and final (red) integrations, which were outside of the transit, show identical signals at the maximum level, while those integrations that occured within the transit (light blue, green, orange) show reduced signal. The F322w2 filter covers a wavelength range of approximately 2.4 - 4.0 microns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multi_spec_plot(data_322, 51, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"accompanying_imaging_data\"></a>\n",
    "### Accompanying Imaging TSO Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the imaging time series data collected on one of the shortwave detectors while the grism data above was being collected in the LW channel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_tso_sw_yaml = os.path.join(output_yaml_dir, 'jw88888001001_01101_00001_nrca1.yaml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "img_tso = ImgSim()\n",
    "img_tso.paramfile = img_tso_sw_yaml\n",
    "img_tso.create()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the raw output file. This is the simulated data before any calibration. Detector A1 covers the same area as the [left half of detector A5](https://jwst-docs.stsci.edu/near-infrared-camera/nircam-instrumentation/nircam-detector-overview) (which was used to collect the Grism TSO above). This is why the source appears towards the right edge of A1. The +8 weak lens was also used in this observation, resulting in the expanded PSF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_seed_file = os.path.join(output_data_dir,\n",
    "                             'jw88888001001_01101_00001_nrca1_uncal_F182M_seg001_part001_seed_image.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_seed = fits.getdata(img_seed_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(img_seed[0, 4, :, :], 'Noiseless Seed Image', min=0, max=3000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, let's look at signal levels in all of the integrations, using the noiseless seed image. Since there is no noise and no background, we'll do a simple sum of the signal in a small box. Here we see that the depth of the transit matches that from the lightcurve we created with the Batman parameters above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signals = []\n",
    "cx, cy = 1857, 151\n",
    "for i in range(img_seed.shape[0]):\n",
    "    signals.append(np.sum(img_seed[i, 4, cy-1: cy+2, cx-1: cx+2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, a = plt.subplots()\n",
    "a.plot(np.arange(len(signals)), signals, color='red')\n",
    "a.set_xlabel('Integration Number')\n",
    "a.set_ylabel('Aperture Signal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now look at the raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uncal_file = os.path.join(output_data_dir, 'jw88888001001_01101_00001-seg001_nrca1_uncal.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = fits.getdata(uncal_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By displaying the difference between the last and first groups of one integration, we can remove bias signal that would otherwise make the object difficult to see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cds = 1.* data[0, 4, :, :] - 1. * data[0, 0, :, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image appears to have 4 512-column wide blocks because each block has been read out by a different amplifier. These differences will be minimized once the data have been run through the calibration pipeline. Note that this observation was taken with one of NIRCam's weak lenses in place, which is why the star's PSF has become so large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(cds, 'Difference Image', min=0, max=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"imaging_data\"></a>\n",
    "### Imaging TSO Observation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The imaging time series observation in this proposal contains an observation of WASP-79 using the SUB400P subarray. This is a 400x400 pixel subarray located in the upper right corner of the B1 and B5 detectors. Run Mirage on both the B1 and B5 exposures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_tso_sw_yaml = os.path.join(output_yaml_dir, 'jw88888003001_01101_00001_nrcb1.yaml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "img_tso = ImgSim()\n",
    "img_tso.paramfile = img_tso_sw_yaml\n",
    "img_tso.create()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_tso_lw_yaml = os.path.join(output_yaml_dir, 'jw88888003001_01101_00002_nrcb5.yaml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "img_tso = ImgSim()\n",
    "img_tso.paramfile = img_tso_lw_yaml\n",
    "img_tso.create()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at the seed image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_image_file = os.path.join(output_data_dir, 'jw88888003001_01101_00001_nrcb1_uncal_F210M_seed_image.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with fits.open(seed_image_file) as hdulist:\n",
    "    sw_img_data = hdulist[1].data\n",
    "    grptime = hdulist[0].header['TGROUP']\n",
    "sw_img = sw_img_data[0, -1, :, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the final group of the initial integration on the B1 detector. The edges of the PSF stamp from the library are visible against the zero-signal background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(sw_img,'Imaging TSO', min=0, max=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the same image for the B5 detector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lw_seed_image_file = os.path.join(output_data_dir, 'jw88888003001_01101_00002_nrcb5_uncal_F470N_seed_image.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with fits.open(lw_seed_image_file) as hdulist:\n",
    "    lw_img_data = hdulist[1].data\n",
    "lw_img = lw_img_data[0, -1, :, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pixel scale in the longwave (B5) detector is is twice as large as in the shortwave detector, so the same 400x400 pixel aperture covers 4 times the area on the sky compared to the B1 case above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(lw_img,'Imaging TSO', min=0, max=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rough photometry and lightcurve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As was done in the [Imaging TSO](#imaging_data) section, we perform rough photometry on each of the integrations within the exposure in an effort to plot the lightcurve. As before, the depth of the measured lightcurve matches that in the initial Batman parameters defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_ints, num_groups, ny, nx = sw_img_data.shape\n",
    "sw_img_phot = np.zeros(num_ints)\n",
    "for integ in range(num_ints):\n",
    "    sw_img_phot[integ] = np.sum(sw_img_data[integ, -1, 100:300, 100:300])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, a = plt.subplots()\n",
    "a.plot(np.arange(len(sw_img_phot)) * grptime * (num_groups + 1), sw_img_phot, color='red')\n",
    "a.set_xlabel('Time (sec)')\n",
    "a.set_ylabel('Total Aperture Signal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"calibration_pipeline\"></a>\n",
    "## Run the Calibration Pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the accompanying notebook for a demonstration of how to run the calibration pipeline on these data."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}