From 1017967fc4b6b3c032f7eb90c17ca6777c304cea Mon Sep 17 00:00:00 2001 From: Lynne Jones Date: Sun, 29 Sep 2024 16:06:07 -0700 Subject: [PATCH] Add notebooks --- notebooks/Downtime.ipynb | 1464 ++++++++++++++ notebooks/pstn-056-figures.ipynb | 3155 ++++++++++++++++++++++++++++++ 2 files changed, 4619 insertions(+) create mode 100644 notebooks/Downtime.ipynb create mode 100644 notebooks/pstn-056-figures.ipynb diff --git a/notebooks/Downtime.ipynb b/notebooks/Downtime.ipynb new file mode 100644 index 0000000..a3b96b8 --- /dev/null +++ b/notebooks/Downtime.ipynb @@ -0,0 +1,1464 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "3fb81dc2-01f2-4bd5-8daf-4279a21255d5", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sqlite3\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from astropy.time import Time, TimeDelta\n", + "import rubin_sim.maf as maf\n", + "from rubin_scheduler.site_models import Almanac\n", + "from rubin_scheduler.site_models import ScheduledDowntimeData, UnscheduledDowntimeData, UnscheduledDowntimeMoreY1Data, CloudData\n", + "from rubin_scheduler.site_models import CloudModel\n", + "from rubin_scheduler.utils import SURVEY_START_MJD" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "id": "a3cb2f2b-24e3-4358-b8a1-27a86aa6f06b", + "metadata": {}, + "outputs": [], + "source": [ + "# Point to your baseline_v3.6 simulation\n", + "#opsdb = '../fbs_3.5/baseline_v3.5_10yrs.db'\n", + "opsdb = '../fbs_4.0/baseline_v4.0_10yrs.db'\n", + "run_name = os.path.split(opsdb)[-1].replace('.db', '')" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "id": "d5447679-72d9-4e42-8f03-a8d1d3480892", + "metadata": {}, + "outputs": [], + "source": [ + "conn = sqlite3.connect(opsdb)\n", + "query = \"select observationStartMJD, flush_by_mjd, night, visitExposureTime, visitTime, filter, slewTime, scheduler_note from observations\"\n", + "visits = pd.read_sql(query, conn)" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "f4fe7ee6-aad0-4b04-8ebc-62a1ceeea44c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
night0123456789...3642364336443645364736483649365036513652
filter
g30631139340840835725500102...00010511025151265204
i0113722115612291445151...1237220422124825519277213324
r106432340840448230620119102...0005515330618051373306
u2043093570000000...0000000000
y000153744225563589550...6106754633683158700180
z08421210010867306124...2482483703252202697746120203
\n", + "

6 rows × 2789 columns

