{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook presents METIS LSS L-band simulation of three models of a young stellar object (YSO)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scopesim as sim\n",
    "sim.bug_report()\n",
    "\n",
    "# Edit this path if you have a custom install directory, otherwise comment it out. [For ReadTheDocs only]\n",
    "sim.rc.__config__[\"!SIM.file.local_packages_path\"] = \"../../../\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from astropy.io import fits\n",
    "from astropy import units as u\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you haven't got the instrument packages yet, uncomment the following cell, which will install the packages into `./inst_pkgs`, a subdirectory of your current working directory. This is the default location where scopesim looks for instrument packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sim.download_package([\"instruments/METIS\", \"telescopes/ELT\", \"locations/Armazones\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have downloaded the packages but to a different location, you can set\n",
    "```python\n",
    "sim.rc.__config__[\"!SIM.file.local_packages_path\"] = <\"/path/to/inst/pkgs\">\n",
    "```\n",
    "We recommend, however, to create a working directory for each simulation project and to use the default installation of packages into a subdirectory. Keeping simulation results and configuration files together that way makes it easy to reconstruct later the exact conditions under which a simulation was run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation of source cubes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input data are cubes of three different models of the same YSO, HD100546. We keep the names of FITS files, the `Source` objects and the results of the Scopesim simulations in dictionaries, indexed by short names for the models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitsfiles = {}\n",
    "fitsfiles['cav'] = \"models_Lband_HD100546_cav_f100PAH.cube_3.0mas.fits\"\n",
    "fitsfiles['emptycav'] = \"models_Lband_HD100546_empytcav.cube_3.0mas.fits\"\n",
    "fitsfiles['gap'] = \"models_Lband_HD100546_gap100.cube_3.0mas.fits\"\n",
    "models = list(fitsfiles.keys())\n",
    "print(\"Model names:\", models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FITS files can be downloaded from the Scopesim server. If you already have them, make sure that the files are in the current working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.download_example_data(list(fitsfiles.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file headers are not yet in the form that Scopesim understands and we make two minor modifications: \n",
    "- Set CRVAL to 0, because Scopesim cannot look elsewhere\n",
    "- Set BUNIT keyword (files have UNITS, which is non-standard)\n",
    "- The cubes contain the occasional negative value. We replace these with 0.\n",
    "- We introduce a factor `scale_delt` to increase the pixel size, which makes features more visible. If you want to simulate the original source pixel scale, set `scale_delt` to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources = {}\n",
    "scale_cdelt = 1\n",
    "for model, fitsfile in fitsfiles.items():\n",
    "    with fits.open(fitsfile) as hdul:\n",
    "        hdul[0].header['CRVAL1'] = 0.\n",
    "        hdul[0].header['CRVAL2'] = 0.\n",
    "        hdul[0].header['CDELT1'] *= scale_cdelt\n",
    "        hdul[0].header['CDELT2'] *= scale_cdelt\n",
    "        hdul[0].header['BUNIT'] = hdul[0].header['UNITS']\n",
    "        hdul[0].data[hdul[0].data < 0] = 0\n",
    "        sources[model] = sim.Source(cube=hdul)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get an impression of what the data look like we display the cubes summed along the wavelength direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine plot limits in arcsec from header keywords\n",
    "hdr = sources['cav'].cube_fields[0].header\n",
    "i_lim = np.array([0, hdr['NAXIS1']])\n",
    "x_lim = hdr['CRVAL1'] + hdr['CDELT1'] * (i_lim + 1 - hdr['CRPIX1']) * 3600\n",
    "j_lim = np.array([0, hdr['NAXIS2']])\n",
    "y_lim = hdr['CRVAL2'] + hdr['CDELT2'] * (j_lim + 1 - hdr['CRPIX2']) * 3600\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 4))\n",
    "for i, (model, src) in enumerate(sources.items()):\n",
    "    im = axes[i].imshow(src.cube_fields[0].data.sum(axis=0) + 1e-14,   # add small positive value to avoid 0 in LogNorm\n",
    "                        origin='lower', norm=LogNorm(vmin=1e-4, vmax=10),\n",
    "                        extent=(x_lim[0], x_lim[-1], y_lim[0], y_lim[-1]))\n",
    "    axes[i].set_xlabel(\"arcsec\")\n",
    "    axes[i].set_ylabel(\"arcsec\")\n",
    "    axes[i].set_title(model);\n",
    "fig.subplots_adjust(right=0.9)\n",
    "cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75])\n",
    "fig.colorbar(im, cax=cbar_ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In wavelength, the cubes are sampled on a linear wavelength grid from 3.1 to 4 $\\mu\\mathrm{m}$, with a step size of $0.2\\,\\mu\\mathrm{m}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources['cav'].cube_fields[0].wave"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation with Scopesim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cubes are observed in the long-slit spectroscopic mode in the L band. As usual, there are four steps: `UserCommands` -> `OpticalTrain` -> `observe` -> `readout` to arrive at a detector image. The optical train can be reused for observation of different source when `update=True` is set in `observe()`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exptime = 3600.   # seconds\n",
    "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"lss_l\"],\n",
    "                      properties={\"!OBS.exptime\": exptime,\n",
    "                                 \"!SIM.spectral.spectral_resolution\": 20000,\n",
    "                                 \"!SIM.spectral.spectral_bin_width\": 2e-4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis = sim.OpticalTrain(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {}\n",
    "for model, src in sources.items():\n",
    "    print(f'Observing model \"{model}...\"')\n",
    "    metis.observe(src, update=True)\n",
    "    results[model] = metis.readout(detector_readout_mode='auto')[0]\n",
    "    print(\"-----\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "for i, (model, result) in enumerate(results.items()):\n",
    "    plt.subplot(1, 3, i+1)\n",
    "    plt.imshow(result[1].data, origin='lower', norm=LogNorm())\n",
    "    plt.xlabel(\"pixel\")\n",
    "    plt.ylabel(\"pixel\")\n",
    "    plt.title(model);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is a simulation of a blank-sky observation, which will be used for background subtraction. Given that the source is fairly compact, one could estimate the background from the science observation as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sky = sim.source.source_templates.empty_sky()\n",
    "metis.observe(sky, update=True)\n",
    "bgresult = metis.readout(detector_readout_mode=\"auto\")[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before showing the background subtracted spectra, we convert the pixel numbers to wavelength and spatial position along the slit. This is possible because the current version of Scopesim/METIS by default uses perfectly linear mapping of the spectra onto the detector. The WCS keywords to use are in `metis['spectral_traces'].meta`. The plots will be restricted to the area covered by the spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# meta = metis['spectral_traces'].meta\n",
    "\n",
    "# det_xi = meta['CRVAL1'] + meta['CDELT1'] * (np.arange(2048) + 1 - meta['CRPIX1'])  * u.Unit(meta['CUNIT1'])\n",
    "# det_xi = det_xi.to(u.arcsec)\n",
    "\n",
    "# det_wave = (meta['CRVAL2'] + meta['CDELT2'] * (np.arange(2048) + 1 - meta['CRPIX2'])) * u.Unit(meta['CUNIT2'])\n",
    "# det_wave = det_wave.to(u.um)          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ilim = np.array([750, 1300])    # pixels along spatial direction\n",
    "# jlim = np.array([300, 1750])    # pixels along wavelength direction\n",
    "# xlim = det_xi[ilim].value\n",
    "# ylim = det_wave[jlim].value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure(figsize=(15, 8))\n",
    "# for i, (model, result) in enumerate(results.items()):\n",
    "#     plt.subplot(1, 3, i+1)\n",
    "#     plt.imshow((result[1].data - bgresult[1].data)[jlim[0]:jlim[1], ilim[0]:ilim[1]], origin='lower', norm=LogNorm(),\n",
    "#               extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
    "#     plt.xlabel(r\"position along slit [arcsec]\")\n",
    "#     plt.ylabel(r\"Wavelength [$\\mu$m]\")\n",
    "#     plt.title(model);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results have to be saved to disk explicitely so they can be analysed with external tools. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from pathlib import Path\n",
    "# for i, (model, result) in enumerate(results.items()):\n",
    "#     outfile = Path(fitsfiles[model]).stem + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n",
    "#     result.writeto(outfile, overwrite=True)\n",
    "#     print(fitsfiles[model], \"--->\", outfile)\n",
    "# bgresult.writeto(\"models_Lband_HD100546-background_simulated_LSS_L.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Note on the spectral mapping\n",
    "To use the true spectral mapping of METIS, the trace file has to be set before building the optical train, for example by setting `\"!OBS.trace_file\"` in the `UserCommands`:  \n",
    "```python\n",
    "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"lss_l\"],\n",
    "                      properties={\"!OBS.exptime\": 3600,\n",
    "                                  \"!OBS.trace_file\": \"TRACE_LSS_L.fits\"})\n",
    "```\n",
    "For the METIS long-slit mode, the non-linear parts of the actual mapping are quite small.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rotating the field\n",
    "In Scopesim, the METIS slit is always aligned with the x-axis. To simulate different slit orientations, the input cube has to be rotated. We define a function to do this (for more explanation of the steps see the notebook `LSS_AGN-01_preparation.ipynb`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.wcs import WCS\n",
    "from scipy.interpolate import RectBivariateSpline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotate_cube(incube, angle):\n",
    "    \"\"\"Rotate input cube by angle\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    incube : ImageHDU\n",
    "    angle : float\n",
    "        angle in degrees, counterclockwise from positive x-axis\n",
    "    \"\"\"\n",
    "    rangle = angle * np.pi / 180     # degrees to radians\n",
    "    \n",
    "    wcs_orig = WCS(incube.header).sub(2)\n",
    "    wcs_orig.wcs.ctype = [\"LINEAR\", \"LINEAR\"]   # avoids discontinuity around RA=0 degrees\n",
    "    wcs_rot = WCS(incube.header).sub(2)\n",
    "    wcs_rot.wcs.ctype = [\"LINEAR\", \"LINEAR\"]\n",
    "    \n",
    "    wcs_rot.wcs.pc = [[np.cos(rangle), np.sin(rangle)],\n",
    "                     [-np.sin(rangle), np.cos(rangle)]]\n",
    "    i_orig = np.arange(incube.header['NAXIS1'])\n",
    "    j_orig = np.arange(incube.header['NAXIS2'])\n",
    "    x_orig, y_orig = wcs_orig.all_pix2world(i_orig, j_orig, 0)\n",
    "    x_orig[x_orig >180] -= 360.     # RA continuous across 0\n",
    "    i_rot, j_rot = np.meshgrid(i_orig, j_orig)\n",
    "    x_rot, y_rot = wcs_rot.all_pix2world(i_rot, j_rot, 0)\n",
    "    \n",
    "    for iplane, plane in enumerate(incube.data):\n",
    "        interp = RectBivariateSpline(y_orig, x_orig, plane, kx=1, ky=1)\n",
    "        incube.data[iplane,] = interp(y_rot, x_rot, grid=False)\n",
    "    incube.header['ANGLE'] = angle, \"slit rotation [deg]\"\n",
    "    return incube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function is applied to the data cubes before creating the `Source` object. We choose an angle of 45 degrees. The simulation the proceeds as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle = 45\n",
    "rotsources = {}\n",
    "for model, fitsfile in fitsfiles.items():\n",
    "    with fits.open(fitsfile) as hdul:     \n",
    "        hdul[0].header['CRVAL1'] = 0.\n",
    "        hdul[0].header['CRVAL2'] = 0.\n",
    "        hdul[0].header['CDELT1'] *= scale_cdelt\n",
    "        hdul[0].header['CDELT2'] *= scale_cdelt\n",
    "        hdul[0].header['BUNIT'] = hdul[0].header['UNITS']\n",
    "        hdul[0].data[hdul[0].data < 0] = 0\n",
    "        hdul[0] = rotate_cube(hdul[0], angle)\n",
    "        rotsources[model] = sim.Source(cube=hdul)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotresults = {}\n",
    "for model, src in rotsources.items():\n",
    "    print(f'Observing model \"{model}\"...')\n",
    "    metis.observe(src, update=True)\n",
    "    rotresults[model] = metis.readout(detector_readout_mode='auto')[0]\n",
    "    print(\"-----\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure(figsize=(15, 8))\n",
    "# for i, (model, result) in enumerate(rotresults.items()):\n",
    "#     plt.subplot(1, 3, i+1)\n",
    "#     plt.imshow((result[1].data - bgresult[1].data)[jlim[0]:jlim[1], ilim[0]:ilim[1]], \n",
    "#                origin='lower', norm=LogNorm(),\n",
    "#                extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
    "#     plt.xlabel(r\"position along slit [arcsec]\")\n",
    "#     plt.ylabel(r\"Wavelength [$\\mu$m]\")\n",
    "#     plt.title(model + \", rotated by \" + str(angle) + \" degrees\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i, (model, result) in enumerate(rotresults.items()):\n",
    "#     outfile = Path(fitsfiles[model]).stem + \"-rot_\" + str(angle) + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n",
    "#     result.writeto(outfile, overwrite=True)\n",
    "#     print(fitsfiles[model], \"--->\", outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare the rotated spectrum for the `emptycav` model to the spectrum without rotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6, 12))\n",
    "\n",
    "# ax1.plot(det_xi, (results['emptycav'][1].data - bgresult[1].data)[500, :], label=\"no rotation\")\n",
    "# ax1.plot(det_xi, (rotresults['emptycav'][1].data - bgresult[1].data)[500, :], label=f\"rotation {angle} deg\")\n",
    "# ax1.set_xlim(xlim[0], xlim[-1])\n",
    "# ax1.set_ylim(2e4, 1e9)\n",
    "# ax1.set_xlabel(\"arcsec\")\n",
    "# ax1.semilogy()\n",
    "# ax1.legend()\n",
    "\n",
    "# fig.subplots_adjust(left=0.1)\n",
    "# ax2.imshow(rotsources['emptycav'].cube_fields[0].data.sum(axis=0) + 1e-14, norm=LogNorm(vmin=1e-4, vmax=10),\n",
    "#           extent=(x_lim[0], x_lim[-1], y_lim[0], y_lim[-1]))\n",
    "# ax2.set_xlabel(\"arcsec\")\n",
    "# ax2.set_ylabel(\"arcsec\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flux calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pixel values in the detector images give electrons accumulated over the exposure time. To get back to physical units, e.g. Jy, one has to perform a flux calibration. As in real observations, we do this here with the observational of a standard star. For simplicity, we use a star with a spectrum that is constant at $f_{\\nu} = 1\\,\\mathrm{Jy}$. We observe it with the identical setup as the science targets above, except that the exposure time is reduced to 1 second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std = sim.source.source_templates.star(flux=1 * u.Jy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis.observe(std, update=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_exptime = 1  # second\n",
    "std_result = metis.readout(exptime=std_exptime, detector_readout_mode='auto')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_bgsub = std_result[1].data - bgresult[1].data / exptime\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(std_bgsub[250:1750, 800:1250], origin='lower', extent=(800, 1250, 250, 1750))\n",
    "plt.xlabel('pixel')\n",
    "plt.ylabel('pixel')\n",
    "plt.title(r\"background-subtracted standard exposure ($T_\\mathrm{exp} = 1\\,\\mathrm{s}$)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# xmin, xmax = 975, 1075\n",
    "# xaxis = np.arange(xmin, xmax)\n",
    "# y_1, y_2 = 200, 1800\n",
    "# lam_1, lam_2 = det_wave[y_1], det_wave[y_2]\n",
    "\n",
    "# plt.figure(figsize=(10, 5))\n",
    "# plt.plot(xaxis, std_bgsub[200, xmin:xmax], label=fr\"$\\lambda = {lam_1:.2f}\\,\\mu\\mathrm{{m}}$\")\n",
    "# plt.plot(xaxis, std_bgsub[1800, xmin:xmax], label=fr\"$\\lambda = {lam_2:.2f}\\,\\mu\\mathrm{{m}}$\")\n",
    "# plt.plot(xaxis, std_bgsub[:, xmin:xmax].mean(axis=0), label='average')\n",
    "# plt.vlines((1024 - 10, 1024 + 10), 1, 1e5, colors='black', linestyles='dashed', label='extraction aperture')\n",
    "# plt.xlabel(\"pixel\")\n",
    "# plt.ylabel(\"e/s\")\n",
    "# plt.semilogy()\n",
    "# plt.legend()\n",
    "# plt.title(\"spatial cuts through standard spectrum\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total flux density of 1 Jy from the star is spread over the point spread function, which extends far from the central position. It is cut by the finite aperture of the spectroscopic slit, and will be cut further by the finite extent of the extraction aperture over which we will sum the two-dimensional spectrum. Most of the flux is contained in the core of the PSF; from the figure, we can take an aperture size in the spatial direction of 20 pixels - outside this aperture, the spectrum is dominated by noise.  \n",
    "The slit that we used is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: This used to return \"C-38_1\". The headers are made ESO\n",
    "#       Compliant, and not all headers have been replaced.\n",
    "#       See https://github.com/AstarVienna/irdb/pull/146\n",
    "# std_result[1].header['SLIT']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with a width of 38.1 mas, corresponding to 7 pixels. As this notebook focuses on the simulator rather than the data reduction, we neglect slit and aperture losses and assume that the signal integrated over 20 pixels in the spatial direction corresponds to the input flux density of 1 Jy. The error is on the order of 2 per cent. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hwidth = 10     # pixels, half width of extraction aperture\n",
    "std_1d = std_bgsub[:, 1024-hwidth:1024+hwidth].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.plot(det_wave, std_1d)\n",
    "# plt.xlabel('Wavelength (um)')\n",
    "# exptime = std_result[1].header['EXPTIME']\n",
    "# plt.ylabel('electrons/s')\n",
    "# plt.title(r\"Extracted standard spectrum ($T_{\\mathrm{exp}} = 1\\,\\mathrm{s}$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct the flux conversion curve that takes the observed spectrum to the calibrated spectrum in Jy. Note that there are strong excursions at the edge of the spectroscopic filter and at deep atmospheric absorption features. The object flux at these wavelengths is not recoverable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# flux_conv = 1 * u.Jy / (std_1d * u.electron)\n",
    "# plt.plot(det_wave, flux_conv)\n",
    "# plt.ylim(0, 1e-5)\n",
    "# plt.xlabel(\"Wavelength (um)\")\n",
    "# plt.ylabel(\"Jy / (e/s)\")\n",
    "# plt.title(\"Flux conversion\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We apply the flux calibration to the two-dimensional spectra from above, taking into account the different exposure time, and display them. Note how the PAH feature at around 3.3 um now becomes apparent in the disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# emptycav_Jy     = flux_conv[:, None] * (results['emptycav'][1].data - bgresult[1].data) / exptime\n",
    "# emptycav_Jy_rot = flux_conv[:, None] * (rotresults['emptycav'][1].data - bgresult[1].data) / exptime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 12))\n",
    "# ax1.imshow(emptycav_Jy[jlim[0]:jlim[1], ilim[0]:ilim[1]].value, origin='lower', \n",
    "#            norm=LogNorm(vmin=1e-5, vmax=1), extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
    "# ax1.set_xlabel(r\"position along slit [arcsec]\")\n",
    "# ax1.set_ylabel(r\"Wavelength [$\\mu$m]\") \n",
    "# ax1.set_title(model + \", no rotation\")\n",
    "\n",
    "# img = ax2.imshow(emptycav_Jy_rot[jlim[0]:jlim[1], ilim[0]:ilim[1]].value, origin='lower',\n",
    "#                  norm=LogNorm(vmin=1e-5, vmax=1), extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
    "# ax2.set_xlabel(r\"position along slit [arcsec]\")\n",
    "# ax2.set_ylabel(r\"Wavelength [$\\mu$m]\") \n",
    "# ax2.set_title(model + \", rotated by \" + str(angle) + \" degrees\");\n",
    "\n",
    "# fig.subplots_adjust(right=0.8)\n",
    "# cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "# fig.colorbar(img, cax=cbar_ax, label='Jy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we perform background subtraction and flux calibration on all simulations. The data are now normalised to an exposure time of 1 second and therefore have units 'electrons per second'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i, (model, result) in enumerate(results.items()):\n",
    "#     result[1].data = flux_conv[:, None].value * (result[1].data - bgresult[1].data) / exptime\n",
    "#     result[1].header['BUNIT'] = \"e/s\"\n",
    "#     outfile = (Path(fitsfiles[model]).stem + \"-simulated_LSS_L-bgsub_fluxcal\" \n",
    "#                + Path(fitsfiles[model]).suffix)\n",
    "#     result.writeto(outfile, overwrite=True)\n",
    "#     print(\"--->\", outfile)\n",
    "\n",
    "# for i, (model, result) in enumerate(rotresults.items()):\n",
    "#     result[1].data = flux_conv[:, None].value * (result[1].data - bgresult[1].data) / exptime\n",
    "#     result[1].header['BUNIT'] = \"e/s\"\n",
    "#     outfile = (Path(fitsfiles[model]).stem + \"-rot_\" + str(angle) \n",
    "#                + \"-simulated_LSS_L-bgsub_fluxcal\" + Path(fitsfiles[model]).suffix)\n",
    "#     result.writeto(outfile, overwrite=True)\n",
    "#     print(\"--->\", outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison to input spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_extract = results['emptycav'][1].data[:, 1018:1030].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_extract = sources['emptycav'].cube_fields[0].data[:, 88:91, 87:92].sum(axis=(1, 2))\n",
    "src_wave = sources['emptycav'].cube_fields[0].wave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure(figsize=(12, 6))\n",
    "# plt.plot(det_wave, result_extract, label=\"flux-calibrated simulated spectrum\")\n",
    "# plt.plot(src_wave, src_extract, label=\"input spectrum\")\n",
    "# plt.legend()\n",
    "# plt.xlim(3.1, 4.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}