{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c4e9807e",
   "metadata": {},
   "source": [
    "# `n_reads` and exposure-time-for-SNR inverses\n",
    "\n",
    "This notebook demonstrates two related `wcc_etc.Simulation` capabilities:\n",
    "\n",
    "- **`n_reads`** — the number of detector reads (frames) that span the *total* exposure `t`,\n",
    "  so each frame integrates `t / n_reads`. This models CDS / Fowler-style multiple sampling:\n",
    "  coadding more reads reduces the effective read-noise penalty per unit time, pushing the SNR\n",
    "  toward the photon-limited ceiling. Because the read is **per frame**, splitting a fixed total\n",
    "  time across more reads also lowers the per-frame peak charge — which can keep a bright star\n",
    "  under full-well and flip `is_saturated` from `True` to `False`.\n",
    "\n",
    "- **Exposure-time inverses** — `get_exptime_for_snr(snr)` (analytic) and\n",
    "  `get_image_exptime_for_snr(snr, ...)` (PSF-aware) answer *how long* to integrate to reach a\n",
    "  target SNR. They invert `get_snr` / `get_image_snr`, and we verify the round trip below.\n",
    "\n",
    "`n_reads` is a settable `Simulation` parameter (default 1, validated `>= 1`), updatable via\n",
    "`sim.update(n_reads=N)` or passable per call as `n_reads=`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cd429ca",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Import the package, build a G5V source scene over a zodiacal background, and construct a\n",
    "`Simulation` from the Sony `r` sensor+filter combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e330481a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-06T21:06:33.493747Z",
     "iopub.status.busy": "2026-06-06T21:06:33.493590Z",
     "iopub.status.idle": "2026-06-06T21:06:37.692472Z",
     "shell.execute_reply": "2026-06-06T21:06:37.692227Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import wcc_etc\n",
    "\n",
    "wcc_etc.set_wcc_style()\n",
    "\n",
    "scene = wcc_etc.get_scene(\n",
    "    name=\"G5V\",\n",
    "    mag=15,\n",
    "    background=\"zodi\",\n",
    "    bandpass=\"johnson_r\",\n",
    "    background_prop={\"bandpass\": \"johnson_r\", \"mag\": 22.5},\n",
    ")\n",
    "sim = wcc_etc.Simulation.from_sensor_and_scene(\"sony:r\", scene)\n",
    "print(\"wcc_etc loaded; sim built. default n_reads =\", sim.meta.get(\"n_reads\", 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "307eb6ef",
   "metadata": {},
   "source": [
    "## 1. Analytic round-trip\n",
    "\n",
    "Pick a target SNR, compute the required time with `get_exptime_for_snr`, then feed it back into\n",
    "`get_snr_airy` (the analytic forward that matches the analytic `get_exptime_for_snr`) and confirm we recover the target (within ~1%)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "505236a5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-06T21:06:37.694016Z",
     "iopub.status.busy": "2026-06-06T21:06:37.693899Z",
     "iopub.status.idle": "2026-06-06T21:06:37.993606Z",
     "shell.execute_reply": "2026-06-06T21:06:37.993350Z"
    }
   },
   "outputs": [],
   "source": [
    "snr_target = 100.0\n",
    "t = sim.get_exptime_for_snr(snr_target)  # Quantity in seconds\n",
    "import warnings\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\", DeprecationWarning)\n",
    "    snr_back = float(sim.get_snr_airy(t).value)\n",
    "\n",
    "print(f\"target SNR        = {snr_target:.2f}\")\n",
    "print(f\"exptime for SNR   = {t.value:.3f} s\")\n",
    "print(f\"get_snr_airy(t)   = {snr_back:.4f}\")\n",
    "\n",
    "rel_err = abs(snr_back - snr_target) / snr_target\n",
    "print(f\"relative error    = {rel_err:.2e}\")\n",
    "assert rel_err < 0.01, \"analytic round-trip disagreement\"\n",
    "print(\"analytic round-trip OK (<1%)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "147f1e17",
   "metadata": {},
   "source": [
    "## 2. PSF-aware round-trip\n",
    "\n",
    "The PSF-aware path renders the PSF and solves per aperture. `get_image_exptime_for_snr` returns\n",
    "a dict; we feed `res['time_s']` back into `get_image_snr` and confirm the recovered SNR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "385c27be",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-06T21:06:37.995731Z",
     "iopub.status.busy": "2026-06-06T21:06:37.995589Z",
     "iopub.status.idle": "2026-06-06T21:06:38.108453Z",
     "shell.execute_reply": "2026-06-06T21:06:38.108177Z"
    }
   },
   "outputs": [],
   "source": [
    "res = sim.get_image_exptime_for_snr(snr_target)\n",
    "snr_img_back = sim.get_image_snr(time=res[\"time_s\"])[\"snr\"]\n",
    "\n",
    "print(f\"target SNR          = {snr_target:.2f}\")\n",
    "print(f\"exptime for SNR     = {res['time_s']:.3f} s\")\n",
    "print(f\"aperture radius     = {res['r_aper_mas']:.1f} mas\")\n",
    "print(f\"enclosed fraction   = {res['enclosed_fraction']:.4f}\")\n",
    "print(f\"n_pix               = {res['n_pix']}\")\n",
    "print(f\"get_image_snr(t)    = {snr_img_back:.4f}\")\n",
    "\n",
    "rel_err_img = abs(snr_img_back - snr_target) / snr_target\n",
    "print(f\"relative error      = {rel_err_img:.2e}\")\n",
    "assert rel_err_img < 0.01, \"PSF-aware round-trip disagreement\"\n",
    "print(\"PSF-aware round-trip OK (<1%)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af34f40d",
   "metadata": {},
   "source": [
    "## 3. SNR vs `n_reads` at fixed exposure time\n",
    "\n",
    "At a fixed *total* exposure time, more reads incur the read noise more times. In the\n",
    "read-noise-affected regime SNR falls slightly with `n_reads`; in the photon-limited regime the\n",
    "curve flattens toward the photon ceiling. We use a fainter source so the read-noise term is\n",
    "visible, and overlay the `n_reads=1` (photon-limit-approaching) value as the ceiling reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf9e59ec",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-06T21:06:38.109693Z",
     "iopub.status.busy": "2026-06-06T21:06:38.109598Z",
     "iopub.status.idle": "2026-06-06T21:06:39.739217Z",
     "shell.execute_reply": "2026-06-06T21:06:39.738959Z"
    }
   },
   "outputs": [],
   "source": [
    "scene_faint = wcc_etc.get_scene(\n",
    "    name=\"G5V\",\n",
    "    mag=19,\n",
    "    background=\"zodi\",\n",
    "    bandpass=\"johnson_r\",\n",
    "    background_prop={\"bandpass\": \"johnson_r\", \"mag\": 22.5},\n",
    ")\n",
    "sim_faint = wcc_etc.Simulation.from_sensor_and_scene(\"sony:r\", scene_faint)\n",
    "\n",
    "fixed_t = 60.0\n",
    "reads = [1, 2, 4, 8, 16]\n",
    "snr_vs_reads = []\n",
    "for N in reads:\n",
    "    sim_faint.update(n_reads=N)\n",
    "    snr_vs_reads.append(float(sim_faint.get_snr(fixed_t)[\"snr\"]))\n",
    "\n",
    "ceiling = snr_vs_reads[0]  # n_reads=1 -> least read-noise penalty\n",
    "fig, ax = plt.subplots(figsize=(8, 4.5))\n",
    "ax.plot(reads, snr_vs_reads, \"o-\", label=f\"get_snr at {fixed_t:.0f} s\")\n",
    "ax.axhline(ceiling, color=\"grey\", ls=\":\", lw=1, label=\"n_reads=1 (lowest read penalty)\")\n",
    "ax.set_xlabel(\"n_reads\")\n",
    "ax.set_ylabel(\"SNR\")\n",
    "ax.set_title(\"SNR vs n_reads at fixed total time (mag 19, G5V)\")\n",
    "ax.legend()\n",
    "plt.show()\n",
    "\n",
    "for N, s in zip(reads, snr_vs_reads):\n",
    "    print(f\"n_reads={N:>3}:  SNR({fixed_t:.0f} s) = {s:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "940467f2",
   "metadata": {},
   "source": [
    "## 4. Required time for a fixed SNR vs `n_reads`\n",
    "\n",
    "For a fixed target SNR, the required exposure time grows with `n_reads` because the read-noise\n",
    "variance is paid once per frame. Conversely, fewer reads minimize the read-noise penalty and the\n",
    "required time. Plotted on the same fainter source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3101ed17",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-06T21:06:39.740452Z",
     "iopub.status.busy": "2026-06-06T21:06:39.740369Z",
     "iopub.status.idle": "2026-06-06T21:06:40.900704Z",
     "shell.execute_reply": "2026-06-06T21:06:40.900441Z"
    }
   },
   "outputs": [],
   "source": [
    "snr_goal = 20.0\n",
    "t_vs_reads = []\n",
    "for N in reads:\n",
    "    sim_faint.update(n_reads=N)\n",
    "    t_vs_reads.append(sim_faint.get_exptime_for_snr(snr_goal).value)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4.5))\n",
    "ax.plot(reads, t_vs_reads, \"o-\", color=\"C1\")\n",
    "ax.set_xlabel(\"n_reads\")\n",
    "ax.set_ylabel(f\"time for SNR={snr_goal:.0f} [s]\")\n",
    "ax.set_title(\"Required exposure time vs n_reads (mag 19, G5V)\")\n",
    "plt.show()\n",
    "\n",
    "for N, tt in zip(reads, t_vs_reads):\n",
    "    print(f\"n_reads={N:>3}:  t(SNR={snr_goal:.0f}) = {tt:8.3f} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19861025",
   "metadata": {},
   "source": [
    "## 5. Per-frame saturation flip\n",
    "\n",
    "Saturation is evaluated on a single `t / n_reads` frame. A bright source that saturates in one\n",
    "long frame (`n_reads=1`) can be brought back under the ADC full scale by spreading the same total\n",
    "time across more reads — the per-frame peak charge scales as `1 / n_reads`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0386e88a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-06T21:06:40.902095Z",
     "iopub.status.busy": "2026-06-06T21:06:40.902009Z",
     "iopub.status.idle": "2026-06-06T21:06:42.326501Z",
     "shell.execute_reply": "2026-06-06T21:06:42.326261Z"
    }
   },
   "outputs": [],
   "source": [
    "scene_bright = wcc_etc.get_scene(\n",
    "    name=\"G5V\",\n",
    "    mag=17,\n",
    "    background=\"zodi\",\n",
    "    bandpass=\"johnson_r\",\n",
    "    background_prop={\"bandpass\": \"johnson_r\", \"mag\": 22.5},\n",
    ")\n",
    "sim_bright = wcc_etc.Simulation.from_sensor_and_scene(\"sony:r\", scene_bright)\n",
    "\n",
    "total_t = 60.0\n",
    "adc_max = sim_bright.sensor.adc_max\n",
    "adc_max_val = float(adc_max.value if hasattr(adc_max, \"value\") else adc_max)\n",
    "print(f\"adc_max (full scale) = {adc_max_val:.0f} ADU\\n\")\n",
    "\n",
    "sat_reads = [1, 2, 5, 10, 25, 50]\n",
    "flip_at = None\n",
    "print(f\"{'n_reads':>8} {'peak [ADU]':>12} {'saturated':>10}\")\n",
    "for N in sat_reads:\n",
    "    sim_bright.update(n_reads=N)\n",
    "    peak = sim_bright.get_peak_pixel(total_t, units=\"adu\").value\n",
    "    sat = bool(sim_bright.is_saturated(total_t))\n",
    "    if flip_at is None and not sat:\n",
    "        flip_at = N\n",
    "    print(f\"{N:>8} {peak:>12.0f} {str(sat):>10}\")\n",
    "\n",
    "print(\n",
    "    f\"\\nsaturation flips True -> False at n_reads = {flip_at}  \"\n",
    "    f\"(total time {total_t:.0f} s, mag 17 G5V)\"\n",
    ")\n",
    "assert flip_at is not None, \"expected a saturation flip across the n_reads range\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef8a125a",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- `get_exptime_for_snr` (analytic) and `get_image_exptime_for_snr` (PSF-aware) invert the SNR\n",
    "  calculation; both round-trip back through `get_snr` / `get_image_snr` to the target within ~1%.\n",
    "- `n_reads` coadds frames over the total exposure: at fixed time, more reads add the read-noise\n",
    "  penalty more times (SNR drops toward / sits at the photon ceiling); for a fixed target SNR the\n",
    "  required time grows with `n_reads`.\n",
    "- Saturation is **per frame** (`t / n_reads`): increasing `n_reads` lowers the per-frame peak\n",
    "  (∝ 1/n_reads) and can flip `is_saturated` from `True` to `False`.\n",
    "- `n_reads=1` reproduces the original single-frame behaviour."
   ]
  }
 ],
 "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.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}