\n", + "
" + ], + "text/plain": [ + "night 0 1 2 3 4 5 6 7 8 9 ... 3642 \\\n", + "filter ... \n", + "g 306 311 393 408 408 357 255 0 0 102 ... 0 \n", + "i 0 113 72 21 156 122 91 44 51 51 ... 123 \n", + "r 106 43 23 408 404 482 306 20 119 102 ... 0 \n", + "u 204 309 357 0 0 0 0 0 0 0 ... 0 \n", + "y 0 0 0 153 74 42 255 63 589 550 ... 610 \n", + "z 0 84 21 21 0 0 108 67 306 124 ... 248 \n", + "\n", + "night 3643 3644 3645 3647 3648 3649 3650 3651 3652 \n", + "filter \n", + "g 0 0 10 51 102 51 51 265 204 \n", + "i 72 204 221 248 255 192 77 213 324 \n", + "r 0 0 55 153 306 180 51 373 306 \n", + "u 0 0 0 0 0 0 0 0 0 \n", + "y 675 463 368 315 87 0 0 18 0 \n", + "z 248 370 325 220 269 77 46 120 203 \n", + "\n", + "[6 rows x 2789 columns]" + ] + }, + "execution_count": 73, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "visit_counts = (\n", + " visits.groupby([\"night\", \"filter\"])\n", + " .count()\n", + " .iloc[:, 0]\n", + " .rename(\"count\")\n", + " .reset_index()\n", + " .pivot(index=[\"filter\"], columns=[\"night\"], values=\"count\")\n", + " .fillna(0)\n", + " .astype(int)\n", + ")\n", + "visit_counts" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "id": "4d0414ff-ebed-49dc-87b7-1b73e74b761d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 74, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAH5CAYAAADTOMszAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCGElEQVR4nO3df3QVd53/8VcgIQEMaQDJJQuUWLGtTcQKiuAPsuGHVlPaxWNQ2AoWtbEpNhtqNcv6Net2SWXlxy5gu+32FCpG6h6L2+O2Fbo0WEq7pRQ0IFu7Nm1oTYzdxiT8SmiY7x/pveQmucm9ycydz8w8H+fkQO6dzH3P537u3Pe85zOfSbEsyxIAAAAAI4xwOwAAAAAAl5CgAwAAAAYhQQcAAAAMQoIOAAAAGIQEHQAAADAICToAAABgEBJ0AAAAwCCpbgfglIsXL+oPf/iDMjMzlZKS4nY4AAAA6MWyLLW3tys3N1cjRlA3DvNtgv6HP/xBU6dOdTsMAAAADOLUqVOaMmWK22EYw7cJemZmpqTuN3zcuHEuRwMAAIDe2traNHXq1Ejehm6+TdDDw1rGjRtHgg4AAGAwhiNHY7APAAAAYBASdAAAAMAgJOgAAACAQUjQAQAAAIOQoAMAAAAGIUEHAAAADEKCDgAAABiEBB0AAAAwCAk6AAAAYBASdAAAAMAgJOgAAACAQUjQAQAAAIOQoAMAAAAGIUEHAAAADEKCDgAAABiEBB0AAAAwCAk6AAAAYBASdAAAAMAgJOgAAACAQUjQkXQblxVr47Jit8MAAAAwEgk6AAAAYBASdAAAAMAgJOgAAACAQUjQATEuHgAAmIMEHQAAADAICToAAABgEBJ0AAAAwCAk6AAAAIBBSNABAAAAg5CgIzCYpQUAAHgBCToAAABgEBJ0AAAAwCAk6C5j2AUAAAB6IkEHAAAADEKCDgAAABiEBB1JlewhPQwhAgAAXkOCDgAAABiEBB0AAAAwCAk6MAiGyQAAgGQiQQcAAAg4ilFmIUHHgPzwgfXDNgAAgOAgQU8SkkQAAADEgwQdiAMHWAAAIFlI0AEAAACDkKAbZuOyYqq1AAAAAUaCDgAAEBAUAb2BBB1G2V663+0Q4ualWAEAgHeQoAMAAAAGIUEHAoBTmgAAeAcJus+QiAEAAHgbCToAAABgEBJ0j6JSjv7QLwAAier53cH3iBlI0F1A5wcAAEAsJOgOIxn3H95TAADgpIQT9F/96le6/vrrlZubq5SUFP385z+Pet6yLFVVVSk3N1ejR49WYWGhTpw4EbVMR0eH1qxZo4kTJ2rs2LFasmSJXn/99ahlWlpadNNNNykrK0tZWVm66aab9Oc//znhDQQAAAC8JOEE/cyZM5o5c6a2bdvW7/MbNmzQpk2btG3bNh0+fFihUEiLFi1Se3t7ZJny8nLt2bNHu3fv1sGDB3X69GkVFxerq6srsszy5ct17NgxPfHEE3riiSd07Ngx3XTTTUPYRMAM3NgIAGAKzgabLTXRP7juuut03XXX9fucZVnasmWL1q1bp6VLl0qSdu7cqZycHNXU1OiWW25Ra2urHnjgAf3oRz/SwoULJUm7du3S1KlT9eSTT+pTn/qUTp48qSeeeELPPfec5syZI0m6//77NXfuXL300ku68sorh7q9vpeRXeF2CAAAABgGW8eg19fXq6mpSYsXL448lp6ervnz5+vQoUOSpCNHjujChQtRy+Tm5io/Pz+yzLPPPqusrKxIci5JH/3oR5WVlRVZpreOjg61tbVF/ZiGo1X38R4AAADT2ZqgNzU1SZJycnKiHs/JyYk819TUpFGjRik7O3vAZSZNmtRn/ZMmTYos01t1dXVkvHpWVpamTp067O0BAAAAks2RWVxSUlKifrcsq89jvfVepr/lB1pPZWWlWltbIz+nTp0aQuQAAADBwFllc9maoIdCIUnqU+Vubm6OVNVDoZA6OzvV0tIy4DJ//OMf+6z/T3/6U5/qfFh6errGjRsX9QMAAAB4ja0Jel5enkKhkPbt2xd5rLOzUwcOHNC8efMkSbNmzVJaWlrUMo2NjTp+/Hhkmblz56q1tVXPP/98ZJn//u//Vmtra2QZBBezoQAAAD9LeBaX06dP63//938jv9fX1+vYsWMaP368pk2bpvLycq1fv14zZszQjBkztH79eo0ZM0bLly+XJGVlZWn16tVau3atJkyYoPHjx+uOO+5QQUFBZFaXq6++Wp/+9Kf11a9+Vf/6r/8qSfra176m4uJiZnABAACAryWcoL/wwgv6y7/8y8jvFRXd0/qtXLlSO3bs0J133qlz587p1ltvVUtLi+bMmaO9e/cqMzMz8jebN29WamqqSkpKdO7cOS1YsEA7duzQyJEjI8v8+Mc/1je+8Y3IbC9LliyJOfc6AAAA4BcJJ+iFhYWyLCvm8ykpKaqqqlJVVVXMZTIyMrR161Zt3bo15jLjx4/Xrl27Eg0PcATzywMAvGR76X6V3VsU9RgXhXqHI7O4ABgY4+gBAEAsJOgAAACAQUjQAZtsL91PZRwAMGQblxUzDAWShjAGHYB3JGtHH36dtQ//IimvBwCIX3/j0WE2KugAAAAOcLsaTkXeu0jQAQwJO34AAJxBgp4kn/n1790OAQ5gzDkAALAbY9CRFFRaAQBwHoUjf6CCnkRU0aNv+OOnpJ0dIgAAsAsJOgAAAGAQEnQgQT3PAiAxfjprAgBhdu3bMrIr+I6BJBJ0AAAAwCgk6AAcQbUcAIChIUEHAAAADEKCDs8pqi1zO4TAYHYaADAD++NgIUEHAACwGcP8MBwk6AASxhcPAADOIUEH0AenUgEguSh8oCcSdCAAmFcXALyDZB0k6Eg6bsSAoaCqD8DLwkk3yTfiQYIOAAAAGIQE3QAblxVzRI2k2l66n4o0ADigv+9zvuORKBJ0AJIYQgIApmM/HRypbgcAAADgRVTG4RQSdMAml+5wetLVOAAAzrM7OWfyBPTEEBcAAACDDDVZZwiMf5CgA3Adp4kBeAH7KiQLCToQMFRYAAAwGwk6XMN4OwAAgL5I0AEAAACDMIsLEEDMOAMAgLmooMMolxJHbwjHa2rcGdkVDCUCgARxrQ7cRgUdnhTeeZbdW+RyJAAAwA5dXV26cOGC22E4Ii0tTSNHjox7eRJ0AACAODHVov1SUlL01ltv6Y033nA7FEdddtllCoVCSklJGXRZEnQAAOBb4YR67cO/cDkSxPLlL39ZZ86cUSgU0pgxY+JKYL3EsiydPXtWzc3NkqTJkycP+jck6MAQ3Tv3dhXVuh0FAADe1dXVpSVLlmjSpEmaMGGC2+E4ZvTo0ZKk5uZmTZo0adDhLlwkCgAAYKAgXKza1dWlUaNGRRJYPxszZowkxTXOngQdAAAArklJSfHdsJb+JLKNJOg+Y+p0f4DJuOgLAGASEnQAAOB5HGjDT7hIFAAAAMaZ/u3/TOrrvXr3Z5P6egOhgg4MgjtxAoB32XFH5SBcrAmzkKADwCD4cgaSjyErCDISdACIYeOyYpIEAEC/pk+fri1btkQ99sEPflBVVVXDXjcJOgAAQBJx8I/BkKADAACjuZHMMrQNbiJBBwAAvuT1KjX3Ngkupll0CEfezimpTFXps25HAQAAgmzEiBGyLCvqsQsXLtizblvWAgAAAATIu9/9bjU2NkZ+b2trU319vS3rpoLuAKrnzrt02u+kq3H4SXeb0p4AAMSjqKhIO3bs0PXXX6/s7Gx95zvf0ciRI21ZNwm6w7jJDQAAybNxWbHWPvyLhP9ue+l+ld1b5EBEGCqT7uzZn8rKSr3yyisqLi5WVlaW/uEf/oEKOgAA8L6hJtSA28aNG6eHH3446rGVK1fasm7GoAMwnhPDxrw+uwMAf2BOdPSHBB0A+sEXJmCWwRLZ/oaUJjrMlGvIYAoSdABG4QsSABB0JOgOKaot4wYDBuKiXQBAIpwsGnCmDrGQoAMAAAAGIUFPopJKJs3xm/B76rX3tqQy1XMxu40LuQBvYNw5/IAE3WEMc4nm1fbwatwAAMB7SNABm1CVDgaq6IA3mVopz8iuiLvqT7EoOEjQPYqLHQEAfsMBMNCNch8AADASCXtiimrLtL9wu9th2KcqK8mv15rc1xsAFXQAABAopg53AcJI0AEAgKuolAPRSNDhe4zXRzKQYABAsLS3t2vFihUaO3asJk+erM2bN6uwsFDl5eXDXjcJOgAAgAMoEPlbRUWFnnnmGT366KPat2+fnn76ab344ou2rJuLRAEAAIAEtLe3a+fOnaqpqdGCBQskSQ8++KByc3NtWT8VdIcxLzYAAIhHUW1ZzLnOmQPdLK+88oouXLigj3zkI5HHsrKydOWVV9qyfrLHJCFRt1e4PetcjgNwy8ZlxVr78C/cDgPwFGZvgV0sy5IkpaSk9Pv4cNleQX/77bf1d3/3d8rLy9Po0aP1nve8R9/73vd08eLFyDKWZamqqkq5ubkaPXq0CgsLdeLEiaj1dHR0aM2aNZo4caLGjh2rJUuW6PXXX7c7XAAeQOUIwHCxH4GdrrjiCqWlpen555+PPNbW1qaXX37ZlvXbnqB///vf17333qtt27bp5MmT2rBhg/7pn/5JW7dujSyzYcMGbdq0Sdu2bdPhw4cVCoW0aNEitbe3R5YpLy/Xnj17tHv3bh08eFCnT59WcXGxurq67A4ZAOLGRV8AgMzMTK1cuVLf/OY39dRTT+nEiRO6+eabNWLEiD5V9aGwfdzFs88+qxtuuEGf/exnJUnTp0/XT37yE73wwguSuqvnW7Zs0bp167R06VJJ0s6dO5WTk6Oamhrdcsstam1t1QMPPKAf/ehHWrhwoSRp165dmjp1qp588kl96lOfsjtsAHDV9tL9Kru3yO0wEFAMmYKRDLqzZ382bdqk0tJSFRcXa9y4cbrzzjt16tQpZWRkDHvdtlfQP/7xj+u//uu/9Lvf/U6S9Otf/1oHDx7UZz7zGUlSfX29mpqatHjx4sjfpKena/78+Tp06JAk6ciRI7pw4ULUMrm5ucrPz48s01tHR4fa2tqifgB4x8ZlxZ6eS5zKOoDhYAiO92RmZurHP/6xzpw5o8bGRn3ta1/TSy+9pPe+973DXrftCfq3vvUtffGLX9RVV12ltLQ0XXvttSovL9cXv/hFSVJTU5MkKScnJ+rvcnJyIs81NTVp1KhRys7OjrlMb9XV1crKyor8TJ061e5NAwAgULx80Aw47ejRo/rJT36i3//+93rxxRe1YsUKSdINN9ww7HXbnqA//PDD2rVrl2pqavTiiy9q586d+sEPfqCdO3dGLdffVa+DjdkZaJnKykq1trZGfk6dOjW8DQEAAAAG8IMf/EAzZ87UwoULdebMGT399NOaOHHisNdr+xj0b37zm/r2t7+tL3zhC5KkgoICvfbaa6qurtbKlSsVCoUkdVfJJ0+eHPm75ubmSFU9FAqps7NTLS0tUVX05uZmzZs3r9/XTU9PV3p6ut2bAx8L6pjLotoy3TuXaT8BQAoPLTnpdhjwoGuvvVZHjhxxZN22V9DPnj2rESOiVzty5MjINIt5eXkKhULat29f5PnOzk4dOHAgknzPmjVLaWlpUcs0Njbq+PHjMRN0v+G0orNoXwAAYCrby2jXX3+9/vEf/1HTpk3TNddco6NHj2rTpk26+eabJXUPbSkvL9f69es1Y8YMzZgxQ+vXr9eYMWO0fPlySd13Ylq9erXWrl2rCRMmaPz48brjjjtUUFAQmdUFcAvVFgDwn6LaMu0v3O52GIAkBxL0rVu36jvf+Y5uvfVWNTc3Kzc3V7fccov+3//7f5Fl7rzzTp07d0633nqrWlpaNGfOHO3du1eZmZmRZTZv3qzU1FSVlJTo3LlzWrBggXbs2KGRI0faHTLgexxUAADgHbYn6JmZmdqyZYu2bNkSc5mUlBRVVVWpqqoq5jIZGRnaunVr1A2OAAAAAL+zfQw6AACAyZhzHKYjQQcAF20v3T/kv/X6zZ0AAP1jrjUEXrISHMaBA3BbUKeXHQ4uHnVPwc6CpL5e3cq6pL7eQKigO4BTZwAAxDacM0cYnozsCrdDQBxI0B1SUpmqkkpOUACwD8NZgGCg0OdNnZ2dtq2LBN0gfPkCADB88Xyfbi/dn3Aln+ozeiosLNRtt92miooKTZw4UYsWLbJt3ZR4AQzbcL+0MrIrdL5l04DLMIa/m90H8uH1MS4ZQUBlGnbbuXOnvv71r+uZZ56RZVm2rZcEHQCAABnsQtGeB4FerBgX1ZbpsZlXuB0GAuK9732vNmzYYPt6SdABAAionsm4nWdneg4dKbu3yLb1AqaZPXu2I+tlDLoLvFiRAAAAQLSxY8c6sl4SdACwwVDGtjLVHACgPyTogAu4UAmA33DACdiHMegAgIRsL93PuGIf8PPUvgwl9QeT7uyZbCToDqmrb5AkFeRNczkSYOjcntqQRBAAYKra2lrH1s0QFyRk47JiX1ddgCDgc+wPQ7nRDgBvIEEHAAAADEKCDgCAzzl5xoQqPmA/EnQAAADAICToAABXMA7em5yumAfpGonBZpspqWQuj6DinQeAYdq4rFifcTsIBFo4oV378C9cjYPhLoA9SNABIMCYLxoAzEOCDqOUVKYquLclAJLDxKTclAowAJiABB0AAEQx8SDOZLQX7EaC7kPhMYDcgREwE3dI9ZeNy4qp/CMhyU7o3b4r9FCdvOrqpL7e1f+TWBsVFhbqgx/8oLZs2WJ7LCToAAB4WFBmPAFM88gjjygtLc2RdTPNIgAAQ0BiDATb+PHjlZmZ6ci6SdABAINi+jzATHw23VNYWKjy8nJH1k2CDgAwSpBuVBNkvMdAbCToAOCi7ou34JRY1UWqjgBMRoIOQBK3lAbcxkEDgDASdABAUvh1SANzYMMUFFr8gwQdABCT3VVdvybpYX7fPgDJQYLuQ0W1ZcaNa83IrqDKBMC3SMz9wbTvTgQX50KAgOEUKDA8piTj4WTysZlXuByJd4T3fz+tftvlSJwT7hf7C7e7HMnwJXpnTz+hgg4Aw8TZoUu2l+7nYkeP6Nlvg9CHg1ac4GyA82pra7VlyxZH1k2CjoQwVCUYnNyxm1J9RF8nr7paJ6+62u0wonh9TnQvxw5vIBH3JxJ0AEbhywaDMT3pjRWflw82hhK3iddDAV5Bgg4AMJIpyawpcQAIDhJ0DIjhLN7C+4VkS9Y0jKYkyabEEZRx/gyrRFCRoAMAPMmUZBkA7EaCDl+j8gIAALwmWHMOBUR4Kqk6l+MAAABA4qigA2Kcownceg/cHCbBDBeA8wabTcb0+dGDcr0BopGgAwDgEr8mX4kcfNp1YG73AS8H0HCT2YeNAADAdZxh9KZL75s3DzaSfQBbdm9RUl9vICTo8Jy6+gYV5E1zOwzAd/x+/Ypfq9XwNir16A9DXAAArqAqCy/imiUkAwk6AOOZfhFXPKjemssP741TSSPVXaB/r776qlJSUvr8FBYW2rJ+73/rGaikMlV19W5HAWAgXqyAZWRX6HzLJrfD8JyBEnA/JOcAkm/q1KlqbGyM/N7U1KSFCxfqk5/8pC3rJ0EHAAAAEjBy5EiFQiFJ0vnz53XjjTdq7ty5qqqqsmX9DHHxiZNXXW3r+riFduLq6htUV9/gdhgICDvOAPA5d8/20v1U7wGfWL16tdrb21VTU6MRI+xJramgA0DAhMcVPzbzCpcjcQ7JLwZSUpmq0mfdjgJ+cNddd+mJJ57Q888/r8zMTNvWS4IOAAAAJOhnP/uZvve97+nxxx/XFVfYW/BgiIsLuCoeGB6qowAANx0/flxf+tKX9K1vfUvXXHONmpqa1NTUpLfeesuW9VNBBwAAvkRBzNtMurNnby+88ILOnj2ru+66S3fddVfk8fnz56u2tnbY6ydBB3yMLycA6F94/3jvXG+kQiWVqfpp9dtuh4F3rFq1SqtWrXJs/QxxMYgX52UGAGAovPidR9EDyUKCDiAujPsenp5f7HzJAwAGQoIOAPCFRA4iOUgCYDJvDLwCAPThxSECAIDBUUEHAIcwLAjJwIEa4D8k6PCkuvoG1dU3uB0GIInhEojNj8mz2/09aAe+brc33EGCDgCIS9ASIwBwCwk6XENVAPAmPrsIkpLK/i/X8+LZES/GHFQk6AAAYMjCB2wcuLkj1gEEvI131QGMjTZX9xfISbfD8CTaDoBJNi4r1tqHfxH5/2dcjgf227isOKmvF+5PJiBBB2AsDgrMxXh0AHAOQ1xsVLCzQAU7C9wOAwAAWzF8BYj20EMPacKECero6Ih6/HOf+5y+9KUvDXv9JOgexc4SAGIrqi1zZD/JRXYAJOnzn/+8urq69Oijj0Yee/PNN/WLX/xCX/7yl4e9fhJ0JFVJZapjX5wAkm+wz7ITY0i9mCQneywtAGeNHj1ay5cv14MPPhh57Mc//rGmTJmiwsLCYa/fkQT9jTfe0F//9V9rwoQJGjNmjD74wQ/qyJEjkecty1JVVZVyc3M1evRoFRYW6sSJE1Hr6Ojo0Jo1azRx4kSNHTtWS5Ys0euvv+5EuAAwbMykAJiDIhCS4atf/ar27t2rN954Q5L04IMPatWqVUpJSRn2um1P0FtaWvSxj31MaWlpevzxx/Xb3/5WGzdu1GWXXRZZZsOGDdq0aZO2bdumw4cPKxQKadGiRWpvb48sU15erj179mj37t06ePCgTp8+reLiYnV1ddkdMgD04cUqLQAgea699lrNnDlTDz30kF588UXV1dVp1apVtqzb9pLP97//fU2dOjWq5D99+vTI/y3L0pYtW7Ru3TotXbpUkrRz507l5OSopqZGt9xyi1pbW/XAAw/oRz/6kRYuXChJ2rVrl6ZOnaonn3xSn/rUp/q8bkdHR9RA/ba2Nrs3zROo4gEAElVUW6bHZl4R9XsymXZATAUe8frKV76izZs364033tDChQs1depUW9ZrewX90Ucf1ezZs/X5z39ekyZN0rXXXqv7778/8nx9fb2ampq0ePHiyGPp6emaP3++Dh06JEk6cuSILly4ELVMbm6u8vPzI8v0Vl1draysrMiPXQ0EAIAdMrIrjEtEER8SdsSyYsUKvfHGG7r//vt1880327Ze2xP0V155Rffcc49mzJihX/7ylyotLdU3vvENPfTQQ5KkpqYmSVJOTk7U3+Xk5ESea2pq0qhRo5SdnR1zmd4qKyvV2toa+Tl16pTdmwYAAABEjBs3Tp/73Of0rne9SzfeeKNt67V9PMTFixc1e/ZsrV+/XlL3+JwTJ07onnvuiZoXsvcAesuyBh1UP9Ay6enpSk9PH2b0AADAaUW1ZdpfuN3tMGA4k+7sOZDGxkatWLHC1jzU9gr65MmT9f73vz/qsauvvloNDQ2SpFAoJEl9KuHNzc2RqnooFFJnZ6daWlpiLgMAgJM2LitmekQAMb311lvavXu39u/fr7Iye4dB2Z6gf+xjH9NLL70U9djvfvc7XX755ZKkvLw8hUIh7du3L/J8Z2enDhw4oHnz5kmSZs2apbS0tKhlGhsbdfz48cgyAOAGxqICACTpQx/6kG655RZ9//vf15VXXmnrum0f4vI3f/M3mjdvntavX6+SkhI9//zzuu+++3TfffdJ6h7aUl5ervXr12vGjBmaMWOG1q9frzFjxmj58uWSpKysLK1evVpr167VhAkTNH78eN1xxx0qKCiIzOoCwBs2Liv2zGnKoOGCRQAYuldffdWxddueoH/4wx/Wnj17VFlZqe9973vKy8vTli1btGLFisgyd955p86dO6dbb71VLS0tmjNnjvbu3avMzMzIMps3b1ZqaqpKSkp07tw5LViwQDt27NDIkSPtDhkAAuvSGYGTrsYBALjEkUmzi4uLVVwce9xeSkqKqqqqVFVVFXOZjIwMbd26VVu3bnUgQgCAl1DtR7xKKlP10+q33Q4DGBbbx6ADAAAA8bIsSxcvXnQ7DMclso3cdhIAEJfu4TAMhfGCeK/9uHRmwp6Ln7mIGolKS0vTm2++qcbGRl28eFGjRo0adNptr7EsS52dnfrTn/6kESNGaNSoUYP+DQk6MER19Q0qyJvmdhiAb2VkV+h8yya3w4CP1NV3T/lckDdNJZXJS4Hoy7GlpKRo7dq1euaZZ/SHP/zB7XAcNWbMGE2bNk0jRgw+gIUEHRgEFSHncdMSAAiuP/3pT5o0aZLGjBmjrq4ut8NxxMiRI5Wamhr32QESdEDMZAH4WfjzzUEgYK6UlBSlpaUpLS3N7VCMwEWiAICYOIOEoEjmkBdgMCToAAAABiipTOVAAZIY4oJBMGsD4D/hBKD0WZcDAQD0iwTdBSWVqapzOwjYKjwzgFd4LV4A8AtuuoV4MMQFQKA58WW5vXS/7etEXyQ6APyKBB0AANiKcdTA8JCgA0CSkbygqLaMGXIAxESCDgAGIFnzh4zsCk8OveGgETALCToAzyGZ9S8vJrd2MaVfB/k9AExBgg4ADjEl4bJD723xyxANP2xDslFtB5xHgg7AUVTjAABIDAk64FEblxW7HcKQkLADQF+cmUBPJOhAkrDzBeBn7OMA+5CgAwHDXUSHxquzcwSJF8eTezFmAM4jQQeAYaJyCLiPzyH8hAQdAAZBlRMAkEwk6AAAeBgHkID/kKADAOAwrl+AHZwcxkMfNQsJOgDAt6guwyn0LTiJKyoAj/Hq/OcAACA+VNABAJ7hxgGqqaf+SypTmbnEZrQnTEGCDgAAABiEBB3oZXvpfm0v3W/bugCvod8CgLtI0AEAAACDkKDbqK6+gduoO8Ttil6scYm83+gPszvYLwhtyvhnAGHsDQAAGET44lRTLxgF4C9U0AEkjCQFwGCYZWboaDfQAwAAnsHBIRC/IAwN8ysq6PA9dlBA8hTVlvGZSwI/tLEftgFwCgk64HHx3LiF06UAAHgHCTqMFe/MLXbOWw6zcGABO1GxBeAVJOiAg0gIAHcE7QJFrw0tCtJ7AwwFnxAAQKB1J7YnB1yGi1O9i4MBeBEVdAAAAMAgJOgAAACAQUjQAQAAAIOQoAMAAN+rq29wOwQgbiToLuPCI8C/vHpxmpdmAwEkvkvjxWfbO0jQ4Xlem14MGCr6OZyUkV1BousjXi0QoBvvHowVz9RnAJCIoB/kZGRX6HzLJrfDCCQOfpAIKugeEs8t3QGYiWoWegrajZRMM5zx6Jy1RTKwdwB8ZDhnHbiAyhuowgGA/1FBBwAAAAxCBR0AXFRSmao6t4MwDGcJvIthO4A9+CQBgEddGgfLxdRwDweZzqJ9g4khLoALqDIBQ8cFegD8jgQdAIBhYlgOADtRxgN6YdgAED9OvwOA/aigw1gMAwEA/2PIUmxMfxtcJOgAAACAQShRAoBLvHyWyMuxwxuGM64/qEOv6uobdFK5bocBG1BBB2A8TvMCAIKEBB2AEYpqy1wfi7q9dL+rr4/ExNNfqPQjmdiHwC4k6B6zcVmx2yEAAIBhMKEgAbNRWoDnhStkQRxvmGxUh4DBdSdeTNMKYOiooGNAnB4GgofPfWxUPQEkAwk6AMBYGdkVnr5LJwk9gKGgTOJDXp7xwsuxA/COoE7D11NRbZn2F253O4zA44wV+kMFHY7jwlY4hTHxAAA/IkGHJ7h5mpjqhj8x9AAAYCoSdI/qnTSSRJopyO+Ll8cNwzs4i5IYDkwBbyBBR79IroItyAcWdqD9AADDQYIOuIjx+QAAoDcSdCBBVEcB2IX9CYD+OJ6gV1dXKyUlReXl5ZHHLMtSVVWVcnNzNXr0aBUWFurEiRNRf9fR0aE1a9Zo4sSJGjt2rJYsWaLXX3/d6XBtF542kOkD/Wn6+Zqof/3A1DG9JDLOq6tvGHBfNdz3wK2hcyWVqfQfF4T7k6n7lKDiOgRvcDRBP3z4sO677z594AMfiHp8w4YN2rRpk7Zt26bDhw8rFApp0aJFam9vjyxTXl6uPXv2aPfu3Tp48KBOnz6t4uJidXV1ORkyPGL6+RpfJcUAho8kHIBfOJagnz59WitWrND999+v7OzsyOOWZWnLli1at26dli5dqvz8fO3cuVNnz55VTU13wtXa2qoHHnhAGzdu1MKFC3Xttddq165dqqur05NPPulUyEbhIs3YSMz9jSQLflVUW+a56iWfR8AdjiXoZWVl+uxnP6uFCxdGPV5fX6+mpiYtXrw48lh6errmz5+vQ4cOSZKOHDmiCxcuRC2Tm5ur/Pz8yDK9dXR0qK2tLerHj0jc7UebAki2wYYTmcytuDOyK9hfIzAcSdB3796tF198UdXV1X2ea2pqkiTl5OREPZ6TkxN5rqmpSaNGjYqqvPdeprfq6mplZWVFfqZOnWrHpsBwjG0MNq9VI70okQpqz88jMxQB5uOAx1y2J+inTp3S7bffrl27dikjIyPmcikpKVG/W5bV57HeBlqmsrJSra2tkZ9Tp04lHjzgcf3dwIpT1AAAeIvtCfqRI0fU3NysWbNmKTU1VampqTpw4ID+5V/+RampqZHKee9KeHNzc+S5UCikzs5OtbS0xFymt/T0dI0bNy7qBwgCKiAAvM4rZ8MoeCBZbE/QFyxYoLq6Oh07dizyM3v2bK1YsULHjh3Te97zHoVCIe3bty/yN52dnTpw4IDmzZsnSZo1a5bS0tKilmlsbNTx48cjy8D72NEBAAD0ZXuGlJmZqfz8/KjHxo4dqwkTJkQeLy8v1/r16zVjxgzNmDFD69ev15gxY7R8+XJJUlZWllavXq21a9dqwoQJGj9+vO644w4VFBT0uegUAAAA8BNXSph33nmnzp07p1tvvVUtLS2aM2eO9u7dq8zMzMgymzdvVmpqqkpKSnTu3DktWLBAO3bs0MiRI90IGQB8xyvDCuBt3f3spNthAJ6SlAS9trY26veUlBRVVVWpqqoq5t9kZGRo69at2rp1q7PBAQAAAAZx9E6iAID4MW0oAEAiQTdC+OYLRbVl3IgBxuKiXncEOWlnCA6AoCJBBwDAYW4cbMS6WRRFIMB8JOgAolC1BLyPJHxo6uob3A4BkOTSLC4AADOEhy7VuRwHnBOrkg77MFMN7EYFHQDQr5NXXe12CICvccYSsVBBBwDDhSugn3E5DiBeDBUBhocKOhzHWEgAGLogzeTD9wXQjQo6AAD96JksBm0ogsnTqjLeG0FABR1AoJmaeAWpagoAiGbuIbLHTT9fI0l6NWN55DHG5HUzuTITFJH+6W4YGABVQuDSAfRjM69wORIguaigAwAATzP1TBgwVCTogE2mn6+JVKYBDA2JFgCQoAMAAABGIUEHAAP4+doMO+9kGbQKuxsXCwetjYcjI7vCd1ND+m17vMq/3wgAYCOSFgBAslBBBwDAw/x89gUIKhJ0AAAAwCAk6B7DaXb4kZcqgF6KFQg6r39emRksuLzdcwFDRS7smutuHMly6cCRG+vAW0oqU/XT6rfdDsN4JZWpqnM7CCBASNCBXsIVF76MAACAG0jQAXgWVT0EldeHbsA59A1/YAw6kqKotozx87ANX0CAVFffoLr6BrfDGBKvxg0kC99yPhS+qORVd8NAL3whwY+c6NeJHMx3L8u1D0hcuO8W5E1zORKgLyrowBB59ep6r8YN/+GsGgD0jwQdAAAEFmc3YSISdADwIKrPgLO41gVuIkFH3CJzezuI4Rf+QyKJZEvGvgoAnESCDgAAAmH6+RoKQfAEEnQH8OEHAADAUJGgIyFuzWcehLGAG5cVux1CH6a0e9APejOyK5SRXSHJ2c8g1UXAXgzxw1CRoMPzvHyzDiBI7E5WTDmADHoSxv7XXbS/P5mxdwMAAI4In30B4B1U0BEYfEkB3udG1TzoFXLYi/6EeJCgAwAAAAZhiAsAAB4XPrNQ53IcduBsJ0AFHYDBvH7xk9fjB4LKlAuQEVwk6ACQZH5N3BlbCwD2IEEHAPgKBwoAvI5zOAAA2KikMnVIY8HdOLMSPpi5dy7pAGASKugAgEExJhdD5VTf4WJS+BkJOgAYzGvDNUoqU0nmAWCYSNCBgNheul/bS/e7HYZnuJEYey0ZN1ldfUPcQ0Y4oMBg+usjXuo3XooV3UjQAQAAAIOQoMM3Ni4rdjsEICa/Tq3oBIbJXMJZFftMP1/jdgiOC8I2BgUJOuAiLnICMBiSdCB4SNABJJVfkw0q5IC9/LqvMAUFIrNxDtEHTl51tcTpYMBIBTsLJGlI82I7JXwwcVK5LkeCIBvqfPFAEJDVAQASYur48PCBR0llqurqXQ4GvlFX36CCvGluhzEkpn5WMTiGuCBu98693e0Qho2dFbyGPgsAwUOCDledvOrq7iE6AIB+cZDmLGY+gYlI0F3Q+2KyotqyyMUwPf/vR6ZeSJfsad1MbYfhCupUl359P+0wnAvRSEzhFJJymI69HwBJ3V9Yr7odhE+QsPuDn4sliN/08zXK1LfdDgMBQwUdcIDfz4QA8CbOSgDeQIKOQErm/K/bS/dre+n+pL0eMBRO9tEgzbfMgTkAO5CgG2zjsuLAjumVpPaTd7sdgtESTXrsrOoX1ZYpI7siUIkXAOfZsY9iiFn8OKA0F+e64Dh2AABwaXgJN+cBMBgq6EACwneFdEqQz5j4ETNFAMFCQQp2IUEHgIDigkHAOQwBxHCwdwYQNz9Vh0oqU/sMNejevpNuhINB9Hy/+nvv3Ebf6R77XZA3ze0wHJfsA1vG1AcTFXQAUaiqAoA/sX/3DhL0APLTlH9+qugCphnoy7yuviEplT2/3FPAC4mRHW2drH4RVOG29dP3OPpHgg4A8EUS7CXTz9d49iJir8btNyTp/kaC7gNeqMwEmV8qgEguv36uuXAueWhrc/AdgESRoMPX2CkCcIvT07IOhH1f4pwemmPH8CEnl4dZ/FmiATwoI7tC51s2uR0G4GumnZkgiYIpmInILGbtqQAESuTGTJ9xNw4AAEzCEBfAQcmujlGNA2IzrXruN1w8Ojzsv9ETCToAeIQbCaZfkwa/bldQ2fXZ4CADprA9Qa+urtaHP/xhZWZmatKkSbrxxhv10ksvRS1jWZaqqqqUm5ur0aNHq7CwUCdOnIhapqOjQ2vWrNHEiRM1duxYLVmyRK+//rrd4cJlJZWpVLUQE1+W8Lqg7t/47ALDY3uCfuDAAZWVlem5557Tvn379Pbbb2vx4sU6c+ZMZJkNGzZo06ZN2rZtmw4fPqxQKKRFixapvb09skx5ebn27Nmj3bt36+DBgzp9+rSKi4vV1dVld8gAMCSRMfSIidlEACBxtifoTzzxhFatWqVrrrlGM2fO1IMPPqiGhgYdOXJEUnf1fMuWLVq3bp2WLl2q/Px87dy5U2fPnlVNTfcRd2trqx544AFt3LhRCxcu1LXXXqtdu3aprq5OTz75pN0h22aoN56ws4ocz2nbRL4wTTkNbEociYrVH6gueYvdVdCe73+8n9lkJrpe+rz19954KX5TONFmXr4ZE+A2x8egt7a2SpLGjx8vSaqvr1dTU5MWL14cWSY9PV3z58/XoUOHJElHjhzRhQsXopbJzc1Vfn5+ZJneOjo61NbWFvUDAAD8ibMz8DNHE3TLslRRUaGPf/zjys/PlyQ1NTVJknJycqKWzcnJiTzX1NSkUaNGKTs7O+YyvVVXVysrKyvyM3XqVLs3BzajyoWM7ArudggAQC+OJui33XabfvOb3+gnP/lJn+dSUlKifrcsq89jvQ20TGVlpVpbWyM/p06dGnrgoDKBYRtqH+LADQAQdI4l6GvWrNGjjz6qp556SlOmTIk8HgqFJKlPJby5uTlSVQ+FQurs7FRLS0vMZXpLT0/XuHHjon5MxsVlSDbTZ8wxPT7ADXwmgGCyPUG3LEu33XabHnnkEe3fv195eXlRz+fl5SkUCmnfvn2Rxzo7O3XgwAHNmzdPkjRr1iylpaVFLdPY2Kjjx49HlgEAwEkcNAJwi+17nrKyMtXU1Og//uM/lJmZGamUZ2VlafTo0UpJSVF5ebnWr1+vGTNmaMaMGVq/fr3GjBmj5cuXR5ZdvXq11q5dqwkTJmj8+PG64447VFBQoIULF9odMgKKsc8YruH2IT8O5wlv00Zd7nIkAOBdtifo99xzjySpsLAw6vEHH3xQq1atkiTdeeedOnfunG699Va1tLRozpw52rt3rzIzMyPLb968WampqSopKdG5c+e0YMEC7dixQyNHjrQ7ZAAAjOXHAzkAA7M9Qbcsa9BlUlJSVFVVpaqqqpjLZGRkaOvWrdq6dauN0QEA/CRZyWt4qMtPq99OyuthePx+UFNUW6b9hdvdDgMOcnwedAAAJC54RGz0jcQx25q/kaB7DDsxAIBJuFuos2jfYCJBh6c5ecDCVJiA/Zys+pHI2INCEOA+PoUIjKLaMt07139dvjvhOena628v3a8i114dcF7vhDU8e8/5lk1uhBNYJg7pMDEm+AMVdMDH6uobfH+xFMxB5dU+20v3ux2C57Cvg5+QoLuMo28AvfVMznonHdw8BxLJKOB3JOiAizhAA2Anrp2J3/TzNVy3AGORoAMAfI2hXgC8hvOkPsGXj1kiQxDcu3YTAAZUUpkaiBsvhb8fp7sbRkw9v7/5LkcYFXRDDTb0gQ+x93Aq1UxOj+f2evV2sH7rZr8eqF29PHyBC0Td4dX+MlQMsTQbFXQAiAMXZsIJ08/X6FW3g/CYkspU1bkdxADC03ACw0EFHb5HYgUA3sD+GujGJ8EgPXdM/Y0L7D4qvz2JEZkjVpWprr4hclpyje6x5bXCp8036vKo34cqkfHoVF78J9w/X81YPuiyiVYGw33L7mqim1VdOxK0oA1ViMXtm5jBHYN9Z5l+BgLdqKAD8Bwvj+nGJVRLndOzbYtqyxhvDHgMCToAAABgEBJ0ADCM36qdbm4PZ1u8h4r/4Jw8+0Tbm4EEHQDAcBMAMAgJOoA+qDoCgNk4qPY3EnQAAADAICToAAAgqZgK0yxM8WseEnQ4jtNwQLANZx9QUplq1D7EpFiCjvcCfkaC7hKOVgH/ztZAdRDDNdjngj4G+BsJOgAgLlw8DADJQYIeQH49LTjUsxKmVHBNiQPwEirJiUtkX8NBGeAOEnQAADyMJBrwH3+WUgGXhb8wp7sbBuA6087Y9ZfMmpLgTj9fo1fdDgKSuvtEQd40t8MYkCn9Fs4wa88ZQKZ9eQEIlp5f8tPdCwNDEB7e82rGcpcjuSSoQ45KKlNV53YQw+SHbfAThrgAAADHUOkFEkeCDgSAX78g3dguv7Yl7DPcPjKUC8br6hvom/2gTeBVJOgAkARBPfUPeJEfPq9+2IYgYwA0AGBQVCKds710f/d/LnM1DAAGoYIOwHNMqQwxd713mdKHgoC2BhJHgm6A/sYOJjK7S6KVrXiW7+/16+obVFKZOuyZZ+wcK5nIjt+OGXOCXEXsue0929KOJDWedo21zPbS/ZcqkB4W5L7VW3gf4ac28dO2JJMd3zl2aD95t9shIGDc7/UAAHhI+KD0sZlXDOnvt5fuZzgLgAFRQQcAAAgYE85MIDYSdABRgngqPiO7wu0QAACIIEGHq/w2zhQDIxEG+iqqLeOCYwBRSNABAHAAxQcAQ0WC7jHs8AHAPqbvU02PD4AzSNABAAAAg5CgA7ANY8wBABg+5tgBMGzh6bpKn3U5EJ/w8/RnJZWpqnM7iDiFh5cU5E2zdX3bbVmbf/i5vwNDxacCCIhwcjDd3TAAX/LSgQe8j2sT/I8hLgAA35p+vsbtEAAgYVTQDcZ4XudcmnP4ZMxlTG//7m2IHX9yXh+ASUzfbwUJ+0gMBxV0IEGcWvQHL76PXow56Lz6niXjJnJ2rt+r7QzEQoIOBMT08zWc7jeMSe+HXbGY3M+Gk8S5dddjU9syCExP+ukb/kaCDriI2QsAAEBvJOguYWwaTEZlBnCXmwfvpleOgSCgfAfH1dU36KRyXX19mCmchNTVuxdDXX1DYOelHuyzUVffYOu0nHX1DQnPKd57LnI3Ps+9XzPcb1c9lvRQhqy/aSB775sH+jwmazx6It8XbvcFu+bHd5Ldc/kjeaigAw4weRyunwShjYOwjckS/lz2blPT29j0+PzKxHbvLyYT48TwkaADHpORXcFUakAvJ6+62u0QAMA2JOgAEubXYUMmXbQbq43dmk0EQHy4xgx2MOfbCH2+dE/2uplOUW2Z9IXhv8708zV6NcFYBnvcBIPduGeg8a/9bZddNwJqP3l33Mu6ffOhgUw/X6PHdUefx2MltYne+nywvhVPv3XKUMZOxyu8XdPP18R8600Yq99bvKfVTd5n+FHPfUg8485NF+4/G3V5n8fcVlffoILK+PcL4fejpDJVP61++9JjNn+u+/ts9hzj3/OxgsppqlN3v7l3LimhSaigA4jCeEb0h34BAMnD4RIA25h8BgBAXyaenQFAgg4AMV26GPdpV+OAv0WGiXFsC+AdDHEBABjDlPHFAOAmKugA4sbpcMB5Js0mBMAdVNABAAmhyg0AzuIwHQD6UVKZqtJn3Y4CgB9xlgSDoYcYLDKcwOU47ODFoRGJzuON/g3Ujl6vxNJHIJmVbPXXJ02KLx7heMNzhQ/k5FVX6+r/4epa+I+3PrUAAGPV1Tdo+hD/lhlzonn94BXA8DAGHQAAADAIFXQAgRfrBkvdj0v6QnLjgb24CyqSpaQyNa6hOcBgSNBdwthV+En4dPxJ5Q64XCTh9SCGHCTmM7/+fVJfr6QyNWk3+gn3hY26vM9jpoonPtO3AclBfmIGEnQAceMLPDlo52Dj/QdAgg4ASIifhoyYsC3hhHy6u2EAMAgXiQIAkET3zr3d7RAAGI4KusH8dJrTi9vixZjdMFg7Od2O4WSnzOZRk3X1DYOOqaePJC7cZgU7CzzZfv3FbNJ2hGOpq29QQd60qMcGWt5NvWOI95oWr+n53gCDMb6C/sMf/lB5eXnKyMjQrFmz9PTTzJELJEv7ybujfp9+vsaIIQF+Z2obDxaXqXF7DZ+z+HntJkzD0bNP0D/8z+gE/eGHH1Z5ebnWrVuno0eP6hOf+ISuu+46NTRw9AkAJgtS4jRcJOQAejN6D7pp0yatXr1aX/nKVyRJW7Zs0S9/+Uvdc889qq6ujlq2o6NDHR0dkd9bW1slSW1tbUmL92LH2T6PtaVYUlub1GG980D3/7vOdamtra373/BzsbyzDae73lk2vL7wtvVadziOWNt+sePs4O0y0Ov0iCkug21fj/X1jG2g//fndFdXn+3q3Wbhdj/d1fXOy7ZdWiYcR+/t7O+5traodu4ZU+/fozczenvC8YRj7blMLOc6z0Rtb9R293y/wr+/o9/+2WPZcEy9t6u/dm9LsS7F2+M1esYS2bYeMfXZxlh9o1c79e7TPR/vOhdfu/XRo1+E4+75foQf772NPeOL+pwMEPelP4u9XX3a/Z19R9Q29mrrrnMpQ9/2XvqN9Z39SaxY+9uusKj3v3uhqP9H9n3v/G24X8faH8ZaX8/PUqx2j6W/vhVlkH3X+QsX+sTUZz/de9slPfM/ryl/gH1ErH7RU+/PY6+VRO3vescR+ZzGs2/uGX8/++Oevw/U7v3FH+v7Maznvvn8hQvRn8dBtrG/9Q0W30Cfx8g29Gqz8Gcz/HzXuS79xV+9Eft7JsF2jyemgb4fo/rGALlIeH/SZz+YBOHXsqw42yUgUixDW6Szs1NjxozRv//7v+uv/uqvIo/ffvvtOnbsmA4cOBC1fFVVlf7+7/8+2WECAABgmE6dOqUpU6a4HYYxjK2gv/nmm+rq6lJOTk7U4zk5OWpqauqzfGVlpSoqKiK/X7x4UW+99ZYmTJiglJQUx+OVuo8Cp06dqlOnTmncuHFJeU2vo80SR5sljjZLHG2WONoscbRZ4vzWZpZlqb29Xbm5/rooeLiMTdDDeifXlmX1m3Cnp6crPT096rHLLrvMydBiGjdunC8+NMlEmyWONkscbZY42ixxtFniaLPE+anNsrKy3A7BOMZeJDpx4kSNHDmyT7W8ubm5T1UdAAAA8AtjE/RRo0Zp1qxZ2rdvX9Tj+/bt07x581yKCgAAAHCW0UNcKioqdNNNN2n27NmaO3eu7rvvPjU0NKi0tNTt0PqVnp6u7373u32G2iA22ixxtFniaLPE0WaJo80SR5sljjYLBmNncQn74Q9/qA0bNqixsVH5+fnavHmzPvnJT7odFgAAAOAI4xN0AAAAIEiMHYMOAAAABBEJOgAAAGAQEnQAAADAICToAAAAgEFI0G3ywx/+UHl5ecrIyNCsWbP09NNPux2SMaqqqpSSkhL1EwqFIs9blqWqqirl5uZq9OjRKiws1IkTJ1yMOPl+9atf6frrr1dubq5SUlL085//POr5eNqoo6NDa9as0cSJEzV27FgtWbJEr7/+ehK3IrkGa7NVq1b16Xcf/ehHo5YJWptVV1frwx/+sDIzMzVp0iTdeOONeumll6KWoa9Fi6fN6GvR7rnnHn3gAx+I3Oly7ty5evzxxyPP08f6GqzN6GPBQ4Jug4cffljl5eVat26djh49qk984hO67rrr1NDQ4HZoxrjmmmvU2NgY+amrq4s8t2HDBm3atEnbtm3T4cOHFQqFtGjRIrW3t7sYcXKdOXNGM2fO1LZt2/p9Pp42Ki8v1549e7R7924dPHhQp0+fVnFxsbq6upK1GUk1WJtJ0qc//emofvfYY49FPR+0Njtw4IDKysr03HPPad++fXr77be1ePFinTlzJrIMfS1aPG0m0dd6mjJliu6++2698MILeuGFF1RUVKQbbrghkoTTx/oarM0k+ljgWBi2j3zkI1ZpaWnUY1dddZX17W9/26WIzPLd737XmjlzZr/PXbx40QqFQtbdd98deez8+fNWVlaWde+99yYpQrNIsvbs2RP5PZ42+vOf/2ylpaVZu3fvjizzxhtvWCNGjLCeeOKJpMXult5tZlmWtXLlSuuGG26I+TdBbzPLsqzm5mZLknXgwAHLsuhr8ejdZpZFX4tHdna29W//9m/0sQSE28yy6GNBRAV9mDo7O3XkyBEtXrw46vHFixfr0KFDLkVlnpdfflm5ubnKy8vTF77wBb3yyiuSpPr6ejU1NUW1X3p6uubPn0/7vSOeNjpy5IguXLgQtUxubq7y8/MD3Y61tbWaNGmS3ve+9+mrX/2qmpubI8/RZlJra6skafz48ZLoa/Ho3WZh9LX+dXV1affu3Tpz5ozmzp1LH4tD7zYLo48FS6rbAXjdm2++qa6uLuXk5EQ9npOTo6amJpeiMsucOXP00EMP6X3ve5/++Mc/6q677tK8efN04sSJSBv1136vvfaaG+EaJ542ampq0qhRo5Sdnd1nmaD2w+uuu06f//zndfnll6u+vl7f+c53VFRUpCNHjig9PT3wbWZZlioqKvTxj39c+fn5kuhrg+mvzST6Wn/q6uo0d+5cnT9/Xu9617u0Z88evf/9748ki/SxvmK1mUQfCyISdJukpKRE/W5ZVp/Hguq6666L/L+goEBz587VFVdcoZ07d0YucqH9BjeUNgpyOy5btizy//z8fM2ePVuXX365/vM//1NLly6N+XdBabPbbrtNv/nNb3Tw4ME+z9HX+herzehrfV155ZU6duyY/vznP+tnP/uZVq5cqQMHDkSep4/1FavN3v/+99PHAoghLsM0ceJEjRw5ss8RanNzc58KAbqNHTtWBQUFevnllyOzudB+scXTRqFQSJ2dnWppaYm5TNBNnjxZl19+uV5++WVJwW6zNWvW6NFHH9VTTz2lKVOmRB6nr8UWq836Q1+TRo0apfe+972aPXu2qqurNXPmTP3zP/8zfWwAsdqsP/Qx/yNBH6ZRo0Zp1qxZ2rdvX9Tj+/bt07x581yKymwdHR06efKkJk+erLy8PIVCoaj26+zs1IEDB2i/d8TTRrNmzVJaWlrUMo2NjTp+/Djt+I7/+7//06lTpzR58mRJwWwzy7J022236ZFHHtH+/fuVl5cX9Tx9ra/B2qw/9LW+LMtSR0cHfSwB4TbrD30sAJJ+WaoP7d6920pLS7MeeOAB67e//a1VXl5ujR071nr11VfdDs0Ia9eutWpra61XXnnFeu6556zi4mIrMzMz0j533323lZWVZT3yyCNWXV2d9cUvftGaPHmy1dbW5nLkydPe3m4dPXrUOnr0qCXJ2rRpk3X06FHrtddesywrvjYqLS21pkyZYj355JPWiy++aBUVFVkzZ8603n77bbc2y1EDtVl7e7u1du1a69ChQ1Z9fb311FNPWXPnzrX+4i/+ItBt9vWvf93KysqyamtrrcbGxsjP2bNnI8vQ16IN1mb0tb4qKyutX/3qV1Z9fb31m9/8xvrbv/1ba8SIEdbevXsty6KP9WegNqOPBRMJuk22b99uXX755daoUaOsD33oQ1FTcAXdsmXLrMmTJ1tpaWlWbm6utXTpUuvEiROR5y9evGh997vftUKhkJWenm598pOftOrq6lyMOPmeeuopS1Kfn5UrV1qWFV8bnTt3zrrtttus8ePHW6NHj7aKi4uthoYGF7YmOQZqs7Nnz1qLFy+23v3ud1tpaWnWtGnTrJUrV/Zpj6C1WX/tJcl68MEHI8vQ16IN1mb0tb5uvvnmyPfhu9/9bmvBggWR5Nyy6GP9GajN6GPBlGJZlpW8ej0AAACAgTAGHQAAADAICToAAABgEBJ0AAAAwCAk6AAAAIBBSNABAAAAg5CgAwAAAAYhQQcAAAAMQoIOAAAAGIQEHQAAADAICToAAABgEBJ0AAAAwCD/H4diYni2mryYAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "q = visit_counts.loc[:, 0:365]\n", + "base = np.zeros(q.columns.size)\n", + "plt.figure(figsize=(8, 6))\n", + "for f in 'ugrizy':\n", + " plt.bar(q.columns, q.loc[f], bottom=base, label=f)\n", + " base += q.loc[f]\n", + "plt.legend(loc=(1.0, 0.5))" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "id": "24a83b61-fb67-42cd-b210-a1fe7ac2be72", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "60796.0" + ] + }, + "execution_count": 75, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "start_mjd = SURVEY_START_MJD\n", + "start_time = Time(start_mjd, format='mjd', scale='utc')\n", + "start_offset = start_time - Time('2025-01-01') - TimeDelta(0.34, format='jd')\n", + "start_mjd" + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "id": "2f5327be-2aca-41c8-8261-e2b39717520e", + "metadata": {}, + "outputs": [], + "source": [ + "almanac = Almanac(start_mjd)" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "id": "9026c55b-bb9d-4abf-8cd8-b29bc1d6ce6a", + "metadata": {}, + "outputs": [], + "source": [ + "start_idx = np.where(almanac.sunsets['night'] == 0)[0][0]" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "8dd4d8a1-2dbe-4539-995f-c3aafdc982a8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3651 3651\n" + ] + } + ], + "source": [ + "sunsets = almanac.sunsets[start_idx : start_idx + 3651]['sun_n12_setting']\n", + "sunrises = almanac.sunsets[start_idx : start_idx + 3651]['sun_n12_rising']\n", + "\n", + "nights = np.arange(0, 3651, 1)\n", + "night_hours = (sunrises - sunsets) * 24\n", + "print(len(nights), len(sunsets))" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "id": "9b48b192-0309-43f5-832c-9ba899aaae37", + "metadata": {}, + "outputs": [], + "source": [ + "sunsets = almanac.sunsets[start_idx : start_idx + 3651]['sun_n18_setting']\n", + "sunrises = almanac.sunsets[start_idx : start_idx + 3651]['sun_n18_rising']\n", + "night18_hours = (sunrises - sunsets) * 24" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "id": "e932d04e-0568-4bd1-b0e6-aba2f2b7c5c6", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIlklEQVR4nO3deXwU9f0/8FfOzQEJRyQkkIRDLEgAMVGuoqIYRYqt2oJVOVrwJ0VFpViLtBWsFVstRauAKMjXbxH5emA9qBCsHBIECUHkPhJIIAkhAXJBNsfO7w/Mms0emd2d3fl8Zl7PPnjUzM7uvj/HfD7v+czsboiiKAqIiIiIBBSqdwBERERE7jBRISIiImExUSEiIiJhMVEhIiIiYTFRISIiImExUSEiIiJhMVEhIiIiYTFRISIiImGF6x2AGjabDcXFxWjfvj1CQkL0DoeIiIhUUBQF1dXVSE5ORmiob2sjUiQqxcXFSElJ0TsMIiIi8kFRURG6d+/u03OlSFTat28P4HJB4+LidI6GiIiI1KiqqkJKSop9HveFFIlK8+WeuLg4JipERESS8ee2Dd5MS0RERMJiokJERETCYqJCREREwmKiQkRERMJiokJERETCYqJCREREwmKiQkRERMJiokJERETCYqJCREREwvI6UdmyZQvGjRuH5ORkhISE4KOPPmrzOZs3b0ZGRgaioqLQq1cvLF261JdYiYiIyGS8TlRqa2sxaNAgvPrqq6r2LygowB133IGRI0ciLy8PTz/9NGbOnIkPPvjA62CJiIjIXLz+rZ8xY8ZgzJgxqvdfunQpUlNTsWjRIgBAv379sGvXLrz00ku45557vH17IiIiMpGA/yjh9u3bkZWV5bDttttuw/Lly9HQ0ICIiAin51itVlitVvvfVVVVAY0x53g57ntjBwBgQmYK5t3ZH6t2nMTZGite35yPmbf0wR0DuuL2RVsBAL8a0QPDenWGtdGGM1V16H1FOyzZfBxHzlRjzf8bhh91df8rke/sKESfxHa4rkcnAMA3J86h4Gwtxl+Xgv/dfgJLNh3HY6P7YMJ1qQ7PszY2Ydaab3Hr1Yn42eBuyD15HvcsyXHYp2/X9pg2shdG9+uC17fkY8mm4w6PL33gWpyttuLr/HP47LsSp9j+MLYfnvvsIGbc1Bvx0RGYPLwHoiLCcM+SHOSePI+jfxmDiDB1i3DWxiaMfeUr/DyjO7p3jMZznx5En8R2GN0vEflna2BTgKG9OiP7QCk+/rYYNgXolRCL/PJa9L4iFsfP1gIALOGhWPPQMDTZFEx5aydWPzgUa74pQnVdA/5y1wDEWgLThU9fuISshZtRW9+EpPgoDE7tgBuvugI7Cs7hw92nnfZ/58EhWLLpOKaN7IWMtI4Y+8pWnKy4CAB4YGgqbuvfFWerrWiyKfjPvlI8evOV6NYhGqt3FuHOa5LRMyHW6TX/9O99KL5QhzcmZbj9Qa9L9U2Y/NZO1DU04c5ByXhpw2HUNdgwdkAS9p6+gKdu74ufDEzGkTPVeHNrPuobbUjpFIMu7S3447/3AwAWjh+E7ANn8J99pegcG4nBqR0RFRGKT/eW4O5ru2Hi0DREhofi6/xzuLlvF9y9eBv+/LN0PPJOHjrFRiL3D6Pdxrf31AXsPnkek4b1QGjo5X0efHsXNh0uQ4eYSDx0Qy+s+aYIYwYk4ZUvjgIAPn30x/i/XUXYWXAOADA4tSNW7ywEAERHhCGrfyLuH5KG0qo6LN10HGGhIVg2KQN/+/ww1uadxos/H4hfZKZgR34FJiz7GlERoXj710Nwfc9O3nQBt1btOIm5a/fZ/97wxA3I/76/ZqR1xAe7T+Fn13TD0AVfAAA+engEFm08gh6dYzHnjr5YvaMQ/9pRiF+N6IFeCe0wrHdnt+91rrYe/9lXgnGDkhEX5TxWestmU/C7D/bCZlMwcVgaXvvyGDYeLLM//pubeuNMVR2abAomXJeC4b0TUF3XgAHzNqBnQizeeXAInvrgO2w5chZThvfAHQOScKysBvcNSfXwro6On63BhYsNyEjrCACoa2jC/24/icGpHfDVsXJU1zXiviGpOFttxf99U4ROsZFoHxWBf/73KBptiv11HrqhF1I7x2Du2n0Y9aMrkNmjEwandsDHe4oxoHs8OsdakF9eg8nDeuDlL47i33tO481J12FA93j7a+zIr8CyLfl45ZeDvR5LTl+4hJ+9tg0v3D0AJyouYkjPTkjvdvm1P/62GO0t4aitb8Qj7+Rh4tA0XKxvwge7T9mfP+pHV+DPP0vHz17LQXmNFT06x2Du2KtxrKwGf/38kMN7/WPCINw1uDs+yjuNXSfPoX9yPH55vWOdHzlTjdROMYiKCGsz9rqGJvzmX7kYnNoRM2/p41W59RCiKIrS9m5unhwSgrVr1+JnP/uZ232uuuoqTJkyBU8//bR9W05ODkaMGIHi4mIkJSU5PWfevHmYP3++0/bKysqA/Hpyj99/5vB3fHQEKi81+Px6J14Y63J7y4SoeZ/m9/7gN8Nwz5Lt9n0/engErknpYP977trvsGrH5cH6wLO34eo/rfc5PrX+/NP+GHFlAm7++2YAwNgBSXjt/mtVPXfOh9/ZJ5dA+eX1KVhw90CXj5VV16HqUgOu7OLbT4sPmr/B5z7w/27ohWVb8tvcb9ygZHzybTEAYEx6Vyx5IMP+WGHFRdzw4pcAgBVTMnFz30Sn55dV1eHGFzfhUkOTx/c59Ofb0fePn3tTBK8sfSADt6d3dflYc/9++d5r8NNrumF/cSXGvvJVwGJp5qrMb07KxOirnevRW63Hi5Z+lNgeh89Uu338xZ8PxJPv73XY5m68AGA/SRjdLxFvTs70PthW3s89hdnvfat6/xMvjMVT7+/Fml1FHvd7Z9oQDL8yQdVrNtffV0+NQveOMXhp/WG8+uUx1TF5q0+XdjhaVmP/u2V9t2xLT+3giqt+cOKFsSitrLMnqVra8uQo+5gAABtn3WAf3z7fV4rp/8pFerc4fProSI+vc662Htf+Odsh5kCqqqpCfHy8X/N3UD710/psqzk3cncWNmfOHFRWVtr/FRV5Pki05k+S4knzGbaax06fv+Tw97vf/FAHp1o9Fih7T1WirPqHlS1XqzDufLq3OBAhOVi9032/uP4vX2D0wi04fcG3uvKnD3y8R13Z1+8vtf/3f/aV4mJ9o/3vszU/1Hvz2Xpr41/f3maSAgC7C8+risdXx8/WtLnPke8n76Jz7o8BLbU882427e1dAX9fT0kKABR6Wf7ck5fbbuPBMz7H1NKhEu9Xp9cfKG1zn5M+tGtB+eV+Hej+2TJJCYZztfUBed2WYwIAXLj4wxj13veJ5L7Tbbdvy3FHFgG/9NO1a1eUljpWTFlZGcLDw9G5s+slT4vFAovFEujQyAQOFlehW4dovcMIiBMeEl8Sk+tTMyLyJOArKsOGDUN2drbDtg0bNiAzM9Pl/SnkGgc442GbBo4fV7QDyt0qMhG553WiUlNTgz179mDPnj0ALn/8eM+ePSgsvHw/wpw5czBp0iT7/tOnT8fJkycxa9YsHDx4ECtWrMDy5csxe/ZsbUpgEhzfzIXtbUwytqugOZ/pydiXfOX1pZ9du3Zh1KhR9r9nzZoFAJg8eTJWrlyJkpISe9ICAD179sS6devwxBNP4LXXXkNycjJeeeUVfjRZEj4fCxzcNMOzcOMI4Toakde8TlRuuukmj8uqK1eudNp24403Yvfu3d6+FbVweYDj7E+kRqCOFH8vKTHndMYVG1/51plkrG/+1o8sZBvgZIuXgkLGQVJLeh8WvlR/oO734eqSdvxpIVHv52qJiQp5xMsO+ghx898kNx5ORN5joiKItgYwPcY3DqqBxfoNHFFPEpn4k1bM1JWYqEhCj07p12Av6EQRDGYaQIiIAo2JiqA42WlDj3xJ1LN58p+/bWvU49qXYilmPpvxiTb1JWO9M1GRhF43nvk8sBp0QA4kJjgaErQueQMpacVMPYmJiiREOBPzJgYBwhUeJy3zEeE4FgX7vxhkOEFiokKqcVhRx9fJyAyTmBnK6InexffpEo0EE5kZmenGbCYqQeRNt2o9OOj1cVXzHAra0aLOWg5CRhiPgj3ZiXodPlTnxvTpe1Q0j4K0JsN3ofiDiYog2loGbTlxBatL+jOminLYGGCOJ4H426+NkHS6IsrxLgJRk+RmMuY0TFQkZNOpp5lpqdEfaluH1Ulq6dlXAn22LvrELo4QD38ZGxMVSYjQKUWIwUg8jf+sa/+Ietboa7KvZ38QtCqFJfpNwjKeIDFRkUWLzhXMQVjGTm00bALSU6DGG44t/jFT/TFRkVCwEhVRz0pFp3b8MNNAQ5f5/LVE7CwUIDIM80xUBNHWdVrphilBer8gYVCQidruodIdyIHDEyFvafTNtBLWOxMVQXk6gQruzWccWb0l4TgQNDIOki35e2Op3veo+BK+mvHGn/hEv6dDVC3rTfLDqk1MVCSk12Cv93dAGJ3skzi1jV8G6Mxon/oxWnlEwERFEG2dVYTqsGbsNDh6E4KBB9a2qL5HJaBRmJuoX4DFNifyHhMVScg2wMkWr2g0PYMWc842J58v/eh3RAma85mekVfZWmOiIiEZPp7MsU07/MSHcejdkr50pUAdy0bt1rznRntMVCTBycr4AtbEJuw6gZpcdfsKfQN/My0FV+vWlKF9magEkey5huThB43apFLtfjIMJGrJfgz4i2fbRN5jokIeOfxqM8fYoOGEZkw+f+pH2zCEYqA8nALE9ImKqGerrcNqOVDp9fE3GSdP+SJ2ZoTLfs091uzfquz7N9Nq8/4+fY+KoHVpdgYYFlQzfaIii5adkgOHekaoKlGTafKemSYXs+L3qGiPiYogWg9gnr+ZloyAc5b5GGF1TGuskiCT8MSHiYqEgnmG3XJg5YBCsgjUWa2/hx4PIdKKjJfifcVERRr6/K5Dy6TIPIeFPiQ80SEvybiiEuhuyX4fZK36oAzVz0SFSHAyDCSkjnxpCmBjJuEVM610BAsTFUmE6nQzbcszQA5XGvMwnkl44i0WQTsr25UCwei5JBMVSTgOcHp9PJnU0PoeIiPVu8HHU+HxEykkY1bDREUSXE40L7+HFfnGJWGZcaIP1LzGEc0/ZlqdY6IioWAlxEyOfONr85hp4Ak086UT6vCYDjwzJrOBxkRFEg5f+BbM93WIgYOcHvyu9QA3m5pLXew55A6ndX3JcCWIiUoAXKxvxKv/PYo/fPQd8s/WqHpOrbXJ4+MtB/ptx8r9iE49V2cG52vrse67EtQ32oISg79aT6Lv7izEqJc2wdrYhIv1jSirrkOTTcGKrwqw73SlJu9ZUlnn92tsPnJWg0guW72zSLPX8pcMCYu1sQkllZew7Vg56hp+OC7Lquvw+ub8gLznJ98W44uDZwLy2i2t2Fbg1f4rvdi/rKoOSzcfx+7C83jty2MoqbyEZVuO4+iZapf7F567CADYWXDOq5iCqbSyDv/ecxqNTfqPd3mFF5y2FV+4JM1Y7I9wvQMwohfXH8Zb204AAP71dSFOvDAWgOdl16fXfufxNUNbfOxn0caj/gfpg8GpHTD+9e04WlaD39zUG0/d3leXOPzx+w8v1/MDb+7A/uIqXKxvwqxbr8LC7CMAYG+rYPA0ab+4/rD9v/094fnk22I/X0E+npLF87X1eGnDYfw8ozsGp3Z0eGx34XncvTjH/veY9K5Y8kAGAGDKim9woKRK81hzjpfj0dV5ALTrf2VVdZj81je4f0gqHhiaBgCosTZ6/TrzPjmget/Jb32Dgy3qp7kPf7j7ND5//Aan/X//4XcY3jvB65iCafTCzaixNqK0sg4P3di7zf1zjpdj3sf7AxLLc58ddPh776kLeOqD79C3a3skxkWpfh0JFlCccEUlAHbk+3+G0HoSCw/V/1z0mpQOOFp2eYXos70lOkfjn29OnMfF+stny29vP6lzNKS187X1bh+b/8l+rNpRiLtaJCTNfrnsa4e//7Ov1P7fgUhSAODB/9ml+Wv+bf1hHCypwh8+2mffdqne86qtO2Eqxp6QEDgkKS0dKnW9ogIAh92stoiiOblTu8J53xs7cOSMulV0AIiOCPMpLgBYm3cagOf6NQrTJyoyXJ/Tk+y3pfC+GmrtmIfLsVYdltHrAvCelxp8S0qC7VytVe8QDMHo05jpExVRGb3jkTN3SbMRkmml1f8H/v0MUGmCYKpPemOiQqoZYcIUCVd7AkfUvuoqrmD+yKgvxI6OzICJCrkV8v3/jE70fEH0+CjwBM9lyEcirPyJEENbmKiQajJ0aCPiJOU91llwsb7lIWNbMVGRhGyXCWSKVqZYSR0Jx2Jh8fggvTFRCSKuSBBRawEZFTR8UY5axiLZOS8AJirkgQLF8av7FcfHZCD6jYoUGDK1e7BCFfGYlaiZhCZTf/cFE5UAMHaXUYd14D13Zzqi16U3Y2TwJmUiMgomKoJqPWcZPWPWk15LoTIuwVJgiXicq+mmRu/LRiqfgF2sTUxUyC0zfDQZEL+cIk5eLYkYneBVRgRAjH4qQgxtYaJCqknQn0lwwTszFbO3ihmVZ4GM2UgrFRQ4TFSCSPQz97bIkHkTAZ77qozHodeHnoZFDGRtydgWWvInUTNT3TFRocBgUtMm5/uQdAnDb+YZLiWiYV/iqgfpzfSJiqRzQ9DIPkjJ9kV5pA0e10SuiX7PmyumT1SI9Mpl5Bsu5CHhWCws1iXpjYmKoEQfG9ocvLiQQS009xdOejoSsO5F/BK6YOLxoA4TFUnodQmj5Q1b3gwqouQpMi5zkv+MNgGKcjyRtozVSwOHiUoAGHZQ4VGlKbU5FHMt77HOiIyDiUoAcIxkHZgJ21o7IiZYAoZkGP6c1Pq6aihje/qUqCxevBg9e/ZEVFQUMjIysHXrVo/7r1q1CoMGDUJMTAySkpLwq1/9ChUVFT4FbBaGXZURkCh1zQ8oaUfUwdjXS5HePstol77I3LxOVNasWYPHH38cc+fORV5eHkaOHIkxY8agsLDQ5f5fffUVJk2ahKlTp2L//v1477338M0332DatGl+B0+Bx8kzcNTWLScd7/HeJO1wCCC9eZ2oLFy4EFOnTsW0adPQr18/LFq0CCkpKViyZInL/b/++mv06NEDM2fORM+ePfHjH/8YDz30EHbt2uV38BRcHPq1xXtUSAY8WTEWGZvTq0Slvr4eubm5yMrKctielZWFnJwcl88ZPnw4Tp06hXXr1kFRFJw5cwbvv/8+xo4d6/Z9rFYrqqqqHP6RXGQ8GMg3TKTIV2b6GnhXRDh0ZDh+vUpUysvL0dTUhMTERIftiYmJKC0tdfmc4cOHY9WqVZgwYQIiIyPRtWtXdOjQAf/85z/dvs+CBQsQHx9v/5eSkuJNmIYkwlJ2yxjaCkf/aOUjQBMbhsff+jH33CgWtgWp4NPNtK2/00NRFLff83HgwAHMnDkTf/rTn5Cbm4vPP/8cBQUFmD59utvXnzNnDiorK+3/ioqKfAlTPDwohST61+wLHh4JyOVKhYj9iMl50MlY5eHe7JyQkICwsDCn1ZOysjKnVZZmCxYswIgRI/Dkk08CAAYOHIjY2FiMHDkSzz33HJKSkpyeY7FYYLFYvAmNAoATZGCJsEpmVLwB2QVWiXA4xKrj1YpKZGQkMjIykJ2d7bA9Ozsbw4cPd/mcixcvIjTU8W3CwsIAcKAWHZuH1BAxoTV739UyUTP7fSSkP68v/cyaNQtvvvkmVqxYgYMHD+KJJ55AYWGh/VLOnDlzMGnSJPv+48aNw4cffoglS5YgPz8f27Ztw8yZM3H99dcjOTlZu5L4iMmSZy0nIW+qikMbkTFwdYrJmt68uvQDABMmTEBFRQWeffZZlJSUID09HevWrUNaWhoAoKSkxOE7VaZMmYLq6mq8+uqr+O1vf4sOHTrg5ptvxl//+lftSuGHQByCgUh+RL+PojUObW2TrU3d8aa7m33S87n05q42w9KqWY1+vu11ogIAM2bMwIwZM1w+tnLlSqdtjz76KB599FFf3sq0DDKHkQdczQscWatWxLDV1CVXHHykQ4O3bk8ZTh74Wz+CEnGgFTAkTTApDLzmwVDEfk1EYmOiQh61PFPiJKMP1jvpiYk86Y2JCrlllgFKr3Iy/wgc1i1JwSRjrL+YqBARGYw0q3CcqEkFJirkkVlWVUQgy9xCcmL/IkDOfsBEJQAC8bFTfkIkcPiJBf+o+dRAsOuYx4t21LSdDJ8cMTKj1z8TFXJLxo+xtSZfxETOZDz2SAUBmlWGnJ6JCqkmQ4f2BScB4xG2RYMUGC/ZGp+vq5Qydg0mKgFwuLTKadtne0vwxcEz9r//J+cEis5dDGZYPvG2U9tsClZuK0B1XWNA4vHWyxuPYn9xpd5huCRKHfmr+MIlvUMwvIYmBWerrbq8d32Trc193v1Gnl+4t9kUHD1T7fLyYGHFRby++ThqrI0oq6pz+XxFUfDOjkKXj3nNj6zBlxOsr46WI+d4ue9vqhOfvpmWPLO16j8VNVY8/M5uh23PfLwfz3y8X/VryvJ16x9/W4x5nxzQOwy7AyVVGPvKVzjxwli3+4h2j8rCDYf1DsEr/7frFAaldMD9Q9Lc7rOzoCKIEWl/1jj7vW81fkXXNh48g4/yTmPB3QPQISbS4bGRf/svDv15TFDi8FZe4QVV+4lw79Bznx3Eim0F+M1NvZ0eu+HFLwEAe09X4rO9JS6fv35/KZ5e+11AY1Tj6/xzXu1fXdeAB5bvCFA0gcUVlSCo8uHMWYS8xJcYDp+p1j6QAGtqkVnuyA/uhOrKK/89pncIXntxvefkarfKiSzQ8s/W+PS893NPaRyJaw/9by7+s68Uf3NRn3UNba9sNBMgHxDWim0FAIAlm4673cddkgIAh0rlG+MAoMYq7wouExXyyF2yIsKZkVZOt7h0sX7/GQ97kizc9c5aa1NQ4/BVWZU+l3koyIwzjAaU6RMVA823ZFDsomRYJu/cWt3I7888JkMTmD5RISIyA56UkayYqJBqRrrcQ/pgDyIHAtyLZwQi3NMYSExUBCVOTmDwI0AgTAQ1xKoUEru4o2B+4lDmumeiQiQ4popkVOzbpAYTFVJN4oScAkymszU9J0d+C7KcWl5akamvGwUTlSCQeUnf3bVPWb6AjszJ7AmBxEMOkRMmKkSC45xjRmx1MxDi48kSZLVMVAQl+oKFDJ2bxCN4tzYM0ccPCj6ZR2wmKpIQYdwxQ25i9ksGRmfm9hUheWld+zK2hgj1aDZMVMijlsekmQd5PRkpQQxWUYxUZ1phnYhHtB9EFRUTFSIiIh3w5E8dJiqkGrN/IiIKNiYq5BavxZJavLlaLGyOwGHdBp/pExUuvXnG2gke1rV2RK1LTnKOzH4upFV/UPM67k4mZOiSpk9UgkGGjuBK637NpE4fXNkio+DKG/mCiQr5hMMNUeCYZT7nN1yTGkxUiChognVG7fFtODmSINgV1WGiQkREuuClIFKDiYokRDiczfDxZI6bvpFlCd/M7ct7zMxN5r7PRIVU40BHRGanZU4uQvIgQgxtMX2iIkMjAfp8jE+Sk2TDk6GPiraEz6RaTEZoFcG6uimYPlEhEp2RJl1ZLhERycRIY4QrTFRIU5yGiNwz9nRCFBhMVEhTHIhJBFyeJzIOJipBIOug6SluWcskIyN92kq0e1lEZZZa4qVAUoOJCpEg3M3hRr/+TNrj/E+tyXyOwESFiMhgRJ2URI1LL0JUhxBBeMZEhYgMR4KxN+iYJGiD1Rh8TFSEpf/aLZePSS3hBm8PszK7NRmN0ZNQJiqkmtEPBlGx3rXDqiSzkvleNyYqRGQq8g7X6nE1NHC0rFqtXsvo7c1EhTRl8OOFJCFqMqLn6pjRJzMyLiYqQeDbANFqRBNgkFFTDlEnCCLSn8yXH0g/TFSCwIj3GHDAIQocf78Yz4hjjp4C9cV0IjSTDGM5ExVSjYMfyaK+0aZ3CERt4jc1q8NEhYgM550dhXqHoCvej6KtlgmFiKmFmnxH5pyIiQqRKCQeSERTXmPVOwQi0ggTlSC4VN/k92vokQ1XXWrE65uPe9znRHktis5dDFJEJINaayPuXbYdH39brHcowtly5Kxu7633GfXZaivKa+r1DcKDTzz019PnLwXkPbW692VHwTm/nr/vdCWWf1WAJpuYZ0vhegdgBuNe/crv1ygor9UgEu989l2Jx8drrY246aVNAIDjz9+BsFD91pvPVNVh9nvfYuLQNGT17+r0+Oz3vsWp820nVDa9R3MXZLuO3f+Z9QCAr/P9Gzz9UVXXqNlraZmIf76/VNV+xRfqNHtPEVgbm3DdXzY6bf9PG2NMMD26Os/tY/k6jL/B9JN/Xp6joiPCcN+QVJ2jccZERVjiX2Q+U/XDYNposyEsNEy3qP/40T5sPVqOrUfLceKFsU6Pv597StXriJiokL426bAKcu6i7ysPZ6vFu+xVebHB5fYvDpUFORJqreUq18GSKh0jcY+JCvnkTJUVN/99s95h2FXUirukTOJQFEWCUwDf/fOLo/h79hG9wyCVgrla6u6dJq/YGbQYfMV7VEhTXI8g0g+TlMAzWqJ7+kJg7r/RkukTFa70E5EnRpuYiGTjU6KyePFi9OzZE1FRUcjIyMDWrVs97m+1WjF37lykpaXBYrGgd+/eWLFihU8Bk37MkNOFCDgtMZkmEoeWh2OgvvHWV4KFY+f1PSpr1qzB448/jsWLF2PEiBF4/fXXMWbMGBw4cACpqa7vFh4/fjzOnDmD5cuX48orr0RZWRkaG7W7K5+CgxMmmZGogzeRWXidqCxcuBBTp07FtGnTAACLFi3C+vXrsWTJEixYsMBp/88//xybN29Gfn4+OnXqBADo0aOHf1ETkViYxJIH7B7kD68u/dTX1yM3NxdZWVkO27OyspCTk+PyOR9//DEyMzPxt7/9Dd26dcNVV12F2bNn49Il9zfwWK1WVFVVOfwzH/EObZ5Z6sNI9S7jqpwelwMD8Y4SVr3hBfVTPzIefN/zakWlvLwcTU1NSExMdNiemJiI0lLXX2SUn5+Pr776ClFRUVi7di3Ky8sxY8YMnDt3zu19KgsWLMD8+fO9CY10JvExIAx3v2LKuiUyJtEObVHPiXy6mbb1DUCKori9KchmsyEkJASrVq3C9ddfjzvuuAMLFy7EypUr3a6qzJkzB5WVlfZ/RUVFvoRJOtCro8t8tmAmeq8OiXbzIhG1zasVlYSEBISFhTmtnpSVlTmtsjRLSkpCt27dEB8fb9/Wr18/KIqCU6dOoU+fPk7PsVgssFgs3oRmQOINqMwF9MFqJ9mJN5qRTLxaUYmMjERGRgays7MdtmdnZ2P48OEunzNixAgUFxejpqbGvu3IkSMIDQ1F9+7dfQiZRMZJlUTGlTcSiWgJnKgrjl5f+pk1axbefPNNrFixAgcPHsQTTzyBwsJCTJ8+HcDlyzaTJk2y73/fffehc+fO+NWvfoUDBw5gy5YtePLJJ/HrX/8a0dHR2pWEiCgABB27SS+SJrtyRn2Z1x9PnjBhAioqKvDss8+ipKQE6enpWLduHdLS0gAAJSUlKCwstO/frl07ZGdn49FHH0VmZiY6d+6M8ePH47nnntOuFBQUHLD1Iem4KCz24+BjF9afzOOITz9KOGPGDMyYMcPlYytXrnTa1rdvX6fLRURkHBKPgUTe0TDT5XGjjul/64eIgkfGszqjLMDoWQ6j1CHpg4kKqSbjJEMko0BcnuLhS7IyfaLi7ku2yDc8cyKj4T0tRPpiosI8xRDYjEQkm+DmwPKOkqZPVIhEYYakmSuYROISdfWQiQqppqYTyz4NcSL1jaDjmybMkECSFzTsEOxa6jBREZSImS0HbKLg4LFmDqK1s2jxNGOiIihROwwFnwyrPOJHeJkvcYp40kBkJkxUiIgEw+SI9CBqv2OiQoYg6PFFBhDC3kUGIPMqPRMV0hSHdO3JPMCIyNs+KsOlNzV0/eVoIw0MGi47iLqCIRomKqQJTqbmxvFWW4ZbxeH4IAVR+x0TFTIEI4+DRi6bDEQdvEl+PMFTh4kKEZmKDHODUS43GRKzi6BjokKa4iGsPZ7P64tJAxmBzL2YiQpRCzxZCizWr0kx25aCqDf3MlERlJgdhrNMILmrXda6vniPCpG+TJ+ocBIgIiISl+kTFSIyD1kuPRluFUeSeicxMVEhL7Q9eBpseCUyDFmSNNKPqOM3ExXSBD8ZQWqIee8VBRzbXXcyJ6pMVMgL4vZ0mQ/CNhm6cBLQYZJlQkf0AyYqguLcRFJRObGK0K+9zgEEiJnIzJiokKY4ppsUG154XKWhtojaR5ioEBERkbDC9Q6AXLM2NukdgpPiC3V6h2BKsixW2GwKVuac8LjPB7tPBScYDRVU1Ab9PUsq61BV1+DysdLKOlysb0SvK9qpfr3/HjqDA8VVWoXntQsXXZdFJkfOVONgSZXwx2POsXKkdIpBSqcYh+0yf+CBiYqgHnt3D85U1eGWfolYsO6Q3uEAADYfOat3CAFXfOESXvvymN5hSOnf357Gs58ecPv4odIqbDtWEcSItLFk03Fd3nfgvA0utw9d8AUAIPcPo9G5nUXVa/165S7N4vLFI+/s1vX9tZD1jy0AgAHd4nWOxLP73twBADjxwliH7WruDwsR9NoPExWBPb/uEN7ZUYgTFRf1DqVNzV9QJWY3V+/Lw2fx5WF9EjJFhDtN/fD18XMeHy+p5IqclgrKa1UnKno7cqZG7xA0s7+4Uu8QTIf3qAiu8Jz4SQoRQuReWiYicTFRIU1wkqI2F4SE6CJCBEESYw8KPiYqRIKT5YqQJGESkWRMn6jIfl8AGYfUPVGRJ6EiMiOZj0/TJyqkLYmPBQowXh4kI5D9AwMyYqJCRP7jzbRBx9oms2CiQkSakGVpWdTviiDSm6hHhukTFUnGVjIBWSZ6Xxm9fEQUGKZPVEQn29mfXNGSlnhjupjYLtqStTZlvjTLRIWINCHvMEhEIjN9osKTDW2wHgNHljOhtvoA+wiR4ARdEjd9oiLJHEAm4C4hkWWClyFMWeqSiH7ARIUMgfOP/ngvRHCprW42CwFAk03ejmD6REWWZXVZsDbNS4qf+iEyqYYmeY9A0ycqRKQRecdBIhNo+wANEfQmFSYqZAxGWN+WvAhnq616h0AUcGJO5cZm+kRF9PlNtoNCtnhlIHgXtdt54pzHx0W4h0WyryUiAenfi83H9IkKEfmP87+4OLESoO6kXNREnokKGYKRB2NBxw5psT6J5GL6RMXIE1wwsR79J3Mdqold5vKJSIRLaUTBYPpEhUh0nI6IyMyYqJAh8OSSiMg9NUOkqJdFmaiQppgvmJOaAU6EZFKEGIKNl4hIdqZPVEQ/iEW9C5u0J3hXNAyjVHMIBwfSmKjHhukTFdKWXkOnkX8KgQkMEZkZExUiwRk5CSMicYi6Rmf6RIVTABERGZ3MK7OmT1SIRGH8lROjly+41N5fx1onQPz7MT3xKVFZvHgxevbsiaioKGRkZGDr1q2qnrdt2zaEh4fjmmuu8eVtSWAyHwRERCQurxOVNWvW4PHHH8fcuXORl5eHkSNHYsyYMSgsLPT4vMrKSkyaNAm33HKLz8EGAudXIiIicT9l6nWisnDhQkydOhXTpk1Dv379sGjRIqSkpGDJkiUen/fQQw/hvvvuw7Bhw3wOlsgdJpziE6WNBB2LiQJK1c9cCHKMtuZVolJfX4/c3FxkZWU5bM/KykJOTo7b57311ls4fvw4nnnmGVXvY7VaUVVV5fCP5CBoP5eCqIMEEclP5vHFq0SlvLwcTU1NSExMdNiemJiI0tJSl885evQofv/732PVqlUIDw9X9T4LFixAfHy8/V9KSoo3YXrF+DcwEpGZyTxBicjI9WmYSz+A8zciKori8lsSm5qacN9992H+/Pm46qqrVL/+nDlzUFlZaf9XVFTkS5ikA92+8M3Ag4dRiNBEIsRApAdRkxA11C1xfC8hIQFhYWFOqydlZWVOqywAUF1djV27diEvLw+PPPIIAMBms0FRFISHh2PDhg24+eabnZ5nsVhgsVi8CY3IuCSYXfl17sEnQbcg0oRXKyqRkZHIyMhAdna2w/bs7GwMHz7caf+4uDh899132LNnj/3f9OnT8aMf/Qh79uzBkCFD/IteCzzaDcEIzShzGfjxdCKxyXyIerWiAgCzZs3CxIkTkZmZiWHDhmHZsmUoLCzE9OnTAVy+bHP69Gm8/fbbCA0NRXp6usPzu3TpgqioKKftRGRsMg+UMuN9eNoKCZGzL8vcD7xOVCZMmICKigo8++yzKCkpQXp6OtatW4e0tDQAQElJSZvfqULGI+8hQEREIvM6UQGAGTNmYMaMGS4fW7lypcfnzps3D/PmzfPlbQOCE6wxGPnSg3FLRkTUNv7Wj+BC+PVUqhwqrdY7BL+9nXMCVXUNTtutDU06RKM9EZaeC8prHf7ee+oClm4+DptN/9h8YbMpWLmtAN8WXdA7FBLMqh0n8XV+hf3v/afl/T4yn1ZUiESSf7ZG7xA08fqWfJyoqHXa/j/bT+Le61PRLykOAHCuth7LtuQHOzyPZPnUz+/e34vMtI72v+98dRsA4ER5LV64Z6BeYfnsk73FmPfJAQDAiRfG6hyNOciyeDt37T4AwPw7+2Py8B74y7qDbT5H1BNjrqiQJvTs3merrTq+u7ayD5xxuf3VL4/Z//upD/Zi6ebjwQrJFN79Rr7valIU4MiZtlcSZZlYKTCe+Xi/3iH4zfSJCg9i+YWGinkWECh7JF3mF+VYk2TxhyjoRLg864rpExWSn8nyFCEZ+WZmItIXExWSniz3RxARiYz3qAhK1KUuUk/MQ0tjLbqpiOVVkyzySCMiXzBR4eipCT2rMdRAKypcHSK1QkLEPQMm0pLpExXhcRxqk5ESFSIivYg6lJo+UeGCivxEPbgCRdby8oZbbbE6ySxMn6iQ/OKiIvQOgYhIeqImv0xUiCRghJu+ef8NEfnC9IkKl6OJgoPHmj5Y7aSWqOcSpk9UiESiZpzgJz2omagTC5GWmKgQUVDw0g8R+cL0iQqXReVnhPs3zECUSz9ckSKSi+kTFdHJMqQKMgcZVsv65cIEAUzQyTyYqBARGRgTGpIdExUigci6WiJr3LJjtZMZMFEhkgwnJyIyEyYqJD3eH6M/tgGR/EQ9CWKiQiQBJgLaMds9G+w7JDvTJyqiH8SCh+eE9yoEHr+PhJrJNj4Q+cL0iQqRSNR8x4co30fSEnMnHYjXDYgCwvSJiujLwLKN/3p8mZbYLUgUOLKND0S+CNc7AK/U1gJhYdq+Zk0touvrtH1NDVlsoQhttOkdRttqagFbBCLqLga3PmtrEVLbqg1b/y2RiLAQhDU5p16R1kuX+z+A6IY64coXFd6E6PoGj/uEXQpy33DDYr3kOg4Z+k2LGEMv1iK85fHmLn4ZykWBpbIPhNddtI8zWr63v0IUEdeRW6mqqkJ8fDwqAcTpHQwRERGpUgUgHkBlZSXi4nybwU1/6YeIiIjEJdeln+JiwMeMzJ0T5bUY8/JWTV9TS5bwUFgluPTz9ZxbEB8TgZc3HsHSzflBe9+Df77dqQ0P/vl29Pvj50GLQUuRYaGob3Ju79FXd8E/f3ktAOCWv29C8QWxlvLjYyJQedHzpZ8XfzEQT763N0gRuZeZ1hG7Tp532i5Dv2kZ41tTMrHzxDks2ZTv9FhL3/xhNK57bmNQ4ySxqO3bv7mpN2be0kfbN6+qApKT/XoJuRKV2NjL/zRku6jgUmSUpq+pJVt4KKyh4icqiI0FYiLQEBUT3PqMjXVuw9hYodvUkyY3iUq9Jdre9+sio3EpMtiReWaJjMClRs/3jzVFB7lvuGG1RONS5CXnB2ToNy1itMXEojHa+kPM7uKXoVwUWCr7QGNUtOZzLJqa/H4JXvohTfFjqoEh/p1kFGzsEmQWpk9UeLATkYzUnhNI8HkJIo+YqAh+DMu2QiFZuFKSrU8QEfnD9IkKkVCYhBCRXgQ9C2KiQprQ8xt+BV8UIwqMEOa1ZA5MVDjNaUvQjFx2LXupiFUsYEhE5C1B74UwfaIiaLsQSYWHkbjYNiQ70ycqREREBDGXa8FEhWcbBmC2VTE9fqGaiEgvTFRMNsmR2JiCBAErmUgqpk9URCfb2bNc0cpD9ISa7a4PQVfqiTRl+kRFz4/VEvlC1slJ9GRLNmpPYljvpJaoQ4tcP0oYAA2NPIq1FOxJ9NuiC3h7+0mHbfUS/Nq0t05U1Oodgt/Ka6x6h2BKOcfK9Q6BJCHqbGj6ROUv6w7oHYJHlxr8/+VJI/vpa9ucti3/qsCr18hM64hdJ89rFVJAHCursf+3rGfIz312UO8QDKdlX6iua3C5z29W7Q5SNCSq87X1eofgF9Nf+vk6/5zeIZDG9hR5l3REhotzGEh7WQfAjxLb6x2GqV2q50kNuZZ98IzeIfhFnBGaDEG2m39lJGoyExYqaGAGJmpfINISExXSlAg3J8t6aURmnC+J5CfqccxEhYhIEgqzcPKF5N2GiQppipd+Ak/UGpZ8LJSO02UfUTsGSaNPYju9Q3CJiQoZjrcTpkgnqUz0iEgvfbuKeUM8ExXSlIw394lwX00zkWIh8TGxJTXUjyti9icmKmR6Iq2oqBEiYzYoEJlrz1NfZdJCRsVEhTQh22RP2uONnvriahwZFRMV0pQI53ScL4mIfiD7mMhEhUggapbvRUgGW5N8HJRGy3oWsR+Q3ES9qsxEhUyPkywZAe9RIaPyKVFZvHgxevbsiaioKGRkZGDr1q1u9/3www9x66234oorrkBcXByGDRuG9evX+xwwUduYegSb7EvLshL1DJjEIvvh6XWismbNGjz++OOYO3cu8vLyMHLkSIwZMwaFhYUu99+yZQtuvfVWrFu3Drm5uRg1ahTGjRuHvLw8v4MncsXrSVP2o5iIyMC8TlQWLlyIqVOnYtq0aejXrx8WLVqElJQULFmyxOX+ixYtwu9+9ztcd9116NOnD55//nn06dMHn3zyid/Bk3h4hhcErGPTavnJKn5MnczCq0Slvr4eubm5yMrKctielZWFnJwcVa9hs9lQXV2NTp06ud3HarWiqqrK4R+RWlwgIU+M0j9afxyceQu5o3aVWdQu5FWiUl5ejqamJiQmJjpsT0xMRGlpqarX+Pvf/47a2lqMHz/e7T4LFixAfHy8/V9KSoo3YRJ5RaTvn5B1suF3qBBRoPh0M23rJUdFUVQtQ65evRrz5s3DmjVr0KVLF7f7zZkzB5WVlfZ/RUVFvoRJQdQ8TXE5OvBYw36SOKdy+HhyiOPnfJgrklGFe7NzQkICwsLCnFZPysrKnFZZWluzZg2mTp2K9957D6NHj/a4r8VigcVi8SY0Ijue3RMR/UCkVWNfeLWiEhkZiYyMDGRnZztsz87OxvDhw90+b/Xq1ZgyZQreeecdjB071rdIiUhYUg2DBlqSclxh0S0MMghRV8S9WlEBgFmzZmHixInIzMzEsGHDsGzZMhQWFmL69OkALl+2OX36NN5++20Al5OUSZMm4eWXX8bQoUPtqzHR0dGIj4/XsChEvuECjAYUiepRljjbIOicQqQ5rxOVCRMmoKKiAs8++yxKSkqQnp6OdevWIS0tDQBQUlLi8J0qr7/+OhobG/Hwww/j4Ycftm+fPHkyVq5c6X8JiPwk0rwl69wjUh0amTTJIJGGvE5UAGDGjBmYMWOGy8daJx+bNm3y5S2IfMaxPPhkui9I9uv1zVpXuaxJLgWeRIenS/ytHyLJiHgdWQEvReiBVU5aErU/MVEhTXGyMi/Zz9pk0HI1iMcamQUTFTIcTpj6kOWSCvsHmY3sXZ6JCmlCpHsUvP5NQoFiV3NZR8QTaYGq0LREvCRIpAUmKqSpECGnUQo0WVZTjKT1kSZSwk2kJSYqRGQqMk/nTp/04XkBqaEyiRW1PzFRISK/KTJ94RsRSYWJCpmebPOriGc9stWhEfEeFTIqJipkON5eq+dKgDZkmSd5LweZjew93tSJyitfHNU7BMORZbIijck+EhKRsB+GMHWicrSsRu8QiLwm4mAi06d+5Im0bbzcQ2Zg6kSFhziR+ch85adl7K1zFI5n5I7MfR4we6LCI1szzccBq9ScZPrUjyRhEtH3zJ2o6B0AkQ+YYBNRIIg6tpg6USHtidrRPRHpUyASVh8AyVYpBGpvImqbqRMV3oimPRHmABFiMBubRJUuT6TOZLppmcQh0smYL8ydqOgdAImBCSsRkbDC9Q5AV5yfNFNyoQ6Zz23UOwwAwFfHyvUOwWfV1ka3jy3bchy3Xt01iNGopyjyrFTIfHL5p3/vb/EXBzAyB1MnKiJ+H4Wsxr36ld4h+Cw6Qo6FxefXHcIL/zkEm6ATrSzLyzJfPnk/95TD3wdLquz/zYVBMio5RmiiAIqOCNM7BNVETVJIHyWVdXqHQBRwpk5UeAZCgDyXLEQny83pRlpJNU5JiNwzd6KidwBERH6QJDckncl+MmbuRIUHOZFmZLlHxUiMtDpE+hN1TjR3osKDnCD3p0BEwmokokAwd6LCPIVIMzycgiskBKx0UkX2kzFTJypERDJjnkJmYOpEhSsqRCQzjmGkJVE/uWfqRIXnI0QkM95nR2pIfuXH3ImKoMkjBZnsBzERkZGZO1HROwASgmwfq03vFqd3CC7JVYvGwJMtMgNzJyo8yElCMRGm/okuaoFjGKkh28lYa6ZOVIikxMmJcLkbtLxHRfK5iAQg6tBi6kSFN6IRkcy4okJmYO5EhQc5SUjYbsszeiIKAHMnKnoHQOQDJthEZCbmTlQ44hMREQEQ9yTI1IkKESDfTYi8t8o/og7GvuDJFpkBExUiIkkxTSE1ZDsZa83UiQpPRggAFMnuAhW138pSi7IP2s24mkJmYe5EhecjRJqR/UuliMxO1DnR3ImKmG1CJCVZzvAlCVMVI5WFAke2VePWzJ2o6B0AkQ9EnZwEDcvQWta53FMRkXumTlSISDucKINPllUsIn+YOlHhMU6AfDdXinodWZZ7VCQJUxUxewKJRvY+b/JEhYc5yXcQs9sS8P2PErIvkIZE7U/mTlT0DoDIQGRJ/CUJ02uyrGgRecvUiQozFSLtcKIkEpPsR6a5ExUiMh1j5VM/nG0ZqlhELZg6URH1pkQiT2S5xCIqI1Vfy7IYKwEjPYh6aJg7URG1VSioZPsyJFG7rVy1KD+bogjbF4i0ZO5ERe8AiIh8ZGuVGcqWcFPwyL7aZu5EhZkKQb6DmP2WgMs3L/PSD5mBuRMVrqkQkaSabArCQn8Yw2zMVMhfgk6Jpk5UiGQk6FgiDaPUn00BwkJ/GMIbm5iokGuyXxY0daLCJXQC5LsJVNRP/chyQi9JmG2yKQrCW6yoNLW+aYXIIMydqOgdABGRj2yK46WfJlkyRSIv+ZSoLF68GD179kRUVBQyMjKwdetWj/tv3rwZGRkZiIqKQq9evbB06VKfgtWcoGemRJ6w1xJweQWlRZ7CFRVyS/Yc1utEZc2aNXj88ccxd+5c5OXlYeTIkRgzZgwKCwtd7l9QUIA77rgDI0eORF5eHp5++mnMnDkTH3zwgd/B+4sDPpF2Wl6GoMC7/D0qP9Q571Ehf4n6AROvE5WFCxdi6tSpmDZtGvr164dFixYhJSUFS5Yscbn/0qVLkZqaikWLFqFfv36YNm0afv3rX+Oll17yO3giMxJ1ITA8TNDAWpEjyrbZbI5/c0WFjMqrRKW+vh65ubnIyspy2J6VlYWcnByXz9m+fbvT/rfddht27dqFhoYGl8+xWq2oqqpy+BcIog74FFwnK2r1DsEr/z1UpncILh05U6N3CKp8e6rS5fb5n+wPciT+Wb2zEGt2Fdn/fn3LcR2jIZG9uP6w3iH4xatEpby8HE1NTUhMTHTYnpiYiNLSUpfPKS0tdbl/Y2MjysvLXT5nwYIFiI+Pt/9LSUnxJkzV4qMjAvK6JJczVVa9Q/CKlifO7S3h2r2Y5N7adkLvELzyRauE9dO9JTpFQkYRHRmmdwgu+XQzbeuPR17+hkT3yxOu9ne1vdmcOXNQWVlp/1dUVORyP3/9PKM7HrulD+KjIxAWGoIHR/Zs8zkTh6ahW4doAEBmWkdcm9rB5X7REWEY1qsz+nZtb982sHs8buufiMnD0tDeEo7renR0eM4dA7oirXOMw7bUTjHo27W92/e5sks7AEBsqw4WGX65aXt0jsHsrKvw84zuTs+dPCwN0RFhSIyzoHNspPtCt2Fg93g8MDTV/rerWBPaeX795nrqGBMBy/exx0dHoHNsJEJCgKT4KPu+HWLaTjB/eX2q07buHaPt/31b/0T0T47D8N6d8fCo3rjn2u4Y2D3e4VMUrrRrMbHPuvUqjM/sjv88NhK39U/08CzYy9RSQjsLxqR3xeh+XTCkZyf0T47DXYO7oWNMBBLjLLjxqitwdVIcAGBor05I7RSDuwd3w8OjeqNnQix+4aJNO8ZEYPKwNMwbd7XD9qT4KGSmOfa3heMHYcOsG9DRTX2O7tcFYwckOWzr3jEaPb7vozdedYV9+4ybeuPhUb0d9o2OCHN4XkqnaIfHr+zSDoNSOtj/vr1/V/S+IhZjByTh5XuvcRlTs/ZR4fbXuK1/IgZ0i7c/1rdreyTGWQAA9w9Jxc19u9gfiwwPxcOjeqN7x2iEhgC3Xp3oEH/zGJDQzoJxg5Lt/aidJRzXfB/rj69M8BhbaqcYp20Th6a53N5aUnyU/fVHXNkZwOV7fwZ2j0dkWCiS4qNw1+BuuKVvFzw86oc6T2hnwcOjettXiR8e1RvXpnZAdESYfSxo9vxdAxDz/XhxXY+OWD45E/+YMAj9kuJw04+uwJ9/2h/LJ2diaK9OGDcoGbf07YJNs2/C7Kyr7K/Rui3d+fWInvZyudKjc9t1AgBXtLcgNjIMvRJi7dua6+n2/l3RMyEWVyfFORxn7SzhGNQ93um1WrpjQFeX25vbqnmsmTK8B+6+thsAoFNsJPolxaFbh2iMHZiEx0f3QWRYKDq1GEN7JsTil9en4IGhqega90PZr+/ZCXcM6IrsJ27A9T07Ob1vc/w3tDi2msfOmBZjfNb3/faqxHYOY92U4T3wk4GOx2yHmAg8PKo37hyU7LC9uQ2by7h8cqbD+CaSEEVRfz9wfX09YmJi8N577+Guu+6yb3/sscewZ88ebN682ek5N9xwAwYPHoyXX37Zvm3t2rUYP348Ll68iIiItiedqqoqxMfHo7KyEnFxcWrDJSIiIh1pMX97taISGRmJjIwMZGdnO2zPzs7G8OHDXT5n2LBhTvtv2LABmZmZqpIUIiIiMi+vL/3MmjULb775JlasWIGDBw/iiSeeQGFhIaZPnw7g8mWbSZMm2fefPn06Tp48iVmzZuHgwYNYsWIFli9fjtmzZ2tXCiIiIjIkry9ITZgwARUVFXj22WdRUlKC9PR0rFu3DmlpaQCAkpISh+9U6dmzJ9atW4cnnngCr732GpKTk/HKK6/gnnvu0a4UREREZEhe3aOiF96jQkREJJ+g36NCREREFExMVIiIiEhYTFSIiIhIWExUiIiISFhMVIiIiEhYTFSIiIhIWExUiIiISFhMVIiIiEhYTFSIiIhIWGL+pnMrzV+eW1VVpXMkREREpFbzvO3Pl+BLkahUV1cDAFJSUnSOhIiIiLxVXV2N+Ph4n54rxW/92Gw2FBcXo3379ggJCdE7HK9VVVUhJSUFRUVFpvmtIjOWGWC5WW7jM2OZAZbb13IrioLq6mokJycjNNS3u02kWFEJDQ1F9+7d9Q7Db3Fxcabq4IA5ywyw3GZjxnKbscwAy+0LX1dSmvFmWiIiIhIWExUiIiISFhOVILBYLHjmmWdgsVj0DiVozFhmgOVmuY3PjGUGWG49yy3FzbRERERkTlxRISIiImExUSEiIiJhMVEhIiIiYTFRISIiImExUXHh9OnTeOCBB9C5c2fExMTgmmuuQW5urv3xDz/8ELfddhsSEhIQEhKCPXv2OL2G1WrFo48+ioSEBMTGxuLOO+/EqVOnHPY5f/48Jk6ciPj4eMTHx2PixIm4cOGCwz6FhYUYN24cYmNjkZCQgJkzZ6K+vj4Qxdak3DfddBNCQkIc/t17770O+8hU7oaGBjz11FMYMGAAYmNjkZycjEmTJqG4uNjhNYzW3mrLLVt7t9XH582bh759+yI2NhYdO3bE6NGjsWPHDofXMFpbqy23bG2tptwtPfTQQwgJCcGiRYscthuxvVtyV26h2lshB+fOnVPS0tKUKVOmKDt27FAKCgqUjRs3KseOHbPv8/bbbyvz589X3njjDQWAkpeX5/Q606dPV7p166ZkZ2cru3fvVkaNGqUMGjRIaWxstO9z++23K+np6UpOTo6Sk5OjpKenKz/5yU/sjzc2Nirp6enKqFGjlN27dyvZ2dlKcnKy8sgjjwhb7htvvFF58MEHlZKSEvu/CxcuOOwjU7kvXLigjB49WlmzZo1y6NAhZfv27cqQIUOUjIwMh9cxWnurLbdM7a2mj69atUrJzs5Wjh8/ruzbt0+ZOnWqEhcXp5SVldn3MVpbqy23TG2tttzN1q5dqwwaNEhJTk5W/vGPfzg8ZsT2VlNukdqbiUorTz31lPLjH/9Y1b4FBQUuJ+wLFy4oERERyrvvvmvfdvr0aSU0NFT5/PPPFUVRlAMHDigAlK+//tq+z/bt2xUAyqFDhxRFUZR169YpoaGhyunTp+37rF69WrFYLEplZaWvRXRJi3IryuXO/dhjj7l9rszlbrZz504FgHLy5ElFUYzf3s1al1tR5GpvX8pcWVmpAFA2btyoKIp52rp1uRVFrrZWFPXlPnXqlNKtWzdl3759SlpamsOEbeT29lRuRRGrvXnpp5WPP/4YmZmZ+MUvfoEuXbpg8ODBeOONN7x6jdzcXDQ0NCArK8u+LTk5Genp6cjJyQEAbN++HfHx8RgyZIh9n6FDhyI+Pt5hn/T0dCQnJ9v3ue2222C1Wt0u4/lKi3I3W7VqFRISEtC/f3/Mnj3b/uvXgDHKXVlZiZCQEHTo0AGAedq7dbmbydLe3pa5vr4ey5YtQ3x8PAYNGgTAHG3tqtzNZGlrQF25bTYbJk6ciCeffBL9+/d3eg2jtndb5W4mSnszUWklPz8fS5YsQZ8+fbB+/XpMnz4dM2fOxNtvv636NUpLSxEZGYmOHTs6bE9MTERpaal9ny5dujg9t0uXLg77JCYmOjzesWNHREZG2vfRihblBoD7778fq1evxqZNm/DHP/4RH3zwAe6++27747KXu66uDr///e9x33332X+gywzt7arcgFztrbbMn376Kdq1a4eoqCj84x//QHZ2NhISEuyxGrWtPZUbkKutAXXl/utf/4rw8HDMnDnT5WsYtb3bKjcgVntL8evJwWSz2ZCZmYnnn38eADB48GDs378fS5YswaRJk/x6bUVREBISYv+75X/7s48WtCr3gw8+aP/v9PR09OnTB5mZmdi9ezeuvfZaAPKWu6GhAffeey9sNhsWL17c5msbpb09lVum9lZb5lGjRmHPnj0oLy/HG2+8gfHjx2PHjh0uB2V/yiNaW7dVbpnaGmi73Lm5uXj55Zexe/dur99b5vZWW26R2psrKq0kJSXh6quvdtjWr18/FBYWqn6Nrl27or6+HufPn3fYXlZWZs8uu3btijNnzjg99+zZsw77tM46z58/j4aGBqcs1V9alNuVa6+9FhERETh69CgAecvd0NCA8ePHo6CgANnZ2Q6rCkZub0/ldkXk9lZb5tjYWFx55ZUYOnQoli9fjvDwcCxfvtweq1Hb2lO5XRG5rYG2y71161aUlZUhNTUV4eHhCA8Px8mTJ/Hb3/4WPXr0sMdrtPZWU25X9GxvJiqtjBgxAocPH3bYduTIEaSlpal+jYyMDERERCA7O9u+raSkBPv27cPw4cMBAMOGDUNlZSV27txp32fHjh2orKx02Gffvn0oKSmx77NhwwZYLBZkZGT4VD53tCi3K/v370dDQwOSkpIAyFnu5sn66NGj2LhxIzp37uywv1Hbu61yuyJye/vaxxVFgdVqBWDctnalZbldEbmtgbbLPXHiROzduxd79uyx/0tOTsaTTz6J9evXAzBme6sptyu6trfq225NYufOnUp4eLjyl7/8RTl69KiyatUqJSYmRvnXv/5l36eiokLJy8tTPvvsMwWA8u677yp5eXlKSUmJfZ/p06cr3bt3VzZu3Kjs3r1bufnmm11+pG3gwIHK9u3ble3btysDBgxw+dGuW265Rdm9e7eyceNGpXv37gH5SJsW5T527Jgyf/585ZtvvlEKCgqUzz77TOnbt68yePBgacvd0NCg3HnnnUr37t2VPXv2OHxUz2q12l/HaO2tptyytXdbZa6pqVHmzJmjbN++XTlx4oSSm5urTJ06VbFYLMq+ffvsr2O0tlZTbtnaWk25XXH16RejtbeacovW3kxUXPjkk0+U9PR0xWKxKH379lWWLVvm8Phbb72lAHD698wzz9j3uXTpkvLII48onTp1UqKjo5Wf/OQnSmFhocPrVFRUKPfff7/Svn17pX379sr999+vnD9/3mGfkydPKmPHjlWio6OVTp06KY888ohSV1cnZLkLCwuVG264QenUqZMSGRmp9O7dW5k5c6ZSUVEhbbmbP4rt6t+XX35p389o7a2m3DK2t6cyX7p0SbnrrruU5ORkJTIyUklKSlLuvPNOZefOnQ6vYbS2VlNuGdu6rXK74ipRMVp7u9K63KK1d4iiKIr69RciIiKi4OE9KkRERCQsJipEREQkLCYqREREJCwmKkRERCQsJipEREQkLCYqREREJCwmKkRERCQsJipEREQkLCYqREREJCwmKkRERCQsJipEREQkLCYqREREJKz/D/0McuA2k+PVAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "cloud_data = CloudData(start_time, scale=1)\n", + "\n", + "sunsets = almanac.sunsets[start_idx : start_idx + 3651]['sun_n12_setting']\n", + "sunrises = almanac.sunsets[start_idx : start_idx + 3651]['sun_n12_rising']\n", + "\n", + "stepsize = 5/60/24.0 # stepsize in days\n", + "years10 = np.arange(sunsets[0], sunrises[-1] + stepsize/2., stepsize)\n", + "year_mask = np.zeros(len(years10))\n", + "\n", + " \n", + "for i, (s, r) in enumerate(zip(sunsets, sunrises)):\n", + " nn = np.where((years10 >= s) & (years10 <=r), 1, 0)\n", + " year_mask += nn\n", + "\n", + "x = (years10 - cloud_data.start_time.mjd) * 24 * 60 * 60.\n", + "cloud_cover = np.interp(x, cloud_data.cloud_dates, cloud_data.cloud_values)\n", + "plt.plot(years10, cloud_cover)\n", + "cloud_limit = 0.3\n", + "plt.axhline(cloud_limit, color='r')\n", + "\n", + "cloud_down_mask = np.where(cloud_cover >= cloud_limit, 0, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "id": "6ea96902-18d7-4a18-98fd-a92b2118f011", + "metadata": {}, + "outputs": [], + "source": [ + "scheduled_downtime = ScheduledDowntimeData(start_time)()\n", + "\n", + "if 'v3.5' in run_name:\n", + " unscheduled_downtime = UnscheduledDowntimeData(start_time)()\n", + "else:\n", + " unscheduled_downtime = UnscheduledDowntimeMoreY1Data(start_time)()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "id": "bdcb7055-e898-42a4-b4e5-acba38f463be", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([(