{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "44ef49d8",
   "metadata": {},
   "source": [
    "# Calibration and error propagation\n",
    "\n",
    "Previously, we compared `opticam` 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, `opticam` will automatically propagate the errors; in this notebook, I will prove that `opticam` does this correctly. Additionally, I will also verify that calibrations applied by `opticam` 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": 1,
   "id": "fdc354ea",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:00.215066Z",
     "iopub.status.busy": "2026-02-16T18:29:00.214974Z",
     "iopub.status.idle": "2026-02-16T18:29:00.217524Z",
     "shell.execute_reply": "2026-02-16T18:29:00.217115Z"
    }
   },
   "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": 2,
   "id": "ecb6071d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:00.218439Z",
     "iopub.status.busy": "2026-02-16T18:29:00.218339Z",
     "iopub.status.idle": "2026-02-16T18:29:00.242771Z",
     "shell.execute_reply": "2026-02-16T18:29:00.242439Z"
    }
   },
   "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) / flat_level**2\n",
    "effective_flat_var = true_flat_var * science_flux**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70542fa0",
   "metadata": {},
   "source": [
    "I have included $\\frac{\\pi}{2 N}$ factors in the dark and flat variances since `opticam` 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": 3,
   "id": "af934725",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:00.244516Z",
     "iopub.status.busy": "2026-02-16T18:29:00.244388Z",
     "iopub.status.idle": "2026-02-16T18:29:00.246846Z",
     "shell.execute_reply": "2026-02-16T18:29:00.246506Z"
    }
   },
   "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": 4,
   "id": "33f7ef94",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:00.247653Z",
     "iopub.status.busy": "2026-02-16T18:29:00.247468Z",
     "iopub.status.idle": "2026-02-16T18:29:00.252556Z",
     "shell.execute_reply": "2026-02-16T18:29:00.252232Z"
    }
   },
   "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": 5,
   "id": "b5186356",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:00.253414Z",
     "iopub.status.busy": "2026-02-16T18:29:00.253294Z",
     "iopub.status.idle": "2026-02-16T18:29:04.128993Z",
     "shell.execute_reply": "2026-02-16T18:29:04.128568Z"
    }
   },
   "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": 6,
   "id": "08f054a1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:04.130409Z",
     "iopub.status.busy": "2026-02-16T18:29:04.130292Z",
     "iopub.status.idle": "2026-02-16T18:29:06.178438Z",
     "shell.execute_reply": "2026-02-16T18:29:06.178042Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:   0%|                                                                                                                                                     |[00:00<?]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  30%|███████████████████████████████████████████▌                                                                                                     |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  60%|███████████████████████████████████████████████████████████████████████████████████████                                                          |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  90%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌              |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] 30 1 bias images.\n"
     ]
    }
   ],
   "source": [
    "import opticam\n",
    "\n",
    "bias_corr = opticam.BiasCorrector(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'biases',\n",
    "    instrument=opticam.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": 7,
   "id": "4c644f00",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:06.179438Z",
     "iopub.status.busy": "2026-02-16T18:29:06.179291Z",
     "iopub.status.idle": "2026-02-16T18:29:14.202991Z",
     "shell.execute_reply": "2026-02-16T18:29:14.202538Z"
    }
   },
   "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": 8,
   "id": "6cecdbf6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:14.204133Z",
     "iopub.status.busy": "2026-02-16T18:29:14.204055Z",
     "iopub.status.idle": "2026-02-16T18:29:15.221456Z",
     "shell.execute_reply": "2026-02-16T18:29:15.221124Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:   0%|                                                                                                                                                     |[00:00<?]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  33%|████████████████████████████████████████████████▎                                                                                                |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  67%|████████████████████████████████████████████████████████████████████████████████████████████████▋                                                |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] 30 1:g dark images.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:   0%|                                                                                                                                                     |[00:00<?]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  30%|███████████████████████████████████████████▌                                                                                                     |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  60%|███████████████████████████████████████████████████████████████████████████████████████                                                          |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  90%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌              |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] 30 1:g dark images.\n"
     ]
    }
   ],
   "source": [
    "# science image corrector\n",
    "dark_corr = opticam.DarkNoiseCorrector(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'darks',\n",
    "    instrument=opticam.OPTICAM_MX(),\n",
    "    bias_corrector=bias_corr,\n",
    "    )\n",
    "\n",
    "# flat-field image corrector\n",
    "flat_dark_corr = opticam.DarkNoiseCorrector(\n",
    "    out_directory=out_dir / 'reduced_flat_dark',\n",
    "    data_directory=out_dir / 'flat_darks',\n",
    "    instrument=opticam.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": 9,
   "id": "1889d203",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:15.222429Z",
     "iopub.status.busy": "2026-02-16T18:29:15.222344Z",
     "iopub.status.idle": "2026-02-16T18:29:19.991357Z",
     "shell.execute_reply": "2026-02-16T18:29:19.991030Z"
    }
   },
   "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": 10,
   "id": "e2d5b953",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:19.992452Z",
     "iopub.status.busy": "2026-02-16T18:29:19.992371Z",
     "iopub.status.idle": "2026-02-16T18:29:20.486780Z",
     "shell.execute_reply": "2026-02-16T18:29:20.486425Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:   0%|                                                                                                                                                     |[00:00<?]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  33%|████████████████████████████████████████████████▎                                                                                                |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory:  67%|████████████████████████████████████████████████████████████████████████████████████████████████▋                                                |[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] 30 1:g flat-field images.\n"
     ]
    }
   ],
   "source": [
    "flat_corr = opticam.FlatFieldCorrector(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'flats',\n",
    "    instrument=opticam.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 `opticam`'s error propagation is. `opticam` uses a helper function, `opticam.utils.helpers.propagate_errors()` to propagate errors, so this is what we will use here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d99176d4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:20.487611Z",
     "iopub.status.busy": "2026-02-16T18:29:20.487535Z",
     "iopub.status.idle": "2026-02-16T18:29:24.638100Z",
     "shell.execute_reply": "2026-02-16T18:29:24.637676Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] 1 master bias image not found. Attempting to create.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master bias image(s) created.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master bias image(s) saved to out/reduced/master_bias.fits.gz.\n",
      "[OPTICAM] 1:g master dark image not found. Attempting to create.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master dark image(s) created.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master dark image(s) saved to out/reduced/master_darks.fits.gz.\n",
      "[OPTICAM] 1:g master flat-field image not found. Attempting to create.\n",
      "[OPTICAM] 1:g master dark image not found. Attempting to create.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master dark image(s) created.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master dark image(s) saved to out/reduced_flat_dark/master_darks.fits.gz.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master flat-field image(s) created.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[OPTICAM] Master flat-field image(s) saved to out/reduced/master_flats.fits.gz.\n"
     ]
    }
   ],
   "source": [
    "from opticam.utils.helpers import propagate_errors\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 opticam\n",
    "    err = propagate_errors(\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(err)\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 `opticam`'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": 12,
   "id": "70eb3028",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.639281Z",
     "iopub.status.busy": "2026-02-16T18:29:24.639205Z",
     "iopub.status.idle": "2026-02-16T18:29:24.673437Z",
     "shell.execute_reply": "2026-02-16T18:29:24.673061Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9930406244538639\n"
     ]
    }
   ],
   "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": 13,
   "id": "25a20e24",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.674230Z",
     "iopub.status.busy": "2026-02-16T18:29:24.674157Z",
     "iopub.status.idle": "2026-02-16T18:29:24.675727Z",
     "shell.execute_reply": "2026-02-16T18:29:24.675478Z"
    },
    "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": 14,
   "id": "f3a44c8f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.676496Z",
     "iopub.status.busy": "2026-02-16T18:29:24.676428Z",
     "iopub.status.idle": "2026-02-16T18:29:24.692696Z",
     "shell.execute_reply": "2026-02-16T18:29:24.692339Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.999480261881968\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(bias_vars / true_bias_var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b65db048",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.693520Z",
     "iopub.status.busy": "2026-02-16T18:29:24.693443Z",
     "iopub.status.idle": "2026-02-16T18:29:24.709638Z",
     "shell.execute_reply": "2026-02-16T18:29:24.709298Z"
    },
    "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": 16,
   "id": "cc5db80e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.710634Z",
     "iopub.status.busy": "2026-02-16T18:29:24.710552Z",
     "iopub.status.idle": "2026-02-16T18:29:24.726792Z",
     "shell.execute_reply": "2026-02-16T18:29:24.726459Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9997573182415519\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(dark_vars / true_dark_var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c9f8070d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.727695Z",
     "iopub.status.busy": "2026-02-16T18:29:24.727613Z",
     "iopub.status.idle": "2026-02-16T18:29:24.743702Z",
     "shell.execute_reply": "2026-02-16T18:29:24.743328Z"
    },
    "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": 18,
   "id": "367b6ffb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.744928Z",
     "iopub.status.busy": "2026-02-16T18:29:24.744855Z",
     "iopub.status.idle": "2026-02-16T18:29:24.761047Z",
     "shell.execute_reply": "2026-02-16T18:29:24.760692Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0004358917515763\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(flat_vars / effective_flat_var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "283832e7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.761938Z",
     "iopub.status.busy": "2026-02-16T18:29:24.761863Z",
     "iopub.status.idle": "2026-02-16T18:29:24.777955Z",
     "shell.execute_reply": "2026-02-16T18:29:24.777594Z"
    },
    "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": 20,
   "id": "1c92eae4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.779235Z",
     "iopub.status.busy": "2026-02-16T18:29:24.779161Z",
     "iopub.status.idle": "2026-02-16T18:29:24.967732Z",
     "shell.execute_reply": "2026-02-16T18:29:24.967376Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCkAAAQvCAYAAAA6r7R6AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAtINJREFUeJzs3Xt8FNXh//93AuQyQBJhJSEfsxgvhUVRSFAMokKJREQrlWqprKKieCFWpFVBJR/ACz9RkYtU6qdVsIFWrYrUSySCiEJEDaIII95dqiZ0RYiw3DO/P/hmmoUEErLJzu6+no/HPh67Z87unhkgObznXOIsy7IEAAAAAAAQZvHhbgAAAAAAAIBESAEAAAAAAByCkAIAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOQEgBoMXExcVp0qRJ4W4GAADAYV199dVq165duJsBxCRCCgBNMm/ePMXFxQU9OnXqpAEDBui1114Ld/MAAECEO7ivkZSUpMzMTBUUFGjWrFn6+eefw91EACHUOtwNABAdpkyZouzsbFmWpcrKSs2bN08XXnih/vWvf+miiy6SJO3cuVOtW/NjBwAANF5NX2Pv3r2qqKjQ8uXLNXbsWE2fPl2LFy/WaaedFu4mAggB/rcAICQGDx6s3r17269HjRql9PR0/f3vf7dDiqSkpHA1DwAARLiD+xoTJkzQsmXLdNFFF+lXv/qVTNNUcnJyk75jx44datu2bVObCqAJmO4BoFmkpaUpOTk5aOTEwWtSfPvtt7r55pvVtWtXJScnq2PHjrrsssv0zTffBH3W3r17NXnyZJ188slKSkpSx44d1a9fP5WWlrbQ2QAAACf65S9/qYkTJ+rbb79VcXGxJOnjjz/W1VdfrRNOOEFJSUnKyMjQtddeqx9//DHovZMmTVJcXJw2bNigK664Qsccc4z69etX73etXbtWxx57rPr376/t27c363kBsYyRFABCYtu2bfL7/bIsS5s3b9bs2bO1fft2eb3eet/z/vvva9WqVRo+fLiOO+44ffPNN3r88cfVv39/bdiwQYZhSDrQiZg6daquu+46nXnmmaqqqtIHH3ygNWvW6Pzzz2+pUwQAAA505ZVX6q677tKSJUt0/fXXq7S0VF999ZWuueYaZWRkaP369XriiSe0fv16vfvuu4qLiwt6/2WXXaaTTz5ZDzzwgCzLqvM73n//fRUUFKh379566aWXmjxiA0D9CCkAhER+fn7Q68TERD355JOHDRGGDBmi3/zmN0FlF198sfLy8vT888/ryiuvlCS98soruvDCC/XEE0+EvuEAACCiHXfccUpNTdWXX34pSbr55pv1hz/8IajOWWedpd/97nd65513dM455wQdO/3007Vw4cJ6P3/lypW68MILdc455+j5559XYmJi6E8CgI3pHgBCYs6cOSotLVVpaamKi4s1YMAAXXfddXrhhRfqfU/tuxB79+7Vjz/+qJNOOklpaWlas2aNfSwtLU3r16/X559/3qznAAAAIlO7du3sXT5q9y927dolv9+vs846S5KC+hc1brzxxno/980331RBQYEGDhyoF154gYACaAGEFABC4swzz1R+fr7y8/M1YsQIvfLKK+revbsKCwu1Z8+eOt+zc+dOFRUVKSsrS4mJiXK5XDr22GO1detWbdu2za43ZcoUbd26Vb/4xS/Uo0cP3X777fr4449b6tQAAIDDbd++Xe3bt5ckbdmyRbfeeqvS09OVnJysY489VtnZ2ZIU1L+oUXPsYLt27dKQIUPUq1cvPfvss0pISGi+EwBgI6QA0Czi4+M1YMAA/fDDD/WOgLjlllt0//336/LLL9ezzz6rJUuWqLS0VB07dlR1dbVd79xzz9WXX36pJ598Uqeeeqr+8pe/KCcnR3/5y19a6nQAAIBD/fvf/9a2bdt00kknSZIuv/xy/d///Z9uvPFGvfDCC1qyZIlKSkokKah/UaO+9SUSExM1ZMgQrV692n4/gObHmhQAms2+ffskqd4VsP/5z39q5MiReuSRR+yyXbt2aevWrYfU7dChg6655hpdc8012r59u84991xNmjRJ1113XbO0HQAARIa//e1vkqSCggL99NNPWrp0qSZPnqyioiK7ztFMGY2Li9OCBQt0ySWX6LLLLtNrr72m/v37h6rZAOrBSAoAzWLv3r1asmSJEhIS5PF46qzTqlWrQ1bRnj17tvbv3x9UdvCWYe3atdNJJ52k3bt3h7bRAAAgoixbtkz33nuvsrOzNWLECLVq1UqSDulfzJgx46g+PyEhQS+88ILOOOMMXXzxxXrvvfea2mQAR8BICgAh8dprr+nTTz+VJG3evFkLFy7U559/rvHjxyslJaXO91x00UX629/+ptTUVHXv3l1lZWV644031LFjx6B63bt3V//+/ZWbm6sOHTrogw8+0D//+U8VFhY2+3kBAABnqOlr7Nu3T5WVlVq2bJlKS0vVpUsXLV68WElJSUpKStK5556radOmae/evfqf//kfLVmyRF9//fVRf29ycrJefvll/fKXv9TgwYP11ltv6dRTTw3hmQGojZACQEjUHlKZlJSkbt266fHHH9cNN9xQ73tmzpypVq1aacGCBdq1a5fOPvtsvfHGGyooKAiq9/vf/16LFy/WkiVLtHv3bnXp0kX33Xefbr/99mY7HwAA4Cw1fY2EhAR16NBBPXr00IwZM3TNNdfYi2ZK0sKFC3XLLbdozpw5sixLgwYN0muvvabMzMyj/u6UlBS9/vrrOvfcc3X++efr7bffttfAABBacdbBY6EAAAAAAADCgDUpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcITW4W5AOFVXV+v7779X+/btFRcXF+7mAADQ4izL0s8//6zMzEzFx3PvornQ5wAAxLqG9jliOqT4/vvvlZWVFe5mAAAQdps2bdJxxx0X7mZELfocAAAccKQ+R0yHFO3bt5d04CKlpKSEuTWIODt2SJmZB55//73Utm142wMAR6GqqkpZWVn270Q0D/ocaBL6HACiQEP7HDEdUtQMt0xJSaHDgMZr1eq/z1NS6DAAiGhMQWhe9DnQJPQ5AESRI/U5mHwKAAAAAAAcgZACAAAAAAA4AiEFAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBEIKQAAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByhdbgbADiFz+eT3++3X7tcLrnd7jC2CAAARAP6GADQcIQUgA50HjwejwKBgF1mGIZM06QTAQAAjhp9DABoHKZ7AJL8fr8CgYCKi4tVXl6u4uJiBQKBoLseAAAAjVVfH+Ptt9/WmjVrtGbNGvl8vnA3EwAcg5EUQC0ej0c5OTnhbgYAAIgyNX0Ml8slwzDk9XrtY4Zh6IUXXtCxxx4riekgAGIbIQUAAADQQtxut0zTtEdr/uc//9Gll16qCy64wK7DdBAAsYyQAgAAAGhBbrc7KICoHVqYpimv1yu/309IASAmEVIgZtVeads0zTrr1C5n6CUAAGgOB4cWABDLCCkQk+pbadvlcklSvfNFGXoJAAAAAM2HkAIx4+CREzUrbXs8HknBIyUOni/K0EsAAAAAaH6EFIgJ9Y2cOOecc+oNHRh6CQAAAAAti5ACMaH2HuV1jZwAAAAAAIRffGPfsGLFCl188cXKzMxUXFycFi1aFHTcsiwVFRWpc+fOSk5OVn5+vj7//POgOlu2bNGIESOUkpKitLQ0jRo1Stu3bw+q8/HHH+ucc85RUlKSsrKyNG3atEPa8txzz6lbt25KSkpSjx499Oqrrzb2dBBjavYoz8nJOaqAwjRNrVmzRmvWrNGmTZuaoYUAAAAAELsaHVLs2LFDp59+uubMmVPn8WnTpmnWrFmaO3euVq9erbZt26qgoEC7du2y64wYMULr169XaWmpXn75Za1YsUKjR4+2j1dVVWnQoEHq0qWLysvL9dBDD2nSpEl64okn7DqrVq3S7373O40aNUoffvihhg4dqqFDh+qTTz5p7CkBR1R7Ic3c3Fzl5uYqJycn3M0CAABRihsjAGJVo6d7DB48WIMHD67zmGVZmjFjhu655x5dcsklkqSnn35a6enpWrRokYYPHy7TNFVSUqL3339fvXv3liTNnj1bF154oR5++GFlZmZqwYIF2rNnj5588kklJCTolFNO0dq1azV9+nQ7zJg5c6YuuOAC3X777ZKke++9V6WlpXrsscc0d+7co7oYQH3qWkhzdK2dPwAAAEKhrh3GXMnJ+k8Y2wQALSmka1J8/fXXqqioUH5+vl2WmpqqPn36qKysTMOHD1dZWZnS0tLsgEKS8vPzFR8fr9WrV+vXv/61ysrKdO655yohIcGuU1BQoAcffFA//fSTjjnmGJWVlWncuHFB319QUHDI9BPEpto7eUgHQoWmYiFNAADQEAfvKNYY3BgBEOsaPd3jcCoqKiRJ6enpQeXp6en2sYqKCnXq1CnoeOvWrdWhQ4egOnV9Ru3vqK9OzfG67N69W1VVVUEPRJ+anTxqpmXk5ubK6/XKMAy5XK5wNw8A0ECsg4VIdHA/5Gj6IG63215Dq2bBbwCIFSENKZxu6tSpSk1NtR9ZWVnhbhKaQe2dPMrLy+2HaZrNNhJi7dq18vl8zfLZABCrWAcLkaiufkhz9kEAINqEdLpHRkaGJKmyslKdO3e2yysrK9WzZ0+7zubNm4Pet2/fPm3ZssV+f0ZGhiorK4Pq1Lw+Up2a43WZMGFC0BSRqqoqgoooVrOTR3NwuVwykpOlnTslSWf36ycZBp0QAAgh1sFCJGvOfggARLOQjqTIzs5WRkaGli5dapdVVVVp9erVysvLkyTl5eVp69atKi8vt+ssW7ZM1dXV6tOnj11nxYoV2rt3r12ntLRUXbt21THHHGPXqf09NXVqvqcuiYmJSklJCXoAR8PtdmvNmjX267/+5S8KBAJB62AAAJrPkdbBknTEdbBq6tS1DtbGjRv1008/2XVqf09NnZrvqQtTTNFcGL0JINo1OqTYvn271q5dq7Vr10o60Emo+WEZFxensWPH6r777tPixYu1bt06XXXVVcrMzNTQoUMlHUiVL7jgAl1//fV67733tHLlShUWFmr48OHKzMyUJF1xxRVKSEjQqFGjtH79ej3zzDOaOXNm0CiIW2+9VSUlJXrkkUf06aefatKkSfrggw9UWFjY9KsCNEDtUTjdunULY0sAIPY4fR0sppgiVOzRm//P2f36yePxEFQAiFqNDik++OAD9erVS7169ZIkjRs3Tr169VJRUZEk6Y477tAtt9yi0aNH64wzztD27dtVUlKipKQk+zMWLFigbt26aeDAgbrwwgvVr1+/oLmfqampWrJkib7++mvl5ubqD3/4g4qKioLmkPbt21cLFy7UE088odNPP13//Oc/tWjRIp166qlHfTEAAABCYcKECdq2bZv92LRpU7ibhAjF6E0AsabRa1L0799flmXVezwuLk5TpkzRlClT6q3ToUMHLVy48LDfc9ppp+ntt98+bJ3LLrtMl1122eEbDAAAoo7T18FKTExUYmLiUZwZcChGbwKIJTG1uwcAAIgOTl8HCwAAHB1CCgAA4EisgwUAQOwJ6RakAAAAofLBBx9owIAB9uua4GDkyJGaN2+e7rjjDu3YsUOjR4/W1q1b1a9fvzrXwSosLNTAgQMVHx+vYcOGadasWfbxmnWwxowZo9zcXLlcrnrXwbrnnnt011136eSTT2YdLAAAmgkhBQAAcCTWwQIAIPYw3QMAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARWDgTAAAAaAKfzye/3y9JMk0zzK0BgMhGSAEAAAAcJZ/PJ4/Ho0AgYJcZhiGXyxXGVgFA5CKkAAAAAI6S3+9XIBBQcXGxPB6PJMnlcsntdjfr99YesdES3wcALYWQAgAAAGgij8ejnJycZv+ejh07yjAMeb1eu8wwDJmmSVABICoQUgAAAAARIisrS6ZpBq2B4fV65ff7CSkARAVCCgAAACCCuN1uAgkAUYuQAggh5ocCABDdau/kIbGbBwCEGiEFokK4t/5ifigAANGvrp08JHbzAIBQIqRAxHPC1l/MDwUAIPrVtZOHxOhJAAglQgpEvHBt/XUw5ocCABAbWmonDwCIRYQUiBpO7DCwRgUAAAAANBwhBdAMXC4Xa1QAAAAAQCMRUgDNwO12s0YFAAAAADQSIQXQTFijAgAAAAAaJz7cDQAAAAAAAJAIKQAAAAAAgEMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByBkAIAAAAAADgCW5ACAAAAEc40Tfu5y+ViG3QAEYuQAgAAAIhQLpdLhmHI6/XaZYZhyDRNggoAEYmQAgAAAIhQbrdbpmnK7/dLOjCiwuv1yu/3E1IAiEiEFAAAAEAEc7vdBBIAogYLZwIAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByBNSkQkXw+X9ACUQAAAACAyEdIgYjj8/nk8XgUCATsMsMw5HK5wtgqAAAQjbgxAgAti5ACEcfv9ysQCKi4uFgej0fSgT3CWdUaAACEEjdGAKDlEVIgYnk8HuXk5IS7GQAAIEpxYwQAWh4hBQAAAHAY3BgBgJbD7h4AAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOwMKZQAuqvb86q4MDAAAAQDBCCqAFuFwuGYYhr9drlxmGIdM0CSoAAAAA4P8hpABagNvtlmma8vv9kg6MqPB6vfL7/YQUAAAg5GqP3pQYwQkgchBSAC3E7XbTOQAAAM2qrtGbEiM4AUQOQgoAAAAgShw8elNiBCeAyEJIAQAAAEQRRm8CiGRsQQoAAAAAAByBkAIAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjsAUpIoLP57P3+zZNM8ytAQAAAAA0B0IKOJ7P55PH41EgELDLDMOQy+UKY6sAAAAAAKFGSAHH8/v9CgQCKi4ulsfjkSS5XC653e4wtwwAAAAAEEqEFIgYHo9HOTk54W4GAAAAAKCZEFIAYXTw+hqMEAEAAAAQywgpgDBwuVwyDENerzeo3DAMmaZJUAEAQJiwWDcAhBchBRAGbrdbpmnanSDpQEfI6/XK7/cTUgAAEAYs1g0A4UdIAYSJ2+0mjAAAwEFYrBsAwo+QAgAAAKiFxboBIHziw90AAAAAAAAAiZACAAAAAAA4BCEFAAAAAABwBEIKAAAAAADgCCycCQAAAMQA0zTt5+xaAsCpCCkAAACAKOZyuWQYhrxer11mGIZM0ySoAOA4hBQAAABAFHO73TJNU36/X9KBERVer1d+v5+QAoDjhHxNiv3792vixInKzs5WcnKyTjzxRN17772yLMuuY1mWioqK1LlzZyUnJys/P1+ff/550Ods2bJFI0aMUEpKitLS0jRq1Cht3749qM7HH3+sc845R0lJScrKytK0adNCfToAAABAxHO73crJyVFOTo48Hk+4mwMA9Qp5SPHggw/q8ccf12OPPSbTNPXggw9q2rRpmj17tl1n2rRpmjVrlubOnavVq1erbdu2Kigo0K5du+w6I0aM0Pr161VaWqqXX35ZK1as0OjRo+3jVVVVGjRokLp06aLy8nI99NBDmjRpkp544olQnxIAAAAAAGgBIZ/usWrVKl1yySUaMmSIJOn444/X3//+d7333nuSDoyimDFjhu655x5dcsklkqSnn35a6enpWrRokYYPHy7TNFVSUqL3339fvXv3liTNnj1bF154oR5++GFlZmZqwYIF2rNnj5588kklJCTolFNO0dq1azV9+vSgMAOINCxqBQAAACBWhXwkRd++fbV06VJ99tlnkqSPPvpI77zzjgYPHixJ+vrrr1VRUaH8/Hz7PampqerTp4/KysokSWVlZUpLS7MDCknKz89XfHy8Vq9ebdc599xzlZCQYNcpKCjQxo0b9dNPP9XZtt27d6uqqiroAThF7UWtcnNzlZubK4/HI5/PF+6mAYAjMcUUAIDoE/KQYvz48Ro+fLi6deumNm3aqFevXho7dqxGjBghSaqoqJAkpaenB70vPT3dPlZRUaFOnToFHW/durU6dOgQVKeuz6j9HQebOnWqUlNT7UdWVlYTzxYInZpFrcrLy1VeXq7i4mIFAgF7kSsAQDCmmAIAEH1CPt3j2Wef1YIFC7Rw4UJ7CsbYsWOVmZmpkSNHhvrrGmXChAkaN26c/bqqqoqgAo7idruZ3gEADcQUUwAAok/IR1Lcfvvt9miKHj166Morr9Rtt92mqVOnSpIyMjIkSZWVlUHvq6ystI9lZGRo8+bNQcf37dunLVu2BNWp6zNqf8fBEhMTlZKSEvQAAACRiSmmAABEn5CHFIFAQPHxwR/bqlUrVVdXS5Kys7OVkZGhpUuX2serqqq0evVq5eXlSZLy8vK0detWlZeX23WWLVum6upq9enTx66zYsUK7d27165TWlqqrl276phjjgn1aQEAAIdhiikAANEn5CHFxRdfrPvvv1+vvPKKvvnmG7344ouaPn26fv3rX0uS4uLiNHbsWN13331avHix1q1bp6uuukqZmZkaOnSoJMnj8eiCCy7Q9ddfr/fee08rV65UYWGhhg8frszMTEnSFVdcoYSEBI0aNUrr16/XM888o5kzZwZN5wAAANGr9hTTNWvWaP78+Xr44Yc1f/78cDdNEyZM0LZt2+zHpk2bwt0kAAAiQsjXpJg9e7YmTpyom2++WZs3b1ZmZqZuuOEGFRUV2XXuuOMO7dixQ6NHj9bWrVvVr18/lZSUKCkpya6zYMECFRYWauDAgYqPj9ewYcM0a9Ys+3hqaqqWLFmiMWPGKDc3Vy6XS0VFRcwNBQAgRtSeYipJPXr00LfffqupU6dq5MiRQVNMO3fubL+vsrJSPXv2lNS8U0wTExObfpIAAMSYkIcU7du314wZMzRjxox668TFxWnKlCmaMmVKvXU6dOighQsXHva7TjvtNL399ttH21QAABDBGjPFtCaUqJlietNNN0kKnmKam5srqe4ppnfffbf27t2rNm3aSGKKKQAAzSXkIQUAAEBLqJli6na7dcopp+jDDz/U9OnTde2110oKnmJ68sknKzs7WxMnTqx3iuncuXO1d+/eOqeYTp48WaNGjdKdd96pTz75RDNnztSjjz4arlMHQsI0Tfu5y+VihzEAjkBIAQAAIhJTTIGj43K5ZBiGvF6vXWYYhkzTJKgAEHaEFAAAICIxxRQ4Om63W6Zpyu/3SzowosLr9crv9xNSAAg7QgoAAAAgxrjdbgIJAI4U8i1IAQAAAAAAjgYhBQAAAAAAcASmewAAACBm+Xy+oLUZAADhRUgBR6LDAAAAmpvP55PH41EgELDLDMOQy+UKY6sAILYRUsBx6DAAAICW4Pf7FQgEVFxcLI/HI+nA9pwsKAkA4UNIAcehwwAAAFqSx+NRTk5OuJsBABAhBRyMDgMAAAAAxBZ29wAAAAAAAI5ASAEAAAAAAByBkAIAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMAR2IIUcDjTNO3nLpdLbrc7jK0BAAAAgOZDSAE4lMvlkmEY8nq9dplhGDJNk6ACAAAAQFQipAAcyu12yzRN+f1+SQdGVHi9Xvn9fkIKAAAAAFGJkAJwMLfbTSABAAAAIGawcCYAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjsHAmAAAAAJmmaT93uVws3g0gLAgpAAAAgBjmcrlkGIa8Xq9dZhiGTNMkqADQ4ggpAAAAgBjmdrtlmqb8fr+kAyMqvF6v/H4/IQWAFkdIAQAAgJjh8/mC/jOOA9xuN4EEAEcgpAAAAEBM8Pl88ng8CgQCdplhGHK5XGFsFQCgNkIKAAAAxAS/369AIKDi4mJ5PB5JLBAJAE5DSAEAAICY4vF4lJOTE+5mAADqEB/uBgAAAAAAAEiEFAAAAAAAwCGY7gFEmNorkTOPFgAAAEA0IaQAIoTL5ZJhGPJ6vXaZYRgyTZOgAgAAAEBUIKSAI7Bn+ZG53W6Zphl0nbxer/x+PyEFAAAAgKhASIGwY8/yhnO73QQSAAAAAKIWIQXCjj3LAQAAAAASIQUchD3LAQAAACC2sQUpAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOELrcDcAAAAAgPOYpmk/d7lccrvdYWwNgFhBSAEAAADA5nK5ZBiGvF6vXWYYhkzTJKgA0OwIKQAAAADY3G63TNOU3++XdGBEhdfrld/vJ6QA0OwIKQAAAAAEcbvdBBIAwoKFMwEAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByBkAIAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIrcPdAABNY5qm/dzlcrGnOQAAAICIRUgBRCiXyyXDMOT1eu0ywzBkmiZBBQAAAICIREgBRCi32y3TNOX3+yUdGFHh9Xrl9/sJKQAAAABEJEIKIIK53W4CCQAADsPn8wUF+gAAZyOkAAAAQFTy+XzyeDwKBAJ2mWEYcrlcYWwVAOBwCCkAAAAQlfx+vwKBgIqLi+XxeCSxyDQAOB0hBQAAAKKax+NRTk5OuJsBAGiA+HA3AAAAAAAAQCKkAAAAAAAADsF0D4QFK20DAAAAAA5GSIEWx0rbAAAAkaf2jSUWIAXQXAgp0OJYaRsAACByuFwuGYYhr9drlxmGIdM06b8BCLlmWZPiu+++k9frVceOHZWcnKwePXrogw8+sI9blqWioiJ17txZycnJys/P1+effx70GVu2bNGIESOUkpKitLQ0jRo1Stu3bw+q8/HHH+ucc85RUlKSsrKyNG3atOY4HTSTmpW2c3Jy+AUHAADgUG63W6Zpqry8XOXl5SouLlYgELCn7gJAKIU8pPjpp5909tlnq02bNnrttde0YcMGPfLIIzrmmGPsOtOmTdOsWbM0d+5crV69Wm3btlVBQYF27dpl1xkxYoTWr1+v0tJSvfzyy1qxYoVGjx5tH6+qqtKgQYPUpUsXlZeX66GHHtKkSZP0xBNPhPqUAAAAgJjmdrvtm0s1I2EBoDmEfLrHgw8+qKysLD311FN2WXZ2tv3csizNmDFD99xzjy655BJJ0tNPP6309HQtWrRIw4cPl2maKikp0fvvv6/evXtLkmbPnq0LL7xQDz/8sDIzM7VgwQLt2bNHTz75pBISEnTKKado7dq1mj59elCYAQAAAAAAIkPIR1IsXrxYvXv31mWXXaZOnTqpV69e+r//+z/7+Ndff62Kigrl5+fbZampqerTp4/KysokSWVlZUpLS7MDCknKz89XfHy8Vq9ebdc599xzlZCQYNcpKCjQxo0b9dNPP4X6tAAAgAMxxRQAgOgS8pDiq6++0uOPP66TTz5Zr7/+um666Sb9/ve/1/z58yVJFRUVkqT09PSg96Wnp9vHKioq1KlTp6DjrVu3VocOHYLq1PUZtb/jYLt371ZVVVXQAwAARCammAIAEH1CPt2jurpavXv31gMPPCBJ6tWrlz755BPNnTtXI0eODPXXNcrUqVM1efLksLYBAACEBlNMAQCIPiEfSdG5c2d17949qMzj8cjn80mSMjIyJEmVlZVBdSorK+1jGRkZ2rx5c9Dxffv2acuWLUF16vqM2t9xsAkTJmjbtm32Y9OmTUdzigAAwAGcPMWU0ZsAABydkIcUZ599tjZu3BhU9tlnn6lLly6SDtzhyMjI0NKlS+3jVVVVWr16tfLy8iRJeXl52rp1q8rLy+06y5YtU3V1tfr06WPXWbFihfbu3WvXKS0tVdeuXYOGedaWmJiolJSUoAcAAIhMTp5iOnXqVKWmptqPrKysJp4tAACxIeQhxW233aZ3331XDzzwgL744gstXLhQTzzxhMaMGSNJiouL09ixY3Xfffdp8eLFWrduna666iplZmZq6NChkg6MvLjgggt0/fXX67333tPKlStVWFio4cOHKzMzU5J0xRVXKCEhQaNGjdL69ev1zDPPaObMmRo3blyoTwkAADhQdXW1cnJy9MADD6hXr14aPXq0rr/+es2dOzfcTWP0JgAARynka1KcccYZevHFFzVhwgRNmTJF2dnZmjFjhkaMGGHXueOOO7Rjxw6NHj1aW7duVb9+/VRSUqKkpCS7zoIFC1RYWKiBAwcqPj5ew4YN06xZs+zjqampWrJkicaMGaPc3Fy5XC4VFRUxNxQxzzRN+7nL5ZLb7Q5jawCg+dQ3xfT555+XFDzFtHPnznadyspK9ezZ067THFNMExMTlZiYeJRnBgBA7Ap5SCFJF110kS666KJ6j8fFxWnKlCmaMmVKvXU6dOighQsXHvZ7TjvtNL399ttH3U4gmrhcLhmGIa/Xa5cZhiHTNAkqAESlxkwxrQklaqaY3nTTTZKCp5jm5uZKqnuK6d133629e/eqTZs2ko48xRQAABydkE/3ABAebrdbpmmqvLxc5eXlKi4uViAQkN/vD3fTAKBZMMUUAIDo0ywjKQCEh9vtZtQEgJjBFFMAAKIPIQUAAIhYTDEFACC6MN0DAAAAAAA4AiEFAAAAAABwBKZ7AAAAIGr4fD570eja23IDACIDIQUAAACigs/nk8fjUSAQsMsMw5DL5QpjqwAAjUFIAQAAgKjg9/sVCARUXFwsj8cjSXK5XOx8BQARhJACAAAAUcXj8SgnJyfczQAAHAUWzgQAAAAAAI7ASAoAAAAAjVZ7YVKm1QAIFUIKAAAAAA3mcrlkGIa8Xq9dZhiGTNMkqADQZIQUAAAAABrM7XbLNM2grV69Xq/8fj8hBYAmI6QAAAAA0Chut5tAAkCzYOFMAAAAAADgCIQUAAAAAADAEZjugRbh8/mC5i0CAAAAAHAwQgo0O5/PJ4/Ho0AgYJcZhiGXyxXGVgEAAAAAnIaQAs3O7/crEAiouLhYHo9HEntpAwAAAAAORUiBFuPxeJSTkxPuZgAAAAAAHIqFMwEAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI7AmhRAlKu95SsLlgIAAABwMkIKIEq5XC4ZhiGv12uXGYYh0zQJKgAAAAA4EiEFEKXcbrdM05Tf75d0YESF1+uV3+8npAAAACHH6E0AoUBIAUQxt9tNBwEAADQrRm8CCCVCCgAAAABHjdGbAEKJkAIAAABAkzB6E0CosAUpAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBEIKQAAAAAAgCO0DncDALQs0zTt5y6XS263O4ytAQAAAID/IqQAYoTL5ZJhGPJ6vXaZYRgyTZOgAgAAAIAjEFIAMcLtdss0Tfn9fkkHRlR4vV75/X5CCgAAAACOQEgBxBC3200gAQAAAMCxWDgTAAAAAAA4AiMpAAAAAIQci3UDOBqEFGgWPp8vaO0DAAAAxAYW6wbQFIQUCDmfzyePx6NAIGCXGYYhl8sVxlYBAIBoxI0R52GxbgBNQUiBkPP7/QoEAiouLpbH45HEED8AABB63BhxLhbrBnC0CCnQbDwej3JycsLdDAAAEKW4MQIA0YeQAgAAABGNGyMAED3YghQAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARWDgTiHG195RnRXQAAAAA4URIAcQol8slwzDk9XrtMsMwZJomQQUAAACAsCCkAGKU2+2WaZry+/2SDoyo8Hq98vv9hBQAAAAAwoKQAohhbrebQAIAAACAY7BwJgAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBFYkwIAAABAs2PbcwANQUgBAAAAoNmw7TmAxiCkAAAAANBs2PYcQGMQUgAAAABoVmx7DqChWDgTAAAAAAA4AiEFAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEVg4EyHh8/mCVmwGAAAAAKCxmn0kxf/3//1/iouL09ixY+2yXbt2acyYMerYsaPatWunYcOGqbKyMuh9Pp9PQ4YMkWEY6tSpk26//Xbt27cvqM7y5cuVk5OjxMREnXTSSZo3b15znw7q4PP55PF4lJubq9zcXHm9XhmGIZfLFe6mAQAAAAAiSLOGFO+//77+/Oc/67TTTgsqv+222/Svf/1Lzz33nN566y19//33uvTSS+3j+/fv15AhQ7Rnzx6tWrVK8+fP17x581RUVGTX+frrrzVkyBANGDBAa9eu1dixY3Xdddfp9ddfb85TQh38fr8CgYCKi4tVXl6u8vJymabJNlMAgBbFjREAACJfs4UU27dv14gRI/R///d/OuaYY+zybdu26a9//aumT5+uX/7yl8rNzdVTTz2lVatW6d1335UkLVmyRBs2bFBxcbF69uypwYMH695779WcOXO0Z88eSdLcuXOVnZ2tRx55RB6PR4WFhfrNb36jRx99tLlOCUfg8XiUk5OjnJwcAgoAQIvixggAANGh2UKKMWPGaMiQIcrPzw8qLy8v1969e4PKu3XrJrfbrbKyMklSWVmZevToofT0dLtOQUGBqqqqtH79ervOwZ9dUFBgf0Zddu/eraqqqqAHAACIbNwYAQAgejRLSPGPf/xDa9as0dSpUw85VlFRoYSEBKWlpQWVp6enq6Kiwq5TO6CoOV5z7HB1qqqqtHPnzjrbNXXqVKWmptqPrKysozo/AADgHE68MQIAAI5OyHf32LRpk2699VaVlpYqKSkp1B/fJBMmTNC4cePs11VVVQQVAABEsJobI++///4hx1rqxkhycvIh3717927t3r3bfs3oTQAAGibkIynKy8u1efNm5eTkqHXr1mrdurXeeustzZo1S61bt1Z6err27NmjrVu3Br2vsrJSGRkZkqSMjIxDFrWqeX2kOikpKXV2FiQpMTFRKSkpQQ8AABCZam6MLFiwwHE3Rhi9CQDA0Ql5SDFw4ECtW7dOa9eutR+9e/fWiBEj7Odt2rTR0qVL7fds3LhRPp9PeXl5kqS8vDytW7dOmzdvtuuUlpYqJSVF3bt3t+vU/oyaOjWfAQAAopuTb4xMmDBB27Ztsx+bNm0KxSlDB3ZjWbNmjdasWSPTNMPdHABAiIV8ukf79u116qmnBpW1bdtWHTt2tMtHjRqlcePGqUOHDkpJSdEtt9yivLw8nXXWWZKkQYMGqXv37rryyis1bdo0VVRU6J577tGYMWOUmJgoSbrxxhv12GOP6Y477tC1116rZcuW6dlnn9Urr7wS6lMCAAAOVHNjpLZrrrlG3bp105133qmsrCz7xsiwYcMk1X1j5P7779fmzZvVqVMnSXXfGHn11VeDvudIN0YSExPtPgtCx+fzyePxKBAI2GWGYcjlcoWxVThaB4dMLpeLHeIAhD6kaIhHH31U8fHxGjZsmHbv3q2CggL96U9/so+3atVKL7/8sm666Sbl5eWpbdu2GjlypKZMmWLXyc7O1iuvvKLbbrtNM2fO1HHHHae//OUvKigoCMcpAQCAFsaNkdjj9/sVCARUXFwsj8cjif/YRiKXyyXDMOT1eoPKDcOQaZr8eQIxrkVCiuXLlwe9TkpK0pw5czRnzpx639OlS5dD7locrH///vrwww9D0UQAABCFuDESnTwej3JycsLdDBwlt9st0zTl9/vtMtM05fV65ff7CSmAGBeWkRQAnKv20EvuTgGINNwYASKD2+2mjwGgToQUACTVPfSSYZcAAAAAWhIhBQBJhw69ZNglAAAAgJZGSAHAxtBLAAAAAOEUH+4GAAAAAAAASIQUAAAAAADAIQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4Aitw90AAM5mmmbQa5fLJbfbHabWAAAAAIhmhBQA6uRyuWQYhrxeb1C5YRgyTZOgAgAAAEDIEVIAqJPb7ZZpmvL7/XaZaZryer3y+/2EFAAAIORqj+Bk9CYQmwgp0Gg+n++Q/7giOrndbjoHAACg2dU1gpPRm0BsIqRAo/h8Pnk8HgUCgaBywzDkcrnC1CoAAABEsoNHcDJ6E4hdhBRoFL/fr0AgoOLiYnk8Hruc4XgAAABoCkZwApAIKXCUPB6PcnJywt0MAAAAAEAUiQ93AwAAAAAAACRCCgAAAAAA4BCEFAAAAAAAwBEIKQAAAAAAgCMQUgAAAAAAAEdgdw8AAAA4ls/nk9/vlySZphnm1gAAmhshBQAAABzJ5/PJ4/EoEAjYZYZhyOVyhbFVAIDmREgBAAAAR/L7/QoEAiouLpbH45EkuVwuud3uMLcMANBcCCkAAADgaB6PRzk5OeFuBgCgBbBwJgAAAAAAcARGUgBotNoLlzHsFgAAAECoEFIAaDCXyyXDMOT1eu0ywzBkmiZBBQAAAIAmI6QA0GBut1umaQZtBef1euX3+wkpAABAyDF6E4g9hBQAGsXtdtNBAAAAzYrRm0DsIqQAAAAA4CiM3gRiFyEFAAAAAMdh9CYQm9iCFAAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBEIKQAAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI7QOtwNABD5TNO0n7tcLvY0BwAAAHBUCCkAHDWXyyXDMOT1eu0ywzBkmiZBBQAAAIBGI6QAcNTcbrdM05Tf75d0YESF1+uV3+8npAAAACHH6E0g+hFSAGgSt9tNBwEAADQrRm8CsYOQAgAAAICjMXoTiB2EFAAAAAAcj9GbQGxgC1IAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjsCYFjsjn8wUtUgQAANAcavc5JPodABCLCClwWD6fTx6PR4FAwC4zDEMulyuMrQIAANGmrj6HRL8DAGINIQUOy+/3KxAIqLi4WB6PR9KBfapZWRkAAIRSXX0OiX4HAMQaQgo0iMfjUU5OTribAQAAohx9DgCIbYQUAEKu9hxi7oABAAAAaChCCgAh43K5ZBiGvF6vXWYYhkzTJKgAAAAhx40RIPoQUgAIGbfbLdM0g3aD8Xq98vv9dBoAAEDIcGMEiF6EFABCyu120zkAAADNihsjQPQipAAAAAAQcbgxAkSn+HA3AAAAAAAAQCKkAAAAAAAADkFIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AhsQQoAAAAgKpimaT93uVxsUQpEIEIKAAAAABHN5XLJMAx5vV67zDAMmaZJUAFEGEIKAAAAABHN7XbLNE35/X5JB0ZUeL1e+f1+QgogwhBSAGh2DL0EAADNze1208cAokDIF86cOnWqzjjjDLVv316dOnXS0KFDtXHjxqA6u3bt0pgxY9SxY0e1a9dOw4YNU2VlZVAdn8+nIUOGyDAMderUSbfffrv27dsXVGf58uXKyclRYmKiTjrpJM2bNy/UpxOTfD6f1qxZozVr1gT95xJorNpDL3Nzc5WbmyuPxyOfzxfupgEAAABwoJCHFG+99ZbGjBmjd999V6Wlpdq7d68GDRqkHTt22HVuu+02/etf/9Jzzz2nt956S99//70uvfRS+/j+/fs1ZMgQ7dmzR6tWrdL8+fM1b948FRUV2XW+/vprDRkyRAMGDNDatWs1duxYXXfddXr99ddDfUoxxefzyePx2P+h9Hq9MgxDLpcr3E1DBKoZelleXq7y8nIVFxcrEAjYQzEBoCm4MRL5uDECADhYyKd7lJSUBL2eN2+eOnXqpPLycp177rnatm2b/vrXv2rhwoX65S9/KUl66qmn5PF49O677+qss87SkiVLtGHDBr3xxhtKT09Xz549de+99+rOO+/UpEmTlJCQoLlz5yo7O1uPPPKIJMnj8eidd97Ro48+qoKCglCfVszw+/0KBAIqLi6Wx+ORxPB8NA1DLwE0l5obI2eccYb27dunu+66S4MGDdKGDRvUtm1bSQdujLzyyit67rnnlJqaqsLCQl166aVauXKlpP/eGMnIyNCqVav0ww8/6KqrrlKbNm30wAMPSPrvjZEbb7xRCxYs0NKlS3Xdddepc+fO9DmaoObGSCAQsMu4MQIAaPY1KbZt2yZJ6tChgySpvLxce/fuVX5+vl2nW7ducrvdKisr01lnnaWysjL16NFD6enpdp2CggLddNNNWr9+vXr16qWysrKgz6ipM3bs2OY+pZjg8XiUk5MT7mYAAFAvboxENm6MAADqEvLpHrVVV1dr7NixOvvss3XqqadKkioqKpSQkKC0tLSguunp6aqoqLDr1A4oao7XHDtcnaqqKu3cubPO9uzevVtVVVVBDwAAEB0ae2NEUr03RqqqqrR+/Xq7Tl03Rmo+oy70ORqu5sZITk4OAQUAoHlDijFjxuiTTz7RP/7xj+b8mgabOnWqUlNT7UdWVla4mwQAAELAaTdG6HMAAHB0mi2kKCws1Msvv6w333xTxx13nF2ekZGhPXv2aOvWrUH1KysrlZGRYdc5eFGrmtdHqpOSkqLk5OQ62zRhwgRt27bNfmzatKlJ5wgAAJzBaTdG6HMAAHB0Qh5SWJalwsJCvfjii1q2bJmys7ODjufm5qpNmzZaunSpXbZx40b5fD7l5eVJkvLy8rRu3Tpt3rzZrlNaWqqUlBR1797drlP7M2rq1HxGXRITE5WSkhL0AAAAkc2JN0bocwAAcHRCHlKMGTNGxcXFWrhwodq3b6+KigpVVFTYwyFTU1M1atQojRs3Tm+++abKy8t1zTXXKC8vT2eddZYkadCgQerevbuuvPJKffTRR3r99dd1zz33aMyYMUpMTJQk3Xjjjfrqq690xx136NNPP9Wf/vQnPfvss7rttttCfUoAAMCBnHxjBAAAHJ2Q7+7x+OOPS5L69+8fVP7UU0/p6quvliQ9+uijio+P17Bhw7R7924VFBToT3/6k123VatWevnll3XTTTcpLy9Pbdu21ciRIzVlyhS7TnZ2tl555RXddtttmjlzpo477jj95S9/YZVtAABixJgxY7Rw4UK99NJL9o0R6cANkeTk5KAbIx06dFBKSopuueWWem+MTJs2TRUVFXXeGHnsscd0xx136Nprr9WyZcv07LPP6pVXXgnbuQMAEK1CHlJYlnXEOklJSZozZ47mzJlTb50uXbro1VdfPezn9O/fXx9++GGj2wggvEzTtJ+z3RyAo8WNEQBHQp8DiDwhDykAoD4ul0uGYcjr9dplhmHINE06DQAajRsjAOpDnwOIXIQUAFqM2+2WaZry+/2SDtzd8Hq98vv9dBgAAEDI0OcAIhchBYAW5Xa76RwAAIBmR58DiEwh390DAAAAAADgaBBSAAAAAAAARyCkAAAAAAAAjsCaFJDP5wtaVAgAAAAAgHAgpIhxPp9PHo9HgUDALjMMQy6XK4ytQqxhD3MAANAS6HMAzkdIEeP8fr8CgYCKi4vl8Xgk8QMbLYc9zAEAQEugzwFEDkIKSJI8Ho9ycnLC3QzEGPYwBwAALYE+BxA5CCkAhBV7mAMAgJZAnwOIDIQUAAAAaBEs1g0AOBJCCgAAADQ7FusGADQEIQUAAACaHYt1AwAagpACAAAALYbFugEAhxMf7gYAAAAAAABIhBQAAAAAAMAhmO4BwHFqr/jOfGUAANBc6HMAzkNIAcAxXC6XDMOQ1+u1ywzDkGmadBoAAEDI0OcAnIuQAoBjuN1umaYpv98v6cDdDa/XK7/fT4cBAACEDH0OwLkIKQA4itvtpnMAAACaHX0OwJlYOBMAAAAAADgCIQUAAAAAAHAEpnsAAAAAgNjtA3ACQooY5PP5ghYJAgAAAGIZu30AzkFIEWN8Pp88Ho8CgYBdZhiGXC5XGFsFAAAAhA+7fQDOQUgRY/x+vwKBgIqLi+XxeCQxlA3Ox9BLAADQ3NjtA3AGQooY5fF4lJOTE+5mAIfF0EsAAAAgthBSAHAshl4CAAAAsYWQAoCjMfQSAACEy8GLzDPtFGh+hBQAAABoFuwohkhV15RTiWmnQEsgpAAAAEDIsaMYItnBU04lpp0CLYWQAgAAACHHjmKIdEw5BcKDkAJAxGFLUgCIHOwoBgBoDEIKABGDLUkBAACA6EZIASBisCUpAAAAEN0IKQBEFOaHAgAAANGLkAJAxGMPcwAAACA6EFIAiFjsYQ4AAABEF0IKABGLPcwBAEBLY5cxoHkRUsQAn88XtNAgEE1YowIAALSE+nYZe+GFF3TsscfadeiXAE1DSBHlfD6fPB6PAoGAXWYYhlwuVxhbBQAAAESWg0dw/uc//9Gll16qCy64wK7DlFOg6Qgpopzf71cgEFBxcbE8Ho8kEl4AAADgaBw8gpOt0YHQI6SIER6PRzk5OeFuBgAAiGJMMUWsYdopEHqEFAAAAGgyppgCAEKBkAJAVGLlbQBoWUwxBQCEAiEFgKjCytsAEF5MMUWs40YJ0DSEFACiCitvAwCAcKjvRgl9DqBxCCkARB1W3gYAAC3t4Bsl9DmAo0NIASDqsfI2AABoCfQ5gKaLD3cDAAAAAAAAJEZSRCX2KAcAAACcgYU0gcYhpIgy7FEONAwdBgAA0JxYSBM4OoQUUYY9yoHDo8MAAABaAgtpAkeHkCJKsUc5UDc6DAAAoKWwkCbQeIQUAGIOHQYAABAuTDkFDo+QAgBEhwEAjgaLdQMNx5RToGEIKQDEtPo6DC+88IKOPfZYuw6dBwAIxmLdQOMw5RRoGEIKADHt4A7Df/7zH1166aW64IIL7Drc5QCAQ7FYN9B4dU05ZTQnEIyQAkDMO7jDwF0OAGg4FusGjg7TP4C6EVJEAeaDAqHFwpoAAKC51Tf94+2332Z0EmIaIUWEYz4o0DIYigkAAEKt9o0R1skCDiCkiHDMBwWaF0MxAQBAS2CdLOAAQooowXxQoHkwFBMA/osppkDzYp0sgJACAI6oIUMxuasBINoxxRRoeewGglhESBGBuIsBhA8jKwDEKqaYAuHFjRLECkKKCMNdDCD8GFkBIJYxxRQIj4bcKJEIDxH5CCkiDHcxAGdhZAWAaMboTcBZjnSjRGJHEEQ+QooIxV0MwDkYWQEgGjF6E3C2g2+USPXvCEJogUhCSBEBuIsBRI6GDsU8GB0GAE7D6E3A+epbWJNtTBHJCCkcjrsYQORpyFDMg3GXA4AT1HVjhNGbQGRpyDamTEuFkxFSOEztzoF04AcJdzGAyFXXUMyDNWRo5sH4OQAg1LgxAkSnhkxLPVyf42D0QdDcIj6kmDNnjh566CFVVFTo9NNP1+zZs3XmmWeGu1lHpa7OgXTgB8c555zDDwMgQtU1FPNgRxqaeTCGagItL5r6HDUOHjnBjREguh1886QhfY6DcSMFzS2iQ4pnnnlG48aN09y5c9WnTx/NmDFDBQUF2rhxozp16hTu5jXIkToHEv/QgVhwuKGZB2vIOhf83ABCKxr6HAerb+QEN0aA6NaYPsfBGnojpTEjMyT6LQgWZ1mWFe5GHK0+ffrojDPO0GOPPSZJqq6uVlZWlm655RaNHz/+iO+vqqpSamqqtm3bppSUlOZurqTgUKLmH/nBnQPujkaIHTukdu0OPN++XWrbNrztQcyob9RVbUfTQWgsOhTRIRy/CyNRJPY5DlbXlFKv18vIiUhAnwMOcvDPktrq+v9NQzDlJDY09HdhxI6k2LNnj8rLyzVhwgS7LD4+Xvn5+SorK6vzPbt379bu3bvt19u2bZN04GKFSkVFhSoqKuo85vf75fV6tXPnTrssOTlZzz//vD3fs2PHjkpLSwtpm9BMduz47/OqKmn//vC1BTElLS1N7733nn788cc6j9f8rGnM0M2jkZycrOLiYuart7CMjAxlZGSE7PNqft9E8D2LZheJfY6D1dUHkQ78O+7Zs6eysrLsMvogDkSfAw6SlpamtLS0Oo+ddNJJh+2j1OVo+i30QVpGuPocERtS+P1+7d+/X+np6UHl6enp+vTTT+t8z9SpUzV58uRDymv/Ym5pO3fu1LBhw8L2/QiRzMxwtwBocfz8ii4///yzUlNTw90MR4qWPkdddu7cqVNPPTXczUBj0OcA6INEuCP1OSI2pDgaEyZM0Lhx4+zX1dXV2rJlizp27Ki4uLgwtqzlVVVVKSsrS5s2bWJ4bxNwHZuOaxgaXMemi9VraFmWfv75Z2XyH5+Qos9xqFj9NxYqXL+m4fo1Ddevabh+BzS0zxGxIYXL5VKrVq1UWVkZVF5ZWVnvkJTExEQlJiYGldU3VClWpKSkxPQ/lFDhOjYd1zA0uI5NF4vXkBEUh0efI7Ri8d9YKHH9mobr1zRcv6bh+jWszxHfAu1oFgkJCcrNzdXSpUvtsurqai1dulR5eXlhbBkAAIgm9DkAAGg5ETuSQpLGjRunkSNHqnfv3jrzzDM1Y8YM7dixQ9dcc024mwYAAKIIfQ4AAFpGRIcUv/3tb/Wf//xHRUVFqqioUM+ePVVSUnLIwlY4VGJiov73f//3kKGoaByuY9NxDUOD69h0XEMcDn2OpuPfWNNw/ZqG69c0XL+m4fo1TpzFnmMAAAAAAMABInZNCgAAAAAAEF0IKQAAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI5ASBHF5syZo+OPP15JSUnq06eP3nvvvXrrzps3T3FxcUGPpKSkFmyt86xYsUIXX3yxMjMzFRcXp0WLFh3xPcuXL1dOTo4SExN10kknad68ec3eTqdr7HVcvnz5IX8X4+LiVFFR0TINdqCpU6fqjDPOUPv27dWpUycNHTpUGzduPOL7nnvuOXXr1k1JSUnq0aOHXn311RZorTMdzTXk5yLQePQ9jh79jqahv3H06Gc0DX2M0COkiFLPPPOMxo0bp//93//VmjVrdPrpp6ugoECbN2+u9z0pKSn64Ycf7Me3337bgi12nh07duj000/XnDlzGlT/66+/1pAhQzRgwACtXbtWY8eO1XXXXafXX3+9mVvqbI29jjU2btwY9PexU6dOzdRC53vrrbc0ZswYvfvuuyotLdXevXs1aNAg7dixo973rFq1Sr/73e80atQoffjhhxo6dKiGDh2qTz75pAVb7hxHcw0lfi4CjUHfo2nodzQN/Y2jRz+jaehjNAMLUenMM8+0xowZY7/ev3+/lZmZaU2dOrXO+k899ZSVmpraQq2LPJKsF1988bB17rjjDuuUU04JKvvtb39rFRQUNGPLIktDruObb75pSbJ++umnFmlTJNq8ebMlyXrrrbfqrXP55ZdbQ4YMCSrr06ePdcMNNzR38yJCQ64hPxeBxqHvETr0O5qG/kbT0M9oGvoYTcdIiii0Z88elZeXKz8/3y6Lj49Xfn6+ysrK6n3f9u3b1aVLF2VlZemSSy7R+vXrW6K5UaOsrCzomktSQUHBYa856tezZ0917txZ559/vlauXBnu5jjKtm3bJEkdOnSotw5/Hw+vIddQ4uci0FD0PVoeP+dDg/7GoehnNA19jKYjpIhCfr9f+/fvV3p6elB5enp6vfPsunbtqieffFIvvfSSiouLVV1drb59++rf//53SzQ5KlRUVNR5zauqqrRz584wtSrydO7cWXPnztXzzz+v559/XllZWerfv7/WrFkT7qY5QnV1tcaOHauzzz5bp556ar316vv7GItzbQ/W0GvIz0Wg4eh7tDz6HU1Df6Nu9DOahj5GaLQOdwPgDHl5ecrLy7Nf9+3bVx6PR3/+85917733hrFliDVdu3ZV165d7dd9+/bVl19+qUcffVR/+9vfwtgyZxgzZow++eQTvfPOO+FuSsRq6DXk5yLQvPg3hnCiv1E3+hlNQx8jNBhJEYVcLpdatWqlysrKoPLKykplZGQ06DPatGmjXr166YsvvmiOJkaljIyMOq95SkqKkpOTw9Sq6HDmmWfyd1FSYWGhXn75Zb355ps67rjjDlu3vr+PDf0ZEK0acw0Pxs9FoH70PVoe/Y7Qi/X+Bv2MpqGPETqEFFEoISFBubm5Wrp0qV1WXV2tpUuXBiV2h7N//36tW7dOnTt3bq5mRp28vLygay5JpaWlDb7mqN/atWtj+u+iZVkqLCzUiy++qGXLlik7O/uI7+HvY7CjuYYH4+ciUD/6Hi2Pn/OhF6v9DfoZTUMfoxmEd91ONJd//OMfVmJiojVv3jxrw4YN1ujRo620tDSroqLCsizLuvLKK63x48fb9SdPnmy9/vrr1pdffmmVl5dbw4cPt5KSkqz169eH6xTC7ueff7Y+/PBD68MPP7QkWdOnT7c+/PBD69tvv7Usy7LGjx9vXXnllXb9r776yjIMw7r99tst0zStOXPmWK1atbJKSkrCdQqO0Njr+Oijj1qLFi2yPv/8c2vdunXWrbfeasXHx1tvvPFGuE4h7G666SYrNTXVWr58ufXDDz/Yj0AgYNc5+N/0ypUrrdatW1sPP/ywZZqm9b//+79WmzZtrHXr1oXjFMLuaK4hPxeBxqHv0TT0O5qG/sbRo5/RNPQxQo+QIorNnj3bcrvdVkJCgnXmmWda7777rn3svPPOs0aOHGm/Hjt2rF03PT3duvDCC601a9aEodXOUbM11cGPmus2cuRI67zzzjvkPT179rQSEhKsE044wXrqqadavN1O09jr+OCDD1onnniilZSUZHXo0MHq37+/tWzZsvA03iHqun6Sgv5+Hfxv2rIs69lnn7V+8YtfWAkJCdYpp5xivfLKKy3bcAc5mmvIz0Wg8eh7HD36HU1Df+Po0c9oGvoYoRdnWZYV+vEZAAAAAAAAjcOaFAAAAAAAwBEIKQAAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByBkAIAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAmlX//v01duzYcDcj4sTFxWnRokXhbgYAABGDPsfRoc8BpyGkAJrZ1Vdfrbi4ON14442HHBszZozi4uJ09dVXt3zDHGb+/Pnq169fuJvRJPySBwCEE32OhqHPATgbIQXQArKysvSPf/xDO3futMt27dqlhQsXyu12h7FlDbNnz55m/46XXnpJv/rVr5r9ewAAiGb0OY6MPgfgbIQUQAvIyclRVlaWXnjhBbvshRdekNvtVq9evYLqVldXa+rUqcrOzlZycrJOP/10/fOf/7SP79+/X6NGjbKPd+3aVTNnzgz6jOXLl+vMM89U27ZtlZaWprPPPlvffvutpAN3WYYOHRpUf+zYserfv7/9un///iosLNTYsWPlcrlUUFAgSfrkk080ePBgtWvXTunp6bryyivl9/vt9+3YsUNXXXWV2rVrp86dO+uRRx5p0PXZtWuXlixZctgOw0svvaScnBwlJSXphBNO0OTJk7Vv3z5J0pQpU5SZmakff/zRrj9kyBANGDBA1dXVkg7ccXj88cc1ePBgJScn64QTTgi6rpK0adMmXX755UpLS1OHDh10ySWX6Jtvvgmq8+STT+qUU05RYmKiOnfurMLCQknS8ccfL0n69a9/rbi4OPv1kdouSZ9//rnOPfdcJSUlqXv37iotLW3QdQMA4GD0OQ6PPgd9DjgfIQXQQq699lo99dRT9usnn3xS11xzzSH1pk6dqqefflpz587V+vXrddttt8nr9eqtt96SdKBDcdxxx+m5557Thg0bVFRUpLvuukvPPvusJGnfvn0aOnSozjvvPH388ccqKyvT6NGjFRcX16j2zp8/XwkJCVq5cqXmzp2rrVu36pe//KV69eqlDz74QCUlJaqsrNTll19uv+f222/XW2+9pZdeeklLlizR8uXLtWbNmiN+19KlS/U///M/6tatW53H3377bV111VW69dZbtWHDBv35z3/WvHnzdP/990uS7r77bh1//PG67rrrJElz5szRqlWrNH/+fMXH//fH3MSJEzVs2DB99NFHGjFihIYPHy7TNCVJe/fuVUFBgdq3b6+3335bK1euVLt27XTBBRfYd3Uef/xxjRkzRqNHj9a6deu0ePFinXTSSZKk999/X5L01FNP6YcffrBfH6nt1dXVuvTSS5WQkKDVq1dr7ty5uvPOOxv+BwUAwEHoc9SPPgd9DkQAC0CzGjlypHXJJZdYmzdvthITE61vvvnG+uabb6ykpCTrP//5j3XJJZdYI0eOtCzLsnbt2mUZhmGtWrUq6DNGjRpl/e53v6v3O8aMGWMNGzbMsizL+vHHHy1J1vLlyw/bntpuvfVW67zzzrNfn3feeVavXr2C6tx7773WoEGDgso2bdpkSbI2btxo/fzzz1ZCQoL17LPP2sd//PFHKzk52br11lvrbbtlWdb1119v/fGPf6z3+MCBA60HHnggqOxvf/ub1blzZ/v1l19+abVv39668847reTkZGvBggVB9SVZN954Y1BZnz59rJtuusn+vK5du1rV1dX28d27d1vJycnW66+/blmWZWVmZlp33313ve2UZL344ouNavvrr79utW7d2vruu+/s46+99lqdnwUAwOHQ56DPQZ8D0aB1uMIRINYce+yxGjJkiObNmyfLsjRkyBC5XK6gOl988YUCgYDOP//8oPI9e/YEDdGcM2eOnnzySfl8Pu3cuVN79uxRz549JUkdOnTQ1VdfrYKCAp1//vnKz8/X5Zdfrs6dOzeqvbm5uUGvP/roI7355ptq167dIXW//PJLux19+vSxyzt06KCuXbse9nssy9K//vUv+65MXT766COtXLnSvhMgHRiCumvXLgUCARmGoRNOOEEPP/ywbrjhBv32t7/VFVdcccjn5OXlHfJ67dq19nd88cUXat++fVCdXbt26csvv9TmzZv1/fffa+DAgYc9n8a23TRNZWVlKTMzs952AgDQGPQ56kafgz4HIgMhBdCCrr32Wns+4Zw5cw45vn37dknSK6+8ov/5n/8JOpaYmChJ+sc//qE//vGPeuSRR5SXl6f27dvroYce0urVq+26Tz31lH7/+9+rpKREzzzzjO655x6VlpbqrLPOUnx8vCzLCvrsvXv3HtKWtm3bHtK2iy++WA8++OAhdTt37qwvvviiIZfgEO+995727dunvn371ltn+/btmjx5si699NJDjiUlJdnPV6xYoVatWumbb77Rvn371Lp1w3/Ebd++Xbm5uVqwYMEhx4499tigIZyN0dC2AwAQSvQ5DkWfA4gMhBRAC6qZaxgXF2cvDFVb9+7dlZiYKJ/Pp/POO6/Oz1i5cqX69u2rm2++2S778ssvD6nXq1cv9erVSxMmTFBeXp4WLlyos846S8cee6w++eSToLpr165VmzZtDtv2nJwcPf/88zr++OPr/EV84oknqk2bNlq9erW9evhPP/2kzz77rN5zkQ4s8DRkyBC1atXqsN+9ceNGey5mXZ555hm98MILWr58uS6//HLde++9mjx5clCdd999V1dddVXQ65q7RTk5OXrmmWfUqVMnpaSk1Pkdxx9/vJYuXaoBAwbUebxNmzbav39/o9ru8Xi0adMm/fDDD/adp3fffbfe8wQAoCHocxyKPgd9DkSIcM41AWLBwfMxt23bZm3bts1+XXt+qGVZ1t1332117NjRmjdvnvXFF19Y5eXl1qxZs6x58+ZZlmVZM2fOtFJSUqySkhJr48aN1j333GOlpKRYp59+umVZlvXVV19Z48ePt1atWmV988031uuvv2517NjR+tOf/mRZlmWVlJRYcXFx1vz5863PPvvMKioqslJSUg6ZH3rwnM7vvvvOOvbYY63f/OY31nvvvWd98cUXVklJiXX11Vdb+/btsyzLsm688UarS5cu1tKlS61169ZZv/rVr6x27doddn7oKaecYj3//POHvYYlJSVW69atrUmTJlmffPKJtWHDBuvvf/+7PVdz06ZN1jHHHGPNmjUrqH5ZWZn9GZIsl8tl/fWvf7U2btxoFRUVWfHx8db69esty7KsHTt2WCeffLLVv39/a8WKFdZXX31lvfnmm9Ytt9xibdq0ybIsy5o3b56VlJRkzZw50/rss8/sP5saJ598snXTTTdZP/zwg7Vly5YGtX3//v1W9+7drfPPP99au3attWLFCis3N5f5oQCARqPPQZ+DPgeiASEF0MzqWjSqtoM7DNXV1daMGTOsrl27Wm3atLGOPfZYq6CgwHrrrbcsyzqw0NXVV19tpaamWmlpadZNN91kjR8/3u4wVFRUWEOHDrU6d+5sJSQkWF26dLGKioqs/fv3299RVFRkpaenW6mpqdZtt91mFRYWHrHDYFmW9dlnn1m//vWvrbS0NCs5Odnq1q2bNXbsWHvhp59//tnyer2WYRhWenq6NW3atHo/y7Is64svvrASExOt7du3H/E6lpSUWH379rWSk5OtlJQU68wzz7SeeOIJq7q62ho4cKBVUFAQtADVLbfcYp144onWzz//bFnWgQ7DnDlzrPPPP99KTEy0jj/+eOuZZ54J+o4ffvjBuuqqqyyXy2UlJiZaJ5xwgnX99dcHdfDmzp1r/9l07tzZuuWWW+xjixcvtk466SSrdevWVpcuXY7Y9hobN260+vXrZyUkJFi/+MUvrJKSEjoMAIBGo89Bn4M+B6JBnGUdNFEMAFrI9OnT9cYbb+jVV19t9u+Ki4vTiy++eMh+7QAAIPrR5wAix9GtygIAIXDcccdpwoQJ4W4GAACIcvQ5gMjBwpkAwubyyy8PdxMAAEAMoM8BRA6mewAAAAAAAEdgugcAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOQEgBAAAAAAAcgZACAAAAAAA4AiEFAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBEIKQAAAAAAgCMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByBkAJAizn++ON19dVXh7sZAAAAjbZ8+XLFxcVp+fLl4W4KENUIKQCExLx58xQXF1fnY/z48U367O+//16TJk3S2rVrQ9NYAAAQEerrWxz8aEhw8MADD2jRokXN3mYATdM63A0AEF2mTJmi7OzsoLJTTz21SZ/5/fffa/LkyTr++OPVs2fPJn0WAACIHH/729+CXj/99NMqLS09pNzj8Rzxsx544AH95je/0dChQ0PZRAAhRkgBIKQGDx6s3r17h7sZAAAgCni93qDX7777rkpLSw8pBxA9mO4BIGy2bNmiP/7xj+rRo4fatWunlJQUDR48WB999JFdZ/ny5TrjjDMkSddcc409rHPevHlhajUAAHCSHTt26A9/+IOysrKUmJiorl276uGHH5ZlWXaduLg47dixQ/Pnz7f7EjXrZH377be6+eab1bVrVyUnJ6tjx4667LLL9M0334TnhIAYx0gKACG1bds2+f3+oDKXy1Vn3a+++kqLFi3SZZddpuzsbFVWVurPf/6zzjvvPG3YsEGZmZnyeDyaMmWKioqKNHr0aJ1zzjmSpL59+zb7uQAAAGezLEu/+tWv9Oabb2rUqFHq2bOnXn/9dd1+++367rvv9Oijj0o6MG3kuuuu05lnnqnRo0dLkk488URJ0vvvv69Vq1Zp+PDhOu644/TNN9/o8ccfV//+/bVhwwYZhhG28wNiESEFgJDKz88/pKz2nYzaevTooc8++0zx8f8d1HXllVeqW7du+utf/6qJEycqPT1dgwcPVlFRkfLy8hjeCQAAbIsXL9ayZct033336e6775YkjRkzRpdddplmzpypwsJCnXjiifJ6vbrxxht1wgknHNKXGDJkiH7zm98ElV188cXKy8vT888/ryuvvLLFzgcA0z0AhNicOXNUWloa9KhPYmKiHVDs379fP/74o9q1a6euXbtqzZo1LdVkAAAQoV599VW1atVKv//974PK//CHP8iyLL322mtH/Izk5GT7+d69e/Xjjz/qpJNOUlpaGv0RIAwYSQEgpM4888wGL5xZXV2tmTNn6k9/+pO+/vpr7d+/3z7WsWPH5moiAACIEt9++60yMzPVvn37oPKa3T6+/fbbI37Gzp07NXXqVD311FP67rvvgkaAbtu2LbQNBnBEhBQAwuaBBx7QxIkTde211+ree+9Vhw4dFB8fr7Fjx6q6ujrczQMAADHglltu0VNPPaWxY8cqLy9PqampiouL0/Dhw+mPAGFASAEgbP75z39qwIAB+utf/xpUvnXr1qDFNuPi4lq6aQAAIAJ06dJFb7zxhn7++eeg0RSffvqpfbxGff2Jf/7znxo5cqQeeeQRu2zXrl3aunVr8zQawGGxJgWAsGnVqtUhi2o+99xz+u6774LK2rZtK0l0FgAAQJALL7xQ+/fv12OPPRZU/uijjyouLk6DBw+2y9q2bVtnX6Ku/sjs2bODpqECaDmMpAAQNhdddJGmTJmia665Rn379tW6deu0YMECnXDCCUH1TjzxRKWlpWnu3Llq37692rZtqz59+ig7OztMLQcAAE5w8cUXa8CAAbr77rv1zTff6PTTT9eSJUv00ksvaezYsfY2o5KUm5urN954Q9OnT1dmZqays7PVp08fXXTRRfrb3/6m1NRUde/eXWVlZXrjjTdYHwsIE0IKAGFz1113aceOHVq4cKGeeeYZ5eTk6JVXXtH48eOD6rVp00bz58/XhAkTdOONN2rfvn166qmnCCkAAIhx8fHxWrx4sYqKivTMM8/oqaee0vHHH6+HHnpIf/jDH4LqTp8+XaNHj9Y999yjnTt3auTIkerTp49mzpypVq1aacGCBdq1a5fOPvtsvfHGGyooKAjTWQGxLc46eGwTAAAAAABAGLAmBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBEIKQAAAAAAgCMQUgAAAAAAAEdodEixYsUKXXzxxcrMzFRcXJwWLVoUdNyyLBUVFalz585KTk5Wfn6+Pv/886A6W7Zs0YgRI5SSkqK0tDSNGjVK27dvD6rz8ccf65xzzlFSUpKysrI0bdq0Q9ry3HPPqVu3bkpKSlKPHj306quvNvZ0AAAAAACAQ7Ru7Bt27Nih008/Xddee60uvfTSQ45PmzZNs2bN0vz585Wdna2JEyeqoKBAGzZsUFJSkiRpxIgR+uGHH1RaWqq9e/fqmmuu0ejRo7Vw4UJJUlVVlQYNGqT8/HzNnTtX69at07XXXqu0tDSNHj1akrRq1Sr97ne/09SpU3XRRRdp4cKFGjp0qNasWaNTTz21QedSXV2t77//Xu3bt1dcXFxjLwUAABHPsiz9/PPPyszMVHw8AyybC30OAECsa3Cfw2oCSdaLL75ov66urrYyMjKshx56yC7bunWrlZiYaP3973+3LMuyNmzYYEmy3n//fbvOa6+9ZsXFxVnfffedZVmW9ac//ck65phjrN27d9t17rzzTqtr167268svv9waMmRIUHv69Olj3XDDDQ1u/6ZNmyxJPHjw4MGDR8w/Nm3a1ODfn2g8+hw8ePDgwYPHgceR+hyNHklxOF9//bUqKiqUn59vl6WmpqpPnz4qKyvT8OHDVVZWprS0NPXu3duuk5+fr/j4eK1evVq//vWvVVZWpnPPPVcJCQl2nYKCAj344IP66aefdMwxx6isrEzjxo0L+v6CgoJDpp/Utnv3bu3evdt+bVmWJGnTpk1KSUlp6ukD9duxQ8rMPPD8+++ltm3D2x4A+H+qqqqUlZWl9u3bh7spUa3m+tLnABqIvhMQdRra5whpSFFRUSFJSk9PDypPT0+3j1VUVKhTp07BjWjdWh06dAiqk52dfchn1Bw75phjVFFRcdjvqcvUqVM1efLkQ8pTUlLoMKB5tWr13+cpKfyiBeA4TEFoXjXXlz4H0ED0nYCodaQ+R0xNPp0wYYK2bdtmPzZt2hTuJgEAAAAAgP8npCFFRkaGJKmysjKovLKy0j6WkZGhzZs3Bx3ft2+ftmzZElSnrs+o/R311ak5XpfExET7DgZ3MgAAAAAAcJaQhhTZ2dnKyMjQ0qVL7bKqqiqtXr1aeXl5kqS8vDxt3bpV5eXldp1ly5apurpaffr0seusWLFCe/futeuUlpaqa9euOuaYY+w6tb+npk7N9wAAAAAAgMjS6JBi+/btWrt2rdauXSvpwGKZa9eulc/nU1xcnMaOHav77rtPixcv1rp163TVVVcpMzNTQ4cOlSR5PB5dcMEFuv766/Xee+9p5cqVKiws1PDhw5X5/xbHueKKK5SQkKBRo0Zp/fr1euaZZzRz5syghTJvvfVWlZSU6JFHHtGnn36qSZMm6YMPPlBhYWHTrwoAAAAAAGhxjV4484MPPtCAAQPs1zXBwciRIzVv3jzdcccd2rFjh0aPHq2tW7eqX79+KikpUVJSkv2eBQsWqLCwUAMHDlR8fLyGDRumWbNm2cdTU1O1ZMkSjRkzRrm5uXK5XCoqKtLo0aPtOn379tXChQt1zz336K677tLJJ5+sRYsW6dRTTz2qCwEAAAAAAMIrzqrZhzMGVVVVKTU1Vdu2bWN9CjSvHTukdu0OPN++nRWqATgGvwtbBtcZaCT6TkDUaejvwpja3QMAAAAAADgXIQUAAHCkFStW6OKLL1ZmZqbi4uK0aNGioOOWZamoqEidO3dWcnKy8vPz9fnnnwfV2bJli0aMGKGUlBSlpaVp1KhR2r59e1Cdjz/+WOecc46SkpKUlZWladOmHdKW5557Tt26dVNSUpJ69OihV199NeTnCwAACCkAAIBD7dixQ6effrrmzJlT5/Fp06Zp1qxZmjt3rlavXq22bduqoKBAu3btsuuMGDFC69evV2lpqV5++WWtWLEiaI2rqqoqDRo0SF26dFF5ebkeeughTZo0SU888YRdZ9WqVfrd736nUaNG6cMPP9TQoUM1dOhQffLJJ8138gAAxCjWpGB+KFoC8yoBOFSk/C6Mi4vTiy++aO8WZlmWMjMz9Yc//EF//OMfJUnbtm1Tenq65s2bp+HDh8s0TXXv3l3vv/++evfuLUkqKSnRhRdeqH//+9/KzMzU448/rrvvvlsVFRVKSEiQJI0fP16LFi3Sp59+Kkn67W9/qx07dujll1+223PWWWepZ8+emjt3boPaHynXGXAM+k5A1GFNCgAAELW+/vprVVRUKD8/3y5LTU1Vnz59VFZWJkkqKytTWlqaHVBIUn5+vuLj47V69Wq7zrnnnmsHFJJUUFCgjRs36qeffrLr1P6emjo131OX3bt3q6qqKugBAACOjJACAABEnIqKCklSenp6UHl6erp9rKKiQp06dQo63rp1a3Xo0CGoTl2fUfs76qtTc7wuU6dOVWpqqv3Iyspq7CkCABCTCCkAAABCbMKECdq2bZv92LRpU7ibBABARCCkAAAAEScjI0OSVFlZGVReWVlpH8vIyNDmzZuDju/bt09btmwJqlPXZ9T+jvrq1ByvS2JiolJSUoIeAADgyAgpAABAxMnOzlZGRoaWLl1ql1VVVWn16tXKy8uTJOXl5Wnr1q0qLy+36yxbtkzV1dXq06ePXWfFihXau3evXae0tFRdu3bVMcccY9ep/T01dWq+BwAAhA4hBQAAcKTt27dr7dq1Wrt2raQDi2WuXbtWPp9PcXFxGjt2rO677z4tXrxY69at01VXXaXMzEx7BxCPx6MLLrhA119/vd577z2tXLlShYWFGj58uDIzMyVJV1xxhRISEjRq1CitX79ezzzzjGbOnKlx48bZ7bj11ltVUlKiRx55RJ9++qkmTZqkDz74QIWFhS19SQAAiHqtw90AAACAunzwwQcaMGCA/bomOBg5cqTmzZunO+64Qzt27NDo0aO1detW9evXTyUlJUpKSrLfs2DBAhUWFmrgwIGKj4/XsGHDNGvWLPt4amqqlixZojFjxig3N1cul0tFRUUaPXq0Xadv375auHCh7rnnHt111106+eSTtWjRIp166qktcBUAAIgtcZZlWeFuRLiwZzkaw+fzye/3269dLpfcbnfD3sxe3wAcit+FLYPrDAQ7uF8lHdS3ou8ERJ2G/i5kJAXQAD6fTx6PR4FAwC4zDEOmaTY8qAAAAIhRtUOJ//znP7r00kuD+lUSfSsABxBSAA3g9/sVCARUXFwsj8cj0zTl9Xrl9/v5RQoAAHAY9d3sKSkp0bHHHitJdt/q7bfflsfjUfzOneoZpvYCCC9CCqARPB6PcnJy6jzWpOkgAAAAUergmz3Sof0kl8slwzDk9XolSYakHeFoLICwI6QAmsA0TUl1D1tkyCIAAMB/He5mj9vtlmma9g2fzz78ULruupZsHgCHIKQAjsLBab8UPGyR6SAAAACN43a77X5T/M6dYW4NgHAhpACOwsFpv8T0DgAAAABoKkIK4CjVTvsBAAAAAE1HSAHUo/ZCmDVrTwAAAAAAmg8hBfD/HGn/bsMw5HK5GvWZNeEG22gBAAAAwJERUgBq2P7djVlzgm20AAAAAKDxCCkANWz/7sZgGy0AAIDQWLt2raqTk1mkHIgRhBRALYfbv7ux2EYLAADg6HTs2NF+fna/fgrowChX0zQJKoAoFx/uBgAAAABAbVlZWfbzle+8o+LiYgUCgaDt3wFEJ0ZSAAAAAHCsnj17qjo5OdzNANBCGEkBAAAAAAAcgZEUAAAAACJCzfbuUtMWOQfgXIQUAAAAABzt4O3dJRbSBKIVIQUAAAAARzt4e3fTNOX1euX3+wkpgChDSAG0MPb6BgAA0cbn89W780btKRpNUXt7dwDRi5ACaAHs9Q0AAKKVz+eTx+NRIBCot45hGHK5XC3YKgCRipACaAEH7/W9/ptvGKIIAACigt/vVyAQUHFxsTweT511GEEKoKEIKYAWxl7fAAAgGnk8HuXk5IS7GQAiXHy4GwAAAAAAACARUgAAAAAAAIcgpAAAAAAAAI5ASAEAAAAAAByBkAIAAAAAADgCu3sAAAAAiEimadrP2eYUiA6EFAAAAAAiisvlkmEY8nq9dplhGDJNk6ACiHCEFAAAAAAiitvtlmma8vv9kg6MqPB6vfL7/YQUQIQjpEDM8vl8Qb/YwuHg72WYIgAAQMO43W76TUAUIqRATPL5fPJ4PAoEAnaZYRhyuVwt8v11DVGsaQPDFAEAAADEKkIKxCS/369AIKDi4mJ5PB5JLTuK4eAhihLDFAEAAACAkAIxzePxKCcnJyzfzRBFAACA0GK3DyDyEVIAAAAAaBQnrO1VG7t9ANGDkAIAAABAg4V7ba+6sNsHED0IKQAAAAA0WLjX9qoPU2mB6EBIAQAAAKDRwrm2F4DoFR/uBgAAAAAAAEiEFAAAAAAAwCGY7oGY4bRVqAEAAAAAwQgpEBOcuAo1AAAAACAYIQViglNXoQYAAAAA/BchBWIKq1ADAAAAgHMRUgAOU3u9DEZ7AAAAAIglhBSAQ7hcLhmGIa/Xa5cZhiHTNAkqAAAAAMQEQgrAIdxut0zTDNqBxOv1yu/3E1IAAAAAiAmEFICDuN1uAgkAAAAAMSs+3A0AAAAAAACQGEkBAAAA4Ah8Pl/QlFQAaC6EFAAAAADq5fP55PF4FAgE7DLDMORyucLYKgDRipACAAAAQL38fr8CgYCKi4vl8XgksU06gOZDSAEAAADgiDwej3JycsLdDABRjoUzAQAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcAS2IAUAAAAQlUzTtJ+7XC653e4wtgZAQxBSAAAAAIgqLpdLhmHI6/XaZYZhyDRNggrA4QgpAIfjDgAAAEDjuN1umaYpv98v6UB/yuv1yu/305cCHI6QAlHJ5/PZv5Sk4P/oRwruAAAAABw9t9tNnwmIQIQUiDo+n08ej0eBQCCo3DAMuVyuMLWq8bgDAAAAACDWEFIg6vj9fgUCARUXF8vj8djlkThVgjsAAAAAAGIJIQWilsfjUU5OTribAQAAAABooPhQf+D+/fs1ceJEZWdnKzk5WSeeeKLuvfdeWZZl17EsS0VFRercubOSk5OVn5+vzz//POhztmzZohEjRiglJUVpaWkaNWqUtm/fHlTn448/1jnnnKOkpCRlZWVp2rRpoT4dAAAAAADQQkIeUjz44IN6/PHH9dhjj8k0TT344IOaNm2aZs+ebdeZNm2aZs2apblz52r16tVq27atCgoKtGvXLrvOiBEjtH79epWWlurll1/WihUrNHr0aPt4VVWVBg0apC5duqi8vFwPPfSQJk2apCeeeCLUpwQAAAAAAFpAyKd7rFq1SpdccomGDBkiSTr++OP197//Xe+9956kA6MoZsyYoXvuuUeXXHKJJOnpp59Wenq6Fi1apOHDh8s0TZWUlOj9999X7969JUmzZ8/WhRdeqIcffliZmZlasGCB9uzZoyeffFIJCQk65ZRTtHbtWk2fPj0ozAAAAAAAAJEh5CMp+vbtq6VLl+qzzz6TJH300Ud65513NHjwYEnS119/rYqKCuXn59vvSU1NVZ8+fVRWViZJKisrU1pamh1QSFJ+fr7i4+O1evVqu865556rhIQEu05BQYE2btyon376qc627d69W1VVVUEPAAAAAADgDCEfSTF+/HhVVVWpW7duatWqlfbv36/7779fI0aMkCRVVFRIktLT04Pel56ebh+rqKhQp06dghvaurU6dOgQVCc7O/uQz6g5dswxxxzStqlTp2ry5MkhOEsgfEzTtJ9H4o4lAAAAAFCfkIcUzz77rBYsWKCFCxfaUzDGjh2rzMxMjRw5MtRf1ygTJkzQuHHj7NdVVVXKysoKY4uAhnO5XDIMQ16v1y4zDEOmaRJUAAAAAIgKIQ8pbr/9do0fP17Dhw+XJPXo0UPffvutpk6dqpEjRyojI0OSVFlZqc6dO9vvq6ysVM+ePSVJGRkZ2rx5c9Dn7tu3T1u2bLHfn5GRocrKyqA6Na9r6hwsMTFRiYmJTT9JIAzcbrdM05Tf75d0YESF1+uV3+8npAAAAAAQFUK+JkUgEFB8fPDHtmrVStXV1ZKk7OxsZWRkaOnSpfbxqqoqrV69Wnl5eZKkvLw8bd26VeXl5XadZcuWqbq6Wn369LHrrFixQnv37rXrlJaWqmvXrnVO9QCigdvtVk5OjnJycuTxeMLdHAAAEKV8Pp/WrFmjNWvWBE01BYDmFvKRFBdffLHuv/9+ud1unXLKKfrwww81ffp0XXvttZKkuLg4jR07Vvfdd59OPvlkZWdna+LEicrMzNTQoUMlSR6PRxdccIGuv/56zZ07V3v37lVhYaGGDx+uzMxMSdIVV1yhyZMna9SoUbrzzjv1ySefaObMmXr00UdDfUoAAABAzPD5fPJ4PAoEAnaZYRhyuVxhbBWAWBHykGL27NmaOHGibr75Zm3evFmZmZm64YYbVFRUZNe54447tGPHDo0ePVpbt25Vv379VFJSoqSkJLvOggULVFhYqIEDByo+Pl7Dhg3TrFmz7OOpqalasmSJxowZo9zcXLlcLhUVFbH9KAAAANAEfr9fgUBAxcXF9shNFusG0FJCPt2jffv2mjFjhr799lvt3LlTX375pe67776grULj4uI0ZcoUVVRUaNeuXXrjjTf0i1/8IuhzOnTooIULF+rnn3/Wtm3b9OSTT6pdu3ZBdU477TS9/fbb2rVrl/7973/rzjvvDPXpAAAAh9q/f78mTpyo7OxsJScn68QTT9S9994ry7LsOpZlqaioSJ07d1ZycrLy8/P1+eefB33Oli1bNGLECKWkpCgtLU2jRo3S9u3bg+p8/PHHOuecc5SUlKSsrCxNmzatRc4RCCePx2NPMyWgANBSQh5SAAAAtIQHH3xQjz/+uB577DGZpqkHH3xQ06ZN0+zZs+0606ZN06xZszR37lytXr1abdu2VUFBgXbt2mXXGTFihNavX6/S0lK9/PLLWrFiRdDIzKqqKg0aNEhdunRReXm5HnroIU2aNElPPPFEi54vgKYzTdNea8Pn84W7OQDqEPLpHgAAAC1h1apVuuSSSzRkyBBJ0vHHH6+///3veu+99yQdGEUxY8YM3XPPPbrkkkskSU8//bTS09O1aNEiDR8+XKZpqqSkRO+//7569+4t6cDU1QsvvFAPP/ywMjMztWDBAu3Zs0dPPvmkEhIS7C3Wp0+fzjRTIEKwlTsQORhJAQAAIlLfvn21dOlSffbZZ5Kkjz76SO+8844GDx4sSfr6669VUVGh/Px8+z2pqanq06ePysrKJEllZWVKS0uzAwpJys/PV3x8vFavXm3XOffcc4OmrhYUFGjjxo366aef6mzb7t27VVVVFfQAED41W7mXl5ervLxcxcXFCgQC9tbuAJyDkRQAACAijR8/XlVVVerWrZtatWql/fv36/7779eIESMkSRUVFZKk9PT0oPelp6fbxyoqKtSpU6eg461bt1aHDh2C6mRnZx/yGTXH6tr6fOrUqZo8eXIIzhJAqLjdbkZNABGAkRQAACAiPfvss1qwYIEWLlyoNWvWaP78+Xr44Yc1f/78cDdNEyZM0LZt2+zHpk2bwt0kAAAiAiMpAABARLr99ts1fvx4DR8+XJLUo0cPffvtt5o6dapGjhypjIwMSVJlZaU6d+5sv6+yslI9e/aUJGVkZGjz5s1Bn7tv3z5t2bLFfn9GRoYqKyuD6tS8rqlzsMTERCUmJjb9JAEAiDGMpEBU8Pl89krNpmmGuzkAgBYQCAQUHx/clWnVqpWqq6slSdnZ2crIyNDSpUvt41VVVVq9erXy8vIkSXl5edq6davKy8vtOsuWLVN1dbX69Olj11mxYoX27t1r1yktLVXXrl3rnOoBAACOHiMpEPF8Pp88Ho8CgYBdZhiGXC5XGFsFAGhuF198se6//3653W6dcsop+vDDDzV9+nRde+21kqS4uDiNHTtW9913n04++WRlZ2dr4sSJyszM1NChQyVJHo9HF1xwga6//nrNnTtXe/fuVWFhoYYPH67MzExJ0hVXXKHJkydr1KhRuvPOO/XJJ59o5syZevTRR8N16gAARC1CCkQ8v9+vQCCg4uJieTweSQe2mWJhJACIbrNnz9bEiRN18803a/PmzcrMzNQNN9ygoqIiu84dd9yhHTt2aPTo0dq6dav69eunkpISJSUl2XUWLFigwsJCDRw4UPHx8Ro2bJhmzZplH09NTdWSJUs0ZswY5ebmyuVyqaioiO1HAQBoBoQUiBoej0c5OTnhbgYAoIW0b99eM2bM0IwZM+qtExcXpylTpmjKlCn11unQoYMWLlx42O867bTT9Pbbbx9tUwEAQAOxJgUAAAAAAHAERlIAAAAAMc7n88nv90sSi5ADCCtCCgAAACCGsQg5ACchpAAAAABiGIuQA3ASQgoAAAAALEIOwBFYOBMAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIhBQAAAAAAMARWDgTiHC19zJnJW4AAAAAkYyQAohQLpdLhmHI6/XaZYZhyDRNggoAAAAAEYmQAohQbrdbpmnK7/dLOjCiwuv1yu/3E1IAAAAAiEiEFEAEc7vdBBIAAAAAogYhBQAAAICYxNpegPMQUgAAAACIKaztBTgXIQUAAACAmMLaXoBzEVIAAAAAiDms7QU4U3y4GwAAAAAAACARUgAAAAAAAIdgugciks/nC5pDCAAAAACIfIQUiDg+n08ej0eBQMAuMwxDLpcrjK0CAAAAADQVIQUijt/vVyAQUHFxsTwejyT2tQYAAGgMRqUCcCpCCkQsj8ejnJyccDcDAAAgojAqFYCTEVIAAAAAMYRRqQCcjJACAAAAiEGMSgXgRGxBCgAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBHY3QOIMqZp2s/ZTgwAAABAJCGkAKKEy+WSYRjyer12mWEYMk2ToAIAAABARCCkAKKE2+2WaZry+/2SDoyo8Hq98vv9hBQAAAAAIgIhBRBF3G43gQQAAACAiMXCmQAAAAAAwBEIKQAAAAAAgCMw3QMAAAAAxC5pgBMQUgAAAACIaeySBjgHIQUAAACAmMYuaYBzEFIAAAAAiHnskgY4AwtnAgAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcAR290BE8Pl8QVtCAQAAAACiDyEFHM/n88nj8SgQCNhlhmHI5XKFsVUAAAAAgFAjpIDj+f1+BQIBFRcXy+PxSJJcLhf7WAMAAABAlCGkQMTweDzKyckJdzMAAAAAAM2EhTMBAAAAAIAjMJICiHK1FxplmgwAAAAAJyOkAKKUy+WSYRjyer12mWEYMk2ToAIAAACAIxFSAFHK7XbLNM2grVu9Xq/8fj8hBQAAAABHIqQAopjb7SaQAAAAOEq1p81KTJ0FWgIhBQAAAADUUte0WYmps0BLIKQAAAAAopzP5wuaAorDO3jarMTUWaClEFIAAAD8/+3deXxU1f3/8XdYsgwQEAYCaAZRWhjUAomCEb4FFY1t2orS4jaIiFJpsAKPqlARBVupK25BtFWgjSj6FdwbRBRQFrEBLMsIAurwswYclwQYCUvO7w++3GZCEpIwk3tn5vV8PObxyNx7MvmcIcw9ed97zgXiWCAQkNfrVSgUsra5XC653W4bq3I+ps0C9iCkAAAAAOJYMBhUKBRSYWGhvF6vJNZWAOBchBQAAABAAvB6vcrKyrK7DACoVRO7CwAAAAAAAJAIKQAAAAAAgEMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAAByBkAIAAAAAADgCIQUAAAAAAHAEQgoAAAAAAOAIzewuAEDj8vv91tdut1sej8fGagAAAADgv6JyJcWXX34pn8+ndu3aKS0tTWeddZb+9a9/WfuNMZoyZYo6deqktLQ0DR48WJ9++mnYa3z77be65pprlJ6erjZt2mjUqFHau3dvWJt///vf+p//+R+lpqYqMzNT999/fzS6A8QFt9stl8sln8+n7OxsZWdny+v1KhAI2F0aAAAAAEiKQkjx3XffqX///mrevLn++c9/avPmzXrooYd00kknWW3uv/9+PfbYY5o1a5Y+/PBDtWjRQrm5udq/f7/V5pprrtGmTZu0ePFivfHGG1q+fLlGjx5t7S8rK9PFF1+sLl26qLi4WA888IDuvvtuPf3005HuEmwQCAS0du1arV27NuzMPxrO4/HI7/eruLhYxcXFKiwsVCgUUjAYtLs0AAAAAJAUheke9913nzIzMzV79mxrW9euXa2vjTF65JFHNHnyZF166aWSpL///e/KyMjQK6+8oiuvvFJ+v19FRUX66KOPdPbZZ0uSHn/8cf385z/Xgw8+qM6dO+u5557TgQMH9Oyzzyo5OVlnnHGG1q9fr4cffjgszEDsCQQC8nq9CoVC1jaXyyW3221jVfHB4/EwvQMAAACAY0X8SorXXntNZ599tn7zm9+oQ4cO6tOnj/76179a+z/77DOVlJRo8ODB1rbWrVurX79+WrVqlSRp1apVatOmjRVQSNLgwYPVpEkTffjhh1abn/70p0pOTrba5ObmasuWLfruu++qra28vFxlZWVhDzhPMBhUKBRSYWGhddbf7/fzxzUAAAAAxLmIhxQ7duzQk08+qR/96EdatGiRxowZo9///veaO3euJKmkpESSlJGREfZ9GRkZ1r6SkhJ16NAhbH+zZs3Utm3bsDbVvUbln1HV9OnT1bp1a+uRmZl5gr1FNHm9XmVlZSkrK4uAAgAAAAASQMRDioqKCmVlZenee+9Vnz59NHr0aN14442aNWtWpH9UvU2aNEmlpaXWY+fOnXaXBAAAAAAA/k/EQ4pOnTqpZ8+eYdsq30GgY8eOkqRdu3aFtdm1a5e1r2PHjtq9e3fY/kOHDunbb78Na1Pda1T+GVWlpKQoPT097AEAAAAAAJwh4iFF//79tWXLlrBtW7duVZcuXSQdWUSzY8eOWrJkibW/rKxMH374oXJyciRJOTk5+v7771VcXGy1effdd1VRUaF+/fpZbZYvX66DBw9abRYvXqzu3buH3UkEAADEL257DgBAfIl4SDF+/HitXr1a9957r7Zt26Z58+bp6aefVn5+viQpKSlJ48aN05/+9Ce99tpr2rBhg6699lp17txZQ4YMkXTkyotLLrlEN954o9asWaMVK1Zo7NixuvLKK9W5c2dJ0tVXX63k5GSNGjVKmzZt0vz58/Xoo49qwoQJke4SAABwIG57DsAOfr9fa9eu1dq1a62rxQFEkImC119/3Zx55pkmJSXF9OjRwzz99NNh+ysqKsydd95pMjIyTEpKirnwwgvNli1bwtp888035qqrrjItW7Y06enpZuTIkWbPnj1hbT7++GMzYMAAk5KSYk4++WTzl7/8pV51lpaWGkmmtLS0YR1FVBQXFxtJpri42O5SImfvXmOkI4+9e+2uxhgTp+8zgHqL5WPh7bffbgYMGFDj/oqKCtOxY0fzwAMPWNu+//57k5KSYp5//nljjDGbN282ksxHH31ktfnnP/9pkpKSzJdffmmMMWbmzJnmpJNOMuXl5WE/u3v37nWuNZbfZ8S+mDzmO3Ds9MUXXxiXy2UkWQ+Xy2W++OILu0sDYkJdj4XNohF8/OIXv9AvfvGLGvcnJSVp2rRpmjZtWo1t2rZtq3nz5tX6c37yk5/o/fffb3CdAAAgdr322mvKzc3Vb37zGy1btkwnn3yyfve73+nGG2+UdPzbnl955ZXHve35ZZddVuNtz++77z599913TDMFEoTH45Hf71cwGJR05IoKn8+nYDDIneiACIpKSAEAABBtR297PmHCBP3xj3/URx99pN///vdKTk7WiBEjInrb865dux7zGkf3VRdSlJeXq7y83HpeVlZ2gr0F6i4QCFh/SEtH/phGZHg8HgIJIMoIKQAAQEyqqKjQ2WefrXvvvVeS1KdPH23cuFGzZs3SiBEjbK1t+vTpmjp1qq01IDEFAgF5vV6FQqGw7S6XS26326aqAKDuIr5wJgAAQGNw8m3PJ02apNLSUuuxc+fOhnQRqLdgMKhQKKTCwkIVFxdbD7/fzxUAAGICIQUAAIhJTr7teUpKitLT08MeQGPyer3KysqyHgQUAGIFIQUAAIhJ3PYcAID4w5oUAAAgJp1zzjlauHChJk2apGnTpqlr16565JFHdM0111htbrvtNu3bt0+jR4/W999/rwEDBqioqEipqalWm+eee05jx47VhRdeqCZNmmjo0KF67LHHrP2tW7fW22+/rfz8fGVnZ8vtdmvKlCkaPXp0o/YXAIBEQEgBAABiFrc9BwAgvjDdAwAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AisSQEkOL/fb33tdru5RRkAAAAA2xBSAAnK7XbL5XLJ5/NZ21wul/x+P0EFAAAAAFsQUgAJyuPxyO/3KxgMSjpyRYXP51MwGCSkAAAAAGALQgoggXk8HgIJAAAAAI7BwpkAAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOQEgBAAAAAAAcgZACAAAAAAA4AiEFAAAAAABwhGZ2FwAAAAAAscrv91tfu91ueTweG6sBYh8hBQAAAADUk9vtlsvlks/ns7a5XC75/X6CCuAEEFIAAAAAQD15PB75/X4Fg0FJR66o8Pl8CgaDhBTACSCkgCMEAoGwD3gAAADUDeMo+3g8HgIJIMIIKWC7QCAgr9erUChkbXO5XHK73TZWBQAA4HyMowDEG0IK2C4YDCoUCqmwsFBer1cSiw4BAADUBeMoAPGGkAKO4fV6lZWVZXcZAAAAMYdxFIB40cTuAgAAAAAAACSupABQBff6BgAAAGAXQgoAkrjXNwAAAAD7EVIAkMS9vgEAAADYj5ACgIV7fQMAAACwEwtnAgAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4AiEFAAAAAAAwBEIKQAAAAAAgCM0s7sAAAAAAIgXfr/f+trtdsvj8dhYDRB7CCkAAAAA4AS53W65XC75fD5rm8vlkt/vJ6gA6oGQAgAAAABOkMfjkd/vVzAYlHTkigqfz6dgMEhIAdQDIQUAAAAARIDH4yGQAE4QC2cCAAAAAABH4EoK2CIQCIRdCgcAAAAAACEFGl0gEJDX61UoFLK2uVwuud1uG6sCAACIDZzsARDPCCnQ6ILBoEKhkAoLC+X1eiVxeyYnqzr44d8KAAD7cLIHQLwjpIBtvF6vsrKy7C4DNajuNloSt9ICAMBOnOwBEO8IKQBUq+pttCRupQUAgFNwsgdAvCKkAFAjbqMFAAAAoDFxC1IAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOQEgBAAAAAAAcgZACAAAAAAA4AiEFAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjtDM7gIAxB6/32997Xa75fF4bKwGAAAAQLwgpABQZ263Wy6XSz6fz9rmcrnk9/sJKgAAAACcMEIKAHXm8Xjk9/sVDAYlHbmiwufzKRgMElIAAAAAOGGEFADqxePxEEgAAAAAiApCCgAAAMDBAoFA2FWMiC2s5QXUDyEFAAAA4FCBQEBer1ehUMja5nK55Ha7bawKdcFaXkDDEFIAAAAADhUMBhUKhVRYWCiv1yuJs/GxgrW8gIYhpAAAAAAczuv1Kisry+4yUE+s5QXUXxO7CwAAAAAAAJAIKQAAAAAAgEMQUgAAAAAAAEcgpAAAAAAAAI5ASAEAAAAAABwh6iHFX/7yFyUlJWncuHHWtv379ys/P1/t2rVTy5YtNXToUO3atSvs+wKBgPLy8uRyudShQwfdeuutOnToUFibpUuXKisrSykpKerWrZvmzJkT7e4AAAAAAIAoiWpI8dFHH+mpp57ST37yk7Dt48eP1+uvv66XXnpJy5Yt03/+8x9dfvnl1v7Dhw8rLy9PBw4c0MqVKzV37lzNmTNHU6ZMsdp89tlnysvL0/nnn6/169dr3LhxuuGGG7Ro0aJodgkNFAgEtHbtWq1du1Z+v9/ucgAAAAAADtQsWi+8d+9eXXPNNfrrX/+qP/3pT9b20tJSPfPMM5o3b54uuOACSdLs2bPl9Xq1evVqnXvuuXr77be1efNmvfPOO8rIyFDv3r11zz336Pbbb9fdd9+t5ORkzZo1S127dtVDDz0k6ci9oz/44APNmDFDubm50eoWGiAQCMjr9SoUClnbXC6X3G63jVUBAAAAAJwmaldS5OfnKy8vT4MHDw7bXlxcrIMHD4Zt79Gjhzwej1atWiVJWrVqlc466yxlZGRYbXJzc1VWVqZNmzZZbaq+dm5urvUa1SkvL1dZWVnYA9EXDAYVCoVUWFio4uJiFRcXy+/3y+Px2F0aAAAAAMBBonIlxQsvvKC1a9fqo48+OmZfSUmJkpOT1aZNm7DtGRkZKikpsdpUDiiO7j+6r7Y2ZWVl+uGHH5SWlnbMz54+fbqmTp3a4H7hxHi9XmVlZdldBgAAAADAoSJ+JcXOnTt1yy236LnnnlNqamqkX/6ETJo0SaWlpdZj586ddpcEAAAihMW6AQCIfREPKYqLi7V7925lZWWpWbNmatasmZYtW6bHHntMzZo1U0ZGhg4cOKDvv/8+7Pt27dqljh07SpI6dux4zADi6PPjtUlPT6/2KgpJSklJUXp6etgDAADEPhbrBgAgPkQ8pLjwwgu1YcMGrV+/3nqcffbZuuaaa6yvmzdvriVLlljfs2XLFgUCAeXk5EiScnJytGHDBu3evdtqs3jxYqWnp6tnz55Wm8qvcbTN0dcA0Hj8fr9195ZAIGB3OQASTOXFuk866SRr+9HFuh9++GFdcMEFys7O1uzZs7Vy5UqtXr1akqzFugsLC9W7d2/97Gc/0z333KOCggIdOHBAksIW6/Z6vRo7dqx+/etfa8aMGbb0FwCAeBbxkKJVq1Y688wzwx4tWrRQu3btdOaZZ6p169YaNWqUJkyYoPfee0/FxcUaOXKkcnJydO6550qSLr74YvXs2VPDhw/Xxx9/rEWLFmny5MnKz89XSkqKJOmmm27Sjh07dNttt+mTTz7RzJkz9eKLL2r8+PGR7hKAGrjdbrlcLvl8PmVnZys7O1ter5egAkCjYrFuAADiR9RuQVqbGTNmqEmTJho6dKjKy8uVm5urmTNnWvubNm2qN954Q2PGjFFOTo5atGihESNGaNq0aVabrl276s0339T48eP16KOP6pRTTtHf/vY3bj8KNCKPxyO/369gMCjpyBUVPp9PwWCQu7cAaBQs1g0AQHxplJBi6dKlYc9TU1NVUFCggoKCGr+nS5cueuutt2p93UGDBmndunWRKBFAA3k8HgIJALY4ulj34sWLHblY94QJE6znZWVlyszMtLEiAABiQ8SnewAAADQGFusGACD+EFIAAICYxGLdAGIRC44DtbNlTQoAAIATdXSx7soqL9YtyVqsu23btkpPT9fNN99c42Ld999/v0pKSqpdrPuJJ57Qbbfdpuuvv17vvvuuXnzxRb355puN22EkjEAgELbeE+JD5QXHj3K5XPL7/UydBSohpAAAAHGLxboRawKBgLxer0KhkLXN5XLJ7XbbWBUigQXHgbohpAAAAHGDxboR64LBoEKhkAoLC+X1eiUdOQPPH7HxgQXHgeMjpAAAAAAcxuv1Kisry+4yAKDRsXAmAAAAAABwBEIKAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAAR+DuHgAizu/3W19z2zQAAAAAdUVIASBi3G63XC6XfD6ftc3lcsnv9xNUAAAAADguQgoAEePxeOT3+xUMBiUduaLC5/MpGAwSUgAAAAA4LkIKABHl8XgIJAAAAAA0CAtnAgAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAILZyIqAoFA2B0eAAAAAByr8ljZ7XazADkSHiEFIi4QCMjr9SoUClnbXC6X3G63jVUBAAAAzuF2u+VyueTz+axtLpdLfr+foAIJjZACERcMBhUKhVRYWCiv1yuJVBgAAACozOPxyO/3h1197PP5FAwGGTcjoRFSIGq8Xq+ysrLsLgMAAABwJI/HQyABVMHCmQAAAAAAwBG4kgIAAACwEQuOA8B/EVIAAAAANmHBcQAIR0gBAAAA2IQFxwEgHCEFgKjj/t8AANSOBccB4AhCCgBRw/2/AQAAANQHIQWAqOH+3wAAAADqg5ACQFRx/28AAAAAddXE7gIAAAAAAAAkQgoAAAAAAOAQhBQAAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOQEgBAAAAAAAcgZACAAAAAAA4QjO7CwAAAAAAHOH3+62v3W63PB6PjdUAjY+QAgAAAABs5na75XK55PP5rG0ul0t+v5+gAgmFkAIAAAAAbObxeOT3+xUMBiUduaLC5/MpGAwSUiChEFIAAAAAgAN4PB4CCSQ8Fs4EAAAAAACOQEgBAAAAAAAcgekeAAAAQCMKBAJh6w4AAP6LkAJAo+PWWgCARBUIBOT1ehUKhaxtLpdLbrfbxqoAwDkIKRARnBFAXXBrLQBAogsGgwqFQiosLJTX65VEYA8AlRFS4IRxRgB1xa21AAA4wuv1Kisry+4yAMBxCClwwjgjgPrg1loAAAAAakJIgYjhjAAAAAAA4ERwC1IAAAAAAOAIhBQAAAAAAMARCCkAAAAAAIAjEFIAAAAAAABHIKQAAAAAAACOQEgBAAAAAAAcgZACAAAAAAA4AiEFAAAAAABwBEIKAAAAAADgCM3sLgAA/H6/9bXb7ZbH47GxGgAAAOdgnIREQ0gBwDZut1sul0s+n8/a5nK55Pf7OQADAICExjgJiYqQAoBtPB6P/H6/gsGgpCNnCnw+n4LBIAdfAACQ0BgnIVERUgCwlcfj4UALAIhrgUAg7A9NoK4YJyEREVIAAAAAURIIBOT1ehUKhaxtLpdLbrfbxqoAwLkIKQAAAIAoCQaDCoVCKiwslNfrlcTihwBQG0IKAAAAIMq8Xq+ysrLsLgMAHK+J3QUAAAAAAABIhBQAAAAAAMAhmO6Bequ8QrXEKtWIvMq/U8zbBQAAABIHIQXqpboVqiVWqUZkuN1uuVwu+Xw+a5vL5ZLf7yeoAAAAABIAIQXqpboVqiXOdiMyPB6P/H5/2L3kfT6fgsEgv18AAABAAiCkQIOwQjWixePxEEgAAAAACYqFMwEAAAAAgCNwJQUAAAAAxAgWGEe8I6QAAAAAAIdjgXEkCkIKAAAAAHA4FhhHoiCkAAAAACIoEAiE/SEJRAoLjCMRRHzhzOnTp+ucc85Rq1at1KFDBw0ZMkRbtmwJa7N//37l5+erXbt2atmypYYOHapdu3aFtQkEAsrLy5PL5VKHDh1066236tChQ2Ftli5dqqysLKWkpKhbt26aM2dOpLsDAAAA1FkgEJDX61V2drays7Pl8/nkcrnkdrvtLg0AYkLEQ4ply5YpPz9fq1ev1uLFi3Xw4EFdfPHF2rdvn9Vm/Pjxev311/XSSy9p2bJl+s9//qPLL7/c2n/48GHl5eXpwIEDWrlypebOnas5c+ZoypQpVpvPPvtMeXl5Ov/887V+/XqNGzdON9xwgxYtWhTpLgEAAAB1EgwGFQqFVFhYqOLiYhUXF7NmAADUQ8SnexQVFYU9nzNnjjp06KDi4mL99Kc/VWlpqZ555hnNmzdPF1xwgSRp9uzZ8nq9Wr16tc4991y9/fbb2rx5s9555x1lZGSod+/euueee3T77bfr7rvvVnJysmbNmqWuXbvqoYcekiR5vV598MEHmjFjhnJzcyPdLQAAAKDOvF6vsrKy7C4DAGJOxK+kqKq0tFSS1LZtW0lScXGxDh48qMGDB1ttevToIY/Ho1WrVkmSVq1apbPOOksZGRlWm9zcXJWVlWnTpk1Wm8qvcbTN0dcAAADxjSmmAADEn6iGFBUVFRo3bpz69++vM888U5JUUlKi5ORktWnTJqxtRkaGSkpKrDaVA4qj+4/uq61NWVmZfvjhh2rrKS8vV1lZWdgDgPP5/X6tXbtWa9euVSAQsLscAA7BFFMAAOJPVO/ukZ+fr40bN+qDDz6I5o+ps+nTp2vq1Kl2lwGgjrgfOIDaMMUUAID4E7UrKcaOHas33nhD7733nk455RRre8eOHXXgwAF9//33Ye137dqljh07Wm2qXop59Pnx2qSnpystLa3amiZNmqTS0lLrsXPnzhPqI4DoOno/8KMLjxUWFioUClm3dQOAypw0xZSrNwEAaJiIhxTGGI0dO1YLFy7Uu+++q65du4btz87OVvPmzbVkyRJr25YtWxQIBJSTkyNJysnJ0YYNG7R7926rzeLFi5Wenq6ePXtabSq/xtE2R1+jOikpKUpPTw97AHA2j8ejrKwsZWVlyev12l0OAIdy2hTT6dOnq3Xr1tYjMzPzhPsIAEAiiHhIkZ+fr8LCQs2bN0+tWrVSSUmJSkpKrIN469atNWrUKE2YMEHvvfeeiouLNXLkSOXk5Ojcc8+VJF188cXq2bOnhg8fro8//liLFi3S5MmTlZ+fr5SUFEnSTTfdpB07dui2227TJ598opkzZ+rFF1/U+PHjI90lAADgcEenmL7wwgt2lyKJqzcBAGioiK9J8eSTT0qSBg0aFLZ99uzZuu666yRJM2bMUJMmTTR06FCVl5crNzdXM2fOtNo2bdpUb7zxhsaMGaOcnBy1aNFCI0aM0LRp06w2Xbt21Ztvvqnx48fr0Ucf1SmnnKK//e1vzA0FACDBHJ1iunz58hqnmFa+mqLqFNM1a9aEvV4kppimpKRYJ1YAAEDdRTykMMYct01qaqoKCgpUUFBQY5suXbrorbfeqvV1Bg0apHXr1tW7RgAAEPuMMbr55pu1cOFCLV26tNYppkOHDpVU/RTTP//5z9q9e7c6dOggqfopplXHJMebYgoAABomqnf3AAAAiJb8/HzNmzdPr776qjXFVDoytTQtLS1simnbtm2Vnp6um2++ucYppvfff79KSkqqnWL6xBNP6LbbbtP111+vd999Vy+++KLefPNN2/oOAEC8IqQAAAAxiSmmACD5/f6w5263m1u1I6YRUgCIOZUPxhyIgcTFFFMAicztdsvlcsnn84Vtd7lc8vv9jI8QswgpcFyBQEDBYFDSsUkt0JiqOxhzIAYAAInI4/HI7/db43TpyFjd5/MpGAwyNkLMIqRArQKBgLxer0KhkLXN5XLJ7XbbWBUSVdWDMQdiAACQyDweD2MgxB1CCtQqGAwqFAqpsLBQXq9XEpfXw14cjAEAAID4RUiBOvF6vcrKyrK7DAAAAMdhaiwARA4hBQAAANBATI0FgMgipAAAAAAaiKmxABBZhBQAAADACWJqLABERhO7CwAAAAAAAJC4kgJAHKi6SBmX2QIAAACxiZACQMxyu91yuVzy+Xxh210ul/x+P0EFAABISJVP4HDyBrGGkAJAzPJ4PPL7/dZt36QjB2Wfz6dgMMgBGQAAJJTqTuBw8gaxhpACQEzzeDwcdAEAAHTsCRxO3iAWEVIAAAAAQJzgBA5iHXf3AAAAAAAAjkBIAQAAAAAAHIHpHgAAAEAdBQKBYxZsBgBEDiEFjlH54MuBF7GKW28BACItEAjI6/UqFAqFbXe5XHK73TZVBQDxhZACYao7+HLgRSzh1lsAgGgJBoMKhUIqLCyU1+u1thOGw+k4eYNYQkiBMNUdfPkgQyzh1lsAgGjzer3KysqyuwzguDh5g1hESIFqcfBFLOPWWwAAAJy8QWwipAAAAACAOMXJG8QabkEKAAAAAAAcgZACAAAAAAA4AiEFAAAAAABwBEIKAAAAAADgCCycCSAhcH9wAAAAwPkIKQDENe4PDgAAAMQOQgoAcY37gwMAAACxg5ACQNzj/uAAgIYKBAJhQTcQD5gGCycjpAAAAACqEQgE5PV6FQqFrG0ul0tut9vGqoCGYxosYgEhBThDgITEGQQAwPEEg0GFQiEVFhbK6/VK4piB2MY0WMQCQooExxkCJBrOIAAA6svr9SorK8vuMoCIYBosnI6QIsFxhgCJhjMIAAAAgHMRUkASZwiQWDiDAAAAADhTE7sLAAAAAAAAkAgpAAAAAACAQxBSAAAAAAAAR2BNCgAAAOD/cGt2ALAXIQUAKHwgyh1uACAxcWt2ALAfIQWAhOZ2u+VyueTz+axtLpdLfr+foAIAEgy3Zkei4mQNnISQAkBC83g88vv9YZf2+nw+BYNBDtAAkKC4NTsSBSdr4ESEFAASnsfj4UAMAAASDidr4ESEFABQDS57BAAAiYCTNXAaQgoAqITLHgEAAAD7EFIAQCVc9ggAABIdV5TCToQUAFAFlz0CAIBExBWlcAJCigQUCATCzhIDAAAAAFeUwgkIKRJMIBCQ1+tVKBSytrlcLrndbhurAgAAsAcnb4BwXFEKuxFSJJhgMKhQKKTCwkJ5vV5JzDMD6oK5mQAQfzh5AwDOQ0iRoLxer7KysuwuA3A85mYCQPzi5A0AOA8hBQDUgrmZABD/OHkDAM5BSAEAx8HcTAAAAKBxEFIAAAAgYbBQJlB/rM2FxkRIAQAAgITAQplA/bA2F+xASJEAOGMAAADAQplAfbE2F+xASBHnOGMARAeXPQJA7GKhTKDuWJsLjY2QIs5xxgCIrJoue1ywYIHat29vteH/GAAAAFB/hBQJgjMGQGRUvezx66+/1uWXX65LLrnEasNcTQBwDqa9AkBsIaQAgHqqetljdXM133//fa5eAgCbMe0ViA6mvSKaCCkA4ARVDi1YBRsAnINpr0BkMc5BYyCkAIAIqmkV7MpXVkgMkgGgMTHtFYgM7vaBxkBIAQARdrwrKyTOOgBAtLAGBRBd3O0D0UZIAQBRVPWMg8RZBwCIFtagAOzBGhWIJEKKOMQZBMBZOOMAAI2DNSiAxsUaFYgGQoo4wxkEAACQ6FiDAmgcdVmLi6AQ9UVIEWc4gwDEDi6NBAAAsY67nCHSCCniFGcQAOeq6QC+YMECtW/f3mrDwRwAAMQS7v6BSCCkAIBGVvUA/vXXX+vyyy/XJZdcYrUhtACAumEtLsBZWIsLJ4qQAgBsUPUAXpfQgkslASAca3EBsaFqgMjJF9SGkCIOcAYBiH21hRZcKgkA/1V13MNaXIBzVTfFVeLkC2pHSBHjOIMAxKfqLpVkoU0Aia6mcc///M//8JkIOFDVKa4SdwDB8RFSxDju5gHEPxbaBIAjGPcAsafqiRfuAILjIaSIE9zNA4hfLLQJIJFVN62VcQ8Qu7gDCI6HkCIGsQYFkHhYaBNAIqg8xpH++/nGtFYgvhxvWmtVnHxJLIQUMYY1KABIdVtok7meAGJJdWMc6cg4p6ioiCvFgDhV0+KalXHyJbHEfEhRUFCgBx54QCUlJerVq5cef/xx9e3b1+6yIopVrAEcT+XQoi5rWNQFny1AuEQYczS2441xJD6LgHhX3eKalXHyJfHEdEgxf/58TZgwQbNmzVK/fv30yCOPKDc3V1u2bFGHDh3sLi8iWMUaQH3VZQ2LumCdC+C/EmHM0RgqhxI1TeVgjAMknuqmfxzVkJMvjFliW5IxxthdREP169dP55xzjp544glJUkVFhTIzM3XzzTdr4sSJx/3+srIytW7dWqWlpUpPT492udWqOveyqqPJIVdOxLh9+6SWLY98vXev1KKFvfUg4Rzvs6aqmv54YEAQf5xwLIwF8TDmiIb6fLbU5XOFzxFYGDuhkuMFnFUxZnGmuh4LY/ZKigMHDqi4uFiTJk2ytjVp0kSDBw/WqlWrbKysdvX9DyZxVgHAiavtDEVN6ns1RkOmlNit6iCl6h9cDGIgxe6YIxJqCyHqOo6pjPUlADREbWtxVdVYYxY+v6InZkOKYDCow4cPKyMjI2x7RkaGPvnkk2q/p7y8XOXl5dbz0tJSSUcSnUgpKSlRSUlJtfuCwaB8Pp9++OEHa1taWppefvnlWhe+bNeundq0aRPROtHI9u3779dlZdLhw/bVAtRRmzZt1KZNG0lSt27dtGbNGn3zzTfVtj36+VbfKSV2S0tLU2Fhodxud42f0Uf3O0XHjh3VsWPHiL3e0WNLDF9YGXWxOOaIhOr+T1RVl3FMZe3atVNmZmbYNsY3qBZjJ9Si8hilqsYaszhxjBBpdo05YjakaIjp06dr6tSpx2yverBsTD/88IOGDh1q28+HDTp3trsCAP/neJ/BifQZvWfPHrVu3druMuKGE8cc0ZBI/0dgI8ZOcCA+/xrueGOOmA0p3G63mjZtql27doVt37VrV41pz6RJkzRhwgTreUVFhb799lu1a9dOSUlJUa03FpSVlSkzM1M7d+6Mq/myTsH7Gz28t9HF+xs9TnhvjTHas2ePOvNHQI2cOuZwwu9PY0vEPkv0O5H6nYh9luh3ovS7rmOOmA0pkpOTlZ2drSVLlmjIkCGSjgwAlixZorFjx1b7PSkpKUpJSQnbVtNlQoksPT09If6T2IX3N3p4b6OL9zd67H5vuYKidk4fc9j9+2OHROyzRL8TSSL2WaLfiaAuY46YDSkkacKECRoxYoTOPvts9e3bV4888oj27dunkSNH2l0aAACII4w5AABoHDEdUlxxxRX6+uuvNWXKFJWUlKh3794qKio6ZmErAACAE8GYAwCAxhHTIYUkjR07tsZLLVE/KSkpuuuuu465PBWRwfsbPby30cX7Gz28t7HFaWOORPz9ScQ+S/Q7kfqdiH2W6Hei9ft4kgz3HAMAAAAAAA7QxO4CAAAAAAAAJEIKAAAAAADgEIQUAAAAAADAEQgpAAAAAACAIxBSJJiCggKdeuqpSk1NVb9+/bRmzZoa286ZM0dJSUlhj9TU1EasNnYsX75cv/zlL9W5c2clJSXplVdeOe73LF26VFlZWUpJSVG3bt00Z86cqNcZq+r7/i5duvSY392kpCSVlJQ0TsExZPr06TrnnHPUqlUrdejQQUOGDNGWLVuO+30vvfSSevToodTUVJ111ll66623GqHa2NKQ95bPXVRVn+P2oEGDqv3sy8vLa8SKT1x9+ixJjzzyiLp37660tDRlZmZq/Pjx2r9/fyNVGzn16ffBgwc1bdo0nX766UpNTVWvXr1UVFTUiNWeuEQdO9W331999ZWuvvpq/fjHP1aTJk00bty4Rqkz0urb7wULFuiiiy5S+/btlZ6erpycHC1atKhxio2Q+vb5gw8+UP/+/dWuXTulpaWpR48emjFjRuMU6zCEFAlk/vz5mjBhgu666y6tXbtWvXr1Um5urnbv3l3j96Snp+urr76yHl988UUjVhw79u3bp169eqmgoKBO7T/77DPl5eXp/PPP1/r16zVu3DjdcMMNMffh21jq+/4etWXLlrDf3w4dOkSpwti1bNky5efna/Xq1Vq8eLEOHjyoiy++WPv27avxe1auXKmrrrpKo0aN0rp16zRkyBANGTJEGzdubMTKna8h763E5y7+q77H7QULFoT97mzcuFFNmzbVb37zm0auvOHq2+d58+Zp4sSJuuuuu+T3+/XMM89o/vz5+uMf/9jIlZ+Y+vZ78uTJeuqpp/T4449r8+bNuummm3TZZZdp3bp1jVx5wyXq2Km+/S4vL1f79u01efJk9erVK8rVRU99+718+XJddNFFeuutt1RcXKzzzz9fv/zlL+P6d7xFixYaO3asli9fLr/fr8mTJ2vy5Ml6+umno1ypAxkkjL59+5r8/Hzr+eHDh03nzp3N9OnTq20/e/Zs07p160aqLn5IMgsXLqy1zW233WbOOOOMsG1XXHGFyc3NjWJl8aEu7+97771nJJnvvvuuUWqKJ7t37zaSzLJly2psM2zYMJOXlxe2rV+/fua3v/1ttMuLaXV5b/ncRWX1PW5XNWPGDNOqVSuzd+/eaJUYcfXtc35+vrngggvCtk2YMMH0798/qnVGWn373alTJ/PEE0+Ebbv88svNNddcE9U6oyVRx0516XdlAwcONLfcckvU6mks9e33UT179jRTp06NfEGNoKF9vuyyy4zP54t8QQ7HlRQJ4sCBAyouLtbgwYOtbU2aNNHgwYO1atWqGr9v79696tKlizIzM3XppZdq06ZNjVFu3Fu1alXYv4Uk5ebm1vpvgfrr3bu3OnXqpIsuukgrVqywu5yYUFpaKklq27ZtjW34/W2Yury3Ep+7OKKhx+3KnnnmGV155ZVq0aJFtMqMqIb0+bzzzlNxcbE1NWLHjh1666239POf/7xRao6EhvS7vLz8mKlgaWlp+uCDD6Jaq5049iS2iooK7dmz57jH0Hiybt06rVy5UgMHDrS7lEZHSJEggsGgDh8+rIyMjLDtGRkZNc7T7969u5599lm9+uqrKiwsVEVFhc477zz9v//3/xqj5LhWUlJS7b9FWVmZfvjhB5uqih+dOnXSrFmz9PLLL+vll19WZmamBg0apLVr19pdmqNVVFRo3Lhx6t+/v84888wa29X0+8uaHzWr63vL5y6Oashxu7I1a9Zo48aNuuGGG6JVYsQ1pM9XX321pk2bpgEDBqh58+Y6/fTTNWjQoJia7tGQfufm5urhhx/Wp59+qoqKCi1evNia7hOvGDsltgcffFB79+7VsGHD7C4l6k455RSlpKTo7LPPVn5+fkx9jkdKM7sLgHPl5OQoJyfHen7eeefJ6/Xqqaee0j333GNjZUDtunfvru7du1vPzzvvPG3fvl0zZszQP/7xDxsrc7b8/Hxt3Lgxrs/E2aWu7y2fu4iUZ555RmeddZb69u1rdylRtXTpUt17772aOXOm+vXrp23btumWW27RPffcozvvvNPu8qLm0Ucf1Y033qgePXooKSlJp59+ukaOHKlnn33W7tKAiJs3b56mTp2qV199NSHWF3v//fe1d+9erV69WhMnTlS3bt101VVX2V1WoyKkSBBut1tNmzbVrl27wrbv2rVLHTt2rNNrNG/eXH369NG2bduiUWJC6dixY7X/Funp6UpLS7OpqvjWt29f/viuxdixY/XGG29o+fLlOuWUU2ptW9Pvb10/SxJNfd7bqvjcTVwnctzet2+fXnjhBU2bNi2aJUZcQ/p85513avjw4daZxrPOOkv79u3T6NGjdccdd6hJE+dfNNyQfrdv316vvPKK9u/fr2+++UadO3fWxIkTddpppzVGybZg7JSYXnjhBd1www166aWXjpnuE6+6du0q6cjn2a5du3T33XcnXEjh/E9uRERycrKys7O1ZMkSa1tFRYWWLFkSdtauNocPH9aGDRvUqVOnaJWZMHJycsL+LSRp8eLFdf63QP2tX7+e391qGGM0duxYLVy4UO+++651YKwNv79105D3tio+dxPXiRy3X3rpJZWXl8vn80W7zIhqSJ9DodAxQUTTpk0lHfk/GAtO5N86NTVVJ598sg4dOqSXX35Zl156abTLtQ3HnsTz/PPPa+TIkXr++edj7lbKkVJRUaHy8nK7y2h8Ni/ciUb0wgsvmJSUFDNnzhyzefNmM3r0aNOmTRtTUlJijDFm+PDhZuLEiVb7qVOnmkWLFpnt27eb4uJic+WVV5rU1FSzadMmu7rgWHv27DHr1q0z69atM5LMww8/bNatW2e++OILY4wxEydONMOHD7fa79ixw7hcLnPrrbcav99vCgoKTNOmTU1RUZFdXXC0+r6/M2bMMK+88or59NNPzYYNG8wtt9ximjRpYt555x27uuBYY8aMMa1btzZLly41X331lfUIhUJWm6qfDStWrDDNmjUzDz74oPH7/eauu+4yzZs3Nxs2bLCjC47VkPeWz11UVt/j9lEDBgwwV1xxRWOXGxH17fNdd91lWrVqZZ5//nmzY8cO8/bbb5vTTz/dDBs2zK4uNEh9+7169Wrz8ssvm+3bt5vly5ebCy64wHTt2jWm7mqVqGOn+vbbGGO1z87ONldffbVZt25dzB0X6tvv5557zjRr1swUFBSEHUO///57u7pQb/Xt8xNPPGFee+01s3XrVrN161bzt7/9zbRq1crccccddnXBNoQUCebxxx83Ho/HJCcnm759+5rVq1db+wYOHGhGjBhhPR83bpzVNiMjw/z85z83a9eutaFq5zt6y8uqj6Pv54gRI8zAgQOP+Z7evXub5ORkc9ppp5nZs2c3et2xor7v73333WdOP/10k5qaatq2bWsGDRpk3n33XXuKd7jq3ldJYb+PVT8bjDHmxRdfND/+8Y9NcnKyOeOMM8ybb77ZuIXHgIa8t3zuoqr6HLeNMeaTTz4xkszbb7/dyJVGTn36fPDgQXP33Xdbn/mZmZnmd7/7XUz9sX5Uffq9dOlS4/V6TUpKimnXrp0ZPny4+fLLL22ouuESdezUkH5X175Lly6NXvuJqG+/Bw4cWGv7WFDfPj/22GPmjDPOMC6Xy6Snp5s+ffqYmTNnmsOHD9vTARslGRMj18IBAAAAAIC4xpoUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgBRNWjQII0bN87uMmJOUlKSXnnlFbvLAAAgZjDmaBjGHHAaQgogyq677jolJSXppptuOmZffn6+kpKSdN111zV+YQ4zd+5cDRgwwO4yTggHeQCAnRhz1A1jDsDZCCmARpCZmakXXnhBP/zwg7Vt//79mjdvnjwej42V1c2BAwei/jNeffVV/epXv4r6zwEAIJ4x5jg+xhyAsxFSAI0gKytLmZmZWrBggbVtwYIF8ng86tOnT1jbiooKTZ8+XV27dlVaWpp69eql//3f/7X2Hz58WKNGjbL2d+/eXY8++mjYayxdulR9+/ZVixYt1KZNG/Xv319ffPGFpCNnWYYMGRLWfty4cRo0aJD1fNCgQRo7dqzGjRsnt9ut3NxcSdLGjRv1s5/9TC1btlRGRoaGDx+uYDBofd++fft07bXXqmXLlurUqZMeeuihOr0/+/fv19tvv13rgOHVV19VVlaWUlNTddppp2nq1Kk6dOiQJGnatGnq3LmzvvnmG6t9Xl6ezj//fFVUVEg6csbhySef1M9+9jOlpaXptNNOC3tfJWnnzp0aNmyY2rRpo7Zt2+rSSy/V559/Htbm2Wef1RlnnKGUlBR16tRJY8eOlSSdeuqpkqTLLrtMSUlJ1vPj1S5Jn376qX76058qNTVVPXv21OLFi+v0vgEAUBVjjtox5mDMAecjpAAayfXXX6/Zs2dbz5999lmNHDnymHbTp0/X3//+d82aNUubNm3S+PHj5fP5tGzZMklHBhSnnHKKXnrpJW3evFlTpkzRH//4R7344ouSpEOHDmnIkCEaOHCg/v3vf2vVqlUaPXq0kpKS6lXv3LlzlZycrBUrVmjWrFn6/vvvdcEFF6hPnz7617/+paKiIu3atUvDhg2zvufWW2/VsmXL9Oqrr+rtt9/W0qVLtXbt2uP+rCVLlujkk09Wjx49qt3//vvv69prr9Utt9yizZs366mnntKcOXP05z//WZJ0xx136NRTT9UNN9wgSSooKNDKlSs1d+5cNWny34+5O++8U0OHDtXHH3+sa665RldeeaX8fr8k6eDBg8rNzVWrVq30/vvva8WKFWrZsqUuueQS66zOk08+qfz8fI0ePVobNmzQa6+9pm7dukmSPvroI0nS7Nmz9dVXX1nPj1d7RUWFLr/8ciUnJ+vDDz/UrFmzdPvtt9f9HwoAgCoYc9SMMQdjDsQAAyCqRowYYS699FKze/duk5KSYj7//HPz+eefm9TUVPP111+bSy+91IwYMcIYY8z+/fuNy+UyK1euDHuNUaNGmauuuqrGn5Gfn2+GDh1qjDHmm2++MZLM0qVLa62nsltuucUMHDjQej5w4EDTp0+fsDb33HOPufjii8O27dy500gyW7ZsMXv27DHJycnmxRdftPZ/8803Ji0tzdxyyy011m6MMTfeeKP5wx/+UOP+Cy+80Nx7771h2/7xj3+YTp06Wc+3b99uWrVqZW6//XaTlpZmnnvuubD2ksxNN90Utq1fv35mzJgx1ut1797dVFRUWPvLy8tNWlqaWbRokTHGmM6dO5s77rijxjolmYULF9ar9kWLFplmzZqZL7/80tr/z3/+s9rXAgCgNow5GHMw5kA8aGZXOAIkmvbt2ysvL09z5syRMUZ5eXlyu91hbbZt26ZQKKSLLroobPuBAwfCLtEsKCjQs88+q0AgoB9++EEHDhxQ7969JUlt27bVddddp9zcXF100UUaPHiwhg0bpk6dOtWr3uzs7LDnH3/8sd577z21bNnymLbbt2+36ujXr5+1vW3bturevXutP8cYo9dff906K1Odjz/+WCtWrLDOBEhHLkHdv3+/QqGQXC6XTjvtND344IP67W9/qyuuuEJXX331Ma+Tk5NzzPP169dbP2Pbtm1q1apVWJv9+/dr+/bt2r17t/7zn//owgsvrLU/9a3d7/crMzNTnTt3rrFOAADqgzFH9RhzMOZAbCCkABrR9ddfb80nLCgoOGb/3r17JUlvvvmmTj755LB9KSkpkqQXXnhBf/jDH/TQQw8pJydHrVq10gMPPKAPP/zQajt79mz9/ve/V1FRkebPn6/Jkydr8eLFOvfcc9WkSRMZY8Je++DBg8fU0qJFi2Nq++Uvf6n77rvvmLadOnXStm3b6vIWHGPNmjU6dOiQzjvvvBrb7N27V1OnTtXll19+zL7U1FTr6+XLl6tp06b6/PPPdejQITVrVvePuL179yo7O1vPPffcMfvat28fdglnfdS1dgAAIokxx7EYcwCxgZACaERH5xomJSVZC0NV1rNnT6WkpCgQCGjgwIHVvsaKFSt03nnn6Xe/+521bfv27ce069Onj/r06aNJkyYpJydH8+bN07nnnqv27dtr48aNYW3Xr1+v5s2b11p7VlaWXn75ZZ166qnVHohPP/10NW/eXB9++KG1evh3332nrVu31tgX6cgCT3l5eWratGmtP3vLli3WXMzqzJ8/XwsWLNDSpUs1bNgw3XPPPZo6dWpYm9WrV+vaa68Ne370bFFWVpbmz5+vDh06KD09vdqfceqpp2rJkiU6//zzq93fvHlzHT58uF61e71e7dy5U1999ZV15mn16tU19hMAgLpgzHEsxhyMORAj7JxrAiSCqvMxS0tLTWlpqfW88vxQY4y54447TLt27cycOXPMtm3bTHFxsXnsscfMnDlzjDHGPProoyY9Pd0UFRWZLVu2mMmTJ5v09HTTq1cvY4wxO3bsMBMnTjQrV640n3/+uVm0aJFp166dmTlzpjHGmKKiIpOUlGTmzp1rtm7daqZMmWLS09OPmR9adU7nl19+adq3b29+/etfmzVr1pht27aZoqIic91115lDhw4ZY4y56aabTJcuXcySJUvMhg0bzK9+9SvTsmXLWueHnnHGGebll1+u9T0sKioyzZo1M3fffbfZuHGj2bx5s3n++eetuZo7d+40J510knnsscfC2q9atcp6DUnG7XabZ555xmzZssVMmTLFNGnSxGzatMkYY8y+ffvMj370IzNo0CCzfPlys2PHDvPee++Zm2++2ezcudMYY8ycOXNMamqqefTRR83WrVutf5ujfvSjH5kxY8aYr776ynz77bd1qv3w4cOmZ8+e5qKLLjLr1683y5cvN9nZ2cwPBQDUG2MOxhyMORAPCCmAKKtu0ajKqg4YKioqzCOPPGK6d+9umjdvbtq3b29yc3PNsmXLjDFHFrq67rrrTOvWrU2bNm3MmDFjzMSJE60BQ0lJiRkyZIjp1KmTSU5ONl26dDFTpkwxhw8ftn7GlClTTEZGhmndurUZP368GTt27HEHDMYYs3XrVnPZZZeZNm3amLS0NNOjRw8zbtw4a+GnPXv2GJ/PZ1wul8nIyDD3339/ja9ljDHbtm0zKSkpZu/evcd9H4uKisx5551n0tLSTHp6uunbt695+umnTUVFhbnwwgtNbm5u2AJUN998szn99NPNnj17jDFHBgwFBQXmoosuMikpKebUU0818+fPD/sZX331lbn22muN2+02KSkp5rTTTjM33nhj2ABv1qxZ1r9Np06dzM0332zte+2110y3bt1Ms2bNTJcuXY5b+1FbtmwxAwYMMMnJyebHP/6xKSoqYsAAAKg3xhyMORhzIB4kGVNlohgANJKHH35Y77zzjt56662o/6ykpCQtXLjwmPu1AwCA+MeYA4gdDVuVBQAi4JRTTtGkSZPsLgMAAMQ5xhxA7GDhTAC2GTZsmN0lAACABMCYA4gdTPcAAAAAAACOwHQPAAAAAADgCIQUAAAAAADAEQgpAAAAAACAIxBSAAAAAAAARyCkAAAAAAAAjkBIAQAAAAAAHIGQAgAAAAAAOAIhBQAAAAAAcARCCgAAAAAA4Aj/H6CSE5ABx/eKAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1280x1280 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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, `opticam` 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": 21,
   "id": "c591828b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.968552Z",
     "iopub.status.busy": "2026-02-16T18:29:24.968464Z",
     "iopub.status.idle": "2026-02-16T18:29:24.984782Z",
     "shell.execute_reply": "2026-02-16T18:29:24.984449Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999779237287508\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(calibrated_images / science_flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f52323de",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:24.985643Z",
     "iopub.status.busy": "2026-02-16T18:29:24.985567Z",
     "iopub.status.idle": "2026-02-16T18:29:25.001766Z",
     "shell.execute_reply": "2026-02-16T18:29:25.001296Z"
    },
    "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": 23,
   "id": "d75e9d7a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T18:29:25.002927Z",
     "iopub.status.busy": "2026-02-16T18:29:25.002850Z",
     "iopub.status.idle": "2026-02-16T18:29:25.164535Z",
     "shell.execute_reply": "2026-02-16T18:29:25.164208Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAHHCAYAAABtF1i4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPxlJREFUeJzt3Xl4U1Xi//FPKXQJUBBTymJDBQWDC1AULOgP0GpFRHCDUcOm4DLUL8jXBVSoqEPHBQYVBDdApzLigitMwal2UKgyA+IaQdYwSgsRZWmASnt+f/AlQ2hKm9L2Nu379Tx5HnJy7r0nx2v76b3nnhNhjDECAACwSAOrGwAAAOo3wggAALAUYQQAAFiKMAIAACxFGAEAAJYijAAAAEsRRgAAgKUIIwAAwFKEEQAAYCnCCFCP9O3bV3379vW/37p1qyIiIrRgwQJ/2ciRI9WkSZOab1wVW7BggSIiIrR169Zy62ZnZ6tr166KiYlRRESEfvvtN40cOVJJSUnV3k4AhBGgVtu0aZNuv/12tW/fXjExMYqLi1Pv3r319NNP68CBA1Y376T5fD49/PDDys3NtawNv/zyi4YMGaLY2FjNnj1bf/3rX9W4cWPL2gPURw2tbgCA4JYsWaIbbrhB0dHRGj58uM455xwVFRXps88+07333qvvvvtOL7zwwkkdo127djpw4IAaNWpURa0Ojc/n09SpUyUp4IpNTfrXv/6lffv26dFHH1VqaqolbQDqO8IIUAtt2bJFf/jDH9SuXTt9/PHHat26tf+zsWPHauPGjVqyZMlJHyciIkIxMTEnvZ+jDh8+rJKSEkVFRVXZPqvbzp07JUnNmze3tiFAPcZtGqAWeuKJJ7R//369/PLLAUHkqDPOOEPjxo3zv58/f74uueQStWzZUtHR0ercubPmzJlT7nGCjRk5avPmzUpLS1Pjxo3Vpk0bPfLIIzp2ke+j2z711FOaOXOmOnTooOjoaH3//fcqKirSlClT1L17dzVr1kyNGzfWxRdfrE8++SRg+/j4eEnS1KlTFRERoYiICD388MP+Oj/88IOuv/56tWjRQjExMTr//PP1/vvvl2rrd999p0suuUSxsbE67bTT9Nhjj6mkpKTc79+3b1+NGDFCknTBBRcoIiJCI0eODFo3NzdXERERpW4pHd+HO3fuVHx8vPr27RvQXxs3blTjxo01dOjQctsF1DdcGQFqoQ8++EDt27dXr169KlR/zpw5Ovvss3X11VerYcOG+uCDD/THP/5RJSUlGjt2bMjHLy4u1hVXXKELL7xQTzzxhLKzs5WRkaHDhw/rkUceCag7f/58HTx4ULfddpuio6PVokUL7d27Vy+99JJuvPFGjRkzRvv27dPLL7+stLQ0rV69Wl27dlV8fLzmzJmjO++8U9dcc42uvfZaSdJ5550n6UjA6N27t9q2bauJEyeqcePGeuONNzR48GC9/fbbuuaaayRJ+fn56tevnw4fPuyv98ILLyg2Nrbc7/nggw+qU6dOeuGFF/TII4/o9NNPV4cOHULur2O1bNlSc+bM0Q033KBnn31W//M//6OSkhKNHDlSTZs21XPPPXdS+wfqJAOgVtmzZ4+RZAYNGlThbXw+X6mytLQ00759+4CyPn36mD59+vjfb9myxUgy8+fP95eNGDHCSDJ33XWXv6ykpMQMGDDAREVFmV27dgVsGxcXZ3bu3BlwnMOHD5tDhw4FlP36668mISHB3HLLLf6yXbt2GUkmIyOjVPsvvfRSc+6555qDBw8GtKNXr17mzDPP9JeNHz/eSDJffPGFv2znzp2mWbNmRpLZsmVLqX0fa/78+UaS+de//hVQPmLECNOuXTv/+08++cRIMp988klAvWB9aIwxN954o7HZbGbDhg3mySefNJLMu+++e8K2APUVt2mAWmbv3r2SpKZNm1Z4m2OvAuzZs0der1d9+vTR5s2btWfPnkq1Iz093f/viIgIpaenq6ioSP/4xz8C6l133XX+2y1HRUZG+seNlJSUaPfu3Tp8+LDOP/98rV27ttxj7969Wx9//LGGDBmiffv2yev1yuv16pdfflFaWpp+/PFH/fTTT5KkpUuX6sILL1SPHj3828fHx+vmm2+u1PeuKrNmzVKzZs10/fXXa/LkyRo2bJgGDRpkaZuA2iqswsiKFSs0cOBAtWnTRhEREXr33XdD3ocxRk899ZQ6duyo6OhotW3bVn/605+qvrFAJcXFxUmS9u3bV+FtVq5cqdTUVDVu3FjNmzdXfHy8HnjgAUmqVBhp0KCB2rdvH1DWsWNHSSo1b8fpp58edB+vvPKKzjvvPMXExOjUU09VfHy8lixZUqH2bNy4UcYYTZ48WfHx8QGvjIwMSf8deLpt2zadeeaZpfbRqVOnco9TnVq0aKFnnnlGX3/9tZo1a6ZnnnnG0vYAtVlYjRkpLCxUly5ddMstt/jvL4dq3LhxWr58uZ566imde+652r17t3bv3l3FLQUqLy4uTm3atNG3335bofqbNm3SpZdeqrPOOkszZsxQYmKioqKitHTpUv3lL3+p0EDOkxFsbEZWVpZGjhypwYMH695771XLli0VGRmpzMxMbdq0qdx9Hm3zPffco7S0tKB1zjjjjJNreIgiIiKClhcXF5e5zbJlyyRJv/76q/7zn//wxA5QhrAKI/3791f//v3L/PzQoUN68MEH9be//U2//fabzjnnHD3++OP++QvcbrfmzJmjb7/91v9XU1l/1QFWuuqqq/TCCy8oLy9PKSkpJ6z7wQcf6NChQ3r//fflcDj85cc+uRKqkpISbd682X81RJI2bNggSRWalfStt95S+/bttXjx4oBf4kevahxV1i/4o1dlGjVqVO7cH+3atdOPP/5Yqnz9+vXltjMUp5xyiiTpt99+Cyjftm1b0PrZ2dl66aWXdN999+m1117TiBEj9MUXX6hhw7D6sQvUiLC6TVOe9PR05eXl6fXXX9fXX3+tG264QVdccYX/B9XRJxQ+/PBDnX766UpKStLo0aO5MoJa57777lPjxo01evRoFRQUlPp806ZNevrppyUdGZ8hKeAx0j179mj+/Pkn1YZZs2b5/22M0axZs9SoUSNdeuml5W4brE1ffPGF8vLyAurZbDZJpX/Bt2zZUn379tXzzz+vHTt2lNr/rl27/P++8sor9fnnn2v16tUBn7/22mvltjMU7dq1U2RkpFasWBFQHuzpmN9++02jR49Wjx49NG3aNL300ktau3atpk2bVqVtAuqKOhPRPR6P5s+fL4/HozZt2kg6cok3Oztb8+fP17Rp07R582Zt27ZNb775pl599VUVFxfr7rvv1vXXX6+PP/7Y4m8A/FeHDh20cOFCDR06VE6nM2AG1lWrVunNN9/0z4dx+eWXKyoqSgMHDtTtt9+u/fv368UXX1TLli2D/iKviJiYGGVnZ2vEiBHq2bOn/v73v2vJkiV64IEHSg1WDeaqq67S4sWLdc0112jAgAHasmWL5s6dq86dO2v//v3+erGxsercubMWLVqkjh07qkWLFjrnnHN0zjnnaPbs2brooot07rnnasyYMWrfvr0KCgqUl5en//znP/rqq68kHQluf/3rX3XFFVdo3Lhx/kd727Vrp6+//rpS3z+YZs2a+R/XjYiIUIcOHfThhx/6x64ca9y4cfrll1/0j3/8Q5GRkbriiis0evRoPfbYYxo0aJC6dOlSZe0C6gRLn+U5CZLMO++843//4YcfGkmmcePGAa+GDRuaIUOGGGOMGTNmjJFk1q9f799uzZo1RpL54YcfavorAOXasGGDGTNmjElKSjJRUVGmadOmpnfv3ubZZ58NeOT1/fffN+edd56JiYkxSUlJ5vHHHzfz5s0r9WhrRR/tbdy4sdm0aZO5/PLLjc1mMwkJCSYjI8MUFxeX2vbJJ58s1e6SkhIzbdo0065dOxMdHW26detmPvzww1KPyxpjzKpVq0z37t1NVFRUqcd8N23aZIYPH25atWplGjVqZNq2bWuuuuoq89ZbbwXs4+uvvzZ9+vQxMTExpm3btubRRx81L7/8cpU+2mvMkUeRr7vuOmOz2cwpp5xibr/9dvPtt98G9OF7771nJJnp06cHbLt3717Trl0706VLF1NUVHTCNgH1TYQxx1xHDSMRERF65513NHjwYEnSokWLdPPNN+u7777zXyI+qkmTJmrVqpUyMjI0bdo0/f777/7PDhw4IJvNpuXLl+uyyy6rya8AAABUh27TdOvWTcXFxdq5c6cuvvjioHV69+6tw4cPa9OmTf5ZFo8OymvXrl2NtRUAAPxXWF0Z2b9/vzZu3CjpSPiYMWOG+vXrpxYtWsjhcMjlcmnlypWaPn26unXrpl27diknJ0fnnXeeBgwYoJKSEl1wwQVq0qSJZs6c6Z8qOy4uTsuXL7f42wEAUD+FVRjJzc1Vv379SpWPGDFCCxYs0O+//67HHntMr776qn766SfZ7XZdeOGFmjp1qs4991xJ0s8//6y77rpLy5cvV+PGjdW/f39Nnz5dLVq0qOmvAwAAFGZhBAAA1D11ap4RAAAQfggjAADAUmHxNE1JSYl+/vlnNW3atMzpowEAQO1ijNG+ffvUpk0bNWhQ9vWPsAgjP//8sxITE61uBgAAqITt27frtNNOK/PzsAgjTZs2lXTkyxxdXh1AHVZYKP3fsg76+WepcWNr2wOgUvbu3avExET/7/GyhEUYOXprJi4ujjAC1AfHzqIcF0cYAcJceUMsGMAKAAAsRRgBAACWIowAAABLEUYAAIClCCMAAMBShBEAAGApwggAALAUYQQAAFiKMAIAACxFGAEAAJYijAAAAEsRRgAAgKUIIwAAwFKEEQAAYKmGVjcAADwej7xer/99gwMH1PX//r1u3TqVxMbKbrfL4XBY0j4A1YswAsBSHo9HTqdTPp/PX2aTVPh//+590UXySbLZbHK73QQSoA4ijACwlNfrlc/nU1ZWlpxOp6QjV0Z00UWSpJWffabvtm6Vy+XSp59+6q/DlRKg7iCMAKgVnE6nkpOTj7wpLPSXd+3aVS0SE2Wz2eRyufzlXCkB6g7CCIBaz+FwyO12+8eVuN1uuVwueb1ewghQBxBGAIQFh8NB8ADqKMIIgBp1/JMzbrfbwtYAqA0IIwBqTLAnZ6Qj4z/sdrtFrQJgNcIIgBoT7MkZiSdjgPqOMAKgxgU8OQOg3iOMAAhbx4834QoLEJ4IIwDCjt1uLzXviMTcI0C4IowACDvHzzsiMfcIEM4IIwDCEvOOAHVHA6sbAAAA6jfCCAAAsBS3aQBUq2NnXGW2VQDBEEYAVJtgM64y2yqA4xFGAFSbYDOuMhcIgOMRRgBUO2ZcBXAiDGAFAACW4soIgDrl2EGy3BICwgNhBECdEGyKeKaHB8IDYQRAnXD8FPFMDw+ED8IIgDqDKeKB8MQAVgAAYCnCCAAAsBRhBAAAWIoxIwCqzLHr0EisRQOgYggjAKpEsHVoJNaiAVA+wgiAKhFsHRqJiccAlI8wAqBKsQ4NgFAxgBUAAFiKMAIAACzFbRoAddrxT/QwhgWofQgjAOqkYAvnSSyeB9RGhBEAddLxC+dJLJ4H1FYhjxlZsWKFBg4cqDZt2igiIkLvvvtuudvk5uYqOTlZ0dHROuOMM7RgwYJKNBUAQuNwOJScnOx/HfvIMYDaI+QwUlhYqC5dumj27NkVqr9lyxYNGDBA/fr107p16zR+/HiNHj1ay5YtC7mxAACg7gn5Nk3//v3Vv3//CtefO3euTj/9dE2fPl3SkTkIPvvsM/3lL39RWlpaqIcHAAB1TLU/2puXl6fU1NSAsrS0NOXl5VX3oQEAQBio9gGs+fn5SkhICChLSEjQ3r17deDAAcXGxpba5tChQzp06JD//d69e6u7mQAAwCK1ctKzzMxMNWvWzP9KTEy0ukkAAKCaVHsYadWqlQoKCgLKCgoKFBcXF/SqiCRNmjRJe/bs8b+2b99e3c0EAAAWqfbbNCkpKVq6dGlA2UcffaSUlJQyt4mOjlZ0dHR1Nw3ASfJ4PP55PI6f6RQAKirkMLJ//35t3LjR/37Lli1at26dWrRoIYfDoUmTJumnn37Sq6++Kkm64447NGvWLN1333265ZZb9PHHH+uNN97QkiVLqu5bAKhxHo9HTqdTPp/PX2az2WS32y1sFYBwFHIY+fe//61+/fr530+YMEGSNGLECC1YsEA7duyQx+Pxf3766adryZIluvvuu/X000/rtNNO00svvcRjvUCY83q98vl8ysrK8k8mxrovACoj5DDSt29fGWPK/DzY7Kp9+/bVl19+GeqhAIQBp9Op5ORkq5sBIIzVyqdpAABA/UEYAQAAliKMAAAAS1X7o70AUNsc+xgyg24B6xFGANQbdrtdNptNLpfLX2az2eR2uwkkgIUIIwDqDYfDIbfbHTBRm8vlktfrJYwAFiKMAKhXHA4HwQOoZRjACgAALEUYAQAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYiknPAFSIx+Pxz1wqBa7vAgAngzACoFwej0dOp1M+ny+g3GazyW63W9QqAHUFYQRAubxer3w+n7KysuR0Ov3lrHgLoCoQRgBUmNPpVHJystXNAFDHMIAVAABYijACAAAsxW0aAPXe8U8GMRYGqFmEEQD1lt1ul81mk8vlCii32Wxyu90EEqCGEEYA1FsOh0Nut7vU/Ckul0ter5cwAtQQwgiAes3hcBA6AIsxgBUAAFiKMAIAACxFGAEAAJYijAAAAEsRRgAAgKUIIwAAwFKEEQAAYCnCCAAAsBRhBAAAWIowAgAALEUYAQAAlmJtGgBBeTwe/wJybrfb4tYAqMsIIwBK8Xg8cjqd8vl8/jKbzSa73W5hqwDUVYQRAKV4vV75fD5lZWXJ6XRKkux2O6vbAqgWhBEAZXI6nUpOTra6GQDqOAawAgAASxFGAACApQgjAADAUowZAYAgjn2cmcG7QPUijADAMex2u2w2m1wul7/MZrPJ7XYTSIBqQhgBgGM4HA653e6ACd9cLpe8Xi9hBKgmhBEAOI7D4SB4ADWIAawAAMBShBEAAGApwggAALAUYQQAAFiKMAIAACxVqTAye/ZsJSUlKSYmRj179tTq1atPWH/mzJnq1KmTYmNjlZiYqLvvvlsHDx6sVIMBAEDdEnIYWbRokSZMmKCMjAytXbtWXbp0UVpamnbu3Bm0/sKFCzVx4kRlZGTI7Xbr5Zdf1qJFi/TAAw+cdOMBAED4CzmMzJgxQ2PGjNGoUaPUuXNnzZ07VzabTfPmzQtaf9WqVerdu7duuukmJSUl6fLLL9eNN95Y7tUUAABQP4QURoqKirRmzRqlpqb+dwcNGig1NVV5eXlBt+nVq5fWrFnjDx+bN2/W0qVLdeWVV5Z5nEOHDmnv3r0BLwAAUDeFNAOr1+tVcXGxEhISAsoTEhL0ww8/BN3mpptuktfr1UUXXSRjjA4fPqw77rjjhLdpMjMzNXXq1FCaBgAAwlS1P02Tm5uradOm6bnnntPatWu1ePFiLVmyRI8++miZ20yaNEl79uzxv7Zv317dzQQAABYJ6cqI3W5XZGSkCgoKAsoLCgrUqlWroNtMnjxZw4YN0+jRoyVJ5557rgoLC3XbbbfpwQcfVIMGpfNQdHS0oqOjQ2kagJPg8Xj8C8NJRxaHA4CaElIYiYqKUvfu3ZWTk6PBgwdLkkpKSpSTk6P09PSg2/h8vlKBIzIyUpJkjKlEkwFUJY/HI6fTKZ/PF1Bus9lkt9stahWA+iTkVXsnTJigESNG6Pzzz1ePHj00c+ZMFRYWatSoUZKk4cOHq23btsrMzJQkDRw4UDNmzFC3bt3Us2dPbdy4UZMnT9bAgQP9oQSAdbxer3w+n7KysuR0Ov3ldrudlWsB1IiQw8jQoUO1a9cuTZkyRfn5+eratauys7P9g1o9Hk/AlZCHHnpIEREReuihh/TTTz8pPj5eAwcO1J/+9Keq+xYATprT6VRycrLVzQBQD4UcRiQpPT29zNsyubm5gQdo2FAZGRnKyMiozKEAAEAdx9o0AADAUoQRAABgqUrdpgGA+ub4x50Z4AtUHcIIAJyA3W6XzWaTy+UKKLfZbHK73QQSoAoQRgDgBBwOh9xud6lJ4Vwul7xeL2EEqAKEEQAoh8PhIHQA1YgBrAAAwFKEEQAAYCnCCAAAsBRhBAAAWIowAgAALEUYAQAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYioXygHrI4/H4V6F1u90WtwZAfUcYAeoZj8cjp9Mpn8/nL7PZbLLb7Ra2CkB9RhgB6hmv1yufz6esrCw5nU5Jkt1ul8PhsLhlAOorwghQTzmdTiUnJ1vdDABgACsAALAWYQQAAFiKMAIAACzFmBEAqKRjH4tmEDBQeYQRAAiR3W6XzWaTy+Xyl9lsNrndbgIJUAmEEQAIkcPhkNvtDpg4zuVyyev1EkaASiCMAEAlOBwOggdQRRjACgAALEUYAQAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYijACAAAsRRgBAACWIowAAABLEUYAAIClCCMAAMBSrNoL1HEej8e/1L10ZLl7AKhNCCNAHebxeOR0OuXz+QLKbTab7Ha7Ra0CgECEEaAO83q98vl8ysrKktPp9Jfb7XY5HA4LWwYA/0UYAeoBp9Op5ORkq5sBAEExgBUAAFiKMAIAACzFbRoAqCLHP6nE2BygYggjAHCS7Ha7bDabXC5XQLnNZpPb7SaQAOUgjADASXI4HHK73aXmc3G5XPJ6vYQRoByVGjMye/ZsJSUlKSYmRj179tTq1atPWP+3337T2LFj1bp1a0VHR6tjx45aunRppRoMALWRw+FQcnKy/3Xso9QATizkKyOLFi3ShAkTNHfuXPXs2VMzZ85UWlqa1q9fr5YtW5aqX1RUpMsuu0wtW7bUW2+9pbZt22rbtm1q3rx5VbQfAACEuZDDyIwZMzRmzBiNGjVKkjR37lwtWbJE8+bN08SJE0vVnzdvnnbv3q1Vq1apUaNGkqSkpKSTazUAAKgzQrpNU1RUpDVr1ig1NfW/O2jQQKmpqcrLywu6zfvvv6+UlBSNHTtWCQkJOuecczRt2jQVFxeXeZxDhw5p7969AS8AAFA3hRRGvF6viouLlZCQEFCekJCg/Pz8oNts3rxZb731loqLi7V06VJNnjxZ06dP12OPPVbmcTIzM9WsWTP/KzExMZRmAgCAMFLtk56VlJSoZcuWeuGFF9S9e3cNHTpUDz74oObOnVvmNpMmTdKePXv8r+3bt1d3MwEAgEVCGjNit9sVGRmpgoKCgPKCggK1atUq6DatW7dWo0aNFBkZ6S9zOp3Kz89XUVGRoqKiSm0THR2t6OjoUJoGAADCVEhXRqKiotS9e3fl5OT4y0pKSpSTk6OUlJSg2/Tu3VsbN25USUmJv2zDhg1q3bp10CACAADql5Bv00yYMEEvvviiXnnlFbndbt15550qLCz0P10zfPhwTZo0yV//zjvv1O7duzVu3Dht2LBBS5Ys0bRp0zR27Niq+xYAACBshfxo79ChQ7Vr1y5NmTJF+fn56tq1q7Kzs/2DWj0ejxo0+G/GSUxM1LJly3T33XfrvPPOU9u2bTVu3Djdf//9VfctAABA2KrUdPDp6elKT08P+llubm6pspSUFH3++eeVORQAAKjjqv1pGgAAgBMhjAAAAEsRRgAAgKUqNWYEQO3l8Xj8S9m73W6LWwMA5SOMAHWIx+OR0+mUz+fzl9lsNtntdgtbBQAnRhgB6hCv1yufz6esrCw5nU5JR2ZOdjgcFrcMAMpGGAHqIKfTqeTkZKubAQAVwgBWAABgKa6MAEA1OnYQMbfMgOAIIwBQDex2u2w2m1wul7/MZrPJ7XYTSIDjEEYAoBo4HA653e6Ax6xdLpe8Xi9hBDgOYQQAqonD4SB4ABXAAFYAAGApwggAALAUYQQAAFiKMAIAACxFGAEAAJYijAAAAEsRRgAAgKUIIwAAwFKEEQAAYCnCCAAAsBRhBAAAWIowAgAALEUYAQAAlmLVXiCMeTwe/xL10pFl6gEg3BBGgDDl8XjkdDrl8/kCym02m+x2u0WtAoDQEUaAMOX1euXz+ZSVlSWn0+kvt9vtcjgcFrYMAEJDGAHCnNPpVHJystXNAIBKI4wAQA06flwPV7IAwggA1Ai73S6bzSaXyxVQbrPZ5Ha7CSSo1wgjAFADHA6H3G53qaefXC6XvF4vYQT1GmEEAGqIw+EgdABBMOkZAACwFGEEAABYijACAAAsRRgBAACWIowAAABLEUYAAIClCCMAAMBShBEAAGApwggAALAUYQQAAFiKMAIAACxFGAEAAJYijAAAAEuxai8QRjwej38JerfbbXFrAKBqEEaAMOHxeOR0OuXz+fxlNptNdrvdwlYBwMkjjABhwuv1yufzKSsrS06nU5Jkt9vlcDgsbhkAnBzCCBBmnE6nkpOTrW4GAFQZwggAWOzY8T9c7UJ9VKmnaWbPnq2kpCTFxMSoZ8+eWr16dYW2e/311xUREaHBgwdX5rAAUKfY7XbZbDa5XC51795d3bt3l9PplMfjsbppQI0KOYwsWrRIEyZMUEZGhtauXasuXbooLS1NO3fuPOF2W7du1T333KOLL7640o0FgLrE4XDI7XZrzZo1WrNmjbKysuTz+fxPTAH1RchhZMaMGRozZoxGjRqlzp07a+7cubLZbJo3b16Z2xQXF+vmm2/W1KlT1b59+5NqMADUJQ6HQ8nJyUpOTvYPTAbqm5DCSFFRkdasWaPU1NT/7qBBA6WmpiovL6/M7R555BG1bNlSt956a4WOc+jQIe3duzfgBQAA6qaQwojX61VxcbESEhICyhMSEpSfnx90m88++0wvv/yyXnzxxQofJzMzU82aNfO/EhMTQ2kmAAAII9U6Hfy+ffs0bNgwvfjiiyFNzDRp0iTt2bPH/9q+fXs1thIAAFgppEd77Xa7IiMjVVBQEFBeUFCgVq1alaq/adMmbd26VQMHDvSXlZSUHDlww4Zav369OnToUGq76OhoRUdHh9I0AAAQpkK6MhIVFaXu3bsrJyfHX1ZSUqKcnBylpKSUqn/WWWfpm2++0bp16/yvq6++Wv369dO6deu4/QIAAEKf9GzChAkaMWKEzj//fPXo0UMzZ85UYWGhRo0aJUkaPny42rZtq8zMTMXExOicc84J2L558+aSVKocAADUTyGHkaFDh2rXrl2aMmWK8vPz1bVrV2VnZ/sHtXo8HjVoUK1DUQAAQB1Sqeng09PTlZ6eHvSz3NzcE267YMGCyhwSAADUUVzCAAAAliKMAAAASxFGAACApSo1ZgRA9fN4PAELph27zDwA1CWEEaAW8ng8cjqd8vl8AeU2my2k2YwBIBwQRoBayOv1yufzKSsrK2AlV7vdLofDYWHLAKDqEUaAWszpdCo5OdnqZgBAtSKMAEAtc/z4IK6Ioa4jjABALWG322Wz2eRyuQLKbTab3G43gQR1FmEEAGoJh8Mht9td6ikql8slr9dLGEGdRRgBgFrE4XAQOlDvMOkZAACwFGEEAABYijACAAAsRRgBAACWIowAAABLEUYAAIClCCMAAMBShBEAAGApwggAALAUYQQAAFiKMAIAACzF2jRALeHxePwLpB2/hDwA1GWEEaAW8Hg8cjqd8vl8/jKbzSa73W5hq1CbHBtQ7XY7i+mhTiGMALWA1+uVz+dTVlaWnE6nJH7h4Ai73S6bzSaXy+Uvs9lscrvdnB+oMwgjQC3idDqVnJxsdTNQizgcDrnd7oBbeC6XS16vlzCCOoMwAgC1nMPhIHigTuNpGgAAYCnCCAAAsBRhBAAAWIowAgAALEUYAQAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYirVpAAt4PB7/wmdS4PLwAFDfEEaAGubxeOR0OuXz+QLKbTab7Ha7Ra1CuDk+wNrtdhbTQ9gijAA1zOv1yufzKSsrS06n01/OLxNUhN1ul81mk8vlCii32Wxyu92cQwhLhBHAIk6nU8nJyVY3A2HG4XDI7XaXus3ncrnk9XoJIwhLhBEACDMOh4PQgTqFp2kAAIClCCMAAMBShBEAAGApwggAALAUYQQAAFiKMAIAACxFGAEAAJYijAAAAEtVKozMnj1bSUlJiomJUc+ePbV69eoy67744ou6+OKLdcopp+iUU05RamrqCesDAID6JeQwsmjRIk2YMEEZGRlau3atunTporS0NO3cuTNo/dzcXN1444365JNPlJeXp8TERF1++eX66aefTrrxQLjweDxau3at1q5dywq9qDZut9t/nnk8HqubA1RYhDHGhLJBz549dcEFF2jWrFmSpJKSEiUmJuquu+7SxIkTy92+uLhYp5xyimbNmqXhw4dX6Jh79+5Vs2bNtGfPHsXFxYXSXMBywVbpZVGzchQWSk2aHPn3/v1S48bWtqeW4xxDbVXR398hrU1TVFSkNWvWaNKkSf6yBg0aKDU1VXl5eRXah8/n0++//64WLVqEcmggbAVbpZcVelGVjl88j4XzEG5CCiNer1fFxcVKSEgIKE9ISNAPP/xQoX3cf//9atOmjVJTU8usc+jQIR06dMj/fu/evaE0E6iVWKUX1YnF8xDOavRpmj//+c96/fXX9c477ygmJqbMepmZmWrWrJn/lZiYWIOtBAAANSmkMGK32xUZGamCgoKA8oKCArVq1eqE2z711FP685//rOXLl+u88847Yd1JkyZpz549/tf27dtDaSYAAAgjIYWRqKgode/eXTk5Of6ykpIS5eTkKCUlpcztnnjiCT366KPKzs7W+eefX+5xoqOjFRcXF/ACAAB1U0hjRiRpwoQJGjFihM4//3z16NFDM2fOVGFhoUaNGiVJGj58uNq2bavMzExJ0uOPP64pU6Zo4cKFSkpKUn5+viSpSZMmanJ0tDwAAKi3Qg4jQ4cO1a5duzRlyhTl5+era9euys7O9g9q9Xg8atDgvxdc5syZo6KiIl1//fUB+8nIyNDDDz98cq0HAABhL+QwIknp6elKT08P+llubm7A+61bt1bmEAAAoJ5gbRoAAGApwggAALBUpW7TACibx+Pxz4QpibVoYJnjzz1m/kVtRRgBqlCwNUKkI+uE2O12i1qF+sZut8tms8nlcgWUs14NaivCCFCFgq1DI/EXKWrW8WvVSKxXg9qNMAJUA9ahgdVYqwbhhAGsAADAUoQRAABgKcIIAACwFGEEAABYijACAAAsRRgBAACW4tFe4CQdO+Mqs62itjv2HGX+G9QWhBHgJASbcZXZVlEbBZuVlRlZUVsQRoCTEGzGVf7aRG10/KyszMiK2oQwAlQBZlxFOGBWVtRWDGAFAACWIowAAABLEUYAAIClCCMAAMBShBEAAGApnqYBQnDsBGcSk5wh/B1/DvNoOqxAGAEqKNgEZxKTnCE8BZsETWIiNFiDMAJUULAJziT+kkR4On4SNImJ0GAdwggQIiY4Q13BJGioLRjACgAALEUYAQAAliKMAAAASzFmBAAQ4NjHfRmgjZpAGAFO4Nh5RZhTBHVdsMd9edQXNYEwApQh2LwizCmCuuz4x3151Bc1hTAClCHYvCJcskZdx+O+sAJhBCgH84oAQPXiaRoAAGAprowA/4dF8IDgWEwP1Y0wAohF8IBgWEwPNYUwAohF8IBgWEwPNYUwAhyDwapAIJ6uQU0gjAAAQsYsrahKhBHUW8yuCoSOWVpRHQgjqJeYXRWoHGZpRXUgjKBeYnZVoPIYR4KqRhhBvcaAVaBqMBcJTgZhBPUCE5oB1YO5SFAVCCOo85jQDKg+zEWCqkAYQZ3HhGZA9SprDAmP/6KiCCOok4I9tsv4EKBm8PgvQkUYQZ3DY7uAtXj8F6EijCDsBRucymO7gLWC3brhiRuUhTCCsHaiwakXX3wxP+iAWoAnblAewgjCSkWugkj8xQXUJid64ubTTz/lCiYqF0Zmz56tJ598Uvn5+erSpYueffZZ9ejRo8z6b775piZPnqytW7fqzDPP1OOPP64rr7yy0o1G/XFs+Ni1a5euvfZaroIAYej42zZlDXJdvHix4uPjA+rx/3bdF3IYWbRokSZMmKC5c+eqZ8+emjlzptLS0rR+/Xq1bNmyVP1Vq1bpxhtvVGZmpq666iotXLhQgwcP1tq1a3XOOedUyZdA3XD8VY9g4cNmsyk7O5sfVkCYO/5qydH/36+44oqAegSU+iHCGGNC2aBnz5664IILNGvWLElSSUmJEhMTddddd2nixIml6g8dOlSFhYX68MMP/WUXXnihunbtqrlz51bomHv37lWzZs20Z88excXFhdJc1GIVvepx7A8ifgjVE4WFUpMmR/69f7/UuLG17UGNqMgfJFLwgHI8flbUDhX9/R3SlZGioiKtWbNGkyZN8pc1aNBAqampysvLC7pNXl6eJkyYEFCWlpamd999N5RDoxY5/gdGZXDVA8DxynoCJ1hAOf4KyvEqElgqgp9BNSOkMOL1elVcXKyEhISA8oSEBP3www9Bt8nPzw9aPz8/v8zjHDp0SIcOHfK/37Nnj6QjCauq5efnn7AtCOT1euVyuXTgwIGT3ldsbKzefvtt//wfp556qhITE0vVq47/7qjlCgv/+++9e6XiYuvaAks1b95czZs3978/44wztHr1av3yyy9lbnP051R5gaUiYmNjlZWVVefnKWrVqpVatWpV5fs9+vO7vJswtfJpmszMTE2dOrVUebBfVAhfBw4c0HXXXWd1M1DbtWljdQtQj/Fzqmrs27dPzZo1K/PzkMKI3W5XZGSkCgoKAsoLCgrKTFStWrUKqb4kTZo0KeDWTklJiXbv3q1TTz1VERERoTQ5bOzdu1eJiYnavn17vR4XQz8cQT/QB0fRD0fQD+HZB8YY7du3T23K+aMipDASFRWl7t27KycnR4MHD5Z0JCjk5OQoPT096DYpKSnKycnR+PHj/WUfffSRUlJSyjxOdHS0oqOjA8qOvUxXl8XFxYXNSVad6Icj6Af64Cj64Qj6Ifz64ERXRI4K+TbNhAkTNGLECJ1//vnq0aOHZs6cqcLCQo0aNUqSNHz4cLVt21aZmZmSpHHjxqlPnz6aPn26BgwYoNdff13//ve/9cILL4R6aAAAUAeFHEaGDh2qXbt2acqUKcrPz1fXrl2VnZ3tH6Tq8XjUoEEDf/1evXpp4cKFeuihh/TAAw/ozDPP1LvvvsscIwAAQFIlB7Cmp6eXeVsmNze3VNkNN9ygG264oTKHqjeio6OVkZFR6vZUfUM/HEE/0AdH0Q9H0A91uw9CnvQMAACgKjUovwoAAED1IYwAAABLEUYAAIClCCMAAMBShJFqNHv2bCUlJSkmJkY9e/bU6tWrT1h/5syZ6tSpk2JjY5WYmKi7775bBw8e9H/+8MMPKyIiIuB11llnVffXOGmh9MPvv/+uRx55RB06dFBMTIy6dOmi7Ozsk9pnbVDVfRCO58KKFSs0cOBAtWnTRhERERVaLDM3N1fJycmKjo7WGWecoQULFpSqE07nQnX0QX04F3bs2KGbbrpJHTt2VIMGDQIm0TzWm2++qbPOOksxMTE699xztXTp0qpvfBWpjj5YsGBBqXMhJiamer5AFSOMVJNFixZpwoQJysjI0Nq1a9WlSxelpaVp586dQesvXLhQEydOVEZGhtxut15++WUtWrRIDzzwQEC9s88+Wzt27PC/Pvvss5r4OpUWaj889NBDev755/Xss8/q+++/1x133KFrrrlGX375ZaX3abXq6AMp/M6FwsJCdenSRbNnz65Q/S1btmjAgAHq16+f1q1bp/Hjx2v06NFatmyZv064nQvV0QdS3T8XDh06pPj4eD300EPq0qVL0DqrVq3SjTfeqFtvvVVffvmlBg8erMGDB+vbb7+tyqZXmeroA+nI7KzHngvbtm2rqiZXL4Nq0aNHDzN27Fj/++LiYtOmTRuTmZkZtP7YsWPNJZdcElA2YcIE07t3b//7jIwM06VLl2ppb3UJtR9at25tZs2aFVB27bXXmptvvrnS+7RadfRBOJ4Lx5Jk3nnnnRPWue+++8zZZ58dUDZ06FCTlpbmfx9u58KxqqoP6sO5cKw+ffqYcePGlSofMmSIGTBgQEBZz549ze23336SLax+VdUH8+fPN82aNauydtUkroxUg6KiIq1Zs0apqan+sgYNGig1NVV5eXlBt+nVq5fWrFnjv8S8efNmLV26VFdeeWVAvR9//FFt2rRR+/btdfPNN8vj8VTfFzlJlemHQ4cOlbqsGBsb6/9LrzL7tFJ19MFR4XQuVEZeXl5Av0lSWlqav9/C7VyojPL64Ki6fi5UREX7qq7bv3+/2rVrp8TERA0aNEjfffed1U2qEMJINfB6vSouLvZPkX9UQkKC8vPzg25z00036ZFHHtFFF12kRo0aqUOHDurbt2/AbZqePXtqwYIFys7O1pw5c7RlyxZdfPHF2rdvX7V+n8qqTD+kpaVpxowZ+vHHH1VSUqKPPvpIixcv1o4dOyq9TytVRx9I4XcuVEZ+fn7Qftu7d68OHDgQdudCZZTXB1L9OBcqoqy+qivnQkV06tRJ8+bN03vvvaesrCyVlJSoV69e+s9//mN108pFGKklcnNzNW3aND333HNau3atFi9erCVLlujRRx/11+nfv79uuOEGnXfeeUpLS9PSpUv122+/6Y033rCw5VXr6aef1plnnqmzzjpLUVFRSk9P16hRowLWO6rrKtIH9eFcQMVwLuColJQUDR8+XF27dlWfPn20ePFixcfH6/nnn7e6aeWqPz/ha5DdbldkZKQKCgoCygsKCtSqVaug20yePFnDhg3T6NGjde655+qaa67RtGnTlJmZqZKSkqDbNG/eXB07dtTGjRur/DtUhcr0Q3x8vN59910VFhZq27Zt+uGHH9SkSRO1b9++0vu0UnX0QTC1/VyojFatWgXtt7i4OMXGxobduVAZ5fVBMHXxXKiIsvqqrpwLldGoUSN169YtLM4Fwkg1iIqKUvfu3ZWTk+MvKykpUU5OjlJSUoJu4/P5Sv31HxkZKUkyZSwftH//fm3atEmtW7euopZXrcr0w1ExMTFq27atDh8+rLfffluDBg066X1aoTr6IJjafi5URkpKSkC/SdJHH33k77dwOxcqo7w+CKYungsVUZm+quuKi4v1zTffhMe5YPUI2rrq9ddfN9HR0WbBggXm+++/N7fddptp3ry5yc/PN8YYM2zYMDNx4kR//YyMDNO0aVPzt7/9zWzevNksX77cdOjQwQwZMsRf53//939Nbm6u2bJli1m5cqVJTU01drvd7Ny5s8a/X0WF2g+ff/65efvtt82mTZvMihUrzCWXXGJOP/108+uvv1Z4n7VNdfRBOJ4L+/btM19++aX58ssvjSQzY8YM8+WXX5pt27YZY4yZOHGiGTZsmL/+5s2bjc1mM/fee69xu91m9uzZJjIy0mRnZ/vrhNu5UB19UB/OBWOMv3737t3NTTfdZL788kvz3Xff+T9fuXKladiwoXnqqaeM2+02GRkZplGjRuabb76p0e9WUdXRB1OnTjXLli0zmzZtMmvWrDF/+MMfTExMTECd2oowUo2effZZ43A4TFRUlOnRo4f5/PPP/Z/16dPHjBgxwv/+999/Nw8//LDp0KGDiYmJMYmJieaPf/xjwC+goUOHmtatW5uoqCjTtm1bM3ToULNx48Ya/EaVE0o/5ObmGqfTaaKjo82pp55qhg0bZn766aeQ9lkbVXUfhOO58MknnxhJpV5Hv/uIESNMnz59Sm3TtWtXExUVZdq3b2/mz59far/hdC5URx/Ul3MhWP127doF1HnjjTdMx44dTVRUlDn77LPNkiVLauYLVUJ19MH48eP9/y8kJCSYK6+80qxdu7bmvtRJiDCmjHsAAAAANYAxIwAAwFKEEQAAYCnCCAAAsBRhBAAAWIowAgAALEUYAQAAliKMAAAASxFGAFS7vn37avz48VY3I+xERETo3XfftboZQLUjjAA1YOTIkYqIiNAdd9xR6rOxY8cqIiJCI0eOrPmG1TKvvPKKLrroIqubcVIIEEDoCCNADUlMTNTrr7+uAwcO+MsOHjyohQsXyuFwWNiyiikqKqr2Y7z33nu6+uqrq/04AGoXwghQQ5KTk5WYmKjFixf7yxYvXiyHw6Fu3boF1C0pKVFmZqZOP/10xcbGqkuXLnrrrbf8nxcXF+vWW2/1f96pUyc9/fTTAfvIzc1Vjx491LhxYzVv3ly9e/fWtm3bJB25UjN48OCA+uPHj1ffvn397/v27av09HSNHz9edrtdaWlpkqRvv/1W/fv3V5MmTZSQkKBhw4bJ6/X6tyssLNTw4cPVpEkTtW7dWtOnT69Q/xw8eFDLly8/YRh57733lJycrJiYGLVv315Tp07V4cOHJUmPPPKI2rRpo19++cVff8CAAerXr59KSkokHblqMWfOHPXv31+xsbFq3759QL9K0vbt2zVkyBA1b95cLVq00KBBg7R169aAOvPmzdPZZ5+t6OhotW7dWunp6ZKkpKQkSdI111yjiIgI//vy2i5JP/74o/7f//t/iomJUefOnfXRRx9VqN+AuoAwAtSgW265RfPnz/e/nzdvnkaNGlWqXmZmpl599VXNnTtX3333ne6++265XC7985//lHQkrJx22ml688039f3332vKlCl64IEH9MYbb0iSDh8+rMGDB6tPnz76+uuvlZeXp9tuu00REREhtfeVV15RVFSUVq5cqblz5+q3337TJZdcom7duunf//63srOzVVBQoCFDhvi3uffee/XPf/5T7733npYvX67c3FytXbu23GPl5OSobdu2Ouuss4J+/umnn2r48OEaN26cvv/+ez3//PNasGCB/vSnP0mSHnzwQSUlJWn06NGSpNmzZ2vVqlV65ZVX1KDBf3/UTZ48Wdddd52++uor3XzzzfrDH/4gt9stSfr999+Vlpampk2b6tNPP9XKlSvVpEkTXXHFFf4rQ3PmzNHYsWN122236ZtvvtH777+vM844Q5L0r3/9S5I0f/587dixw/++vLaXlJTo2muvVVRUlL744gvNnTtX999/f8X/QwHhzuqV+oD6YMSIEWbQoEFm586dJjo62mzdutVs3brVxMTEmF27dplBgwb5V+s8ePCgsdlsZtWqVQH7uPXWW82NN95Y5jHGjh1rrrvuOmOMMb/88ouRZHJzc0/YnmONGzcuYJXQPn36mG7dugXUefTRR83ll18eULZ9+3Yjyaxfv97s27fPREVFmTfeeMP/+S+//GJiY2PNuHHjymy7McaMGTPG3HPPPWV+fumll5pp06YFlP31r381rVu39r/ftGmTadq0qbn//vtNbGysee211wLqSzJ33HFHQFnPnj3NnXfe6d9fp06dTElJif/zQ4cOmdjYWLNs2TJjjDFt2rQxDz74YJntlGTeeeedkNq+bNky07Bhw4DVmf/+978H3RdQFzW0NAkB9Ux8fLwGDBigBQsWyBijAQMGyG63B9TZuHGjfD6fLrvssoDyoqKigNs5s2fP1rx58+TxeHTgwAEVFRWpa9eukqQWLVpo5MiRSktL02WXXabU1FQNGTJErVu3Dqm93bt3D3j/1Vdf6ZNPPlGTJk1K1d20aZO/HT179vSXt2jRQp06dTrhcYwx+uCDD/xXdoL56quvtHLlSv/VBOnI7aqDBw/K5/PJZrOpffv2euqpp3T77bdr6NChuummm0rtJyUlpdT7devW+Y+xceNGNW3aNKDOwYMHtWnTJu3cuVM///yzLr300hN+n1Db7na7lZiYqDZt2pTZTqAuI4wANeyWW27xjzGYPXt2qc/3798vSVqyZInatm0b8Fl0dLQk6fXXX9c999yj6dOnKyUlRU2bNtWTTz6pL774wl93/vz5+p//+R9lZ2dr0aJFeuihh/TRRx/pwgsvVIMGDWSMCdj377//XqotjRs3LtW2gQMH6vHHHy9Vt3Xr1tq4cWNFuqCU1atX6/Dhw+rVq1eZdfbv36+pU6fq2muvLfVZTEyM/98rVqxQZGSktm7dqsOHD6thw4r/mNu/f7+6d++u1157rdRn8fHxAbd7QlHRtgP1FWEEqGFHxx9ERET4B4Ueq3PnzoqOjpbH41GfPn2C7mPlypXq1auX/vjHP/rLNm3aVKpet27d1K1bN02aNEkpKSlauHChLrzwQsXHx+vbb78NqLtu3To1atTohG1PTk7W22+/raSkpKC/5Dt06KBGjRrpiy++8D8h9Ouvv2rDhg1lfhfpyODOAQMGKDIy8oTHXr9+vX98RjCLFi3S4sWLlZubqyFDhujRRx/V1KlTA+p8/vnnGj58eMD7o1eckpOTtWjRIrVs2VJxcXFBj5GUlKScnBz169cv6OeNGjVScXFxSG13Op3avn27duzY4b969fnnn5f5PYE6x+LbREC9cPwYjT179pg9e/b43x87ZsQYYx588EFz6qmnmgULFpiNGzeaNWvWmGeeecYsWLDAGGPM008/beLi4kx2drZZv369eeihh0xcXJzp0qWLMcaYzZs3m4kTJ5pVq1aZrVu3mmXLlplTTz3VPPfcc8YYY7Kzs01ERIR55ZVXzIYNG8yUKVNMXFxcqTEjx4/z+Omnn0x8fLy5/vrrzerVq83GjRtNdna2GTlypDl8+LAxxpg77rjDtGvXzuTk5JhvvvnGXH311aZJkyYnHDNy9tlnm7fffvuEfZidnW0aNmxoHn74YfPtt9+a77//3vztb3/zj9/Yvn27OeWUU8wzzzwTUD8vL8+/D0nGbrebl19+2axfv95MmTLFNGjQwHz33XfGGGMKCwvNmWeeafr27WtWrFhhNm/ebD755BNz1113me3btxtjjFmwYIGJiYkxTz/9tNmwYYP/v81RZ555prnzzjvNjh07zO7duyvU9uLiYtO5c2dz2WWXmXXr1pkVK1aY7t27M2YE9QZhBKgBwQaMHuv4MFJSUmJmzpxpOnXqZBo1amTi4+NNWlqa+ec//2mMOTLIdeTIkaZZs2amefPm5s477zQTJ070h5H8/HwzePBg07p1axMVFWXatWtnpkyZYoqLi/3HmDJliklISDDNmjUzd999t0lPTy83jBhjzIYNG8w111xjmjdvbmJjY81ZZ51lxo8f7x/0uW/fPuNyuYzNZjMJCQnmiSeeKHNfxhizceNGEx0dbfbv319uP2ZnZ5tevXqZ2NhYExcXZ3r06GFeeOEFU1JSYi699FKTlpYWMPj0rrvuMh06dDD79u0zxhwJI7NnzzaXXXaZiY6ONklJSWbRokUBx9ixY4cZPny4sdvtJjo62rRv396MGTMmIDzOnTvX/9+mdevW5q677vJ/9v7775szzjjDNGzY0LRr167cth+1fv16c9FFF5moqCjTsWNHk52dTRhBvRFhzHE3jgGgBs2YMUP/+Mc/tHTp0mo/VkREhN55551Sc6wAsBbzjACw1GmnnaZJkyZZ3QwAFmIAKwBLHTthGoD6iTACoN7grjRQO3GbBgAAWIowAgAALEUYAQAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABY6v8D6DkpKAhgTSMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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 `opticam`'s calibration error propagation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opticam6",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}