pandora-moonlight.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f6e4167",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8499ccde",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandorasat as ps\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import warnings\n",
    "import astropy.units as u\n",
    "plt.style.use(ps.pandorastyle)\n",
    "import os\n",
    "phoenixpath = f\"/Users/chedges/repos/pandora-psf/src/pandorapsf/data/phoenix\"\n",
    "os.environ[\"PYSYN_CDBS\"] = phoenixpath\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.filterwarnings(\"ignore\", message=\"Extinction files not found in \")\n",
    "    # Third-party\n",
    "    import pysynphot\n",
    "\n",
    "def get_phoenix_model(teff, logg, vmag):\n",
    "    logg1 = logg.value if isinstance(logg, u.Quantity) else logg\n",
    "    star = pysynphot.Icat(\n",
    "        \"phoenix\",\n",
    "        teff.value if isinstance(teff, u.Quantity) else teff,\n",
    "        0,\n",
    "        logg1 if np.isfinite(logg1) else 5,\n",
    "    )\n",
    "    star_norm = star.renorm(\n",
    "        vmag, \"vegamag\", pysynphot.ObsBandpass(\"johnson,V\")\n",
    "    )\n",
    "    star_norm.convert(\"Micron\")\n",
    "    star_norm.convert(\"flam\")\n",
    "\n",
    "    mask = (star_norm.wave >= 0.1) * (star_norm.wave <= 3)\n",
    "    wavelength = star_norm.wave[mask] * u.micron\n",
    "    wavelength = wavelength.to(u.angstrom)\n",
    "\n",
    "    sed = star_norm.flux[mask] * u.erg / u.s / u.cm**2 / u.angstrom\n",
    "    return wavelength, sed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c0b1326e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solar spectrum renormalized to moon visible magnitude\n",
    "wavelength, moon_spectrum = get_phoenix_model(5777, 4.43, -12.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "15f16523",
   "metadata": {},
   "outputs": [],
   "source": [
    "VISDA = ps.visibledetector.VisibleDetector()\n",
    "NIRDA = ps.irdetector.NIRDetector()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5c090016",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x136e607c0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7EAAALBCAYAAACOdiYsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAABcSAAAXEgFnn9JSAAB2tklEQVR4nO3deXhU5d3/8c8kk20CCUsQElllEUVBwAVExLoA1kQEt8elgFtdW/uzdjOiVI21Kra2auvjY9VqXVqrKC64oCC7CsoOLiBb2AIBQlYmOb8/jhNAJvuZc+6Zeb+uK9chc+ac8wXPNeaT733fx2dZliUAAAAAAKJAgtcFAAAAAADQWIRYAAAAAEDUIMQCAAAAAKIGIRYAAAAAEDUIsQAAAACAqEGIBQAAAABEDUIsAAAAACBqEGIBAAAAAFGDEAsAAAAAiBqEWAAAAABA1CDEAgAAAACiBiEWAAAAABA1CLEAAAAAgKhBiAUAAAAARA1CrEM+//xzjR8/Xr169ZLP59Odd95Z7/uXLVsmv9+vzp07u1QhAAAAAEQ/QqxD5s6dqwULFui0005TZmZmg+//xS9+ofbt27tQGQAAAADEDkKsQ372s5/pq6++0rPPPqs2bdrU+96pU6dq7dq1uvrqq90pDgAAAABiBCHWIQkJjfunrKqq0u23364HHnhAKSkpEa4KAAAAAGJLzIfYRYsW6YEHHtC4cePUuXNn+Xw++Xy+Bo8rLy/XXXfdpT59+ig1NVU5OTm6+uqrtXnz5hbV8+c//1kdOnTQpZde2qLzAAAAAEA88ntdQKTde++9euONN5p0TEVFhc4880wtWLBA2dnZGjNmjL777js988wzeuutt7RgwQIdddRRTa5l27ZtKigo0PTp05t8LAAAAAAgDjqxQ4cO1aRJk/Tmm29qy5YtjRrCe99992nBggUaOnSovvrqK73yyitauHChpkyZoh07djR7Lusdd9yh0aNHa+jQoc06HgAAAADinfGd2Llz52rLli266KKL6nxPSUmJHnnkEd15551KTEw8ZN9vfvObJl2vqqpKjz32mCTp8ccfV6tWrWr33XbbbXruuec0a9YsLVq0SIMHD270eZcvX64XXnhBCxYs0O7duyXZHV/LsrR7924FAgElJyc3qVYAAAAAiDdGh9j9+/fr8ssvV2Fhofx+vy644ILD3lNaWqof//jHmjNnjnJycnTddde16Jpz587Vnj171LNnTw0cOPCw/RdddJGWLl2qadOmNSnEfvPNN6qqqtKgQYMO29e2bVv97W9/0w033NCi2gEAAAAg1hkdYpOSkvTyyy9r5MiRuvTSS/Xaa6/pvPPOq91fXl6u3NxczZkzR5deeqkjj6xZsmSJJIUNmwe/vnTp0iad97TTTtPHH398yGvPPvus3n77bf3nP/9Rnz59mlEtAAAAAMQXo0OsZM9pffvtt3Xuuefqwgsv1BtvvKFRo0apsrJSY8aM0cyZMzV27Fi98MILhw0lbo4NGzZIkjp37hx2f+j19evXH/L6jh07NGvWLElSWVmZVq9erVdffVXp6ek699xzlZWVpTPOOOOQY2bOnKmUlJTDXgcAAAAAhGd8iJWk008/XW+++aZyc3N1wQUX6L///a8ef/xxffDBBzrvvPP0yiuvyO935q+yb98+SVIgEAi7Pz09XZI9D/dgK1as0MUXX1z7/X//+1/997//Vbdu3fTdd985UhsAAAAAxLuoWZ34rLPO0muvvaaamhqdd955eueddzRy5Ej997//VVJSktfl6YwzzpBlWYd91RdgJ0+erE2bNrlXJAAAAABEuagJsZI0cuTIQx5P86tf/apRj8xpitBqxGVlZWH3l5aWSpJat27t6HUBAAAAAA2LmhBbU1Oj8ePHa9asWerfv7/8fr8uvPBCffbZZ45ep2vXrpJUZ4c09Hq3bt0cvS4AAAAAoGFREWIty9I111yjF198UcOGDdPcuXP1r3/9S6WlpRo1apS++OILx641YMAASdLixYvD7g+93r9/f8euCQAAAABoHONDrGVZuuGGG/Tss8/q5JNP1jvvvKNWrVrpkksu0bPPPqs9e/bonHPO0bJlyxy53rBhw5SZmalvv/1WX3755WH7X331VUlSXl6eI9cDAAAAADSe8SH25z//uf73f/9XgwYN0nvvvaeMjIzafVdeeaWeeuop7dq1S2effbZWrVrV4uslJyfrlltukSTdfPPNtXNgJemRRx7R0qVLNWLECA0ePLjF1wIAAAAANI3Rj9gJBoP65ptvNGDAAH3wwQdq06bNYe+5+uqrVVVVpdtuu01btmzRMcccc8j+t99+W/fee2/t91VVVZKkIUOG1L42adIknXfeebXf33nnnfrwww81b9489e7dW8OHD9f69eu1cOFCdejQQf/4xz8c/psCAAAAABrDZ1mW5XUR9SkvL1dZWZnat29f7/s2btyoLl26HPb6s88+q6uuuqreY5955hlNnDjxsOv+4Q9/0IsvvqiNGzeqXbt2Gj16tO6991517ty5yX8PAAAAAEDLGR9iAQAAAAAIMX5OLAAAAAAAIYRYAAAAAEDUIMQCAAAAAKIGIRYAAAAAEDUIsQAAAACAqGFUiL3iiis0ePBgXXHFFV6XAgAAAAAwkN/rAg62evVqLV682OsyAAAAAACGMqoTCwAAAABAfQixAAAAAICoQYgFAAAAAEQNQiwAAAAAIGoQYgEAAAAAUYMQCwAAAACIGoRYAAAAAEDUIMQCAAAAAKIGIRYAAAAAEDUIsQAAAACAqEGIBQAAAABEDUIsAAAAACBq+L0uIJxgMKiioqKw+wKBgAKBgMsVAQAAAABMYGSILSwsVH5+fth9ubm5ysvLc7kiAAAAAIAJjAyxOTk5KigoCLuPLiwAAAAAxC8jQ6zf71dWVpbXZQAAAAAADMPCTgAAAACAqEGIBQAAAABEDUIsAAAAACBqEGIBAAAAAFGDEAsAAAAAiBqEWAAAAABA1CDEAgAAAACiBiEWAAAAABA1CLEAAAAAgKhBiAUAAAAARA1CLAAACK+mRioq8roKAAAOQYgFAADhXX+91KGDtGWL15UAAFDL73UB4QSDQRXV8ZvfQCCgQCDgckUAAMSh//s/e7tli5Sd7W0tAAB8z8gQW1hYqPz8/LD7cnNzlZeX53JFAADEmV27Dvy5qsq7OgAA+AEjQ2xOTo4KCgrC7qMLCwCACz766MCfy8u9qwMAgB8wMsT6/X5lZWV5XQYAAPHr/fcP/JkQCwAwCAs7AQCAQ9XUSO+8c+B7QiwAwCCEWAAAcKjPPpM2b5b69LG/J8QCAAxCiAUAAId6/XV7e8UV9pYQCwAwCCEWAAAcYFnSa69JOTnSiBH2a4RYAIBBHAuxM2fOlM/na/DrnnvuceqSAADAaatWSV9/LV1wgRR6IgAhFgBgEMdWJ+7UqZMmTJgQdl91dbVeeOEFSdLw4cOduiQAAHDatGn29vzzpbQ0+8+EWACAQRwLsX379tWzzz4bdt+7776rF154QV26dNEZZ5zh1CUBAIDTpk2TWrWSzjhD2rTJfo0QCwAwiCtzYkNd2CuuuEI+n8+NSwIAgKbasUOaN08aNUpKSaETCwAwUsRDbGlpqd544w1J0k9+8pNIXw4AADTXO+/YCzvl5dnfE2IBAAaKeIh97bXXVFpaqoEDB+rYY4+N9OUAAEBzvfmmlJAg/fjH9veEWACAgSIeYkNDienCAgBgsMpK6f33paFDpQ4d7NdSUiSfjxALADCKYws7hbNlyxbNmDFDiYmJuuyyy5p03OTJk5t93ZYcCwBAXJo5U9q3z16VOMTnk1JTCbEAAKNEtBP70ksvqbq6Wuecc446deoUyUsBAICWeOstexuaDxuSlkaIBQAYJaKd2OYOJc7OzqabCgCAmz74QOraVerb99DXCbEAAMNErBO7atUqffHFF2rVqpUuuOCCSF0GAAC01KZN0po10lln2UOID+b3S8GgN3UBABBGxELs888/L0kaN26cAoFApC4DAABaasYMe3v22YfvS0iwH7sDAIAhIhJiLcvSiy++KIlViQEAMN6HH9rbM888fJ/PJ9XUuFsPAAD1iEiInT17ttavX68jjzxSZ4b7HyIAADCDZdmd2OOOk8ItwkgnFgBgmIiE2NCCTpdffrkSEiL+KFoAANBcq1ZJW7bY82HDoRMLADCM4wmzsrJSr776qiTpyiuvdPr0AADASfXNh5XoxAIAjOP4I3ZSUlK0a9cup08LAAAi4cMPpcRE6fTTw++nEwsAMAxjfQEAiFfBoDRzpnTKKVJGRvj30IkFABiGEAsAQLz6/HNp796658NKdGIBAMYhxAIAEK8amg8r0YkFABiHEAsAQLz68EMpEJCGDKn7PXRiAQCGIcQCABCPSkuluXOlESOk5OS630cnFgBgGEIsAADxaNYsaf9+6Zxz6n8fnVgAgGEcf8SOE4LBoIqKisLuCwQCCgQCLlcEAECM+eADe9tQiKUTCwAwjJEhtrCwUPn5+WH35ebmKi8vz+WKAACIMR98IGVnS/361f8+OrEAAMMYGWJzcnJUUFAQdh9dWAAAWmjzZmnFCuknP7FDan3oxAIADGNkiPX7/crKyvK6DAAAYtN779nbUaMafi+dWACAYVjYCQCAePPee3Y4HTmy4ffSiQUAGIYQCwBAPKmutufDDhokdejQ8PvpxAIADEOIBQAgnnz+uVRc3LihxBKdWACAcQixAADEk6bMh5XoxAIAjEOIBQAgnrz3ntS6tTR0aOPeTycWAGAYQiwAAPGiuFhasEA680wpKalxx9CJBQAYhhALAEC8mDHDDqSNHUos0YkFABiHEAsAQLxo6nxYiU4sAMA4hFgAAFrKsuxH15jMsuwQ27u3dNRRjT+OTiwAwDCEWAAAmqusTPrDH+xQ2Lq1NHasNGuW11WFt2KFtHFj07qwEp1YAIBxCLEAADTH/v3SBRdId9whVVRIgwZJb74p/ehH0j/+4XV1h3v1VXs7dmzTjqMTCwAwDCEWAIDm+MUvpA8+kK65RtqwQZozR1q1Sure3X7tqae8rvBQr74qZWVJp5/etOPoxAIADOP3uoBwgsGgioqKwu4LBAIKBAIuVwQAwEH+/nfpiSfsruvf/nbgcTV9+tjDiX/0I+mnP7Xnyd5wg7e1Sna4XrFCuu46yd/E//XTiQUAGMbIEFtYWKj8/Pyw+3Jzc5WXl+dyRQAAfG/WLOlnP7Pnwf7nP4c/b7VLF2nmTPtZrDfeKCUnS1df7UmptUJDiS++uOnH0okFABjGZ1nm/Hp18ODBWrx4sfr3768ZM2aEfQ+dWACAZzZvlgYMkCorpQULpH796n5vYaF06qlSUZG0cqXUtat7df7QgAHSpk3S1q2Hh+6GjBkjvfuuVFUVmdoAAGgiIzuxfr9fWVlZXpcBAMChHn9c2rlTeuWV+gOsJOXkSE8+KY0eLf3859LUqa6UeJhly6SlS+15uk0NsBKdWACAcVjYCQCAxqiulp57TurWTbroosYdM2qUdOml0htv2F9eeO45eztxYvOOZ04sAMAwhFgAABrj/fftIcITJ9rBrrH+9CcpI8OeR1tWFrHywtq/X3rhBalnT2nYsOadg04sAMAwhFgAABoj9OzXCROadlx2tnTvvdLGjXagddN770nbttk1+3zNO0cosNONBQAYghALAEBDioulN9+0H53To0fTj7/hBql3b+mBB+zFldwSGko8fnzzzxEKv4RYAIAhCLEAADRk6lR7dd4rrmje8cnJ0oMPSvv2SXff7Whpddq50w7eZ55pz+NtLjqxAADDEGIBAGjIK69Ifr80dmzzzzFmjHT66dL//Z+0fLlztdXl5Zft4N3cBZ1CQp1Y5sUCAAxBiAUAoD5FRdKHH0ojR0rt2jX/PD6fNGWKHQZ/9Svn6qvLs89KrVpJ48a17Dx0YgEAhiHEAgBQn9desx+vc+mlLT/XiSdKV14pTZ9uL7oUKYsXS59/Ll1yiZSe3rJz0YkFABiGEAsAQH1eeUVKSbGHAzvh/vul1FTp9tvtcBwJoVWQb7ml5eeiEwsAMAwhFgCAumzdKs2cKZ17rpSZ6cw5u3SRbrvNnhcbemyPkzZvtufD/uhH0sCBLT8fnVgAgGEIsQAA1OXVV+3w5sRQ4oP99rfSEUdIkyZJJSXOnvuxx6Rg0A7KTqATCwAwjJEhNhgMqqioKOxXWVmZ1+UBAOLFK69IaWlSbq6z523dWrrnHmnbNvvRO07Zt0/6+9+lPn2kH//YmXPSiQUAGMbvdQHhFBYWKj8/P+y+3Nxc5eXluVwRACDubNokzZkjXXyxvcqv0665RvrLX+wVi6+/XurcueXnfO45afdu6Q9/ONBBbSk6sQAAwxgZYnNyclRQUBB2XyAQcLkaAEBceuMNe3vJJZE5v98vPfyw3THNz7cDaEtUV0t//rP9GKDx4x0pURKdWACAcYwMsX6/X1lZWV6XAQCIZ++8IyUl2c+HjZTRo6VzzpH++U/p1lulQYOaf65p06RvvrEDsZO/8KUTCwAwjJFzYgEA8FRZmfTRR9Lw4VJGRuSu4/PZ3VifT/rlL1sWFB95REpOduaxOgejEwsAMAwhFgCAH5o5U6qocG5xpPr072/Pj505U3r66ead49NPpdmzpcsvlzp1crS82hBLJxYAYAhCLAAAP/T22/b2vPPcud6DD9rPj731VmnZsqYda1l2Fzchwd46LTScmE4sAMAQhFgAAA5mWfZ82B49pKOPdueabdtK//qXtH+/NGqUtG5d4499+WV7FeXrr5eOO8752ujEAgAMQ4gFAOBgq1ZJ331nd2FDAc4Nw4dLL7wgbd1qLya1bVvDx5SUSL/6lR2C7703MnXRiQUAGIYQCwDAwd591966MR/2hy65RHriCXuV4dGjpT176n5vRYU9B3bzZun++6X27SNTE51YAIBhCLEAABxs5kz7Ga4jRnhz/RtukO65R/ryS3tocWHh4e/ZuFE64wzprbekn/zEHkocKXRiAQCGIcQCABBSU2PPLz3xRGeftdpUd94p/e530sKF0sCB0osvSpWV0qJFdmDt1cve99vfSs8+G9lhz3RiAQCG8XtdAAAAxli+XNq9256f6iWfzx4iPHCgdPXV0hVXSFddJVVV2fuHDpXuvtvu1EZaqBNLiAUAGIIQCwBAyOzZ9tbrEBty8cX2sOHnnpPefFMaNEi67DLp5JPdW3QqdB2GEwMADEGIBQAgZM4ceztsmLd1HKxDB+n22+0vL9CJBQAYhjmxAABIdkibPdt+1mq7dl5XYw46sQAAwxjZiQ0GgyoqKgq7LxAIKODlYhsAgNi0fr39uJrzz/e6ErPQiQUAGMbIEFtYWKj8/Pyw+3Jzc5WXl+dyRQCAmBeaD3vaad7WYRo6sQAAwxgZYnNyclRQUBB2H11YAEBEhObDmrKokynoxAIADGNkiPX7/crKyvK6DABAPJk9W+raVerSxetKzEInFgBgGBZ2AgCgqEhatYoubDh0YgEAhiHEAgAwd669ZT7s4ejEAgAMQ4gFAID5sHWjEwsAMAwhFgCAOXOktm2lY47xuhLz0IkFABiGEAsAiG9lZdKiRdKwYQe6jjiATiwAwDCO/996x44duv3223X00UcrLS1N7dq106BBg/SrX/3K6UsBANByn30m7d/PfNi60IkFABjG0RC7aNEiHXPMMZoyZYqSkpI0ZswYDRkyRLt27dKf/vQnJy8FAIAzQvNhCbHh0YkFABjGsefE7tixQ6NHj1Z5ebneeOMNnX/++Yfs//TTT526FAAAzpk7V0pOlgYP9roSM9GJBQAYxrEQe/fdd6uoqEiPP/74YQFWkk4++WSnLgUAgDOqq6V586STTpJSU72uxkx0YgEAhnFkOHF5ebleeOEFpaen66qrrnLilAAARN6KFdKePQwlrg+dWACAYRzpxH7++ecqKSnRaaedprS0NL377rv64IMPVFFRoT59+uiSSy5RTk6OE5cCAMA5c+fa22HDvK3DZHRiAQCGcSTErly5UpJ0xBFH6IILLtAbb7xxyP477rhDTz/9tC677DInLgcAgDNCizqdeqq3dZiMTiwAwDCOhNji4mJJ0ptvvqnExEQ9/vjjuvjii1VWVqbHHntMDz/8sCZMmKBjjjlGJ5xwQoPn27JliyZPntzselpyLAAgjsyZIx17rNS+vdeVmItOLADAMI7Mia35/rezwWBQ99xzj2666SZ16NBB3bp100MPPaSLL75Y+/fv10MPPeTE5QAAaLmNG6UNG5gP2xA6sQAAwzjSiW3VqlXtn8Mt7HTVVVfpP//5j2bNmtWo82VnZ9NNBQBEVmg+LCG2fnRiAQCGcaQT261bN0lSIBBQhw4dDtvfvXt3SdL27duduBwAAC0Xmg/Lok71oxMLADCMIyF24MCBkuxH7VRWVh62f9euXZIO7dgCAOCpOXOk7GypRw+vKzEbnVgAgGEcCbFdu3bVgAEDZFlW2CHDoddCYRcAAE/t2SMtXWoPJQ51GhEenVgAgGEcCbGS9Otf/1qSdPvtt2vLli21r3/55ZeaMmWKJOmGG25w6nIAADTf/Pl2Z5H5sA2jEwsAMIwjCztJ0uWXX673339fzz33nI499lideuqpKi8v17x581RZWanrrrtOF198sVOXAwCg+ULzYQmxDaMTCwAwjGMhVpKeeeYZDRs2TE8++aRmzpwpn8+nQYMG6frrr9eECROcvBQAAM03Z47UqpXUv7/XlZiPTiwAwDCOhlifz6frrrtO1113nZOnBQDAOVVV0qef2qsS+x3932BsohMLADCMY3NiAQCICgsXSuXl0vDhXlcSHejEAgAMQ4gFAMSX6dPt7ejR3tYRLejEAgAMQ4gFAMSX6dOl9u2lwYO9riQ60IkFABiGEAsAiB/btkmLF0vnnCMlJnpdTXSgEwsAMAwhFgAQP95/394ylLjx6MQCAAxj5LKMwWBQRUVFYfcFAgEFAgGXKwIAxITQfNiRI72tI5rQiQUAGMbIEFtYWKj8/Pyw+3Jzc5WXl+dyRQCAqGdZ0owZ0vHHS9nZXlcTPejEAgAMY2SIzcnJUUFBQdh9dGEBAM2yapU9J/ayy7yuJLrQiQUAGMbIEOv3+5WVleV1GQCAWPLRR/b2zDO9rSPa0IkFABiGhZ0AAPHh44/tQHb66V5XEl3oxAIADEOIBQDEvpoaO8SeeKKUmel1NdGFTiwAwDCEWABA7FuyRCouZihxc9CJBQAYhhALAIh9ofmwP/qRt3VEIzqxAADDEGIBALHv44+lpCRp2DCvK4k+dGIBAIYhxAIAYltFhTRzpjR0qJSe7nU10YdOLADAMIRYAEBsmzFDKi2V8vK8riQ60YkFABiGEAsAiG1vvmlvx4zxto5oRScWAGAYQiwAIHbV1Nghtm9fqXdvr6uJTnRiAQCGIcQCAGLXZ59JW7fShW0JOrEAAMP4vS4AABCH9u+XZs+Wvv1W6tBBOuUUKTvb+eu88Ya9JcQ2H51YAIBhCLEAAHe9+aZ07bXSjh0HXktNlaZMkW688UBocsIbb0gdO9ohGc1DJxYAYBgjhxMHg0EVFRWF/SorK/O6PABAc1iW9Nvf2l3R6mrpnnukd96RnnrK7sLefLN0+eXOdfxmz5ZWrpTGjj0QxNB0dGIBAIYxshNbWFio/Pz8sPtyc3OVx2MSACD6/Oc/0h//aHdF//1vqWvXA/suuUQaP156+WXp+OOlO+5o+fXuv98Or7fd1vJzxTM6sQAAwxgZYnNyclRQUBB2XyAQcLkaAECLVVVJv/udlJEhvfWWlJV16P6MDOnFF+2AO2mSNHy4/dVcixZJ06dL//M/rErcUnRiAQCGMTLE+v1+Zf3wBxwAQPT6+9+ltWulP/zh8AAbEgjYHdoTT5Quu0z68su639uQ+++3t7/7XfOOxwF0YgEAhmGSEAAgskpL7fmvnTtLt95a/3uPOUZ64glp82Z7eHFzun8rV0qvvSbl5Un9+zevZhwQCrF0YgEAhiDEAgAi6913pZ07pdtvl9LSGn7/hAn217vvSg8+2PTrPfCAvXViXi0IsQAA4xBiAQCR9dpr9vbCCxt/zOOPS/362UH0vfcaf9y6dfbc2jPPlIYMaVqdCI8QCwAwDCEWABA5lZX2Qk6nnGIPJ26s9HRp6lQpM9OeH/vtt407rqDAfnwPXVjnEGIBAIYhxAIAImfmTKmkRLrggqYf26uX3VXdvVs64wxpxYr63//hh9LTT9vvPfPMpl8P4RFiAQCGIcQCACJn2jR7e/75zTv+3HOlF16Qtm61H7kzb1749+3YYc+jTU+X/vGPA4+FQcsRYgEAhiHEAgAiw7KkN9+Ueva0Vx1urssvt4ckV1VJZ59tP//1YBs32iG5sNB+lE+PHi2rG4fiObEAAMMQYgEAkbF0qR0w8/Ja3hkdNUr66CN7deMxY6Sbb5b+8x/pN7+RBgyQFiyQfvtb6cornakdB/CcWACAYfxeFwAAiFFvvmlvmzuU+IdOPln6+GPpJz+xnyX7xBP26126SP/8p5Sb68x1cCiGEwMADEOIBQBExrRp9urCp53m3Dn795e+/NL+WrRIGjzYfi0x0blr4FCEWACAYQixAADnFRZKn31mPx4nKcnZc/t80sCB9hcijxALADCMkSE2GAyqqKgo7L5AIKBAIOByRQCAJnn7bXubl+dtHWg5QiwAwDBGhtjCwkLl5+eH3Zebm6s8figCALO9/77dMR01yutK0FKEWACAYYwMsTk5OSooKAi7jy4sABiupkaaOdMe7tuundfVoKV4xA4AwDBGhli/36+srCyvywAANMeKFVJRkTRhgteVwAk8YgcAYBieEwsAcNbHH9vbH/3I2zrgDIYTAwAMQ4gFADjr44/tR94MH+51JXACIRYAYBhCLADAOTU10qxZ0oknShkZXlcDJxBiAQCGIcQCAJyzZIlUXMxQ4lhCiAUAGIYQCwBwzsyZ9pYQGzsIsQAAwxBiAQDOmTvXDj1Dh3pdCZxCiAUAGIYQCwBwhmXZIbZ/f6l1a6+rgVN4TiwAwDCEWACAM9avl7ZulU491etK4CSeEwsAMAwhFgDgjPnz7S0hNrYwnBgAYBhCLADAGfPm2Vvmw8YWQiwAwDCEWACAM+bNkzp2lHr08LoSOIkQCwAwDCEWANBy+/bZz4g99dQDCwEhNhBiAQCG8XtdQDjBYFBFRUVh9wUCAQUCAZcrAgDU67PPpOpq5sPGIkIsAMAwRobYwsJC5efnh92Xm5urvLw8lysCANQrNB+WEBt7CLEAAMMYGWJzcnJUUFAQdh9dWAAw0Lx5UnKyNGiQ15XAaTwnFgBgGCNDrN/vV1ZWltdlAAAao6bGfrzOoEFSaqrX1cBpPCcWAGAYFnYCALTMmjVScTFDiWMVnVgAgGEIsQCAlpk/394SYmNXQgIhFgBgDEIsAKBlQos6DR3qbR2IHEIsAMAghFgAQMvMmyd17y7l5HhdCSKFEAsAMAghFgDQfLt2SatWMZQ41hFiAQAGcTTEnnHGGfL5fHV+TZ8+3cnLAQC8tmCBvWUocWzz+QixAABjROQROxdeeKFatWp12OtHHnlkJC4HAPBKaD4sndjYlpDAI3YAAMaISIh9+OGH1b1790icGgBgknnzpEBA6t/f60oQSQwnBgAYhDmxAIDmCQalhQulU06R/BH5nShMQYgFABiEEAsAaJ6lS6WyMoYSxwNCLADAIBH51fnTTz+tnTt3KiEhQX369NEFF1ygrl27RuJSAACvzJ1rbwmxsY8QCwAwSERC7H333XfI97fffrsmTZqkSZMmReJyAAAvzJljr1rLysSxjxALADCIoyH29NNP17XXXqtTTz1V2dnZ2rhxo1599VXdd999uuuuu5SRkaFbb721wfNs2bJFkydPbnYdLTkWANAIlmWH2H79pLZtva4GkUaIBQAYxNE5sffcc4+uvPJKHXXUUUpLS1OfPn10xx13aOrUqZLscFleXu7kJQEAXli/XioslE47zetK4AaeEwsAMIgry0mOHDlSJ554oj7//HMtXLhQZ5xxRr3vz87OppsKACYLzYclxMYHnhMLADCIa6sT9+7dW5I9VBgAEOXmzLG3w4Z5WwfcwXBiAIBBXAuxxcXFkqT09HS3LgkAiJQ5c6Qjj5S6dfO6EriBEAsAMIgrIXbHjh2aPXu2JGnQoEFuXBIAECnFxdKKFXYX1ufzuhq4gRALADCIYyF23rx5mjp1qqqrqw95/bvvvtPYsWNVWlqq888/X507d3bqkgAAL8yfb8+PZD5s/CDEAgAM4tjCTl999ZWuuuoqderUSYMGDVKbNm20fv16LVq0SBUVFerXr5+eeuoppy4HAPAKizrFH0IsAMAgjoXYU045RTfeeKMWLlyozz77TMXFxUpPT9cJJ5ygiy++WDfeeKPS0tKcuhwAwCtz5kitWknHH+91JXALj9gBABjEsRB7zDHH6IknnnDqdAAAE1VVSZ9+Kg0fLvldeUobTEAnFgBgEH4CAYBYVFQkvfGG1LGjdOqpUrt2zpx39mypokIaMcKZ8yE68JxYAIBBCLEAEGs++EAaP17autX+vk0baeZMacCAlp972jR7m5vb8nMheiQkSD9YuBEAAK+49pxYAIAL5s6VRo2SSkulv/5VevhhqaREuvHGlg8HtSw7xHbpIvXv70y9iA4MJwYAGIROLADEkj/9yQ6bc+YcCJrr19uB9l//kn7yk+afe/Vqae1a6aabeD5svCHEAgAMQicWAGLFpk3S1KnSyJGHdkp//3spK0v69a/trmxzMZQ4fhFiAQAGMTLEBoNBFRUVhf0qKyvzujwAMNP//q89b/Hmmw99vW1b6f777Tmy997b/PNPmyalp0s/+lHL6kT0IcQCAAzisyxzlhscPHiwFi9erKysLI0bNy7se3Jzc5WXl+dyZQBguGDQnquakiJ9+62UmHjo/upq6eSTpWXL7K+jj27a+QsL7fOPGSO99ppzdSM6nHSStH27PTQdAACPGTknNicnRwUFBWH3BQIBl6sBgCjw+ed2pzU///AAK9mv/eUv0mmnSb/4hfTOO02b1/rCC3Yn7oorHCsZUYRH7AAADGJkiPX7/crKyvK6DACIHh98YG9Hjqz7PcOGSVdeaQfSt99u/NxWy5Kee84elsx82PjEcGIAgEGMnBMLAGiiDz6w56sOGVL/+/74R/t9v/iFVFnZuHMvWiStXCldfrk9XBnxhxALADAIIRYAol1JiTR/vr3gUnJy/e/NyZEmTbLnzf71r407/3PP2dsJE1pWJ6IXIRYAYBBCLABEu1mz7IWdzjmnce+/9VZ7kaYpU6SqqvrfW1JiDz8+9ljpxBNbXiuiEyEWAGAQQiwARLv337e3jQ2xqanSTTfZC0FNnVr/e596Stq92x5+3JSFoBBbCLEAAIMQYgEg2s2caQ8T7tu38cdcc4099PiJJ+p+T2Wl3a3NzpbGj29xmYhihFgAgEEIsQAQzUpKpOXLpVNPbVqntEMH6ZJL7KHIM2aEf8+//mU/H/YXv2BBp3jn8xFiAQDGIMQCQDT77DP7ETgNrUoczu9/b69UfNVV0s6dh+7bulW6+24pM1O6/npnakX04jmxAACDEGIBIJotWGBvmxNijzpKevRRaeNG6eKLpeJi+/W1a6URI6RNm+z9mZnO1YvoxHBiAIBB/F4XAABogQULJL9fGjSoecdffbW0eLE9N3bAAOnMM6XXXpP27ZMee4zH6sBGiAUAGIROLABEK8uyQ+wJJ0hpac07h89nh9XQKsTPPSf16CG9/bZ0881OVotoRogFABiETiwARKt166QdO+wFmlrC55OuvVYaN84+X58+PE4HhyLEAgAMQogFgGi1cKG9bc582HDatbO/gB8ixAIADGJkiA0GgyoqKgq7LxAIKBAIuFwRABjos8/s7ckne1sHYh+P2AEAGMTIEFtYWKj8/Pyw+3Jzc5WXl+dyRQBgoGXLpEBA6tXL60oQ6+jEAgAMYmSIzcnJUUFBQdh9dGEB4HtLl0rHHWcHDCCSQveYZTFfGgDgOSNDrN/vV1ZWltdlAIC5tm2Ttm+XGJkCNxBiAQAG4df3ABCNli2zt/37e1sH4kMoxDKkGABgAEIsAESjUIg9/nhv60B8IMQCAAxCiAWAaLR0qb0lxMINhFgAgEEIsQAQjZYvl7KzJdYPgBsIsQAAgxBiASDaWJa0erV0zDFeV4J4EVrMiRALADAAIRYAok1hobRvn3T00V5Xgnhx8OrEAAB4jBALANFmzRp7S4iFWxhODAAwCCEWAKJNKMT27ettHYgfhFgAgEEIsQAQbejEwm2hEFtd7W0dAACIEAsA0Wf1aik1Vera1etKEC/oxAIADEKIBYBos2aN1Lv3gWABRFpior0lxAIADMBPQAAQTcrLpfXrmQ8Ld9GJBQAYxO91AeEEg0EVFRWF3RcIBBQIBFyuCAAM8fXX9mNOmA8LNzEnFgBgECNDbGFhofLz88Puy83NVV5enssVAYAhWNQJXmA4MQDAIEaG2JycHBUUFITdRxcWQFwjxMILDCcGABjEyBDr9/uVlZXldRkAYB5CLLzAcGIAgEFY2AkAosnq1VJ2tpSR4XUliCcMJwYAGIQQCwDRwrLsTixdWLiN4cQAAIMQYgEgWmzdKpWUEGLhPoYTAwAMQogFgGixerW95RmxcBvDiQEABiHEAkC0YFEneIXhxAAAgxBiASBaEGLhFYYTAwAMQogFgGixZo2UkiJ16+Z1JYg3DCcGABiEEAsA0WL1aql37wOBAnALw4kBAAYhxAJANKiokL77jqHE8EboFycMJwYAGIAQCwDR4Jtv7OfEEmLhBTqxAACDEGIBIBqwqBO8RIgFABjE73UB4QSDQRUVFYXdFwgEFAgEXK4IADzGM2LhJYYTAwAMYmSILSwsVH5+fth9ubm5ysvLc7kiAPAYnVh4iU4sAMAgRobYnJwcFRQUhN1HFxZAXFqzRurYUcrM9LoSxCNCLADAIEaGWL/fr6ysLK/LAAAzWJY9nPiEE7yuBPGK4cQAAIOwsBMAtNTq1dIZZ0ijR0v79zt//q1bpb17GUoM79CJBQAYhBALAC3x1lvS4MHSrFnSe+9Jdcznb5Hly+3tccc5f26gMQixAACDEGIBoLkqKqRrrpHS0qSPPpKGDJEeesj+s5NCIfb44509L9BYDCcGABiEEAsAzfXii9L27Xb39Uc/sr9PSpLuvNOex+qUZcvsLZ1YeIVOLADAIIRYAGgOy5L+9CepdWu7GytJPXpIEydK8+dLc+Y4d63ly+2ViTt0cO6cQFMQYgEABolYiN25c6eOOOII+Xw+9erVK1KXAQBvzJxph8trr5UyMg68fttt9vaxx5y5Tk2NtGIFQ4nhLYYTAwAMErEQ+8tf/lJFRUWROj0AeOvf/7a311576Ot9+0pnny3997/S7t0tv866dVJZGUOJ4S06sQAAg0QkxM6YMUPPPfecrrvuukicHgC8VVMjTZ1qP/Lm2GMP33/BBXbHatasll8rNB+WTiy8RIgFABjE8RBbXl6u66+/Xscee6xuv/12p08PAN5bsMB+duu4ceH3n322vf3ww5Zfa+lSe0uIhZcYTgwAMIjf6RP+/ve/19q1azVr1iwlJSU5fXoA8N7rr9vbsWPD7+/TR+rc2ZkQu2SJ3QVjODG8RCcWAGAQRzuxS5cu1ZQpU3TVVVdp+PDhTp4aAMxgWdJrr9kh9cQTw7/H57O7satXS5s3t+x6X35pD1tOS2vZeYCWIMQCAAziWIitqanRtddeqzZt2ujBBx906rQAYJbly6W1a+0urM9X9/tCQ4pnzGj+tfbuta81YEDzzwE4geHEAACDODac+K9//as+++wzPfPMM2rfvn2LzrVlyxZNnjy52ce35FgAqNfMmfZ29Oj633fWWfb2ww+l8eObd63Qok4nnNC84wGn0IkFABjEkU7shg0bdOedd2rEiBGaOHGiE6cEADPNmWN3YIcOrf99nTpJ/frZIdaymnetL7+0t3Ri4bVQJ5YQCwAwgCOd2JtvvllVVVX6+9//7sTplJ2dTTcVgHksS5o71w6nbds2/P6zz5YefVRatSr8o3gasmSJvSXEwmuhTizDiQEABnAkxL711ltq06aNbrjhhkNer6iokCRt3rxZZ5xxhiTp5ZdfVqdOnZy4LAC4a8MGe6Gm3NzGvT8UYmfMaH6I7dDB7uoCXmI4MQDAII7Nid29e7dmzZoVdl9FRUXtvlCwBYCoM3euvT3ttMa9f9gwe+jx3LnSz37WtGtVV9tzYkPnALzEcGIAgEEcmRNrWVbYr3Xr1kmSevbsWfta9+7dnbgkALhvzhx7O2xY497ftq39fNfZs5s+L/brr6XychZ1ghkYTgwAMIijz4kFgJg2d66UkyM15Zdxp50mFRZK333XtGsxHxYmYTgxAMAghFgAaIw9e5o3vPf79QD0/vtNux4hFiZhODEAwCCEWABojIUL7SHBp57atONGjZL8fmnatKYd9+WXUnKy1Ldv044DIoHhxAAAgzi2sFM43bt3l9Xc5yMCgEkWLLC3DT0f9ocyM6XTT7dXKC4rkwKBho+xLOmLL+xH+SQlNb1WwGkMJwYAGIROLAA0xoIFdme0OQst5eZKFRXSRx817v3ffitt3dr0ri8QKQwnBgAYhBALAA2xLDvEDhokpaQ0/fjQc2Xfeqtx7//kE3t7+ulNvxYQCQwnBgAYhBALAA35+mupuFgaMqR5x/fuLfXpY4fYxkyxmD3b3g4f3rzrAU5jODEAwCCEWABoSGg+bHNDrGR3Yzdvthdsasgnn0i9eknZ2c2/HuAkhhMDAAxCiAWAhjgRYvPy7G1DQ4o3bZLWrqULC7MwnBgAYBBCLAA0ZMECqVMnqWvX5p9j2DB7peKGHrUzY4a9PfPM5l8LcBrDiQEABonoI3aaKxgMqqioKOy+QCCgQGMeUQEATigtlZYutTupPl/zz5OUJI0eLb3yir3ycKdO4d/34Yf29uyzm38twGkMJwYAGMTIEFtYWKj8/Pyw+3Jzc5UXGpYHAJG2aJE9hLIlQ4lDcnPtEDttmnTddYfvr66WPvhAOv74ukMu4AWGEwMADGJkiM3JyVFBQUHYfXRhAbhq/nx7e8opLT9Xbq79iJ7nnw8fYt97T9q2TfrpT1t+LcBJDCcGABjEyBDr9/uVlZXldRkAIH38sZSc7EyIbdNGGjdOeukl6Ztv7BWID/bkk3ZYuPball8LcBLDiQEABmFhJwCoy/790pw50tChUlqaM+e86ip7+/TTh76+caO9cvG557ZsASkgEhhODAAwCCEWAOry2Wf2wk4/+pFz5zzrLOnoo6VHH5W+++7A6088YXe5rr/euWsBTmE4MQDAIIRYALGtulqaPl264grpz39u2rEff2xvnQyxCQnSY49J5eXSLbdIliW9/7704INSv352JxYwDcOJAQAGMXJOLAA4Zvx46cUX7T+/+KLUo4c0Zkzjjv3oI3sYsRPzYQ929tnS5Zfb9QwebD/Cp3Vr6fXXJT8fyzAQw4kBAAahEwsgdi1caAfFs86SZs+W2reXJkyQvv224WMrK6V586Rhw+wVhZ32179KP/mJtHy53eV66SWpd2/nrwM4gU4sAMAghFgAsSs/3+4g/eUv0mmn2YF2717p4ovtRZvqs2CBVFHh7FDig7VrJ/3zn9KmTdKaNQwjhtmYEwsAMAghFkBs+vhjacYMu9t57LH2ayNHSnfcIX3xhfTwww0fL0UuxIYccYTUvXtkrwG0FMOJAQAGIcQCiE3PPWdv77zz0NcnTZL69pV+/3upsLDu499+W2rVSjrxxMjVCEQLhhMDAAxCiAUQe2pq7BWJ+/WTevU6dF9Kir0ScGWl9NBD4Y//4gvp88+lSy+VkpIiXy9gOoYTAwAMQogFEHuWL5e2bZNGjQq/PzdXOuEE6ckn7ff90JNP2tsbbohYiUBUYTgxAMAghFgAsWfuXHtb13xWn0+66y77Wa0//7n9rNaQkhLpX/+yH33DUGLARicWAGAQIx9IGAwGVVRUFHZfIBBQIBBwuSIAUWXePHs7ZEjd77ngAmncOOnf/5YGDpR++1v79aeflvbtk66/PuJlAlElIYEQCwAwgpEhtrCwUPn5+WH35ebmKi8vz+WKAESV+fPtZ65mZdX9Hp/PXvzpq6/sFYuLiqRu3aRf/lLq1Em67DL36gWiQUICw4kBAEYwMsTm5OSooKAg7D66sADqtX279O230oQJDb+3VSvpzTeliy6SpkyxX+vRQ3rvPXsfgAMSE+nEAgCMYGSI9fv9yqqvgwIAdZk/394OHdq49/foIX32mT0P9uOPpfvvtzuxAA7FcGIAgCGMDLEA0Gyh+bCnntr4YxISpJ/8xP4CEB7DiQEAhmB1YgCxZf58qXVr6dhjva4EiC2JiYRYAIARCLEAYkdVlT00eMgQ+wduAM4hxAIADEGIBRA7liyRKioaPx8WQOMRYgEAhiDEAogdixfb25NP9rYOIBYRYgEAhiDEAogdy5fb2+OP97YOIBYRYgEAhiDEAogdy5fbizp16eJ1JUDsIcQCAAxBiAUQO1askPr1k3w+rysBYg8hFgBgCEIsgNiwfbu0Y4d03HFeVwLEJp4TCwAwBCEWQGwIzYclxAKRkZgo1dR4XQUAAIRYADFixQp726+ft3UAsYrhxAAAQxBiAcQGOrFAZBFiAQCG8HtdQDjBYFBFRUVh9wUCAQUCAZcrAmC85cul9u2ljh29rgSITYmJUlWV11UAAGBmiC0sLFR+fn7Yfbm5ucrLy3O5IgBGsyw7xA4cyMrEQKTQiQUAGMLIEJuTk6OCgoKw++jCAjjM5s3S3r0MJQYiiRALADCEkSHW7/crKyvL6zIARIvQfFgWdQIihxALADAECzsBiH4s6gREHiEWAGAIQiyA6MfjdYDII8QCAAxBiAUQ/VaskLKzpXbtvK4EiF2EWACAIQixAKLfN99Iffp4XQUQ2xISCLEAACMQYgFEt127pOJiqVcvrysBYhudWACAIQixAKLbN9/YW0IsEFmEWACAIQixAKIbIRZwByEWAGAIQiyA6BYKsT17elsHEOsIsQAAQxBiAUS3r7+2t4RYILISE6WaGq+rAACAEAsgyn3zjXTEEVJGhteVALEtFGIty+tKAABxjhALIHpVV0vLlkn9+nldCRD7EhPtLd1YAIDH/F4XEE4wGFRRUVHYfYFAQIFAwOWKABhp9WqptFQ66SSvKwFiXyjEVlcf+DMAAB4wMsQWFhYqPz8/7L7c3Fzl5eW5XBEAI33+ub098URv6wDiwcEhFgAADxkZYnNyclRQUBB2H11YALUWLbK3gwd7WwcQDwixAABDGBli/X6/srKyvC4DgOkWLZLatpV69PC6EiD2EWIBAIZgYScA0SkYlL74wu7C+nxeVwPEvoTvf2QgxAIAPEaIBRCd1qyRyssZSgy4hU4sAMAQhFgA0enLL+3twIGelgHEDUIsAMAQhFgA0WnJEnvbv7+3dQDxghALADAEIRZAdFq6VEpNlXr39roSID4QYgEAhnAsxD7yyCMaN26cevfurczMTKWkpKhbt24aP368li1b5tRlAMC2ZIl03HGS38hF1oHYQ4gFABjCsRB7//33691331W7du101lln6bzzzlNqaqqef/55DR48WG+99ZZTlwIQ77Zvl7ZuZSgx4CZCLADAEI61MN544w0NHjxYqamph7z+xBNP6Oabb9a1116rTZs2yU/XBEBLhRZ1GjDA0zKAuEKIBQAYwrFO7LBhww4LsJJ00003qWfPntq2bZtWrlzp1OUAxLP58+3tKad4WwcQT0IhtqbG2zoAAHHPlYWdkpKSJEnJycluXA5ArJs/X0pJ4fE6gJvoxAIADBHxEPv8889rzZo16t27t3qziiiAlqqpkRYskAYPlvjFGOAeQiwAwBCOT1B96KGHtGLFCpWWlmrVqlVasWKFcnJy9NJLLykx9D9AAGiu1aulPXukU0/1uhIgvhBiAQCGcDzEvvfee5oxY0bt9926ddM///lPDR48uNHn2LJliyZPntzsGlpyLADDhebDDh3qbR1AvEn4fvAWIRYA4DHHhxN/+OGHsixLxcXF+uSTT9S7d2+NGDFCBQUFTl8KQDwixALeoBMLADBExJ5306ZNGw0fPlzvvPOOhg4dqkmTJmnkyJE66aSTGjw2OzubbiqA8ObPl7p1k7Kzva4EiC+EWACAISK+sFNSUpIuvfRSWZaladOmRfpyAKLJV19JP/2pNGaMPc+1Ibt3SytX0oUFvECIBQAYImKd2INlZWVJknbs2OHG5QBEg48+kkaOPPAD8fnnS++9J4V53nStuXPtLSEWcB8hFgBgCFeeEztr1ixJUs+ePd24HIBo8MQT9uNyPvhAuvVW6ZNPpNtuq/+YN96wt6NHR74+AIcixAIADOFIiJ07d66mT5+umpqaQ17fv3+//vrXv+r5559XWlqaLr30UicuByDaWZY0Z440YIB09tnSI49Iw4ZJTz0lffdd+GOqq6WpU6V+/aQ+fdysFoBEiAUAGMOR4cRff/21rrrqKmVlZWnw4MFq3769ioqKtGzZMm3ZskWpqal69tln1aVLFycuByDaffuttG2bdMkl9vcJCdLkydI550gPPCD9/e+HHzNnjrRjh3TDDa6WCuB7hFgAgCEcCbEjRozQHXfcoVmzZmnp0qUqKipScnKyunfvrosuukg///nP1atXLycuBSAWzJ5tb0877cBrZ50lDRkiPfOMNGmSdOSRhx7z2mv2dtw4d2oEcChCLADAEI6E2B49evAcWACNt2CBvR027MBrPp90xx32Ak9TpthDjEP277dDbI8e9hBkAO4jxAIADOHKwk4AcIjFi6WOHQ/vtubmSv37S08+KW3ffuD1P/1J2rRJuuYaO+wCcB8hFgBgCEIsAHft3y8tWyYNHHj4Pp9PuusuqaxMmjjRXr14/Xrp97+XevWSfvlL18sF8L1QiP3BIo4AALjNlefEAkCtr76SKivDh1jJnvM6caL07LPSxRdLq1fbofbxx+t/hiyAyEr4/vfedGIBAB4jxAJw17p19rauxd58PvsZsl9+ac+DTUyU7r5bGjnStRIBhBHqxAaD3tYBAIh7hFgA7tq40d7W98ittDRp1ixpwwapd28pJcWd2gDUzf/9jwx0YgEAHmNOLICmsSxp8+bmH9+YECtJGRnScccRYAFTEGIBAIYwMsQGg0EVFRWF/SorK/O6PCC+/f73UufO0nvvNe/4xoZYAGZhODEAwBBGDicuLCxUfn5+2H25ubnKy8tzuSIAtaZMsbcffCCNGtX04zdskNq1k9LTna0LQGSFOrGEWACAx4wMsTk5OSooKAi7LxAIuFwNgFqVlVJFhf3nTz9t3jk2bqQLC0QjhhMDAAxhZIj1+/3KysryugwAP7Ry5YEuzBdfNP34mhpp06bmdXABeItOLADAEEbOiQVgqE2bDvx53z6ppKRpx2/fLu3fTycWiEbMiQUAGIIQC6Dxtm61tyecYG+3bGna8SzqBEQvhhMDAAxBiAXQeKEQO3Dgod83FiEWiF4MJwYAGIIQCxgsWBNUdY1BXY89e+xtz572tri4aceHQmzXrs7VBMAdDCcGABiCEAsYqCJYobs+vkut7m8l/71+BQoC+n/T/5/2Vu71trB9++xt5872dvfuph2/YYO9pRMLRB+GEwMADGHk6sRAPKsMVurUp0/VF1u/0LEdjtVxRxynNUVr9OeFf9b7a9/Xl9d/qaTEJG+K+2GIbU4n1ueTjjzS2boARB7DiQEAhqATC7hsR+kOTZk3RaNeGKXef+2tuz++W/ur99fuf+KzJ/TF1i/0s5N/pi+v/1KvXPSKFl+/WLcNuU0rd6zUs18+613x+/bZP8h27Gh/39RO7MaN9rHJyY6XBiDCGE4MADAEIRZwSWFJoa6fdr2OfORI3f7B7Zq9frZ2V+zWPZ/coyteu0KStL96vwpmF6hrZlf98ew/1nZcE3wJmnzGZGWkZOipxU9595fYt09q1Upq29b+vjmdWIYSA9GJ4cQAAEMQYoEI21OxR/kz8tXrL730v4v/VycfebJevvBlFf+mWJtv26yzepyl/676r3aU7tCCTQu0s3ynrj7haqUlpR1yntYprXVur3P1WeFn2rZvmzd/mVCIbdPG/r4pndhg0H4kDyEWiE4MJwYAGIIQC0SIZVl68vMn1fMvPXX/nPt1VNujNO2yaZp91WxdetylSvGnKDkxWRMGTFCNVaM317ypd795V5J0bu9zw54zt0+uJOmdr99x7e9xiFCIDQSkpKSmdWILC6WaGlYmBqIVw4kBAIYgxAIR8rfP/6Yb3r5Bqf5U/eP8f2jJDUuU2ydXPp/vkPed1+c8JfgSNHXNVE1dPVXt09prcPbgsOc8+6izJUmfbPgk4vWHFQqxPp/djW1KJ5ZnxALRjeHEAABDGLk6cTAYVFFRUdh9gUBAgUDA5YqApvlq51e6/f3b1b1Ndy3+6WK1TWtb53vbpbXTCZ1O0FtfvSVJ+uXQXyoxITHsezu16qSebXtqzoY5Eam7Qfv2ST162H9u27ZpnVhCLBDdGE4MADCEkSG2sLBQ+fn5Yffl5uYqLy/P5YqAxgvWBDVh6gRVBCv03AXP1RtgQ8b3H6/FWxbr9qG364/n/LHe95505El6efnLKq0qVXpyulNlN05pqZT+/TXbtLGHCDcWz4gFohvDiQEAhjAyxObk5KigoCDsPrqwMN1Dcx/Sgk0LdNuQ23R6t9MbdcytQ27Vpcddqk6tOjX43qPbHy3J7vYOzB7YolqbJBiUKirs4cSS3YldsaLxx9OJBaIbw4kBAIYwMsT6/X5lZWV5XQbQZEu2LtHdM+/WMVnHqOCs8L+IqUtjAqwk9c3qK0las3ONuyG2tNTehkJsmzb2a/v324s8NWTjRvuH4E6N+3sCMAzDiQEAhmBhJ8AhlcFKjZ86XjVWjf459p9K9adG5DqhTuyaojUROX+d9u2ztwd3YqXGL+60davUseOBIYkAokvC9z8yEGIBAB4jxAIOuXvm3Vq6banuPP1OnZhzYsSu06d9H0nS6p2rI3aNsH4YYkPPim3s4k5FRRIjLIDo5fPZv4RiODEAwGOEWMABn6z/RA/OfVAn5Zyk/OHhFyVzSnpyujpndPa+ExsKsXv2NO74nTsJsUC08/vpxAIAPEeIBVpob+VejX99vFL9qXp+7PNKSmzE/NAWOrr90fpq51cRv84h6gqxjRlOHAza72vfPgKFAXBNYiIhFgDgOUIs0EK3vHOL1u9Zrykjp+jorKNdueYR6UeodH+pKoOVrlxP0oEQe/AjdqTGdWJ37bK3dGKB6Ob3M5wYAOA5QizQAs99+ZyeX/q8cvvk6oYTb3DtuhkpGZKkkqoS1655WCc2M9PeNqYTW1Rkb+nEAtGN4cQAAAMQYoFmWl20Wje9c5OObH2knhnzjHw+n2vXzkyxA+SeikbOR3VCS0Lszp32lhALRDeGEwMADECIBZqhfH+5Ln31UlUEK/TShS8pK+DuMNlQJ3Zv5V73LvrDEBvahp4fW59QJ5bhxEB0YzgxAMAAfq8LAKLRL9//pZZuW6p7f3Svhncb7vr1QyF2T6WLndhQWA2F19Dc2MaEWDqxQGxgODEAwAB0YoEmen3V6/rb53/TmT3O1O9O+50nNWSm2kN56cQCcBXDiQEABiDEAk2wo3SHrn/rerVPa68Xxr6gxIRET+owYjgxnVgg/jCcGABgAIYTA01wy7u3aEfZDr104UvKbp3tWR2eLuwUCq9pafa2KZ1YQiwQ3RhODAAwgJEhNhgMqij0Q+8PBAIBBQIBlysCpFdXvqp/r/i3xvYdq0v7XeppLZ52YkMhNiFBCgQaH2KTkqTWrSNXH4DIYzgxAMAARobYwsJC5efnh92Xm5urvLw8lytCvNtRukM3vX2T2qe119/O+5urj9MJx7MQm5pq/xAb0rq1tLcRNWzaJB15pOTxvxuAFvL7pcpKr6sAAMQ5I0NsTk6OCgoKwu6jCwsvHDyMuGOrjl6XU7uwk6urE1dUHBhCHJKVdWCocH02bpSOPTYydQFwD8OJAQAGMDLE+v1+ZbGKKQwRGkY87phxng8jDvGkE1tRYXdiD5aVJS1b1vBxO3dKnTtHrjYA7mA4MQDAAKxODNTj4GHET/z4Cc+HEYek+lOVnJjsfoj9YSe2Qwdp1676f6gNrUzML6aA6MfqxAAAAxBigXrcOv1W7Sjbocd//LgRw4gPlpGS4f5w4h92YkOrDe/eXfdxhFggdjCcGABgAEIsUIeP1n2kl5a/pPOPPl+X9LvE63IOk5mS6e4jdsrLDw+xGfawZu2ppw4erwPEDoYTAwAMQIgFwqiqrtIt79yiVH+qHh39qDHDiA+WmZrpfSc2015gqt4QG+rEEmKB6MdwYgCAAQixQBh/WfgXrSpapd+d9jt1b9Pd63LCykjJcLcTG25ObCjE1veYHUIsEDsYTgwAMAAhFviBXeW7dM+se3RU26P062G/9rqcOmWmZGpv5V5ZluXOBZvbiWU4MRA7GE4MADCAkY/YAbz0l4V/UUlVif523t+U6k9t+ACPZKZmqtqqVtn+MqUnp0f2YpYVPsQ2Zk4sCzsBsYPhxAAAA9CJBQ6yt3KvHl34qHq166VLjzPjmbB1yUyxu6CuzIutrLS3zIkF4lsoxLo1AgQAgDAIscBBHv/0ce2u2K3fnfY7+RPMHqhQG2LdmBdbUWFv6wqx9c2J3bpVatVKSo9wtxhA5Pm//1xkSDEAwEOEWOB7pVWlemTBI+qa2VVX9r/S63IalJnqYic2FGLrWtipvk5sYaGUkxOZugC4KynJ3hJiAQAeIsQC38v/KF9FZUX6zbDfKDkx2etyGmREJ7Yxc2IJsUDsCHVi9+/3tg4AQFwjxAKSXl35qh5d+KiGdx2u6wZd53U5jZKRYgdIVzqx5eX2tqnDicvLpeJiQiwQK0KdWEIsAMBDRk76CwaDKgo9luMHAoGAAoGAyxUhln2y/hP95PWf6Ij0I/TyRS8rKTHJ65IapXY4sZed2KQke4hxXZ3YLVvsLSEWiA0MJwYAGMDIEFtYWKj8/Pyw+3Jzc5WXl+dyRYhVi7csVt5LeUpKSNI7l7+jnNbRE7ZCw4n3VtazqJJT6poTK9lDiusKsYWF9pYQC8QGhhMDAAxgZIjNyclRQUFB2H10YeGUdcXrNOqFUaoMVmr6ldM1OGew1yU1iScLO/2wEyvZQ4oJsUB8oBMLADCAkSHW7/crKyvL6zIQw2qsGk2YOkFFZUWaeulUndH9DK9LajIjFnaS7BC7dWv44wixQGyhEwsAMAALOyEu/XnBnzV7w2zdctItGtN3jNflNIurndi6FnaS6MQC8YSFnQAABiDEIu6s3LFSd8y4Q73b9dYfz/mj1+U0W3pSuhJ8Cd4+J1aS2re3VyeurDx8XyjEZmdHrjYA7mE4MQDAAIRYxJX91fs1/vXx2l+zX/8c+08FkqJ3jrXP51NGSob3w4k7drS327cfvq+wUGrTRmIuOxAbGE4MADAAIRZx5f7Z92vRlkX6zbDfaEjnIV6X02JtUttod8XuyF+ovhDbqZO93bbt8H2FhQwlBmIJnVgAgAEcC7FlZWWaOnWqrrnmGh199NFKTU1Venq6BgwYoHvuuUf79u1z6lJAsywqXKT7Zt+n/h376+4Rd3tdjiPaprZVcUVx5C9UX4jt0MHe7thx+D5CLBBbmBMLADCAYyH2xRdf1NixY/WPf/xDiYmJOv/88zV8+HCtW7dOd999t0466SRtDzfcEHBBRbBC46eOl08+/fOCfyrFn+J1SY5om9ZWxeUuhNjQwk7h5sRm2gtMqaTk0Nf37bNfI8QCsYPhxAAAAzgWYpOSkvTTn/5UK1eu1MqVK/Xvf/9b06dP15o1azRw4ECtXr1av/jFL5y6HNAkd350p1buWKnJZ0zWgE4DvC7HMW1T22pP5R5V11RH9kL1dWIzMuzt3r2Hvr5li70lxAKxg+HEAAADOBZiJ0yYoCeffFLHHHPMIa9nZ2fr8ccflyS99tprqqqqcuqSQKN8sv4TPTL/EQ3pPES/HvZrr8txVNvUtpJceMxOc0IsKxMDsYdOLADAAK4s7DRggN35qqys1M6dO924JCBJKqks0cSpE5XqT9VzFzwnf4Lf65Ic1TbNDrERH1Lckk4sIRaIHXRiAQAGcOUn+rVr10qyhxy3a9fOjUsCkqRfffArrdu9Tn8Z/Rf1ad/H63IcF+rERnxxp/qeExsKsXt+0A0mxAKxh4WdAAAGcKUT++ijj0qSRo8erZSU2FhQB+ab/s10PbnoSZ3Z40zdfPLNXpcTEe3S7F8KRbwTG1rYqb7ViUOhNYQQC8QehhMDAAwQ8U7sO++8o6efflpJSUm69957G3XMli1bNHny5GZfsyXHIjYUlxfrmjevUUZKhp4Z84wSfLH5SOTa4cRudWKTkw/fl5pqPyt2/fpDXyfEArGH4cQAAANENMSuXr1aV155pSzL0kMPPVQ7NxaItFvevUWFJYV6Zswz6prZ1etyIqZ2OLEbc2JTUyWfL/z+7t2l77479LUtW+yhxoFAZGsD4B6GEwMADBCxELt582aNHj1axcXFuu2223Trrbc2+tjs7Gy6qWi2md/N1IvLXlRenzxNGDDB63IiKtSJ3VW+K7IXqqgIPx82pFs3acGCA2FXkjZvpgsLxJrQcGI6sQAAD0VkjOWuXbs0cuRIrV+/XldddZUefvjhSFwGOIxlWbpjxh3yJ/j159F/lq+uzmGMcG1hp/Ly8PNhQ7p3t7ehIcVlZdJXX0n9+kW2LgDuohMLADCA4yF23759Ovfcc7Vy5UqNGzdOTz31VMwHCZhj2lfTNH/TfP100E91VNujvC4n4lx9xE59IbZXL3v71Vf2dvlyqaZGGjgwsnUBcBcLOwEADOBoiK2srNSYMWP06aefatSoUXrppZeUmJjo5CWAOtVYNcr/KF9p/jTdefqdXpfjisyUTEkuLexUX4g95hh7u2qVvV250t7SiQViCws7AQAM4FiIra6u1mWXXaaPPvpIw4cP12uvvabkcCuZAhHy0rKXtHz7cv38lJ8ru3V8zMVMTEhUZkqm9yG2b197+8MQe+yxka0LgLsYTgwAMIBjCzs99thjev311yVJWVlZuummm8K+7+GHH1ZWVpZTlwUkSVXVVbpr5l3KTMnUb4b9xutyXNUurZ07z4nt1Knu/e3bS0ccYQ8jluwwm5Qk9ewZ2boAuIuFnQAABnAsxBYXH/ghOhRmw5k8eTIhFo57evHTWlu8VgVnFtTOE40XbdPaRn514rKy+lcnlqSTT5amT5dKS6XPP7eHEvsj/ihqAG6iEwsAMIBjw4knT54sy7Ia/OoeWsUUcEjZ/jLd+8m9OiL9CN16SuMf5RQr2qa2dacT29DzXkeMsLszDz4obd0qnXZaZGsC4D5CLADAABF5xA7gpv9b/H/asm+L8ofnKz053etyXNc2ra32VO5RdU11ZC5QXS1VVTXciT33XHt7zz1SYqL0859Hph4A3mE4MQDAAIRYRLWq6io9NO8hdUzvqOsGXed1OZ4IPSt2d8XuyFygvNzeNtSJ7ddPmjjR/vO110q9e0emHgDeoRMLADAAE9YQ1f655J/atHeT/nj2H5WW1ECnMEaFQmxxRbHaB9o7f4GyMnvbUCdWkh5+WDrySOn//T/n6wDgPTqxAAADEGIRtYI1QT0w5wG1TW2rG0+80etyPBNayCpi82Ib24mV7FWK77svMnUA8B6dWACAARhOjKj17xX/1rfF3+rnp/xcrVNae12OZw7uxEZEUzqxAGJbKMTSiQUAeIgQi6hUY9Xo/tn3q1VyK/38lPheQMioTiyA2BYaTkwnFgDgISOHEweDQRUVFYXdFwgEFOCH6bj35po3tWLHCv3q1F+pXVo7r8vxVOjvTycWQMQxnBgAYAAjQ2xhYaHy8/PD7svNzVVeXp7LFcEklmXpD3P+oJTEFN029Davy/Fc7XDiSHViQyGWXx4BoBMLADCAkSE2JydHBQUFYffRhcVH6z7Sp5s/1Y0n3qhOrTp5XY7naocTR6oTGxpOTCcWgM8nJSfbz44GAMAjRoZYv9+vrKwsr8uAof4w5w9K9CXqV6f+yutSjBDqxO4q3xWZC9CJBXCwpCRCLADAUyzshKjy6eZPNWPdDF12/GXq0baH1+UYITM1Uz756MQCcAedWACAxwixiCp/mPMHSdJvh/3W40rMkeBLUGZqJnNiAbgjOZk5sQAATxFiETVWbF+hqaunaszRY9TviH5el2OUtqlt6cQCcAedWACAxwixiBpT5k+RJP32NLqwP9Q2rW3kOrEVFfY2NTUy5wcQXZgTCwDwGCEWUWF76Xb9a9m/dGqXUzWk8xCvyzFORDuxlZX2NiUlMucHEF3oxAIAPEaIRVT4++d/V1V1lf7fkP/ndSlGapfWTnsr96q6ptr5kxNiARyMObEAAI8RYmG8/dX79cRnT6hbZjdd0PcCr8sxUugxO7srdjt/8lDHJTnZ+XMDiD50YgEAHiPEwnjvfP2OtpVu008H/1T+BCMfbey5tml2iI3IkGI6sQAOxpxYAIDHCLEw3nNLnpNPPo0fMN7rUowV6sTuKt/l/MkJsQAORicWAOAxQiyMtqN0h6Z9NU3n9DxHnTM6e12OsWo7sZFYoZjhxAAOxpxYAIDHCLEw2kvLX1KwJqiJAyZ6XYrRQp3YiA4nJsQCkOjEAgA8Z+QEw2AwqKKiorD7AoGAAoGAyxXBK89++awyUjJY0KkBEe3EVlbaP7T6fM6fG0D0YU4sAMBjRobYwsJC5efnh92Xm5urvLw8lyuCF5ZtW6Yvtn6h6wZdp7SkNK/LMVpEO7FVVXRhARyQnCxVV0s1NVICA7oAAO4zMsTm5OSooKAg7D66sPHjtVWvSZIuO+4yjysxX8Q7sSzqBCAk9Eut/fv5bAAAeMLIEOv3+5WVleV1GfDY1DVT1S6tnYZ3G+51KcZrl9ZOEp1YAC4IfR5UVRFiAQCeYBwQjPTd7u/05dYvldcnj2fDNkJGSoZ88kVuYSd+UAUQkpRkb5kXCwDwCCEWRpq6eqokaWzfsd4WEiUSfAlqk9qG4cQAIu/gTiwAAB4gxMJIU1dPVZo/Tef0PMfrUqJG27S22lW+y/kTM5wYwMEOnhMLAIAHCLEwTlFZkWZvmK1RvUYpkMRCXo3VNrUtw4kBRB6dWACAxwixMM60NdNUY9XogqMv8LqUqNI2ra12V+x2/sSEWAAHY04sAMBjhFgY5+UVLyspIUl5R/M84KbISMlQSWWJaqwaZ0/McGIAB6MTCwDwGCEWRtm2b5tmrJ2h0b1G1z42Bo2TkZIhS5ZKq0qdPTGdWAAHY04sAMBjhFgY5T8r/6Nqq1qXH3+516VEncyUTEnSnso9zp6YTiyAg9GJBQB4jBALo7y0/CUFkgLK68NQ4qbKSMmQJO2t3OvcSWtqpGCQTiyAA5gTCwDwGCEWxthVvkvzN87Xub3OVXpyutflRJ2IhNjKSntLiAUQQicWAOAxQiyMMfO7mbJk6eyjzva6lKhUO5y4wsHhxKEQy3BiACGEWACAxwixMMaMtTMkSWf1OMvjSqJTRDqx5eX2NsDzegF8LzQyI/RLLgAAXOb3uoBwgsGgioqKwu4LBAIK8AN1TJqxboa6ZHRRr3a9vC4lKmWm2p1YR0Ns6fcrHaczvBvA91JT7W1Fhbd1AADilpEhtrCwUPn5+WH35ebmKi+PRX9izea9m7Vm5xpNPGGifD6f1+VEpVAn1tHVicvK7C2/OAIQQicWAOAxI0NsTk6OCgoKwu6jCxubPlr3kSTpzO5nelxJ9IrIcOJQiKUTCyCEEAsA8JiRIdbv9ysrK8vrMuCiGeu+nw97FPNhmyu0sFNEhhPzyyMAIaEQy3BiAIBHWNgJnrMsSzPWzVDfrL7KaZ3jdTlRq3Y4sZOrEzOcGMAPhebE0okFAHiEEAvPrdm5Rpv2bmJV4hZqldxKkrS3ioWdAEQQw4kBAB4jxMJz7379riTp3F7nelxJdEtMSFTr5NaRmRNLJxZACCEWAOAxQiw8N/3b6UpOTNYZ3c/wupSol5GSwXBiAJHFnFgAgMcIsfBU2f4yzfpulkZ0G6H0ZIastlRGSgbPiQUQWcyJBQB4jBALT838bqYqqys1utdor0uJCZmpmQwnBhBZDCcGAHiMEAtPMR/WWRkpGdpTyXBiABFEiAUAeIwQC09N/3a6umV2U9+svl6XEhMyUjK0r2qfqmuqnTkhw4kB/FBysr1lTiwAwCOEWHhmbfFafbPrG43uNVo+n8/rcmJCZkqmJGlf1T5nTkgnFsAP+Xx2N5ZOLADAI4RYeGb2+tmSxKrEDspIyZAk54YU7/s+DKelOXM+ALGBEAsA8BAhFp6Zu3GuJGlYl2EeVxI7QiHWscWdSkrsocSJic6cD0BsIMQCADxEiIVn5m6cqy4ZXdQls4vXpcSM0HBix54Vu2+f1Lq1M+cCEDtSU5kTCwDwjN/rAsIJBoMqKioKuy8QCCjA/Lyot6t8l1buWKn/Oe5/vC4lpkSkE9uqlTPnAhA76MQCADxkZIgtLCxUfn5+2H25ubnKy8tzuSI4bf7G+ZIYSuy0iITYtm2dOReA2EGIBQB4yMgQm5OTo4KCgrD76MLGhvmbCLGRkJn6/XBipxZ2KimRunZ15lwAYkdKilRc7HUVAIA4ZWSI9fv9ysrK8roMRNAXW79QcmKyjjviOK9LiSmOd2KZEwsgHObEAgA8xMJO8MSSrUvUr0M/JSUmeV1KTAkt7ORIiN2/3x4uyJxYAD/EcGIAgIcIsXBdUVmRNpds1oBOA7wuJebUPifWidWJS0rsLZ1YAD9EJxYA4CFCLFy3ZOsSSdKAjoRYp9UOJ65yoBO7b5+9JcQC+KG0NKm8XLIsrysBAMQhQixct2QbITZSWiW3kk8+Z4YT04kFUJe0NDvAVlV5XQkAIA45urDTokWL9MEHH+jTTz/Vp59+qs2bN0uSLH5Ti4PUhliGEzvO5/MpIyXD2eHEzIkF8ENpafa2vNyeHwsAgIscDbH33nuv3njjDSdPiRi0ZOsSdc7orHZp7bwuJSZlpGTQiQUQWaEQW1YmtWnjaSkAgPjjaIgdOnSo+vfvr5NOOkknnXSSunfvrkpWL8RBqqqrtHLHSo3qNcrrUmJWZmqmMyGWObEA6nJwJxYAAJc5GmJ/85vfOHk6xKBVO1Zpf81+5sNGUEZKhtYWr235iejEAqgLIRYA4CEWdoKrWNQp8hwfTsycWAA/RIgFAHiIEAtXLd22VBKLOkVSZkqmyvaXKVgTbNmJ6MQCqAshFgDgIUIsXPXNrm+U6EvUUW2P8rqUmFX7rNiWdmOZEwugLoRYAICHHJ0T65QtW7Zo8uTJzT6+JccistYWr1W3Nt3kTzDy1osJoRBbUlnSshWgGU4MoC6EWACAh+jEwjWWZWlt8Vq6sBHWOtnunLa4E8twYgB1IcQCADxkZDssOzubbmoM2lG2Q6X7S3VUG0JsJNV2YqtKWnai0HBiOrEAfogQCwDwEJ1YuCb02Bc6sZHl2JzYkhIpEJASEx2oCkBMIcQCADxEiIVrCLHuaJ3i4HBiurAAwgmF2LIyb+sAAMQlQixcQ4h1h6OdWObDAgiHTiwAwEOEWLiGEOuOg1cnbpF9+wixAMIjxAIAPOTowk5vv/227r333trvq6qqJElDhgypfW3SpEk677zznLwsosS63evUJrWN2qa19bqUmObo6sRHHulARQBiDiEWAOAhR0Psjh07tHDhwsNeP/i1HTt2OHlJRJF1xevUo00Pr8uIeY4OJ2ZOLIBwAgF7y5xYAIAHHB1OPHHiRFmWVe/XxIkTnbwkokjp/tLagIXIceQRO8GgVFHBcGIA4aWn29vSUm/rAADEJebEwjXVNdVK8HHLRZojqxOHnhFLiAUQTkqK/fgtQiwAwAMkCrimxqpRYgLPHI00f4Jfaf60loXYXbvsbZs2jtQEIMb4fHY3NvQLLwAAXESIhWuqLTqxbslIyWhZiN240d527uxMQQBiT6tWdGIBAJ4gUcA1NVYNIdYlGSkZLZsTGwqxXbs6UxCA2JOeTogFAHiCRAHX1Fg1SvQxnNgNrVNat6wTu2GDve3SxZmCAMSeVq0YTgwA8AQhFq5hYSf3NGk48b332gs4PfrogdfoxAJoCJ1YAIBHHH1OrFOCwaCKiorC7gsEAgqEnk+HqMLCTu7JSMlQSWWJLMuSz+er/8133WVvf/EL6cYbpeRkuxObnCx16BDxWgFEKTqxAACPGBliCwsLlZ+fH3Zfbm6u8vLyXK4ITmBhJ/dkpGRof81+VVZXKtWfWvcby8sP/X7zZqlHD7sT26WLlMB/LwB1SE+Xysqkmho+KwAArjIyxObk5KigoCDsPrqw0cmyLEliTqxLWicfeFZsvSE2NPc1ZONGO8Ru2CANHBjBCgFEvVat7G15uR1oAQBwiZEh1u/3Kysry+sy4KAaq0aS6MS6JCMlQ5JUUlmiI9KPqPuN33136PcbNkglJdKePcyHBVC/UHDdt48QCwBwFYkCrqi2qiURYt0SCrENLu4UCrH33WdvN2w4sKgTKxMDqE+oE8viTgAAl5Eo4IpQJ5aFndxx8HDieoVC7PDh9nbDhgNDjOnEAqjPwZ1YAABcRIiFK6pr6MS6qcmd2BNPlFJS6MQCaLxQJ5YQCwBwmZFzYhF7ajuxLOzkito5sVUl9b/xu++kjh2lQMBe0GnZMqlzZ3sfIRZAfVrbIz4IsQAAt9EWgytY2MldTerE9uhh//mii6RNm6SnnpIGDZKOOSayRQKIbhn254z27PG2DgBA3CFRwBUs7OSu1imNmBNbUyNt3y516mR/f/31kt9vf//GG1IiXXMA9QiF2L0N/LIMAACHMZwYrmA4sbsOfsROnUpK7CDbtq39fefO0ty5Unb2gSHFAFCXzEx7S4gFALiMEAtXsLCTuxo1nLi42N6GQqwknXxyBKsCEFMYTgwA8AiJAq7gETvuqn3ETlU9IXb3bnt7cIgFgMZiODEAwCOEWLiChZ3clZ6cLp98Te/EAkBjMZwYAOAREgVcwcJO7krwJah1Suv658SGQmybNq7UBCDGhJ4Ty3BiAIDLSBRwBQs7uS8jJYNOLIDISUy0gyydWACAy4xc2CkYDKqoqCjsvkAgoEAg4HJFaCkWdnJf6+TWhFgAkZWRQYgFALjOyBBbWFio/Pz8sPtyc3OVl5fnckVoKRZ2cl9GSoY2l2yu+w0s7ASgpTIzGU4MAHCdkSE2JydHBQUFYffRhY1OzIl1X0ZKhlYVrar7DXRiAbRURoa0caPXVQAA4oyRIdbv9ysrK8vrMuAgVid2X2hhJ8uy5PP5Dn8DCzsBaKl27aSlS72uAgAQZ0gUcAULO7kvIyVDliyV7i8N/4biYiktTUpJcbcwALGjXTupvNz+AgDAJYRYuIKFndyXkZwhSXUv7rR7N0OJAbRMu3b2NjSyAwAAF5Ao4AoWdnJf65TWkhoIsQwlBtASoRC7a5e3dQAA4gohFq5gYSf3ZaQ00Indu9delAUAmosQCwDwAIkCrmBhJ/eFQmxJZUn4N+zZYz8eAwCaixALAPAAiQKuYGEn99XbiQ0GpdJSQiyAliHEAgA8QIiFK1jYyX2tk+uZE1vyfXeW4cQAWoIQCwDwAIkCrmBhJ/fVDieuCjOceM8ee0snFkBLhELszp3e1gEAiCuEWLiChZ3cV+9w4lCIpRMLoCWOOMLe7tjhbR0AgLhCooArWNjJffU+Ymfv96/RiQXQEpmZUnKytG2b15UAAOIIiQKuYGEn9zWqE0uIBdASPp/djd2+3etKAABxxO91AeEEg0EVFRWF3RcIBBQIBFyuCC3Fwk7ua9ScWIYTA2ipI46gEwsAcJWRIbawsFD5+flh9+Xm5iovL8/litBSLOzkvpTEFCUlJDGcGEBkdeworVolWZbdmQUAIMKMDLE5OTkqKCgIu48ubHRiYSf3+Xw+tU5pzXBiAJHVsaNUXi7t2ye1bu11NQCAOGBkiPX7/crKyvK6DDiIhZ28kZGSoZLKMMOJQ51YhhMDaKnQCsXbthFiAQCuIFHAFSzs5I2MlAw6sQAiKzvb3m7Z4m0dAIC4QYiFK1jYyRutkxsYTkwnFkBLde5sbzdu9LYOAEDcIFHAFSzs5I06O7HFxVJKipSW5n5RAGJLly72lhALAHAJIRauYGEnb2SkZKg8WK5gTfDQHcXFUtu23hQFILaEOrGbNnlbBwAgbpAo4ArmxHqj9lmxP1zciRALwCmdOkmJiXRiAQCuIcTCFaxO7I3WyfZKoYcNKSbEAnBKYqKUk0OIBQC4hkQBV7CwkzdqO7FVP+jE7t5NiAXgnO7dpXXrvK4CABAnSBRwBQs7eSMUYg/pxJaXS5WVhFgAzunVyx7hsWuX15UAAOIAIRauYGEnb7ROCTOcuLjY3hJiATilVy97++233tYBAIgLJAq4goWdvBF2YSdCLACnhULsN994WwcAIC4QYuEKFnbyRtjhxIRYAE4LhdivvvK2DgBAXCBRwBUs7OSNsKsTh0JsmzbuFwQgNvXtK/l80vLlXlcCAIgDfq8LCCcYDKqoqCjsvkAgoEAg4HJFaCkWdvIGnVgArggE7G7s0qVeVwIAiANGhtjCwkLl5+eH3Zebm6u8vDyXK0JLsbCTN8I+YocQCyASjj9emjrVXgE9Lc3ragAAMczIEJuTk6OCgoKw++jCRicWdvJG2E5sqFPSo4cHFQGIWQMGSK+9Ji1ZIg0Z4nU1AIAYZmSI9fv9ysrK8roMOIiFnbwR9hE7K1dKHTtKXbp4VBWAmBQKrvPnE2IBABFFooArWNjJG/4Ev9L8aYcOJ961S+rQwbuiAMSmU06xF3eaP9/rSgAAMY5EAVewsJN3MlIyDu3E7twptWvnXUEAYlNmpj0vdtYsqabG62oAADGMEAtXsLCTd1qntD4QYmtq7IWd2rf3tigAsWnUKGn7dnteLAAAEUKigCtY2Mk7h3Rit22zg2ynTt4WBSA2jR5tb996y9s6AAAxjRALV7Cwk3cyUjJUUvn9nNjvvrO33bt7VQ6AWHb66fac+xdflCzL62oAADHK0URRXl6uu+66S3369FFqaqpycnJ09dVXa/PmzU5eBlGIhZ28E+rEWpYlrV9vv9itm7dFAYhNfr/0P/8jrV4tzZvndTUAgBjlWKKoqKjQmWeeqXvvvVf79u3TmDFj1KVLFz3zzDMaOHCg1q5d69SlEIVY2Mk7rZNba3/NflVWVxJiAUTeTTfZ2wcf9LYOAEDMcizE3nfffVqwYIGGDh2qr776Sq+88ooWLlyoKVOmaMeOHbr66quduhSiEAs7eScjJUOS7CHFX39tv9ijh4cVAYhpfftKF14ovfmmNHOm19UAAGKQI4miqqpKjz32mCTp8ccfV6tWrWr33Xbbberfv79mzZqlRYsWOXE5RCEWdvJOKMTurdwrzZ0rde4sdezocVUAYtof/yilpUkTJ0pbtnhdDQAgxjgSYufOnas9e/aoZ8+eGjhw4GH7L7roIknStGnTnLgcohALO3mndXJrSVLV0i/seWrnnedxRQBiXs+e0mOP2VMYzjhD+vZbrysCAMQQRxLFku+fBzdo0KCw+0OvL1261InLIQqxsJN3Qp3Ytk/8w37hf/7Hw2oAxI2rr7aD7FdfSccdJ/3mN9KqVaxaDABoMb8TJ9mwYYMkqXPnzmH3h15fH1pUJop98a8pqty5LfzO+v7HXMcuX73HNO9/9M0+Zz27ml+nve/I7z7RldukVq+8LqW1bUQtzdjn9PliZN+wVbP09nSp0zfvavspx2nOETulVa/VfS4AcMqZ2cp69vc64Z7/VcaDD0oPPqjKNq21+5gequjYThVZbbW/dUA1yUmqTklSTVKSapKTZCX47ON99tb6flsr3Ou+Q/cBQLw7+sLr1fbInl6XERGOhNh9+/ZJkgKBQNj96enpkqSSkpJGnW/Lli2aPHlys+tpybENSb1jkgZuKI/Y+WPVSaE/vH6rl2XEpUGSaiS93le67vTl2vmfi7wuCUCc8V0pnbVOylsjnbK5RMd9vlQd93tdFQDEtuU9jibEwlY6OV9zd2yt+w31/Ab4sN8kH3JcQ1eOwHnrrbWhcuo9cdhXO6R30DFHHNu8czb0m/XmHhuJfYbVUx1I0wf+9apIrtFf6z4DALhiraS1liV/WYXSduxWUmm5Eqr2K7EqqISq/UrYv18+SwdGmHy/sUcEhXvtwHt9jFQGgFpHDzjN6xIixpEQG1qNuKysLOz+0tJSSVLr1q0bdb7s7OyIdlNb4sSr8r0uAWiSREmjNcTrMgAAAABHOLLKTteuXSVJmzZtCrs/9Hq3bt2cuBwAAAAAIE45EmIHDBggSVq8eHHY/aHX+/fv78TlAAAAAABxypEQO2zYMGVmZurbb7/Vl19+edj+V199VZKUl5fnxOUAAAAAAHHKkRCbnJysW265RZJ08803186BlaRHHnlES5cu1YgRIzR48GAnLgcAAAAAiFOOrU5855136sMPP9S8efPUu3dvDR8+XOvXr9fChQvVoUMH/eMf/3DqUgAAAACAOOVIJ1aSUlNT9fHHH2vSpEkKBAKaOnWq1q9fr4kTJ2rx4sU66qijnLoUAAAAACBO+SzLMuapaoMHD9bixYs1aNAgLVq0yOtyAAAAAACGcawTCwAAAABApBFiAQAAAABRI+pD7OTJkzV58mSvy0CU4H5BU3C/oKm4Z9AU3C9oCu4XNEWs3y9RH2IBAAAAAPGDEAsAAAAAiBqOPSfWScFgUEVFRWH3BQIBBQIBlysCAAAAAJjAyBBbWFio/Pz8sPtyc3OVl5fnckUAAAAAABMYGWJzcnJUUFAQdh9dWAAAAACIX0aGWL/fr6ysLK/LAAAAAAAYhoWdAAAAAABRgxALAAAAAIgahNgmKisr07Rp01RWVsb1ovB6bov1f89Yv57bYv3f04v/ftwzXM/k67kt1v89Y/16bov1f89Yv57b3P77EWKbqKysTG+99ZarNzzXi16x/u8Z69dzW6z/e3rx3497huuZfD23xfq/Z6xfz22x/u8Z69dzm9t/P0IsAAAAACBqEGIBAAAAAFGDEAsAAAAAiBqEWAAAAABA1CDEAgAAAACihs+yLMvrIkLatWun4uJipaWl6ZhjjmnUMVu2bJEkZWdnR7K0WsFgUIWFhcrJyZHf7+d6UXY97heu1xTcL9F9PS+uyT3D9ZqC+4XrNQX3C9drimi5X/r27at//etfTb6eUSE2EAiovLzc6zIAAAAAABE2aNAgLVq0qMnHufOr8EY64ogjtH37dqWmpqpHjx5elwMAAAAAiJC+ffs26zijOrEAAAAAANSHhZ0AAAAAAFGDEAsAAAAAiBqEWAAAAABA1CDEAgAAAACiRtSG2PLyct11113q06ePUlNTlZOTo6uvvlqbN2/2ujS00BlnnCGfz1fn1/Tp08Me15x7wq1j0HKLFi3SAw88oHHjxqlz586190NDTL8vuJ8iozn3i+mfPS05DnUrKyvT1KlTdc011+joo49Wamqq0tPTNWDAAN1zzz3at29fncfy+RKfmnvP8BkTvx555BGNGzdOvXv3VmZmplJSUtStWzeNHz9ey5Ytq/M4PmPqYUWh8vJya8iQIZYkKzs727rkkkusk08+2ZJkdejQwfr222+9LhEtMGLECEuSdeGFF1oTJkw47Gvp0qWHHdOce8KtY+CMMWPGWJIO+6qP6fcF91PkNOd+MfmzpyXHoX5PPfVU7f1xzDHHWBdffLE1atQoq3Xr1pYkq2/fvta2bdsOO47Pl/jV3HuGz5j41b59eys1NdU6+eSTrbFjx1pjx461+vTpY0mykpKSrGnTph12DJ8x9YvKEJufn29JsoYOHWqVlJTUvj5lyhRLkjVixAjvikOLhT7k161b1+hjmnNPuHUMnPHAAw9YkyZNst58801ry5YtVkpKSoOhxPT7gvspcppzv5j82dOS41C/Z5991vrpT39qrVy58pDXCwsLrYEDB1qSrMsuu+yw4/h8iV/NvWf4jIlfc+bMscrLyw97/fHHH7ckWR07drT2799/yD4+Y+oXdSG2srLSyszMtCRZixcvPmx///79LUnW559/7kF1cEJTP+Sbc0+4dQwip6FQYvp9wf3krkiEWO6X2Ddv3jxLkpWSkmJVVlbWvs7nC+pS1z1jWXzGILyePXtakqwlS5bUvsZnTMOibk7s3LlztWfPHvXs2VMDBw48bP9FF10kSZo2bZrbpcEjzbkn3DoG3jH9vuB+in7cL7FvwIABkqTKykrt3Lmz9nU+X1CXuu6Z5uCeiQ9JSUmSpOTk5NrX+IxpmD9iZ46QJUuWSJIGDRoUdn/o9aVLl7pWEyLj6aef1s6dO5WQkKA+ffroggsuUNeuXQ97X3PuCbeOgXdMvy+4n8xl2mdPS45Dy6xdu1aS/UNmu3btal/n8wV1qeueORifMQh5/vnntWbNGvXu3Vu9e/eufZ3PmIZFXYjdsGGDJKlz585h94deX79+vWs1ITLuu+++Q76//fbbNWnSJE2aNOmQ15tzT7h1DLxj+n3B/WQu0z57WnIcWubRRx+VJI0ePVopKSm1r/P5grrUdc8cjM+Y+PXQQw9pxYoVKi0t1apVq7RixQrl5OTopZdeUmJiYu37+IxpWNQNJw4tWx4IBMLuT09PlySVlJS4VhOcdfrpp+v555/Xt99+q7KyMq1Zs0YFBQXy+/266667av8HEdKce8KtY+Ad0+8L7ifzmPrZ05Lj0HzvvPOOnn76aSUlJenee+89ZB+fLwinvntG4jMG0nvvvafnnntOr776qlasWKFu3brppZde0uDBgw95H58xDYu6EIvYd8899+jKK6/UUUcdpbS0NPXp00d33HGHpk6dKkmaPHmyysvLvS0SQMzhswchq1ev1pVXXinLsvTQQw/VznME6tKYe4bPGHz44YeyLEvFxcX65JNP1Lt3b40YMUIFBQVelxZ1oi7EtmrVSpL9oOlwSktLJUmtW7d2rSa4Y+TIkTrxxBO1e/duLVy4sPb15twTbh0D75h+X3A/RQ+vP3tachyabvPmzRo9erSKi4t122236dZbbz3sPXy+4GCNuWfqw2dM/GnTpo2GDx+ud955R4MHD9akSZP02Wef1e7nM6ZhURdiQxPfN23aFHZ/6PVu3bq5VhPcE5r0vmXLltrXmnNPuHUMvGP6fcH9FF28/OxpyXFoml27dmnkyJFav369rrrqKj388MNh38fnC0Iae880hM+Y+JSUlKRLL71UlmUdspIvnzENi7oQGxqesXjx4rD7Q6/379/ftZrgnuLiYkkHxtpLzbsn3DoG3jH9vuB+ii5efva05Dg03r59+3Tuuedq5cqVGjdunJ566in5fL6w7+XzBVLT7pmG8BkTv7KysiRJO3bsqH2Nz5hGiNgTaCPk4IfrfvHFF4ft52HMsWv79u1Wenq6JcnauHFj7evNuSfcOgaRk5KSYtX3EWb6fcH95K6G7pf6eP3Z05Lj0DgVFRXWmWeeaUmyRo0aZVVWVtb7fj5f0NR7pj58xsS3CRMmWJKshx56qPY1PmMaFnUh1rIsKz8/35JknXrqqda+fftqX58yZYolyRoxYoR3xaFF5s6da73++utWMBg85PV169ZZw4YNsyRZ559//mHHNeeecOsYREZjQonp9wX3k3saul9M/+xpyXGoXzAYtMaOHWtJsoYPH26VlpY26jg+X+JXc+4ZPmPi15w5c6x3333Xqq6uPuT1qqoq6y9/+YuVkJBgpaWlWRs2bDhkP58x9YvKEFteXm6dcsopliQrOzvbuuSSS2q/79Chg/Xtt996XSKa6ZlnnrEkWZ06dbJ+/OMfW5dffrk1bNgwKzU11ZJk9evXz9q2bdthxzXnnnDrGDjjrbfesk455ZTaL5/PZ0k65LW33nrrkGNMvy+4nyKnqfeL6Z89LTkO9fvzn/9sSbIkWWPHjrUmTJgQ9mvHjh2HHMfnS/xqzj3DZ0z8Cv23z8rKskaNGmVdfvnl1siRI63s7GxLkpWammq98sorhx3HZ0z9ojLEWpZllZWVWZMmTbJ69uxpJScnW506dbImTpx4yDAMRJ+VK1daN954ozVo0CCrQ4cOlt/vtzIzM60hQ4ZYU6ZMscrKyuo8tjn3hFvHoOVC/xOo7+uZZ5457DjT7wvup8ho6v0SDZ89LTkOdbv77rsbvFckWevWrTvsWD5f4lNz7hk+Y+LX2rVrrTvuuMMaNmyYlZ2dbSUlJVnp6elWv379rJ/97GfW119/XeexfMbUzWdZliUAAAAAAKJA1K1ODAAAAACIX4RYAAAAAEDUIMQCAAAAAKIGIRYAAAAAEDUIsQAAAACAqEGIBQAAAABEDUIsAAAAACBqEGIBAAAAAFGDEAsAAAAAiBqEWAAAAABA1CDEAgAAAACiBiEWAAAAABA1CLEAAAAAgKhBiAUAAAAARA1CLAAAAAAgahBiAQAAAABRgxALAAAAAIgahFgAAAAAQNT4/3n0DcP62LYSAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(wavelength, VISDA.sensitivity(wavelength), c='g')\n",
    "plt.plot(wavelength, NIRDA.sensitivity(wavelength), c='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d73679c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Area of the moon on the sky\n",
    "moon_area = 0.2*u.deg**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "13e5ed29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pointed at the moon, each NIRDA pixel recieves 119635894.68893522 electron / (s pix2)\n",
      "Off-pointed, with 1e-8 scattered light, NIRDA pixel recieves 1.1963589468893523 electron / (s pix2)\n"
     ]
    }
   ],
   "source": [
    "moon_brightness = np.trapz(NIRDA.sensitivity(wavelength) * moon_spectrum, wavelength) \n",
    "moon_brightness = NIRDA.apply_gain(moon_brightness)\n",
    "\n",
    "moon_brightness_per_area = moon_brightness / moon_area\n",
    "moon_brightness_per_pixel = (moon_brightness_per_area * NIRDA.pixel_scale**2).to((u.electron/u.second) * (1/u.pixel**2))\n",
    "\n",
    "print(\"Pointed at the moon, each NIRDA pixel recieves\", moon_brightness_per_pixel)\n",
    "print(\"Off-pointed, with 1e-8 scattered light, NIRDA pixel recieves\", 1e-8 * moon_brightness_per_pixel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "49dfe388",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pointed at the moon, each VISDA pixel recieves 31530361.13844697 electron / (s pix2)\n",
      "Off-pointed, with 1e-8 scattered light, VISDA pixel recieves 0.3153036113844697 electron / (s pix2)\n"
     ]
    }
   ],
   "source": [
    "moon_brightness = np.trapz(VISDA.sensitivity(wavelength) * moon_spectrum, wavelength) \n",
    "# Apply VISDA gain\n",
    "moon_brightness = 0.52 * u.electron/u.DN * moon_brightness\n",
    "\n",
    "moon_brightness_per_area = moon_brightness / moon_area\n",
    "moon_brightness_per_pixel = (moon_brightness_per_area * VISDA.pixel_scale**2).to((u.electron/u.second) * (1/u.pixel**2))\n",
    "\n",
    "print(\"Pointed at the moon, each VISDA pixel recieves\", moon_brightness_per_pixel)\n",
    "print(\"Off-pointed, with 1e-8 scattered light, VISDA pixel recieves\", 1e-8 * moon_brightness_per_pixel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20c38bf5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}