{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Auto Exposure\n",
    "\n",
    "This is a setup/test/demonstration notebook for the `AutoExposure` effect in Scopesim. This effect splits the requested total exposure time into NDIT subexposures of integration time DIT such that the maximum counts in a single subexposure does not exceed a certain fill fraction of the detector full well. The final readout is the sum over the NDIT subexposures, i.e. corresponds to the total requested exposure time.\n",
    "\n",
    "The notebook uses the `irdb/METIS` configuration. The observed source is blank sky, except for the last example where a star of 0 mag is used (Vega)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "import scopesim as sim\n",
    "sim.bug_report()\n",
    "\n",
    "# Edit this path if you have a custom install directory, otherwise comment it out.\n",
    "sim.link_irdb(\"../../../../\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you haven't got the instrument packages yet, uncomment the following cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sim.download_packages([\"METIS\", \"ELT\", \"Armazones\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imaging LM-band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"img_lm\"])\n",
    "metis = sim.OpticalTrain(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis.observe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For `AutoExposure` to work the exposure time has to be given explicitely as a parameter to the `readout` method. If this is not done, the default values for `DIT` and `NDIT` will be used. The following is for an exposure time of 1 second. The resulting readout is divided by `NDIT` to produce the average over the `NDIT` subexposures. After application of the gain to convert from ADU to electrons this allows direct comparison to the detector full well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul = metis.readout(exptime=1)[0]\n",
    "gain = outhdul[1].header[\"ESO DET1 CHIP GAIN\"] * u.electron / u.adu\n",
    "full_well = outhdul[1].header[\"ESO DET1 CHIP FULLWELL\"] * u.electron\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same with a much larger exposure time of 1000 seconds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul = metis.readout(exptime = 1000)[0]\n",
    "gain = outhdul[1].header[\"ESO DET1 CHIP GAIN\"] * u.electron / u.adu\n",
    "full_well = outhdul[1].header[\"ESO DET1 CHIP FULLWELL\"] * u.electron\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The desired fill fraction can be changed with the argument `fill_frac`. The default value of 75 per cent is a typical good value that keeps the detector counts within the linear regime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul = metis.readout(exptime = 1000, fill_frac=0.9)[0]\n",
    "gain = outhdul[1].header[\"ESO DET1 CHIP GAIN\"] * u.electron / u.adu\n",
    "full_well = outhdul[1].header[\"ESO DET1 CHIP FULLWELL\"] * u.electron\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imaging N-band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"img_n\"])\n",
    "metis = sim.OpticalTrain(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis.observe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul = metis.readout(exptime=1)[0]\n",
    "gain = outhdul[1].header[\"ESO DET2 CHIP GAIN\"] * u.electron / u.adu\n",
    "full_well = outhdul[1].header[\"ESO DET2 CHIP FULLWELL\"] * u.electron\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Long-slit spectroscopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"lss_l\"],\n",
    "                       properties={\"!OBS.interp_psf\": False})\n",
    "metis = sim.OpticalTrain(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis.observe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul = metis.readout(exptime=3600.)[0]\n",
    "gain = outhdul[1].header[\"ESO DET1 CHIP GAIN\"] * u.electron / u.adu\n",
    "full_well = outhdul[1].header[\"ESO DET1 CHIP FULLWELL\"] * u.electron\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What happens when the source saturates the detector?\n",
    "We take an N-band image of Vega. `DIT` is automatically set to the minimum value supported by the detector, but the centre of the star still saturates the detector. In the final image, the star's profile is capped at the full well of the detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"img_n\"])\n",
    "metis = sim.OpticalTrain(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = sim.source.source_templates.star()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis.observe(src)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul = metis.readout(exptime=1)[0]\n",
    "gain = outhdul[1].header[\"ESO DET2 CHIP GAIN\"] * u.electron / u.adu\n",
    "full_well = outhdul[1].header[\"ESO DET2 CHIP FULLWELL\"] * u.electron\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a cut through the star to show how its peak saturates the detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(outimg[950:1100, 1024])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npix = (outimg >= full_well).sum()\n",
    "print(\"Number of saturated pixels:\", npix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default values for the detector full well in the various modes reflects our current best knowledge of the properties of the actual METIS detectors. These values can be changed as in the following example, but be aware that this makes the simulations unrealistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NB: There's something wrong with this example, please ignore for the time being.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_well = 1000 * metis.cmds[\"!DET.full_well\"] * u.electron\n",
    "outhdul = metis.readout(exptime=1, full_well=full_well)[0]\n",
    "gain = outhdul[1].header[\"ESO DET2 CHIP GAIN\"] * u.electron / u.adu\n",
    "outimg = outhdul[1].data * u.adu * gain\n",
    "fill_frac = outimg.max() / full_well << u.percent\n",
    "\n",
    "print(\"\\nResult\\n======\")\n",
    "print(f\"Maximum value in readout: {outimg.max():7.1f} (per DIT)\")\n",
    "print(f\"Detector full well: {full_well:13.0f}\")\n",
    "print(f\"Fill fraction: {fill_frac:18.1f}\")"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}