{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "44ef49d8",
   "metadata": {},
   "source": [
    "# Calibration and error propagation\n",
    "\n",
    "Previously, we compared `phoptic` to [`Source-Extractor`](https://sextractor.readthedocs.io/en/latest/index.html) and showed that the two programs produce [very similar results](sextractor_comparison.ipynb). However, we neglected image calibrations in that example to keep things simple.\n",
    "\n",
    "When performing image calibrations, `phoptic` will automatically propagate the errors; in this notebook, I will prove that `phoptic` does this correctly. Additionally, I will also verify that calibrations applied by `phoptic` are correct. \n",
    "\n",
    "## Generating data\n",
    "\n",
    "First, we need to generate some data that contain only noise. We will include typical values for the bias and dark current, and simulate long exposures to allow for a significant dark noise contribution. We will also generate flat-field images, which introduce some additional noise. These images will be generated as if they were produced by OPTICAM-MX, but we'll only simulate data from a single camera for simplicity.\n",
    "\n",
    "Let's define the properties of these generated data: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdc354ea",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:07.689821Z",
     "iopub.status.busy": "2026-01-23T11:19:07.689681Z",
     "iopub.status.idle": "2026-01-23T11:19:07.691953Z",
     "shell.execute_reply": "2026-01-23T11:19:07.691742Z"
    }
   },
   "outputs": [],
   "source": [
    "X, Y = 512, 512  # image dimensions\n",
    "\n",
    "gain = 1.  # gain of OPTICAM-MX\n",
    "read_noise = 1.1  # nominal read noise of OPTICAM-MX\n",
    "dark_curr = 0.1  # nominal dark current of OPTICAM-MX\n",
    "texp = 30.  # exposure time of science images (in seconds)\n",
    "flat_texp = 3.  # exposure time of flat-field images (in seconds)\n",
    "fltr = 'g'  # image filter\n",
    "\n",
    "binning = f'{2048 // X}x{2048 // Y}'  # binning mode\n",
    "\n",
    "bias_value = 400.  # typical bias according to Handbook of CCD Astronomy, Howell\n",
    "\n",
    "flat_level = 40000.  # mean flux of flats\n",
    "\n",
    "# dark noise contributions\n",
    "dark_signal = dark_curr * texp\n",
    "flat_dark_signal = dark_curr * flat_texp\n",
    "\n",
    "n_calibration_images = 30  # number of images to generate for each calibration step\n",
    "\n",
    "science_flux = 1000.  # mean flux in science images (before noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1eec252",
   "metadata": {},
   "source": [
    "We can compute the analytical variances for each calibration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecb6071d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:07.692727Z",
     "iopub.status.busy": "2026-01-23T11:19:07.692651Z",
     "iopub.status.idle": "2026-01-23T11:19:07.719467Z",
     "shell.execute_reply": "2026-01-23T11:19:07.719253Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "true_bias_var = read_noise**2 / n_calibration_images\n",
    "\n",
    "true_dark_var = np.pi / (2 * n_calibration_images) * (dark_signal + read_noise**2 + true_bias_var)\n",
    "\n",
    "true_flat_var = np.pi / (2 * n_calibration_images) * (flat_level + read_noise**2 + true_bias_var + flat_dark_signal)\n",
    "effective_flat_var = true_flat_var * science_flux**2 / flat_level**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70542fa0",
   "metadata": {},
   "source": [
    "I have included $\\frac{\\pi}{2 N}$ factors in the dark and flat variances since `phoptic` computes master darks and master flats using the median (to remove outliers due to, e.g., cosmic rays), while the master bias is computed using the mean. Note that the dark and flat-field variances include a bias contribution, and the flat-field variance further includes a (small) dark-noise contribution. The flat variance also has an additional term because it acts multiplicatively, whereas the bias and dark corrections are additive.\n",
    "\n",
    "Now let's create the directories for storing our simulated data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af934725",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:07.720632Z",
     "iopub.status.busy": "2026-01-23T11:19:07.720510Z",
     "iopub.status.idle": "2026-01-23T11:19:07.723094Z",
     "shell.execute_reply": "2026-01-23T11:19:07.722871Z"
    }
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "out_dir = Path('out')\n",
    "\n",
    "if not out_dir.joinpath('biases').is_dir():\n",
    "    out_dir.joinpath('biases').mkdir(parents=True)\n",
    "\n",
    "if not out_dir.joinpath('darks').is_dir():\n",
    "    out_dir.joinpath('darks').mkdir(parents=True)\n",
    "\n",
    "if not out_dir.joinpath('flats').is_dir():\n",
    "    out_dir.joinpath('flats').mkdir(parents=True)\n",
    "\n",
    "if not out_dir.joinpath('flat_darks').is_dir():\n",
    "    out_dir.joinpath('flat_darks').mkdir(parents=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97a40484",
   "metadata": {},
   "source": [
    "We'll use `numpy`'s `default_rng()` Generator to generate random noise. For reproducibility, I'll pass a fixed seed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33f7ef94",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:07.724473Z",
     "iopub.status.busy": "2026-01-23T11:19:07.724100Z",
     "iopub.status.idle": "2026-01-23T11:19:07.728829Z",
     "shell.execute_reply": "2026-01-23T11:19:07.728629Z"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8aeb4644",
   "metadata": {},
   "source": [
    "Let's now generate our calibration images.\n",
    "\n",
    "### Biases\n",
    "\n",
    "Bias images should be taken with the telescope shutter closed and have an exposure time of 0.0 seconds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5186356",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:07.730024Z",
     "iopub.status.busy": "2026-01-23T11:19:07.729907Z",
     "iopub.status.idle": "2026-01-23T11:19:11.968489Z",
     "shell.execute_reply": "2026-01-23T11:19:11.968236Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "for i in range(n_calibration_images):\n",
    "    bias = rng.normal(bias_value, read_noise, size=(X, Y))\n",
    "    \n",
    "    hdu = fits.PrimaryHDU(bias)\n",
    "    \n",
    "    hdu.header['EXPOSURE'] = 0.0\n",
    "    hdu.header[\"BINNING\"] = binning\n",
    "    hdu.header[\"GAIN\"] = gain\n",
    "    hdu.header['FILTER'] = fltr\n",
    "    hdu.header['INSTRUME'] = 'OPTICAM'\n",
    "    \n",
    "    path = out_dir / 'biases' / f'image_{i}.fits.gz'\n",
    "    if not path.is_file():\n",
    "        hdu.writeto(path, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a7bbb48",
   "metadata": {},
   "source": [
    "Let's pass these bias images to a `BiasCorrector`: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08f054a1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:11.969743Z",
     "iopub.status.busy": "2026-01-23T11:19:11.969625Z",
     "iopub.status.idle": "2026-01-23T11:19:13.662467Z",
     "shell.execute_reply": "2026-01-23T11:19:13.662244Z"
    }
   },
   "outputs": [],
   "source": [
    "import phoptic\n",
    "\n",
    "bias_corr = phoptic.BiasCorrector(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'biases',\n",
    "    instrument=phoptic.OPTICAM_MX(),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a70369c",
   "metadata": {},
   "source": [
    "### Darks\n",
    "\n",
    "Dark images should again be taken with the telescope shutter closed, but now have an exposure time equal to that of the science images. Like all images, dark images also include a bias offset. We'll generate one set of darks for the science images, and another set of darks, using a different exposure time, for the flat-field images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c644f00",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:13.663470Z",
     "iopub.status.busy": "2026-01-23T11:19:13.663299Z",
     "iopub.status.idle": "2026-01-23T11:19:22.597695Z",
     "shell.execute_reply": "2026-01-23T11:19:22.597445Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "for i in range(n_calibration_images):\n",
    "    dark = rng.poisson(dark_signal, size=(X, Y)).astype(np.float64)  # dark noise\n",
    "    dark += rng.normal(bias_value, read_noise, size=(X, Y))  # add bias\n",
    "    \n",
    "    hdu = fits.PrimaryHDU(dark)\n",
    "    \n",
    "    hdu.header['EXPOSURE'] = texp\n",
    "    hdu.header[\"BINNING\"] = binning\n",
    "    hdu.header[\"GAIN\"] = gain\n",
    "    hdu.header['FILTER'] = fltr\n",
    "    hdu.header['INSTRUME'] = 'OPTICAM'\n",
    "    \n",
    "    path = out_dir / 'darks' / f'image_{i}.fits.gz'\n",
    "    if not path.is_file():\n",
    "        hdu.writeto(path, overwrite=True)\n",
    "    \n",
    "    dark = rng.poisson(flat_dark_signal, size=(X, Y)).astype(np.float64)  # flat-field dark noise\n",
    "    dark += rng.normal(bias_value, read_noise, size=(X, Y))  # add bias\n",
    "    \n",
    "    hdu = fits.PrimaryHDU(dark)\n",
    "    \n",
    "    hdu.header['EXPOSURE'] = flat_texp\n",
    "    hdu.header[\"BINNING\"] = binning\n",
    "    hdu.header[\"GAIN\"] = gain\n",
    "    hdu.header['FILTER'] = fltr\n",
    "    hdu.header['INSTRUME'] = 'OPTICAM'\n",
    "    \n",
    "    path = out_dir / 'flat_darks' / f'image_{i}.fits.gz'\n",
    "    if not path.is_file():\n",
    "        hdu.writeto(path, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9b6728e",
   "metadata": {},
   "source": [
    "Let's pass these dark images to `DarkNoiseCorrector`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cecdbf6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:22.598739Z",
     "iopub.status.busy": "2026-01-23T11:19:22.598657Z",
     "iopub.status.idle": "2026-01-23T11:19:23.022459Z",
     "shell.execute_reply": "2026-01-23T11:19:23.022245Z"
    }
   },
   "outputs": [],
   "source": [
    "# science image corrector\n",
    "dark_corr = phoptic.DarkNoiseCorrector(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'darks',\n",
    "    instrument=phoptic.OPTICAM_MX(),\n",
    "    bias_corrector=bias_corr,\n",
    "    )\n",
    "\n",
    "# flat-field image corrector\n",
    "flat_dark_corr = phoptic.DarkNoiseCorrector(\n",
    "    out_directory=out_dir / 'reduced_flat_dark',\n",
    "    data_directory=out_dir / 'flat_darks',\n",
    "    instrument=phoptic.OPTICAM_MX(),\n",
    "    bias_corrector=bias_corr,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bed2f55",
   "metadata": {},
   "source": [
    "Note that I have also passed our previously defined `BiasCorrector` instance to both dark noise correctors to automatically bias-correct the dark images.\n",
    "\n",
    "### Flats\n",
    "\n",
    "Flat-field images should be taken using a bright, uniform field with an exposure time that prevents the detector from saturating. Flat-field images contain a bias offset **and** a dark noise contribution: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1889d203",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:23.023281Z",
     "iopub.status.busy": "2026-01-23T11:19:23.023190Z",
     "iopub.status.idle": "2026-01-23T11:19:28.450169Z",
     "shell.execute_reply": "2026-01-23T11:19:28.449893Z"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(n_calibration_images):\n",
    "    flat = rng.poisson(flat_level, size=(X, Y)).astype(np.float64)\n",
    "    flat += rng.normal(bias_value, read_noise, size=(X, Y))  # add bias\n",
    "    flat += rng.poisson(flat_dark_signal, size=(X, Y)).astype(np.float64)  # add dark noise\n",
    "    \n",
    "    hdu = fits.PrimaryHDU(flat)\n",
    "    \n",
    "    hdu.header['EXPOSURE'] = flat_texp\n",
    "    hdu.header[\"BINNING\"] = binning\n",
    "    hdu.header[\"GAIN\"] = gain\n",
    "    hdu.header['FILTER'] = fltr\n",
    "    hdu.header['INSTRUME'] = 'OPTICAM'\n",
    "    \n",
    "    path = out_dir / 'flats' / f'image_{i}.fits.gz'\n",
    "    if not path.is_file():\n",
    "        hdu.writeto(path, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbe97b5f",
   "metadata": {},
   "source": [
    "Let's pass these flats to a `FlatFieldCorrector`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2d5b953",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:28.451299Z",
     "iopub.status.busy": "2026-01-23T11:19:28.451191Z",
     "iopub.status.idle": "2026-01-23T11:19:28.662995Z",
     "shell.execute_reply": "2026-01-23T11:19:28.662692Z"
    }
   },
   "outputs": [],
   "source": [
    "flat_corr = phoptic.FlatFieldCorrector(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'flats',\n",
    "    instrument=phoptic.OPTICAM_MX(),\n",
    "    bias_corrector=bias_corr,\n",
    "    dark_corrector=flat_dark_corr,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d65e95a",
   "metadata": {},
   "source": [
    "Note that I have passed our previously define `BiasCorrector` instance **and** the `DarkNoiseCorrector` instance that I created to subtract the dark noise from the flat-field images.\n",
    "\n",
    "## Propagation\n",
    "\n",
    "Now that we have defined our calibrators, we can use them to calibrate some synthetic images and evaluate how accurate `phoptic`'s error propagation is. `phoptic` uses a helper function, `phoptic.utils.helpers.propagate_errors()` to propagate errors, so this is what we will use here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d99176d4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:28.664156Z",
     "iopub.status.busy": "2026-01-23T11:19:28.663961Z",
     "iopub.status.idle": "2026-01-23T11:19:33.379905Z",
     "shell.execute_reply": "2026-01-23T11:19:33.379642Z"
    }
   },
   "outputs": [],
   "source": [
    "from phoptic.utils.helpers import combine_variances\n",
    "\n",
    "calibrated_images = []\n",
    "bias_vars = []\n",
    "dark_vars = []\n",
    "flat_vars = []\n",
    "errs = []\n",
    "\n",
    "noiseless_image = np.ones((X, Y)) * science_flux\n",
    "\n",
    "for i in range(100):\n",
    "    image = rng.poisson(noiseless_image).astype(np.float64)\n",
    "    image += rng.normal(bias_value, read_noise, size=(X, Y))  # add bias\n",
    "    image += rng.poisson(dark_signal)  # add dark noise\n",
    "    \n",
    "    # apply bias correction\n",
    "    bias_corrected_image, bias_var = bias_corr.correct(image, camera='1')\n",
    "    bias_vars.append(bias_var)\n",
    "    \n",
    "    # apply dark noise correction\n",
    "    dark_corrected_image, dark_var = dark_corr.correct(bias_corrected_image, key='1:g')\n",
    "    dark_vars.append(dark_var)\n",
    "    \n",
    "    # apply flat correction\n",
    "    flat_corrected_image, flat_var = flat_corr.correct(dark_corrected_image, key='1:g')\n",
    "    flat_vars.append(flat_var)\n",
    "    calibrated_images.append(flat_corrected_image)\n",
    "    \n",
    "    # error propagation function used by phoptic\n",
    "    var = combine_variances(\n",
    "        flat_corrected_image,\n",
    "        bias_var=bias_var,\n",
    "        dark_var=dark_var,\n",
    "        flat_var=flat_var,\n",
    "        background_rms=0.,\n",
    "        read_noise=read_noise,\n",
    "        )\n",
    "    errs.append(np.sqrt(var))\n",
    "\n",
    "calibrated_images = np.asarray(calibrated_images)\n",
    "bias_vars = np.asarray(bias_vars)\n",
    "dark_vars = np.asarray(dark_vars)\n",
    "flat_vars = np.asarray(flat_vars)\n",
    "errs = np.asarray(errs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0aaf722",
   "metadata": {},
   "source": [
    "If `phoptic`'s error propagation is correct, then the average propagated error should be close to the empirical error. We can compute the empirical error by taking the standard deviation of the calibrated images. Let's check that the average ratio between the propagated and empirical errors is close to 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70eb3028",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.381286Z",
     "iopub.status.busy": "2026-01-23T11:19:33.381165Z",
     "iopub.status.idle": "2026-01-23T11:19:33.421178Z",
     "shell.execute_reply": "2026-01-23T11:19:33.420936Z"
    }
   },
   "outputs": [],
   "source": [
    "empirical_err = np.std(calibrated_images, axis=0)\n",
    "\n",
    "measured_err = np.mean(errs, axis=0)\n",
    "\n",
    "print(np.mean(empirical_err / measured_err))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25a20e24",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.422284Z",
     "iopub.status.busy": "2026-01-23T11:19:33.422046Z",
     "iopub.status.idle": "2026-01-23T11:19:33.423832Z",
     "shell.execute_reply": "2026-01-23T11:19:33.423634Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# check errors are within 1%\n",
    "assert np.isclose(np.mean(empirical_err / measured_err), 1., rtol=0.01, atol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e1fae7c",
   "metadata": {},
   "source": [
    "That's reassuring! We can also compare each calibrator's variance to its analytical value:\n",
    "\n",
    "#### Biases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a44c8f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.424682Z",
     "iopub.status.busy": "2026-01-23T11:19:33.424533Z",
     "iopub.status.idle": "2026-01-23T11:19:33.443908Z",
     "shell.execute_reply": "2026-01-23T11:19:33.443686Z"
    }
   },
   "outputs": [],
   "source": [
    "print(np.mean(bias_vars / true_bias_var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b65db048",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.444804Z",
     "iopub.status.busy": "2026-01-23T11:19:33.444721Z",
     "iopub.status.idle": "2026-01-23T11:19:33.463217Z",
     "shell.execute_reply": "2026-01-23T11:19:33.462960Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "assert np.isclose(np.mean(bias_vars / true_bias_var), 1., rtol=0.01, atol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56a144f2",
   "metadata": {},
   "source": [
    "#### Darks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc5db80e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.464313Z",
     "iopub.status.busy": "2026-01-23T11:19:33.464223Z",
     "iopub.status.idle": "2026-01-23T11:19:33.484246Z",
     "shell.execute_reply": "2026-01-23T11:19:33.483976Z"
    }
   },
   "outputs": [],
   "source": [
    "print(np.mean(dark_vars / true_dark_var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9f8070d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.485019Z",
     "iopub.status.busy": "2026-01-23T11:19:33.484935Z",
     "iopub.status.idle": "2026-01-23T11:19:33.503412Z",
     "shell.execute_reply": "2026-01-23T11:19:33.503126Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "assert np.isclose(np.mean(dark_vars / true_dark_var), 1., rtol=0.01, atol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d9d61e0",
   "metadata": {},
   "source": [
    "#### Flats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "367b6ffb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.514395Z",
     "iopub.status.busy": "2026-01-23T11:19:33.514289Z",
     "iopub.status.idle": "2026-01-23T11:19:33.532687Z",
     "shell.execute_reply": "2026-01-23T11:19:33.532465Z"
    }
   },
   "outputs": [],
   "source": [
    "print(np.mean(flat_vars / effective_flat_var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "283832e7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.533532Z",
     "iopub.status.busy": "2026-01-23T11:19:33.533454Z",
     "iopub.status.idle": "2026-01-23T11:19:33.552462Z",
     "shell.execute_reply": "2026-01-23T11:19:33.552191Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "assert np.isclose(np.mean(flat_vars / effective_flat_var), 1., rtol=0.01, atol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8280df48",
   "metadata": {},
   "source": [
    "All variances are well within 1% of their analytical values!\n",
    "\n",
    "Let's compare the propagated errors to their analytical/empirical values graphically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c92eae4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.553562Z",
     "iopub.status.busy": "2026-01-23T11:19:33.553481Z",
     "iopub.status.idle": "2026-01-23T11:19:33.760238Z",
     "shell.execute_reply": "2026-01-23T11:19:33.760033Z"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(\n",
    "    nrows=2,\n",
    "    ncols=2,\n",
    "    figsize=(12.8, 12.8),\n",
    "    )\n",
    "\n",
    "axes[0, 0].set_title('Bias')\n",
    "axes[0, 0].hist(np.mean(bias_vars / true_bias_var, axis=0).flatten(),\n",
    "        bins=100,\n",
    "        histtype='step',\n",
    "        color='k',\n",
    "        )\n",
    "\n",
    "axes[0, 1].set_title('Dark')\n",
    "axes[0, 1].hist(np.mean(dark_vars / true_dark_var, axis=0).flatten(),\n",
    "        bins=100,\n",
    "        histtype='step',\n",
    "        color='k',\n",
    "        )\n",
    "\n",
    "axes[1, 0].set_title('Flat')\n",
    "axes[1, 0].hist(np.mean(flat_vars / effective_flat_var, axis=0).flatten(),\n",
    "        bins=100,\n",
    "        histtype='step',\n",
    "        color='k',\n",
    "        )\n",
    "\n",
    "axes[1, 1].set_title('Total')\n",
    "axes[1, 1].hist((empirical_err / measured_err).flatten(),\n",
    "        bins=100,\n",
    "        histtype='step',\n",
    "        color='k',\n",
    "        )\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.axvline(1, c='r', ls='-')\n",
    "    ax.set_xlabel('Measured / expected')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ba1508f",
   "metadata": {},
   "source": [
    "Clearly, `phoptic` is correctly propagating the errors from each calibration. However, this doesn't prove that the calibrations have actually been applied correctly. To verify this, we can check that the mean flux of the calibrated images is close to the expected value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c591828b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.761150Z",
     "iopub.status.busy": "2026-01-23T11:19:33.761063Z",
     "iopub.status.idle": "2026-01-23T11:19:33.778261Z",
     "shell.execute_reply": "2026-01-23T11:19:33.778066Z"
    }
   },
   "outputs": [],
   "source": [
    "print(np.mean(calibrated_images / science_flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f52323de",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.779092Z",
     "iopub.status.busy": "2026-01-23T11:19:33.779014Z",
     "iopub.status.idle": "2026-01-23T11:19:33.796245Z",
     "shell.execute_reply": "2026-01-23T11:19:33.796007Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "assert np.isclose(np.mean(calibrated_images / science_flux), 1., rtol=0.01, atol=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03ec43d8",
   "metadata": {},
   "source": [
    "Another reassuring result! Comparing the calibrated fluxes graphically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75e9d7a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:33.797283Z",
     "iopub.status.busy": "2026-01-23T11:19:33.797197Z",
     "iopub.status.idle": "2026-01-23T11:19:33.980439Z",
     "shell.execute_reply": "2026-01-23T11:19:33.980146Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.hist(\n",
    "    calibrated_images.flatten() / science_flux,\n",
    "    bins=100,\n",
    "    histtype='step',\n",
    "    color='k',\n",
    "    )\n",
    "\n",
    "ax.axvline(1, c='r')\n",
    "ax.set_xlabel('Measured / expected')\n",
    "ax.set_title('Calibrated flux')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbabc8bd",
   "metadata": {},
   "source": [
    "We can see that, unsurprisingly, the calibrated fluxes follow a Gaussian distribution centred on the expected value.\n",
    "\n",
    "That concludes this demonstration of `phoptic`'s calibration error propagation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opticam7",
   "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.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}