ddf_slopes_tests.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0d4e5edf-445b-4d1a-8b6a-156cb4432797",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'MJD')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAG2CAYAAACqHSU4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+JElEQVR4nO3df1xUdaL/8ff4g1ERJsFgYEXW3dRU1DVNxdrUVNANqbUbbharZWppGhfcNvXbil2TdPeqraapa2L+osc+0lZuReLd1Ap/km7+yrVNTVoQa3Eww8HwfP/oerYRFNCBOcDr+XjM4+Gc+cw5n3Pa5cVhzszYDMMwBAAALKGRrycAAAD+jTADAGAhhBkAAAshzAAAWAhhBgDAQggzAAAWQpgBALAQwgwAgIUQZgAALIQwAwBgIYQZAIAq2LFjh4YPH67w8HDZbDa99dZbHo8bhqHU1FSFh4erefPmGjBggA4fPlzt7RBmAACq4MKFC+revbsWL15c4ePz5s3T/PnztXjxYu3du1dOp1NDhgzR+fPnq7Udy4Q5LS1NNptNSUlJ5rKq/Pbhdrs1efJktW7dWv7+/oqPj1deXp7HmKKiIiUmJsrhcMjhcCgxMVHnzp2rhb0CANQXw4YN0+zZszVixIhyjxmGoYULF2rGjBkaMWKEoqKitHr1an377bdav359tbbTxFsTvhl79+7V8uXL1a1bN4/lV377SE9PV4cOHTR79mwNGTJEx44dU0BAgCQpKSlJmZmZysjIUHBwsFJSUhQXF6fc3Fw1btxYkjRq1Cjl5eUpKytLkjR+/HglJiYqMzOzynO8fPmy/vnPfyogIEA2m81Lew4A9YthGDp//rzCw8PVqFHNnftdvHhRpaWlN70ewzDK/Uy32+2y2+3VWs+JEydUUFCgmJgYj/X0799fOTk5mjBhQrUm5VPnz5832rdvb2RnZxv9+/c3nnnmGcMwDOPy5cuG0+k0XnrpJXPsxYsXDYfDYbz66quGYRjGuXPnjKZNmxoZGRnmmC+//NJo1KiRkZWVZRiGYRw5csSQZOzatcscs3PnTkOS8emnn1Z5nqdPnzYkcePGjRu3KtxOnz59M2m4rpKSEqNVq1ZemWfLli3LLZs5c2alc5BkbNq0ybz/0UcfGZKML7/80mPcuHHjjJiYmGrtn8/PmCdNmqT77rtPgwcP1uzZs83lVfntIzc3V5cuXfIYEx4erqioKOXk5Cg2NlY7d+6Uw+FQnz59zDF9+/aVw+FQTk6OOnbsWOG83G633G63ed/4v6+tXrVqlVq0aOG1/QeA+uTbb7/VY489Zv5VsyaUlpaqqKjopn8eX5nr6dOnFRgYaC6v7tnyD1199m1UcEZeGZ+GOSMjQx9//LH27t1b7rGCggJJUmhoqMfy0NBQnTp1yhzj5+enVq1alRtz5fkFBQUKCQkpt/6QkBBzTEXS0tI0a9ascstbtGhBmAGgErXxkp+3fh4HBgZ6hPlGOJ1OSd83JywszFxeWFhYrmOV8dnFX6dPn9YzzzyjtWvXqlmzZtccdyO/fVw9pqLxla1n2rRpcrlc5u306dPX3SYAoOFq166dnE6nsrOzzWWlpaXavn27+vXrV611+eyMOTc3V4WFherZs6e5rKysTDt27NDixYt17NgxSdf/7cPpdJp/0vjhWXNhYaF5IJxOp86cOVNu+2fPnr3ubzE38uI/AKD++uabb/TZZ5+Z90+cOKEDBw4oKChIbdu2VVJSkubMmaP27durffv2mjNnjlq0aKFRo0ZVazs+O2MeNGiQDh48qAMHDpi3Xr166ZFHHtGBAwf0k5/8pNLfPnr27KmmTZt6jMnPz9ehQ4fMMdHR0XK5XNqzZ485Zvfu3XK5XNX+LQYA0HDt27dPPXr0UI8ePSRJycnJ6tGjh373u99Jkp599lklJSVp4sSJ6tWrl7788ktt2bKl2q+3++yMOSAgQFFRUR7L/P39FRwcbC6v7LcPh8OhsWPHKiUlRcHBwQoKCtLUqVPVtWtXDR48WJLUqVMnDR06VOPGjdOyZcskff92qbi4uGte+AUAwNUGDBhgXghcEZvNptTUVKWmpt7Udnx+Vfb1PPvssyopKdHEiRNVVFSkPn36lPvtY8GCBWrSpIkSEhJUUlKiQYMGKT093XwPsyStW7dOU6ZMMa/ejo+Pv+YntwAA4Es243r5h6m4uFgOh0NvvPEGV2UDwDV8++23GjlypFwu101f6Xwt3vp5XBtzvRGW+UhOAABAmAEAsBTCDACAhRBmAAAshDADAGAhhBkAAAshzAAAWAhhBgDAQggzAAAWQpgBALAQwgwAgIUQZgAALIQwAwBgIYQZAAALIcwAAFgIYQYAwEIIMwAAFkKYAQCwEMIMAICFEGYAACyEMAMAYCGEGQAACyHMAABYCGEGAMBCCHMddEerU76eAgCghjTx9QTwbwQXAECYawHBBQBUFWG+SUQXAOBNhLmaut/yhQL8m/t6GgBgGT88QTnvV+LDmdQPhBkAUCH+IugbhBkAGhiCa22EGQDqCYJbPxBmALA4gtuwEGYA8BGCi4oQZgCoAUQXN8qnH8m5dOlSdevWTYGBgQoMDFR0dLTeffdd8/ExY8bIZrN53Pr27euxDrfbrcmTJ6t169by9/dXfHy88vLyPMYUFRUpMTFRDodDDodDiYmJOnfuXG3sIoB65o5Wp6p0A26UT8+Y27Rpo5deekm33XabJGn16tW6//77tX//fnXp0kWSNHToUK1atcp8jp+fn8c6kpKSlJmZqYyMDAUHByslJUVxcXHKzc1V48aNJUmjRo1SXl6esrKyJEnjx49XYmKiMjMza2M3AdQBxBRW4dMwDx8+3OP+iy++qKVLl2rXrl1mmO12u5xOZ4XPd7lcWrlypdasWaPBgwdLktauXauIiAht3bpVsbGxOnr0qLKysrRr1y716dNHkrRixQpFR0fr2LFj6tixYw3uIQBfI7j1181+4JNVPwzFMq8xl5WV6c9//rMuXLig6Ohoc/m2bdsUEhKiW265Rf3799eLL76okJAQSVJubq4uXbqkmJgYc3x4eLiioqKUk5Oj2NhY7dy5Uw6Hw4yyJPXt21cOh0M5OTnXDLPb7Zbb7TbvFxcXe3uXAdwEgov6yudhPnjwoKKjo3Xx4kW1bNlSmzZtUufOnSVJw4YN00MPPaTIyEidOHFCzz//vO69917l5ubKbreroKBAfn5+atWqlcc6Q0NDVVBQIEkqKCgwQ/5DISEh5piKpKWladasWV7cUwBVQXDR0Pk8zB07dtSBAwd07tw5vfnmmxo9erS2b9+uzp07a+TIkea4qKgo9erVS5GRkXr77bc1YsSIa67TMAzZbDbz/g//fa0xV5s2bZqSk5PN+8XFxYqIiKju7gH4AaILVM7nYfbz8zMv/urVq5f27t2rl19+WcuWLSs3NiwsTJGRkTp+/Lgkyel0qrS0VEVFRR5nzYWFherXr5855syZM+XWdfbsWYWGhl5zXna7XXa7/ab2DWgoCC7gPT4P89UMw/B4bfeHvv76a50+fVphYWGSpJ49e6pp06bKzs5WQkKCJCk/P1+HDh3SvHnzJEnR0dFyuVzas2ePevfuLUnavXu3XC6XGW8AFSO4QO3zaZinT5+uYcOGKSIiQufPn1dGRoa2bdumrKwsffPNN0pNTdWDDz6osLAwnTx5UtOnT1fr1q31y1/+UpLkcDg0duxYpaSkKDg4WEFBQZo6daq6du1qXqXdqVMnDR06VOPGjTPPwsePH6+4uDiuyEaDRXAB6/JpmM+cOaPExETl5+fL4XCoW7duysrK0pAhQ1RSUqKDBw/q9ddf17lz5xQWFqaBAwfqjTfeUEBAgLmOBQsWqEmTJkpISFBJSYkGDRqk9PR08z3MkrRu3TpNmTLFvHo7Pj5eixcvrvX9BWoawQXqPpthGIavJ1EXFBcXy+Fw6NP3fn9T75vzlo+LIn09BdQigou64vyFEt0e+xu5XC4FBgbWyDa89fO4NuZ6Iyz3GjPQkBBcAFcjzEANIboAbgRhBqqJ4AKoSYQZ+D8EF4AVEGbUewQXQF1CmFFnEVwA9RFhhuUQXAANGWFGrSG4AFA5woybRnABNATfffedUlNTtW7dOhUUFCgsLExjxozR//t//0+NGjXy2nYIM66J4ALAv82dO1evvvqqVq9erS5dumjfvn167LHH5HA49Mwzz3htO4S5ASK4AFB9O3fu1P3336/77rtPkvTjH/9YGzZs0L59+7y6HcJcjxBcAKi+4uJij/t2u112u73cuLvvvluvvvqq/v73v6tDhw7629/+pg8//FALFy706nwIcx0UfusKSeN8PQ0AqBciIiI87s+cOVOpqanlxv32t7+Vy+XS7bffrsaNG6usrEwvvviiHn74Ya/OhzBbyPfBBQDUptOnT3t8u1RFZ8uS9MYbb2jt2rVav369unTpogMHDigpKUnh4eEaPXq01+ZDmGsJ0QUAawoMDKzS1z7+5je/0XPPPadf/epXkqSuXbvq1KlTSktLI8xWQnABoGH49ttvy70tqnHjxrp8+bJXt0OYqyms9WoFtmzs62kAgGX88ASluHmZD2dSs4YPH64XX3xRbdu2VZcuXbR//37Nnz9fjz/+uFe3Q5gBABXiL4KeFi1apOeff14TJ05UYWGhwsPDNWHCBP3ud7/z6nYIMwA0MAT3xgQEBGjhwoVef3vU1QgzANQTBLd+IMwAUAcQ3YaDMAOADxFcXI0wA0ANILi4UYQZAKqB4KKmEWYAEMGFdRBmAPUawUVdQ5gB1EkEF/UVYQZgOUQXDRlhBlBrCC5QOcIM4KYRXMB7CDOAayK4QO0jzEADRHAB6yLMQD1CcIG6jzADdQTRBRoGwgz4GMEF8EONfLnxpUuXqlu3bgoMDFRgYKCio6P17rvvmo8bhqHU1FSFh4erefPmGjBggA4fPuyxDrfbrcmTJ6t169by9/dXfHy88vLyPMYUFRUpMTFRDodDDodDiYmJOnfuXG3sIhqw8FtXVOkGAD/k0zPmNm3a6KWXXtJtt90mSVq9erXuv/9+7d+/X126dNG8efM0f/58paenq0OHDpo9e7aGDBmiY8eOKSAgQJKUlJSkzMxMZWRkKDg4WCkpKYqLi1Nubq4aN24sSRo1apTy8vKUlZUlSRo/frwSExOVmZnpmx1HnUZMAdQkm2EYhq8n8UNBQUH6/e9/r8cff1zh4eFKSkrSb3/7W0nfnx2HhoZq7ty5mjBhglwul2699VatWbNGI0eOlCT985//VEREhN555x3Fxsbq6NGj6ty5s3bt2qU+ffpIknbt2qXo6Gh9+umn6tixY5XmVVxcLIfDIVdulAJbNq6Zna+Gf54d5+sp1DsEF7h5xd+UydHzkFwulwIDA2tmG//38/jT936vAP/mN7ye8xdKdHvsb2p0rjfCMq8xl5WV6c9//rMuXLig6OhonThxQgUFBYqJiTHH2O129e/fXzk5OZowYYJyc3N16dIljzHh4eGKiopSTk6OYmNjtXPnTjkcDjPKktS3b185HA7l5ORcM8xut1tut9u8X1xcXAN7jdpAcIH6Kaz16ps6USpuXubF2XiPz8N88OBBRUdH6+LFi2rZsqU2bdqkzp07KycnR5IUGhrqMT40NFSnTp2SJBUUFMjPz0+tWrUqN6agoMAcExISUm67ISEh5piKpKWladasWTe1b6hZBBdAfeTzMHfs2FEHDhzQuXPn9Oabb2r06NHavn27+bjNZvMYbxhGuWVXu3pMReMrW8+0adOUnJxs3i8uLlZERESl+4ObR3ABNGQ+D7Ofn5958VevXr20d+9evfzyy+brygUFBQoLCzPHFxYWmmfRTqdTpaWlKioq8jhrLiwsVL9+/cwxZ86cKbfds2fPljsb/yG73S673X7zOwgTwQWAyvk8zFczDENut1vt2rWT0+lUdna2evToIUkqLS3V9u3bNXfuXElSz5491bRpU2VnZyshIUGSlJ+fr0OHDmnevHmSpOjoaLlcLu3Zs0e9e/eWJO3evVsul8uMN24OwQUA7/FpmKdPn65hw4YpIiJC58+fV0ZGhrZt26asrCzZbDYlJSVpzpw5at++vdq3b685c+aoRYsWGjVqlCTJ4XBo7NixSklJUXBwsIKCgjR16lR17dpVgwcPliR16tRJQ4cO1bhx47Rs2TJJ379dKi4urspXZDdUBBcAap9Pw3zmzBklJiYqPz9fDodD3bp1U1ZWloYMGSJJevbZZ1VSUqKJEyeqqKhIffr00ZYtW8z3MEvSggUL1KRJEyUkJKikpESDBg1Senq6+R5mSVq3bp2mTJliXr0dHx+vxYsX1+7OWgjBBQDrstz7mK2qLryPmeAC8LXafB/zzf48ro253gjLvcaMqiHCAFA/+fSzsgEAgCfCDACAhRBmAAAshDADAGAhhBkAAAshzAAAWAhhBgDAQggzAAAWQpgBALAQwgwAgIUQZgAALIQwAwBgIYQZAAALIcwAAFgIYQYAwEIIMwAAFkKYAQCwEMIMAICFEGYAACyEMAMAYCGEGQAACyHMAABYCGEGAMBCCDMAABZCmAEAsBDCDACAhRBmAAAshDADAGAhhBkAgCr68ssv9eijjyo4OFgtWrTQz372M+Xm5np1G028ujYAAOqpoqIi3XXXXRo4cKDeffddhYSE6B//+IduueUWr26HMAMAUAVz585VRESEVq1aZS778Y9/7PXt8KdsAACqYPPmzerVq5ceeughhYSEqEePHlqxYoXXt8MZcx2wwWhfbtnDtuM+mAkA1D/FxcUe9+12u+x2e7lxn3/+uZYuXark5GRNnz5de/bs0ZQpU2S32/XrX//aa/MhzD5UUXABALUrIiLC4/7MmTOVmppabtzly5fVq1cvzZkzR5LUo0cPHT58WEuXLvVqmH36p+y0tDTdeeedCggIUEhIiB544AEdO3bMY8yYMWNks9k8bn379vUY43a7NXnyZLVu3Vr+/v6Kj49XXl6ex5iioiIlJibK4XDI4XAoMTFR586dq5H92mC0r9INAOB7p0+flsvlMm/Tpk2rcFxYWJg6d+7ssaxTp0764osvvDofn54xb9++XZMmTdKdd96p7777TjNmzFBMTIyOHDkif39/c9zQoUM9Xmz38/PzWE9SUpIyMzOVkZGh4OBgpaSkKC4uTrm5uWrcuLEkadSoUcrLy1NWVpYkafz48UpMTFRmZma15vxn46dqYTS90V0GAFhMYGCgAgMDKx131113lTt5/Pvf/67IyEivzsenYb4SyStWrVqlkJAQ5ebm6p577jGX2+12OZ3OCtfhcrm0cuVKrVmzRoMHD5YkrV27VhEREdq6datiY2N19OhRZWVladeuXerTp48kacWKFYqOjtaxY8fUsWPHGtpDAKh/rvcXv2+NS5IO1d5katF//ud/ql+/fpozZ44SEhK0Z88eLV++XMuXL/fqdiz1GrPL5ZIkBQUFeSzftm2bQkJCdMstt6h///568cUXFRISIknKzc3VpUuXFBMTY44PDw9XVFSUcnJyFBsbq507d8rhcJhRlqS+ffvK4XAoJyenwjC73W653W7z/tUXBwBAfVPXXmK72b9gVveXiDvvvFObNm3StGnT9MILL6hdu3ZauHChHnnkkRueQ0UsE2bDMJScnKy7775bUVFR5vJhw4bpoYceUmRkpE6cOKHnn39e9957r3Jzc2W321VQUCA/Pz+1atXKY32hoaEqKCiQJBUUFJgh/6GQkBBzzNXS0tI0a9YsL+4hAPhGXQuulcXFxSkuLq5Gt2GZMD/99NP65JNP9OGHH3osHzlypPnvqKgo9erVS5GRkXr77bc1YsSIa67PMAzZbDbz/g//fa0xPzRt2jQlJyeb94uLi8tduQcAvkRw6ydLhHny5MnavHmzduzYoTZt2lx3bFhYmCIjI3X8+Pfv43U6nSotLVVRUZHHWXNhYaH69etnjjlz5ky5dZ09e1ahoaEVbuda72MDgJpGcBs2n4bZMAxNnjxZmzZt0rZt29SuXbtKn/P111/r9OnTCgsLkyT17NlTTZs2VXZ2thISEiRJ+fn5OnTokObNmydJio6Olsvl0p49e9S7d29J0u7du+Vyucx4A0BNI7ioCp+GedKkSVq/fr3+8pe/KCAgwHy91+FwqHnz5vrmm2+UmpqqBx98UGFhYTp58qSmT5+u1q1b65e//KU5duzYsUpJSVFwcLCCgoI0depUde3a1bxKu1OnTho6dKjGjRunZcuWSfr+7VJxcXFckQ3gphFceJNPw7x06VJJ0oABAzyWr1q1SmPGjFHjxo118OBBvf766zp37pzCwsI0cOBAvfHGGwoICDDHL1iwQE2aNFFCQoJKSko0aNAgpaenm+9hlqR169ZpypQp5tXb8fHxWrx4cc3vJIA6jeiitvn8T9nX07x5c7333nuVrqdZs2ZatGiRFi1adM0xQUFBWrt2bbXnCKB+IriwKktc/AUA3kJwUdcRZgB1AsFFQ0GYAfgUwQU8EWYANYLgAjeGMAOoNqIL1BzCDMBEcAHfI8xAA0BwgbqDMAN1GMEF6h/CDFgQwQUaLsIM1CKCC6AyhBnwAoILwFsIM1AJogugNhFmNFgEF4AVEWbUOwQXQF1GmFFnEFwADQFhhs8RXAD4N8KMGkNwAaD6CDNuCNEFgJpBmOGB4AKAbxHmBoLgAkDdQJjrOIILAPULYa4jBm5+xuP+hviXfTQTAEBNIsw+dnVwAQANG2GuIQQXAHAjCPMNILoAgJpCmKvpnneekpr5ehYA4FvXOkE5f/GCpL/U7mTqGcIMADDxF0HfI8wA0AAQ3LqDMANAHUZw659qh/ny5ctKT0/Xxo0bdfLkSdlsNrVr107/8R//ocTERNlstpqYJwA0KAS34apWmA3DUHx8vN555x11795dXbt2lWEYOnr0qMaMGaONGzfqrbfeqqGpAkD9QHRxPdUKc3p6unbs2KH//d//1cCBAz0e++tf/6oHHnhAr7/+un796197dZIAUBcQXHhDtcK8YcMGTZ8+vVyUJenee+/Vc889p3Xr1hFmAPUKwUVtqlaYP/nkE82bN++ajw8bNkx//OMfb3pSAFAbCC6sqFph/te//qXQ0NBrPh4aGqqioqKbnhQA3Chii7quWmEuKytTkybXfkrjxo313XffVXl9aWlp2rhxoz799FM1b95c/fr109y5c9WxY0dzjGEYmjVrlpYvX66ioiL16dNHr7zyirp06WKOcbvdmjp1qjZs2KCSkhINGjRIS5YsUZs2bcwxRUVFmjJlijZv3ixJio+P16JFi3TLLbdU4wgA8CWii4ag2ldljxkzRna7vcLH3W53tTa+fft2TZo0SXfeeae+++47zZgxQzExMTpy5Ij8/f0lSfPmzdP8+fOVnp6uDh06aPbs2RoyZIiOHTumgIAASVJSUpIyMzOVkZGh4OBgpaSkKC4uTrm5uWrcuLEkadSoUcrLy1NWVpYkafz48UpMTFRmZma15gzA+wgu8G/VCvPo0aMrHVOdC7+uRPKKVatWKSQkRLm5ubrnnntkGIYWLlyoGTNmaMSIEZKk1atXKzQ0VOvXr9eECRPkcrm0cuVKrVmzRoMHD5YkrV27VhEREdq6datiY2N19OhRZWVladeuXerTp48kacWKFYqOjtaxY8c8ztABeA/BBaqvWmFetWpVTc1DkuRyuSRJQUFBkqQTJ06ooKBAMTEx5hi73a7+/fsrJydHEyZMUG5uri5duuQxJjw8XFFRUcrJyVFsbKx27twph8NhRlmS+vbtK4fDoZycnArD7Ha7Pf4CUFxc7PX9BeoqggsruOedpxTQzP+Gn2/VL9ywzEdyGoah5ORk3X333YqKipIkFRQUSFK5C85CQ0N16tQpc4yfn59atWpVbsyV5xcUFCgkJKTcNkNCQswxV0tLS9OsWbNubqeAOobgAr5XrTBf+XPydVfYpImcTqeGDBmi4cOHV3ndTz/9tD755BN9+OGH5R67+mM+DcOo9KM/rx5T0fjrrWfatGlKTk427xcXFysiIuK62wSsiuACdUe1wuxwOCodc/nyZR0/flx/+tOfNHXqVL3wwguVPmfy5MnavHmzduzY4XEltdPplPT9GW9YWJi5vLCw0DyLdjqdKi0tVVFRkcdZc2Fhofr162eOOXPmTLntnj179ppv/7Lb7de8yA2wCoIL1D819hrz22+/raeeeuq6YTYMQ5MnT9amTZu0bds2tWvXzuPxdu3ayel0Kjs7Wz169JAklZaWavv27Zo7d64kqWfPnmratKmys7OVkJAgScrPz9ehQ4fMD0OJjo6Wy+XSnj171Lt3b0nS7t275XK5zHgDVkJwgYarWmF+/PHHKx1js9m0cuVK3XXXXerVq9d1x06aNEnr16/XX/7yFwUEBJiv9zocDjVv3lw2m01JSUmaM2eO2rdvr/bt22vOnDlq0aKFRo0aZY4dO3asUlJSFBwcrKCgIE2dOlVdu3Y1r9Lu1KmThg4dqnHjxmnZsmWSvn+7VFxcHFdko9YRXQDXU+0vsYiMjFSPHj1kGMZ1x95yyy3auHHjdccsXbpUkjRgwACP5atWrdKYMWMkSc8++6xKSko0ceJE8wNGtmzZYr6HWZIWLFigJk2aKCEhwfyAkfT0dPM9zJK0bt06TZkyxbx6Oz4+XosXL67qrgOVIrgAvMFmVFbYH5g4caIyMjLUtm1bPf7443r00UfNtzbVd8XFxXI4HPr7C1k3dXm+t7wf/7Kvp9BgEFyg6s5fvKAOvxsql8ulwMDAGtmGt34e18Zcb0S1zpiXLFmiBQsWaOPGjXrttdc0bdo03XfffRo7dqxiYmIqvVIasBKCC8CKqv0+ZrvdrocfflgPP/ywTp06pfT0dE2cOFGXLl3SkSNH1LJly5qYJ1BlBBdAXXZTHzBis9lks9lkGIYuX77srTkBFSK4ABqCaofZ7Xabf8r+8MMPFRcXp8WLF2vo0KFq1KhRTcwR9RzBBYB/q1aYf3jx12OPPWZ+mxNwLUQXAKqnWmF+9dVX1bZtW7Vr107bt2/X9u3bKxxX2dukUPcRXAANXVpamqZPn65nnnlGCxcu9Np6qxXmX//611x5Xc8RXACo3N69e7V8+XJ169bN6+uu9geMoG4iuADgHd98840eeeQRrVixQrNnz/b6+i3ztY+4MQQXAG5OcXGxx/3KvsRo0qRJuu+++zR48GDC3JCtD2nmcX8UQQYAr7j6K31nzpyp1NTUCsdmZGTo448/1t69e2tsPoTZAq6OLgCg9pw+fdrjIzmvdbZ8+vRpPfPMM9qyZYuaNau5n9uEuQYRXACwvsDAwCp9VnZubq4KCwvVs2dPc1lZWZl27NihxYsXy+12e3x50o0izDeA4AJAwzNo0CAdPHjQY9ljjz2m22+/Xb/97W+9EmWJMFfbn2+1iywDaOiudYJyseS7Wp5J7QkICFBUVJTHMn9/fwUHB5dbfjMIMwDAxF8EfY8wA0ADQHBrxrZt27y+TsIMAHUYwa1/CDMAWBTRbZgIMwDUMoKL6yHMAOAlBBfeQJgBoBIEF7WJMANosAgurIgwA6h3CC7qMsIMoE4huqjvCDMASyC4wPcIM4AaRXCB6iHMAG4IwQVqBmEG4IHgAr5FmIEGguACdQNhBuo4ggvUL4QZsDCiCzQ8hBnwAYIL4FoIM+BFBBfAzSLMQBUQXAC1hTCjQSO4AKymkS83vmPHDg0fPlzh4eGy2Wx66623PB4fM2aMbDabx61v374eY9xutyZPnqzWrVvL399f8fHxysvL8xhTVFSkxMREORwOORwOJSYm6ty5czW8d/Cl9SHNqnQDAKvx6RnzhQsX1L17dz322GN68MEHKxwzdOhQrVq1yrzv5+fn8XhSUpIyMzOVkZGh4OBgpaSkKC4uTrm5uWrcuLEkadSoUcrLy1NWVpYkafz48UpMTFRmZmYN7RlqEkEFUJ/5NMzDhg3TsGHDrjvGbrfL6XRW+JjL5dLKlSu1Zs0aDR48WJK0du1aRUREaOvWrYqNjdXRo0eVlZWlXbt2qU+fPpKkFStWKDo6WseOHVPHjh29u1O4YQQXAOrAa8zbtm1TSEiIbrnlFvXv318vvviiQkJCJEm5ubm6dOmSYmJizPHh4eGKiopSTk6OYmNjtXPnTjkcDjPKktS3b185HA7l5ORcM8xut1tut9u8X1xcXEN7WP8RXAA14c+32tWs+Y3/fLlY8p0XZ+M9lg7zsGHD9NBDDykyMlInTpzQ888/r3vvvVe5ubmy2+0qKCiQn5+fWrVq5fG80NBQFRQUSJIKCgrMkP9QSEiIOaYiaWlpmjVrlnd3qB4htgBQMywd5pEjR5r/joqKUq9evRQZGam3335bI0aMuObzDMOQzWYz7//w39cac7Vp06YpOTnZvF9cXKyIiIjq7kKdQ3ABwLcsHearhYWFKTIyUsePH5ckOZ1OlZaWqqioyOOsubCwUP369TPHnDlzpty6zp49q9DQ0Gtuy263y263e3kPfIfgAkDdUKfC/PXXX+v06dMKCwuTJPXs2VNNmzZVdna2EhISJEn5+fk6dOiQ5s2bJ0mKjo6Wy+XSnj171Lt3b0nS7t275XK5zHjXdUQXAOoPn4b5m2++0WeffWbeP3HihA4cOKCgoCAFBQUpNTVVDz74oMLCwnTy5ElNnz5drVu31i9/+UtJksPh0NixY5WSkqLg4GAFBQVp6tSp6tq1q3mVdqdOnTR06FCNGzdOy5Ytk/T926Xi4uLq9BXZxBgA6iefhnnfvn0aOHCgef/Ka7qjR4/W0qVLdfDgQb3++us6d+6cwsLCNHDgQL3xxhsKCAgwn7NgwQI1adJECQkJKikp0aBBg5Senm6+h1mS1q1bpylTpphXb8fHx2vx4sW1tJcAAFSdT8M8YMAAGYZxzcffe++9StfRrFkzLVq0SIsWLbrmmKCgIK1du/aG5ggAQG3y6UdyAgAAT4QZAAALIcwAAFgIYQYAwEIIMwAAFkKYAQCwEMIMAICFEGYAACyEMAMAYCGEGQAACyHMAABYCGEGAMBCCDMAABZCmAEAsBDCDACAhRBmAAAshDADAGAhhBkAAAshzAAAWAhhBgDAQggzAAAWQpgBALAQwgwAgIUQZgAALIQwAwBgIYQZAAALIcwAAFRBWlqa7rzzTgUEBCgkJEQPPPCAjh075vXtEGYAAKpg+/btmjRpknbt2qXs7Gx99913iomJ0YULF7y6nSZeXRsAAPVUVlaWx/1Vq1YpJCREubm5uueee7y2HcIMAGjQiouLPe7b7XbZ7fZKn+dyuSRJQUFBXp0PYbawZ75ees3HXg5+qhZnAgD1V0REhMf9mTNnKjU19brPMQxDycnJuvvuuxUVFeXV+RBmH7hecAEAtev06dMKDAw071flbPnpp5/WJ598og8//NDr8yHMXkJsAaBuCgwM9AhzZSZPnqzNmzdrx44datOmjdfnQ5ir6al/rVRgMz9fTwMALOeZr5eq+GKpZvh6IjXEMAxNnjxZmzZt0rZt29SuXbsa2Q5hBgBUir8KSpMmTdL69ev1l7/8RQEBASooKJAkORwONW/e3Gvb8WmYd+zYod///vfKzc1Vfn6+Nm3apAceeMB83DAMzZo1S8uXL1dRUZH69OmjV155RV26dDHHuN1uTZ06VRs2bFBJSYkGDRqkJUuWePx5oaioSFOmTNHmzZslSfHx8Vq0aJFuueWW2tpVALAkglt1S5d+f6wGDBjgsXzVqlUaM2aM17bj0zBfuHBB3bt312OPPaYHH3yw3OPz5s3T/PnzlZ6erg4dOmj27NkaMmSIjh07poCAAElSUlKSMjMzlZGRoeDgYKWkpCguLk65ublq3LixJGnUqFHKy8sz34M2fvx4JSYmKjMzs/Z2FgBqEcH1PsMwamU7Pg3zsGHDNGzYsAofMwxDCxcu1IwZMzRixAhJ0urVqxUaGqr169drwoQJcrlcWrlypdasWaPBgwdLktauXauIiAht3bpVsbGxOnr0qLKysrRr1y716dNHkrRixQpFR0fr2LFj6tixY+3sLAB4AcGt/yz7GvOJEydUUFCgmJgYc5ndblf//v2Vk5OjCRMmKDc3V5cuXfIYEx4erqioKOXk5Cg2NlY7d+6Uw+EwoyxJffv2lcPhUE5OzjXD7Ha75Xa7zftXvwEdALyJ4FbfzV6Ma9UL1Swb5isvqoeGhnosDw0N1alTp8wxfn5+atWqVbkxV55fUFCgkJCQcusPCQkxx1QkLS1Ns2bNuql9AACCi+qybJivsNlsHvcNwyi37GpXj6lofGXrmTZtmpKTk837xcXF5T4dBkDDRXBRUywbZqfTKen7M96wsDBzeWFhoXkW7XQ6VVpaqqKiIo+z5sLCQvXr188cc+bMmXLrP3v2bLmz8R+q6melAqh/iC58ybJhbteunZxOp7Kzs9WjRw9JUmlpqbZv3665c+dKknr27KmmTZsqOztbCQkJkqT8/HwdOnRI8+bNkyRFR0fL5XJpz5496t27tyRp9+7dcrlcZrwBNAwEF3WBT8P8zTff6LPPPjPvnzhxQgcOHFBQUJDatm2rpKQkzZkzR+3bt1f79u01Z84ctWjRQqNGjZL0/Zu6x44dq5SUFAUHBysoKEhTp05V165dzau0O3XqpKFDh2rcuHFatmyZpO/fLhUXF8cV2UA9QXBRn/g0zPv27dPAgQPN+1de0x09erTS09P17LPPqqSkRBMnTjQ/YGTLli3me5glacGCBWrSpIkSEhLMDxhJT08338MsSevWrdOUKVPMq7fj4+O1ePHiWtpLADeK4KIhshm19Y7pOq64uFgOh0P/evFRS3xWNl/7iLqM4NZfxRdLFTRjrVwuV7W+GKJa2/DSz+PamOuNsOxrzADqHoIL3DzCDKBKiC5QOwgz0MARXMBaCDNQTxFcoG4izEAdQ3CB+o0wAxZBcAFIhBmocQQXQHUQZuAmEF0A3kaYgQoQXAC+QpjRoBBcAFZHmFEvEFwA9QVhhqURXAANDWGGTxBcAKgYYYZXEVwAuDmEGVVGdAGg5hFmEFwAsBDCXI8RXACoewhzHUV0AaB+IswWNvHOQdd8bMne/63FmQAAagth9oHrBRcA0LARZi8htgAAbyDM1ZR8R3/5+bfw9TQAwJKS7+gvaa2vp1GnEWYAQKWq/FfBC9/W7EQaAMIMAA0YL8NZD2EGgHqI4NZdhBkA6hCCW/8RZgCwAIKLKwgzANQwoovqIMwAcIMILmoCYQaAqxBc+BJhBtBgEFzUBYQZQJ1HcFGfEGYAlkVw0RARZgA+QXSBihFmAF5FcIGbY+kwp6amatasWR7LQkNDVVBQIEkyDEOzZs3S8uXLVVRUpD59+uiVV15Rly5dzPFut1tTp07Vhg0bVFJSokGDBmnJkiVq06ZNre4LUNcRXKB2WDrMktSlSxdt3brVvN+4cWPz3/PmzdP8+fOVnp6uDh06aPbs2RoyZIiOHTumgIAASVJSUpIyMzOVkZGh4OBgpaSkKC4uTrm5uR7rAhoqggtYi+XD3KRJEzmdznLLDcPQwoULNWPGDI0YMUKStHr1aoWGhmr9+vWaMGGCXC6XVq5cqTVr1mjw4MGSpLVr1yoiIkJbt25VbGxsre4LUJsILlA3WT7Mx48fV3h4uOx2u/r06aM5c+boJz/5iU6cOKGCggLFxMSYY+12u/r376+cnBxNmDBBubm5unTpkseY8PBwRUVFKScn57phdrvdcrvd5v3i4uKa2UGgmgguUL9ZOsx9+vTR66+/rg4dOujMmTOaPXu2+vXrp8OHD5uvM4eGhno8JzQ0VKdOnZIkFRQUyM/PT61atSo35srzryUtLa3c69tATSK4QPUk39Fffv4tbvj5pRe+lbTWexPyEkuHediwYea/u3btqujoaP30pz/V6tWr1bdvX0mSzWbzeI5hGOWWXa0qY6ZNm6bk5GTzfnFxsSIiIqq7C4Akogug6iwd5qv5+/ura9euOn78uB544AFJ358Vh4WFmWMKCwvNs2in06nS0lIVFRV5nDUXFhaqX79+192W3W6X3W73/k6gXiG4ALytToXZ7Xbr6NGj+vnPf6527drJ6XQqOztbPXr0kCSVlpZq+/btmjt3riSpZ8+eatq0qbKzs5WQkCBJys/P16FDhzRv3jyf7Qesj+AC8BVLh3nq1KkaPny42rZtq8LCQs2ePVvFxcUaPXq0bDabkpKSNGfOHLVv317t27fXnDlz1KJFC40aNUqS5HA4NHbsWKWkpCg4OFhBQUGaOnWqunbtal6ljYaF4AKwOkuHOS8vTw8//LC++uor3Xrrrerbt6927dqlyMhISdKzzz6rkpISTZw40fyAkS1btpjvYZakBQsWqEmTJkpISDA/YCQ9PZ33MNczBBdAfWHpMGdkZFz3cZvNptTUVKWmpl5zTLNmzbRo0SItWrTIy7NDbSC4ABoaS4cZ9RvRBYDyCDO8juACqM+WLFmi3//+98rPz1eXLl20cOFC/fznP/fa+gkzqozgAmjo3njjDSUlJWnJkiW66667tGzZMg0bNkxHjhxR27ZtvbINwgyCCwBVNH/+fI0dO1ZPPPGEJGnhwoV67733tHTpUqWlpXllG4S5HiO4AFC5q78L4VofMFVaWqrc3Fw999xzHstjYmKUk5PjtfkQ5jqK6AKAd1z9ccszZ86s8N0+X331lcrKyir8jobKvn+hOggzAKBBO336tAIDA837lX0c8418R0N1EGYAQIMWGBjoEeZrad26tRo3blzu7PiH39HgDY28tiYAAOoxPz8/9ezZU9nZ2R7Ls7OzK/1ipOrgjBkAgCpKTk5WYmKievXqpejoaC1fvlxffPGFnnzySa9tgzADAFBFI0eO1Ndff60XXnhB+fn5ioqK0jvvvGN+h4M3EGYAAKph4sSJmjhxYo2tn9eYAQCwEMIMAICFEGYAACyEMAMAYCGEGQAACyHMAABYCGEGAMBCCDMAABZCmAEAsBDCDACAhRBmAAAshDADAGAhhBkAAAshzAAAWAhhBgDAQggzAAAWQpgBALAQwgwAgIUQZgAALIQwAwBgIYQZAAALIcwAAFgIYQYAwEIaVJiXLFmidu3aqVmzZurZs6c++OADX08JAAAPDSbMb7zxhpKSkjRjxgzt379fP//5zzVs2DB98cUXvp4aAACmBhPm+fPna+zYsXriiSfUqVMnLVy4UBEREVq6dKmvpwYAgKmJrydQG0pLS5Wbm6vnnnvOY3lMTIxycnIqfI7b7Zbb7Tbvu1yu79d1oaTmJgoAddyVn5GGYdTatnz1/JrSIML81VdfqaysTKGhoR7LQ0NDVVBQUOFz0tLSNGvWrHLL1//HlBqZIwDUJ+fPn5fD4aiRdfv5+cnpdHrl57HT6ZSfn58XZuU9DSLMV9hsNo/7hmGUW3bFtGnTlJycbN4/d+6cIiMj9cUXX9TY/9jwb8XFxYqIiNDp06cVGBjo6+nUexzv2lWfj7dhGDp//rzCw8NrbBvNmjXTiRMnVFpaetPr8vPzU7NmzbwwK+9pEGFu3bq1GjduXO7suLCwsNxZ9BV2u112u73ccofDUe/+j2RlgYGBHO9axPGuXfX1eNfGyUuzZs0sF1RvaRAXf/n5+alnz57Kzs72WJ6dna1+/fr5aFYAAJTXIM6YJSk5OVmJiYnq1auXoqOjtXz5cn3xxRd68sknfT01AABMDSbMI0eO1Ndff60XXnhB+fn5ioqK0jvvvKPIyMgqPd9ut2vmzJkV/nkb3sfxrl0c79rF8cb12IzauKYdAABUSYN4jRkAgLqCMAMAYCGEGQAACyHMAABYCGGuIr4y8ubt2LFDw4cPV3h4uGw2m9566y2Pxw3DUGpqqsLDw9W8eXMNGDBAhw8f9hjjdrs1efJktW7dWv7+/oqPj1deXl4t7kXdkZaWpjvvvFMBAQEKCQnRAw88oGPHjnmM4Zh7z9KlS9WtWzfzQ0Oio6P17rvvmo9zrFFVhLkK+MpI77hw4YK6d++uxYsXV/j4vHnzNH/+fC1evFh79+6V0+nUkCFDdP78eXNMUlKSNm3apIyMDH344Yf65ptvFBcXp7KystrajTpj+/btmjRpknbt2qXs7Gx99913iomJ0YULF8wxHHPvadOmjV566SXt27dP+/bt07333qv777/fjC/HGlVmoFK9e/c2nnzySY9lt99+u/Hcc8/5aEZ1nyRj06ZN5v3Lly8bTqfTeOmll8xlFy9eNBwOh/Hqq68ahmEY586dM5o2bWpkZGSYY7788kujUaNGRlZWVq3Nva4qLCw0JBnbt283DINjXhtatWpl/OlPf+JYo1o4Y67Ela+MjImJ8Vh+va+MRPWdOHFCBQUFHsfZbrerf//+5nHOzc3VpUuXPMaEh4crKiqK/xZVcOWrS4OCgiRxzGtSWVmZMjIydOHCBUVHR3OsUS2EuRI38pWRqL4rx/J6x7mgoEB+fn5q1arVNcegYoZhKDk5WXfffbeioqIkccxrwsGDB9WyZUvZ7XY9+eST2rRpkzp37syxRrU0mI/kvFnV+cpI3LgbOc78t6jc008/rU8++UQffvhhucc45t7TsWNHHThwQOfOndObb76p0aNHa/v27ebjHGtUBWfMlbiRr4xE9TmdTkm67nF2Op0qLS1VUVHRNcegvMmTJ2vz5s16//331aZNG3M5x9z7/Pz8dNttt6lXr15KS0tT9+7d9fLLL3OsUS2EuRJ8ZWTtaNeunZxOp8dxLi0t1fbt283j3LNnTzVt2tRjTH5+vg4dOsR/iwoYhqGnn35aGzdu1F//+le1a9fO43GOec0zDENut5tjjerx1VVndUlGRobRtGlTY+XKlcaRI0eMpKQkw9/f3zh58qSvp1annD9/3ti/f7+xf/9+Q5Ixf/58Y//+/capU6cMwzCMl156yXA4HMbGjRuNgwcPGg8//LARFhZmFBcXm+t48sknjTZt2hhbt241Pv74Y+Pee+81unfvbnz33Xe+2i3LeuqppwyHw2Fs27bNyM/PN2/ffvutOYZj7j3Tpk0zduzYYZw4ccL45JNPjOnTpxuNGjUytmzZYhgGxxpVR5ir6JVXXjEiIyMNPz8/44477jDfcoKqe//99w1J5W6jR482DOP7t+/MnDnTcDqdht1uN+655x7j4MGDHusoKSkxnn76aSMoKMho3ry5ERcXZ3zxxRc+2Bvrq+hYSzJWrVpljuGYe8/jjz9u/oy49dZbjUGDBplRNgyONaqOr30EAMBCeI0ZAAALIcwAAFgIYQYAwEIIMwAAFkKYAQCwEMIMAICFEGYAACyEMAN1RGlpqW677TZ99NFH1xxz8uRJ2Ww2HThwwKvbvvPOO7Vx40avrhNAxQgz4CVjxoyRzWaTzWZTkyZN1LZtWz311FPlvpRAkkpKStSqVSsFBQWppKSkSutfvny5IiMjddddd3l76pV6/vnn9dxzz+ny5cu1vm2goSHMgBcNHTpU+fn5OnnypP70pz8pMzNTEydOLDfuzTffVFRUlDp37lzlM9FFixbpiSee8PaUq+S+++6Ty+XSe++955PtAw0JYQa8yG63y+l0qk2bNoqJidHIkSO1ZcuWcuNWrlypRx99VI8++qhWrlxZ6Xo//vhjffbZZ7rvvvs8lu/Zs0c9evRQs2bN1KtXL+3fv7/cc48cOaJf/OIXatmypUJDQ5WYmKivvvrKfPz8+fN65JFH5O/vr7CwMC1YsEADBgxQUlKSOaZx48b6xS9+oQ0bNlTjaAC4EYQZqCGff/65srKy1LRpU4/l//jHP7Rz504lJCQoISFBOTk5+vzzz6+7rh07dqhDhw4KDAw0l124cEFxcXHq2LGjcnNzlZqaqqlTp3o8Lz8/X/3799fPfvYz7du3T1lZWTpz5owSEhLMMcnJyfroo4+0efNmZWdn64MPPtDHH39cbg69e/fWBx98cCOHAkA1NPH1BID65H/+53/UsmVLlZWV6eLFi5Kk+fPne4x57bXXNGzYMLVq1UrS93/+fu211zR79uxrrvfkyZMKDw/3WLZu3TqVlZXptddeU4sWLdSlSxfl5eXpqaeeMscsXbpUd9xxh+bMmeOx/YiICP39739XWFiYVq9erfXr12vQoEGSpFWrVpXbliT96Ec/0hdffKHLly+rUSN+pwdqCv/vArxo4MCBOnDggHbv3q3JkycrNjZWkydPNh8vKyvT6tWr9eijj5rLHn30Ua1evVplZWXXXG9JSYmaNWvmsezo0aPq3r27WrRoYS6Ljo72GJObm6v3339fLVu2NG+33367pO/P3D///HNdunRJvXv3Np/jcDjUsWPHcnNo3ry5Ll++LLfbXcWjAeBGEGbAi/z9/XXbbbepW7du+uMf/yi3261Zs2aZj7/33nv68ssvNXLkSDVp0kRNmjTRr371K+Xl5VX4WvQVrVu3Lnd1d1W+sfXy5csaPny4Dhw44HE7fvy47rnnHnMdNput0nX/61//UosWLdS8efNKtwvgxhFmoAbNnDlTf/jDH/TPf/5T0vcXff3qV78qF8pHHnnkuheB9ejRQ59++qlHMDt37qy//e1vHm+32rVrl8fz7rjjDh0+fFg//vGPddttt3nc/P399dOf/lRNmzbVnj17zOcUFxfr+PHj5eZw6NAh3XHHHTd8LABUDWEGatCAAQPUpUsXzZkzR2fPnlVmZqZGjx6tqKgoj9vo0aO1efNmnT17tsL1DBw4UBcuXNDhw4fNZaNGjVKjRo00duxYHTlyRO+8847+8Ic/eDxv0qRJ+te//qWHH35Ye/bs0eeff64tW7bo8ccfV1lZmQICAjR69Gj95je/0fvvv6/Dhw/r8ccfV6NGjcqdRX/wwQeKiYnx/kEC4IEwAzUsOTlZK1as0JIlS+Tv729eZPVDAwcOVEBAgNasWVPhOoKDgzVixAitW7fOXNayZUtlZmbqyJEj6tGjh2bMmKG5c+d6PC88PFwfffSRysrKFBsbq6ioKD3zzDNyOBzmBVzz589XdHS04uLiNHjwYN11113q1KmTx2vaX375pXJycvTYY49545AAuA6bUZUXqgD43MGDBzV48GB99tlnCggIqLHtXLhwQT/60Y/03//93xo7dqwk6Te/+Y1cLpeWL19eY9sF8D3OmIE6omvXrpo3b55Onjzp1fXu379fGzZs0D/+8Q99/PHHeuSRRyRJ999/vzkmJCRE//Vf/+XV7QKoGGfMQAO3f/9+PfHEEzp27Jj8/PzUs2dPzZ8/X127dvX11IAGiTADAGAh/CkbAAALIcwAAFgIYQYAwEIIMwAAFkKYAQCwEMIMAICFEGYAACyEMAMAYCGEGQAAC/n/jmGTokZ7YFoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from rubin_scheduler.utils import survey_start_mjd, calc_season\n",
    "\n",
    "from rubin_scheduler.data import get_data_dir\n",
    "from rubin_scheduler.scheduler.utils import scheduled_observation\n",
    "from rubin_scheduler.utils import  ddf_locations\n",
    "from rubin_scheduler.scheduler.surveys import generate_ddf_scheduled_obs\n",
    "\n",
    "ra = np.arange(0, 360, 1)\n",
    "mjds = np.arange(survey_start_mjd(), survey_start_mjd()+11*365.25, 1)\n",
    "seasons = calc_season(ra, mjds)\n",
    "seasons = np.swapaxes(seasons, 0, 1)\n",
    "seasons = np.floor(seasons)\n",
    "\n",
    "plt.imshow(seasons, aspect=.1, origin='lower', cmap='Set2')\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"RA (deg)\")\n",
    "plt.ylabel(\"MJD\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "7301c77e-ee31-44b3-9769-45a966a6b579",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ddf_slopes_new(ddf_name, raw_obs, night_season, min_season=60 / 365.25):\n",
    "    \"\"\"\n",
    "    Let's make custom slopes for each DDF\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ddf_name : `str`\n",
    "       The DDF name to use\n",
    "    raw_obs : `np.array`, (N,)\n",
    "        An array with values of 1 or zero. One element per night, value of\n",
    "        1 indicates the night is during an active observing season.\n",
    "    night_season : `np.array`, (N,)\n",
    "        An array of floats with the fractional season value.\n",
    "        Season values range from 0-1 for the first season;\n",
    "        0.5 would be the \"peak\" of seasonal visiblity.\n",
    "        These should ONLY be the night the field should be 'active'\n",
    "        (take out the season ends first).\n",
    "    min_season : `float`, optional\n",
    "        Minimal season length that is useful.\n",
    "        In units of \"season\" (fractions of a year).\n",
    "    \"\"\"\n",
    "\n",
    "    # OK, so 258 sequences is ~1% of the survey\n",
    "    # so a 25.8 sequences is a 0.1% season\n",
    "    # COSMOS is going to be 0.7% for 3 years, then 0.175 for the rest.\n",
    "\n",
    "    ss = 30  # standard season, was 45 .. note this may be affected by season length\n",
    "\n",
    "    int_season = np.floor(night_season)\n",
    "\n",
    "    n_season = night_season[np.where(raw_obs)]\n",
    "    r_season = int_season[np.where(raw_obs)]\n",
    "    season_list = np.unique(r_season)\n",
    "\n",
    "    # Calculate season length for each season\n",
    "    # In general this should be the same each season except for first and last\n",
    "    season_length = np.zeros(len(season_list), float)\n",
    "    for i, season in enumerate(season_list):\n",
    "        match = np.where(r_season == season)\n",
    "        season_length[i] = n_season[match].max() - n_season[match].min()\n",
    "\n",
    "    # Determine goal number of sequences in each season.\n",
    "    season_vals = np.ones(len(season_list), float) * ss\n",
    "    # Throw out seasons which are too short\n",
    "    too_short = np.where(season_length < min_season)\n",
    "    season_vals[too_short] = 0\n",
    "    # Adjust other seasons, relative to the max season length.\n",
    "    season_vals = season_vals * season_length / np.max(season_length)\n",
    "\n",
    "    # Add extra adjustment for COSMOS to boost visits in seasons 0-3 \n",
    "    # Note this means boost in first incomplete season too\n",
    "    if ddf_name == \"COSMOS\":\n",
    "        season_vals[0:3] = season_vals[0:3] * 5\n",
    "        season_vals[3] = season_vals[3] * 2\n",
    "\n",
    "    if ddf_name == \"EDFS_b\":\n",
    "        # EDFS_b ddf visits are allocated some other way\n",
    "        season_vals = season_vals * 0\n",
    "\n",
    "    # Round the season_vals so that we're looking for integer sequences\n",
    "    season_vals = np.round(season_vals)\n",
    "\n",
    "    #print(list(zip(season_vals, season_length)))\n",
    "          \n",
    "    # assign cumulative values\n",
    "    cumulative_desired = np.zeros(raw_obs.size, dtype=float)\n",
    "    for i, season in enumerate(season_list):\n",
    "        in_season = np.where(int_season == season)\n",
    "        cumulative = np.cumsum(raw_obs[in_season])\n",
    "        if cumulative.max() > 0:\n",
    "            cumulative = cumulative / cumulative.max() * season_vals[i]\n",
    "            cumulative_desired[in_season] = cumulative + np.max(cumulative_desired)\n",
    "\n",
    "    return cumulative_desired"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "b1d7edec-3987-4079-b290-79672861feda",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ddf_slopes_old(ddf_name, raw_obs, night_season):\n",
    "    \"\"\"\n",
    "    Let's make custom slopes for each DDF\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ddf_name : str\n",
    "       The DDF name to use\n",
    "    raw_obs : np.array\n",
    "        An array with values of 1 or zero. One element per night, value of\n",
    "        1 indicates the night is during an active observing season.\n",
    "    night_season : np.array\n",
    "        An array of floats with the fractional season value\n",
    "        (e.g., 0.5 would be half way through the first season)\n",
    "    \"\"\"\n",
    "\n",
    "    # OK, so 258 sequences is ~1% of the survey\n",
    "    # so a 25.8 sequences is a 0.1% season\n",
    "    # COSMOS is going to be 0.7% for 3 years, then 0.175 for the rest.\n",
    "\n",
    "    ss = 30  # standard season, was 45\n",
    "\n",
    "    if (ddf_name == \"ELAISS1\") | (ddf_name == \"XMM_LSS\") | (ddf_name == \"ECDFS\"):\n",
    "        # Dict with keys for each season and values of the number of\n",
    "        # sequences to attempt.\n",
    "        season_vals = {\n",
    "            0: 10,\n",
    "            1: ss,\n",
    "            2: ss,\n",
    "            3: ss,\n",
    "            4: ss,\n",
    "            5: ss,\n",
    "            6: ss,\n",
    "            7: ss,\n",
    "            8: ss,\n",
    "            9: ss,\n",
    "            10: 10,\n",
    "        }\n",
    "\n",
    "        round_season = np.floor(night_season)\n",
    "\n",
    "        cumulative_desired = np.zeros(raw_obs.size, dtype=float)\n",
    "        for season in np.unique(round_season):\n",
    "            in_season = np.where(round_season == season)\n",
    "            cumulative = np.cumsum(raw_obs[in_season])\n",
    "            if cumulative.max() > 0:\n",
    "                cumulative = cumulative / cumulative.max() * season_vals[season]\n",
    "                cumulative_desired[in_season] = cumulative + np.max(cumulative_desired)\n",
    "\n",
    "    if ddf_name == \"EDFS_a\":\n",
    "        season_vals = {\n",
    "            0: 10,\n",
    "            1: ss,\n",
    "            2: ss,\n",
    "            3: ss,\n",
    "            4: ss,\n",
    "            5: ss,\n",
    "            6: ss,\n",
    "            7: ss,\n",
    "            8: ss,\n",
    "            9: ss,\n",
    "            10: 10,\n",
    "        }\n",
    "\n",
    "        round_season = np.floor(night_season)\n",
    "\n",
    "        cumulative_desired = np.zeros(raw_obs.size, dtype=float)\n",
    "        for season in np.unique(round_season):\n",
    "            in_season = np.where(round_season == season)\n",
    "            cumulative = np.cumsum(raw_obs[in_season])\n",
    "            if cumulative.max() > 0:\n",
    "                cumulative = cumulative / cumulative.max() * season_vals[season]\n",
    "                cumulative_desired[in_season] = cumulative + np.max(cumulative_desired)\n",
    "\n",
    "    if ddf_name == \"COSMOS\":\n",
    "        # looks like COSMOS has no in-season time for 10 at the\n",
    "        # current start mjd.\n",
    "        season_vals = {\n",
    "            0: 10,\n",
    "            1: ss * 5,\n",
    "            2: ss * 5,\n",
    "            3: ss * 2,\n",
    "            4: ss,\n",
    "            5: ss,\n",
    "            6: ss,\n",
    "            7: ss,\n",
    "            8: ss,\n",
    "            9: ss,\n",
    "            10: 10,\n",
    "        }\n",
    "\n",
    "        round_season = np.floor(night_season)\n",
    "\n",
    "        cumulative_desired = np.zeros(raw_obs.size, dtype=float)\n",
    "        for season in np.unique(round_season):\n",
    "            in_season = np.where(round_season == season)[0]\n",
    "            cumulative = np.cumsum(raw_obs[in_season])\n",
    "            if cumulative.max() > 0:\n",
    "                cumulative = cumulative / cumulative.max() * season_vals[season]\n",
    "                cumulative_desired[in_season] = cumulative + np.max(cumulative_desired)\n",
    "\n",
    "    return cumulative_desired"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "c7348f62-233e-404a-8e58-5b786c9ea36b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimize_ddf_times(\n",
    "    ddf_name,\n",
    "    ddf_RA,\n",
    "    ddf_grid,\n",
    "    sun_limit=-18,\n",
    "    sequence_time=60.0,\n",
    "    airmass_limit=2.5,\n",
    "    sky_limit=None,\n",
    "    g_depth_limit=23.5,\n",
    "    season_unobs_frac=0.1,\n",
    "    old=False\n",
    "):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    ddf : `str`\n",
    "        The name of the DDF\n",
    "    ddf_grid : `np.array`\n",
    "        An array with info for the DDFs. Generated by the\n",
    "        rubin_scheduler.scheduler/surveys/generate_ddf_grid.py` script\n",
    "    season_unobs_frac : `float`\n",
    "        7.2 month observing season if season_unobs_frac = 0.2\n",
    "        (shaves 20% off each end of the full year)\n",
    "    sequence_time : `float`\n",
    "        How long a sequence is expected to be (minutes). Used to make\n",
    "        sure things are not scheduled too close to twilight.\n",
    "    \"\"\"\n",
    "    sun_limit = np.radians(sun_limit)\n",
    "    sequence_time = sequence_time / 60.0 / 24.0  # to days\n",
    "\n",
    "    # XXX-- double check that I got this right\n",
    "    ack = ddf_grid[\"sun_alt\"][0:-1] * ddf_grid[\"sun_alt\"][1:]\n",
    "    night = np.zeros(ddf_grid.size, dtype=int)\n",
    "    night[np.where((ddf_grid[\"sun_alt\"][1:] >= 0) & (ack < 0))] += 1\n",
    "    night = np.cumsum(night)\n",
    "    ngrid = ddf_grid[\"mjd\"].size\n",
    "\n",
    "    # set a sun, airmass, sky masks\n",
    "    sun_mask = np.ones(ngrid, dtype=int)\n",
    "    sun_mask[np.where(ddf_grid[\"sun_alt\"] >= sun_limit)] = 0\n",
    "\n",
    "    # expand sun mask backwards by the sequence time.\n",
    "    n_back = np.ceil(sequence_time / (ddf_grid[\"mjd\"][1] - ddf_grid[\"mjd\"][0])).astype(int)\n",
    "    shadow_indx = np.where(sun_mask == 0)[0] - n_back\n",
    "    shadow_indx = shadow_indx[np.where(shadow_indx >= 0)]\n",
    "\n",
    "    sun_mask[shadow_indx] = 0\n",
    "\n",
    "    airmass_mask = np.ones(ngrid, dtype=int)\n",
    "    airmass_mask[np.where(ddf_grid[\"%s_airmass\" % ddf_name] >= airmass_limit)] = 0\n",
    "\n",
    "    sky_mask = np.ones(ngrid, dtype=int)\n",
    "    if sky_limit is not None:\n",
    "        sky_mask[np.where(ddf_grid[\"%s_sky_g\" % ddf_name] <= sky_limit)] = 0\n",
    "        sky_mask[np.where(np.isnan(ddf_grid[\"%s_sky_g\" % ddf_name]) == True)] = 0\n",
    "\n",
    "    m5_mask = np.zeros(ngrid, dtype=bool)\n",
    "    m5_mask[np.isfinite(ddf_grid[\"%s_m5_g\" % ddf_name])] = 1\n",
    "\n",
    "    if g_depth_limit is not None:\n",
    "        m5_mask[np.where(ddf_grid[\"%s_m5_g\" % ddf_name] < g_depth_limit)] = 0\n",
    "\n",
    "    big_mask = sun_mask * airmass_mask * sky_mask * m5_mask\n",
    "\n",
    "    potential_nights = np.unique(night[np.where(big_mask > 0)])\n",
    "\n",
    "    # prevent a repeat sequence in a night\n",
    "    unights, indx = np.unique(night, return_index=True)\n",
    "    night_mjd = ddf_grid[\"mjd\"][indx]\n",
    "    # The season of each night\n",
    "    night_season = calc_season(ddf_RA, night_mjd)\n",
    "\n",
    "    # Mod by 1 to turn the season value in each night a simple 0-1 value\n",
    "    season_mod = night_season % 1\n",
    "    # Remove the tails of the seasons which are within \"season_unobs_fract\"\n",
    "    # Note that the \"season\" does include times where the field wouldn't\n",
    "    # really have been observable, due to being too close to the sun.\n",
    "    # Small values of season_unobs_frac may not reduce the usable season.\n",
    "    out_season = np.where((season_mod < season_unobs_frac) | (season_mod > (1.0 - season_unobs_frac)))\n",
    "\n",
    "    raw_obs = np.ones(unights.size)\n",
    "    raw_obs[out_season] = 0\n",
    "\n",
    "    if old:\n",
    "        night_season = night_season - np.floor(night_season).min()\n",
    "        cumulative_desired = ddf_slopes_old(ddf_name, raw_obs, night_season)\n",
    "    else:\n",
    "        cumulative_desired = ddf_slopes_new(ddf_name, raw_obs, night_season)\n",
    "\n",
    "    night_mask = unights * 0\n",
    "    night_mask[potential_nights] = 1\n",
    "\n",
    "    unight_sched = match_cumulative(cumulative_desired, mask=night_mask)\n",
    "    cumulative_sched = np.cumsum(unight_sched)\n",
    "\n",
    "    nights_to_use = unights[np.where(unight_sched == 1)]\n",
    "\n",
    "    # For each night, find the best time in the night.\n",
    "    # XXX--probably need to expand this part to resolve the times when\n",
    "    # multiple things get scheduled\n",
    "    mjds = []\n",
    "    for night_check in nights_to_use:\n",
    "        in_night = np.where((night == night_check) & (np.isfinite(ddf_grid[\"%s_m5_g\" % ddf_name])))[0]\n",
    "        m5s = ddf_grid[\"%s_m5_g\" % ddf_name][in_night]\n",
    "        # we could intorpolate this to get even better than 15 min\n",
    "        # resolution on when to observe\n",
    "        max_indx = np.where(m5s == m5s.max())[0].min()\n",
    "        mjds.append(ddf_grid[\"mjd\"][in_night[max_indx]])\n",
    "\n",
    "    return mjds, night_mjd, cumulative_desired, cumulative_sched\n",
    "\n",
    "def match_cumulative(cumulative_desired, mask=None, no_duplicate=True):\n",
    "    \"\"\"Generate a schedule that tries to match the desired cumulative\n",
    "    distribution given a mask\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    cumulative_desired : `np.array`, float\n",
    "        An array with the cumulative number of desired observations.\n",
    "        Elements  are assumed to be evenly spaced.\n",
    "    mask : `np.array`, bool or int (None)\n",
    "        Set to zero for indices that cannot be scheduled\n",
    "    no_duplicate : `bool` (True)\n",
    "        If True, only 1 event can be scheduled per element\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    schedule : `np.array`\n",
    "        The resulting schedule, with values marking number of events\n",
    "        in that cell.\n",
    "    \"\"\"\n",
    "\n",
    "    rounded_desired = np.round(cumulative_desired)\n",
    "    sched = cumulative_desired * 0\n",
    "    if mask is None:\n",
    "        mask = np.ones(sched.size)\n",
    "\n",
    "    valid = np.where(mask > 0)[0].tolist()\n",
    "    x = np.arange(sched.size)\n",
    "\n",
    "    drd = np.diff(rounded_desired)\n",
    "    step_points = np.where(drd > 0)[0] + 1\n",
    "\n",
    "    # would be nice to eliminate this loop, but it's not too bad.\n",
    "    # can't just use searchsorted on the whole array, because then there\n",
    "    # can be duplicate values, and array[[n,n]] = 1 means that extra\n",
    "    # match gets lost.\n",
    "    for indx in step_points:\n",
    "        left = np.searchsorted(x[valid], indx)\n",
    "        right = np.searchsorted(x[valid], indx, side=\"right\")\n",
    "        d1 = indx - left\n",
    "        d2 = right - indx\n",
    "        if d1 < d2:\n",
    "            sched_at = left\n",
    "        else:\n",
    "            sched_at = right\n",
    "\n",
    "        # If we are off the end\n",
    "        if sched_at >= len(valid):\n",
    "            sched_at -= 1\n",
    "\n",
    "        sched[valid[sched_at]] += 1\n",
    "        if no_duplicate:\n",
    "            valid.pop(sched_at)\n",
    "\n",
    "    return sched"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "32f7381c-f2a3-44d4-9964-799fc96e65d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_file = os.path.join(get_data_dir(), \"scheduler\", \"ddf_grid.npz\")\n",
    "mjd_start = survey_start_mjd()\n",
    "ddfs = ddf_locations()\n",
    "ddf_data = np.load(data_file)\n",
    "ddf_grid = ddf_data[\"ddf_grid\"].copy()\n",
    "\n",
    "mjd_max = mjd_start + 10. * 365.25\n",
    "\n",
    "in_range = np.where((ddf_grid[\"mjd\"] >= mjd_start) & (ddf_grid[\"mjd\"] <= mjd_max))\n",
    "ddf_grid = ddf_grid[in_range]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "c0b8e74d-a597-4070-92ee-e257eef0237c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddf_name = \"XMM_LSS\"\n",
    "mjdso, night_mjdo, cumulative_desiredo, cumulative_schedo = optimize_ddf_times(\n",
    "            ddf_name,\n",
    "            ddfs[ddf_name][0],\n",
    "            ddf_grid,\n",
    "            season_unobs_frac=0.1,\n",
    "            sequence_time=60,\n",
    "            old = True\n",
    "        )\n",
    "\n",
    "mjds, night_mjd, cumulative_desired, cumulative_sched = optimize_ddf_times(\n",
    "            ddf_name,\n",
    "            ddfs[ddf_name][0],\n",
    "            ddf_grid,\n",
    "            season_unobs_frac=0.1,\n",
    "            sequence_time=60,\n",
    "            old = False\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "50d54ffc-dc62-4c98-9b76-0f159ec084c0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "season #mjd_new  #mjd_old\n",
      "-1 27 9\n",
      "0 31 29\n",
      "1 30 30\n",
      "2 29 29\n",
      "3 31 31\n",
      "4 30 30\n",
      "5 30 30\n",
      "6 30 30\n",
      "7 29 29\n",
      "8 33 33\n",
      "9 0 0\n",
      "10 0 0\n",
      "300 280\n"
     ]
    }
   ],
   "source": [
    "#plt.plot(night_mjd, cumulative_sched)\n",
    "#print(cumulative_sched[-1])\n",
    "\n",
    "season_mjd = calc_season(ddfs[ddf_name][0], mjds)\n",
    "season_mjdo = calc_season(ddfs[ddf_name][0], mjdso)\n",
    "#plt.scatter(season_mjd, cumulative_sched)\n",
    "\n",
    "print(\"season #mjd_new  #mjd_old\")\n",
    "total = 0\n",
    "totalo = 0\n",
    "for season in np.arange(-1, 11, 1):\n",
    "    match = np.where(np.floor(season_mjd) == season)[0]\n",
    "    total += len(match)\n",
    "    matcho = np.where(np.floor(season_mjdo) == season)[0]\n",
    "    totalo += len(matcho)\n",
    "    print(season, len(match), len(matcho))\n",
    "print(total, totalo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "192e91f3-342e-4ce8-8817-80d34bfd9219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x19fc7e890>]"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWqklEQVR4nO3deVhUZePG8S+yqQgoLiyJioorLqjlWmpumWZWLmn1s7I3fV1JrbRVW7SstMWy9VVb1HK30lJbKLNSARXQXHJXEBdk3zm/P840BmaKM8gA9+e6uN73mTk8nDkXMbfnPs8ZJ8MwDEREREQcSIWS3gERERGRwhRQRERExOEooIiIiIjDUUARERERh6OAIiIiIg5HAUVEREQcjgKKiIiIOBwFFBEREXE4LiW9A1cjPz+fkydP4unpiZOTU0nvjoiIiFwBwzBISUkhICCAChX+/RxJqQwoJ0+eJDAwsKR3Q0RERK7CsWPHqF279r9uUyoDiqenJ2C+QC8vrxLeGxEREbkSycnJBAYGWt/H/02pDCh/1TpeXl4KKCIiIqXMlVyeoYtkRURExOEooIiIiIjDUUARERERh6OAIiIiIg5HAUVEREQcjgKKiIiIOBwFFBEREXE4CigiIiLicBRQRERExOEUOaCcOHGCe++9l+rVq1O5cmVat25NRESE9XnDMJg+fToBAQFUqlSJbt26ERsbW2COrKwsxo8fT40aNfDw8GDAgAEcP37c9lcjIiIiZUKRAkpiYiKdO3fG1dWV9evXs3v3bl577TWqVq1q3Wb27NnMmTOHefPmsW3bNvz8/OjVqxcpKSnWbcLCwli1ahVLly5l8+bNpKam0r9/f/Ly8uz2wkRERKT0cjIMw7jSjadOncovv/zCzz///I/PG4ZBQEAAYWFhPP7444B5tsTX15eXX36ZUaNGkZSURM2aNfnkk08YOnQocOHTidetW0efPn0uux/Jycl4e3uTlJSkz+IREREpJYry/l2kMyhr166lXbt2DB48mFq1ahEaGsoHH3xgff7QoUPEx8fTu3dv62Pu7u507dqVLVu2ABAREUFOTk6BbQICAggJCbFuU1hWVhbJyckFvkRERMSOko7DxwMhYU9J7wlQxIBy8OBB5s+fT3BwMN9++y2jR49mwoQJfPzxxwDEx8cD4OvrW+D7fH19rc/Fx8fj5uZGtWrVLrlNYbNmzcLb29v6FRgYWJTdFhERkX+z9xt4twsc/AG+nAhXXq4UmyIFlPz8fNq0acPMmTMJDQ1l1KhR/Oc//2H+/PkFtiv8McqGYVz2o5X/bZtp06aRlJRk/Tp27FhRdltERET+SV4ObHgalgyFjEQICIU73oPLvGdfC0UKKP7+/jRr1qzAY02bNuXo0aMA+Pn5AVx0JiQhIcF6VsXPz4/s7GwSExMvuU1h7u7ueHl5FfgSERERG5w/BgtuhS1vmuP2o+HBb8EnqGT3y6JIAaVz587s3bu3wGP79u2jbt26AAQFBeHn58fGjRutz2dnZxMeHk6nTp0AaNu2La6urgW2iYuLIyYmxrqNiIiIFKO938B7N8LxreDuDUM+gb4vg4t7Se+ZlUtRNn7kkUfo1KkTM2fOZMiQIWzdupX333+f999/HzCrnbCwMGbOnElwcDDBwcHMnDmTypUrM3z4cAC8vb0ZOXIkkydPpnr16vj4+DBlyhRatGhBz5497f8KRURExJSXA9/NgC1vmeOANjB4AVSrV6K79U+KFFCuv/56Vq1axbRp03juuecICgri9ddf55577rFu89hjj5GRkcGYMWNITEykffv2bNiwAU9PT+s2c+fOxcXFhSFDhpCRkUGPHj1YuHAhzs7O9ntlIiIicsH5o7D8QTi+zRx3GAM9Z4CLW8nu1yUU6T4ojkL3QRERESmCveth1WjIPG9WOgPfgab9r/luFOX9u0hnUERERKQUycuBTdPh13nm2IErncIUUERERMqiUlbpFKaAIiIiUtb8sQ5W/9esdCp6w+0lU+nYQgFFRESkrMjNNlfp/FXpXNcWBi2AanVLdr+uggKKiIhIWZB4BJY/ACcizHHHcdDj2VJT6RSmgCIiIlLa/fG1pdJJMiudgfOhSb+S3iubKKCIiIiUVrnZ5iqd3942x9e1g0H/K5WVTmEKKCIiIqVRGat0ClNAERERKW32fAVrxlgqnaqWSufWkt4ru1JAERERKS1ys2HjM/D7fHNc+3qz0qlap2T3qxgooIiIiJQGiYdh2QNwMtIcdxpvVjrOriW6W8VFAUVERMTR7fkSVo+FLEulc8e70LhvSe9VsVJAERERcVS5WZZK511zXPt688ZrVQNLdr+uAQUUERERR3TukLlK52SUOe40AXo8U2YrncIUUERERBzN7rWwZpxZ6VSqBne8B436FPuPzcrNI3zvab7aFYePhxvTBzQv9p95KQooIiIijiI3CzY8DVvfM8eB7c1VOt61i+1H5uTls/nAGb7ceZKNsadIycoFwKuiC0/c2hQ3lwrF9rP/jQKKiIiIIzh3CJbdD3E7zHHniXDz08VS6eTm5fPbwXN8tesk38TGcz49x/qcn1dFbm3hz22t/HF1drL7z75SCigiIiIlLXY1rB0PWcmWSud9aNTbrj8iP99g62FLKImJ50xqtvW5GlXcuLWFP/1bBtCubjUqVCi5YPIXBRQREZGSkpsFG56Cre+bYztXOoZhEHn0PF/tOsm66DhOJWdZn6ta2ZW+If7c1tKf9vWr4+wAoeTvFFBERERKwrmD5o3XrJVOGNz8lF0qnZgTSazdeZKvd8Vx4nyG9XHPii70ae5H/5b+dG5YA1fnkrm+5EoooIiIiFxrBSodH8sqHdsrncycPJ7/ajef/X7U+piHmzM9m/lyW8sAbmxUA3cXZ5t/zrWggCIiInKt5GSalc62D8xxYAdLpXOdzVMfPJ3K2MVR7IlLBqBviB8DWgXQvUktKrqWjlDydwooIiIi18LZP81VOvG7zHGXSdD9SXC2/a14zY4TPLEymrTsPHw83Jg7tDVdG9W0ed6SpIAiIiJS3GJWwtoJkJ1iVjp3fgDBPW2eNjMnjxlf7mbJVrPSuSHIhzfvDsXPu6LNc5c0BRQREZHikpMJG56EbR+a4zod4a6P7FLp/Hk6lbGfRfJHfApOTjCue0Mm9gjGxYEvfC0KBRQREZHicI0qneqWSuemUl7pFKaAIiIiYm8xK2DtRLPSqVzdvPGa3SqdWJZsPQZA+yAf3hwWiq9X6a90ClNAERERsZecTPj2Cdj+kTmu0wkGfQReATZPfSAhlXGLL1Q647s3ZEIZqnQKU0ARERGxh7N/wrIREB8NOMGNk6HbNLtUOqujTvDEqmjSs/OoUcWN14eG0iW4hu377MAUUERERGwVvRy+nAjZqWalc+f70ND2SicjO4/pa2P5fLtZ6XSob67SqVUGK53CFFBERESuVk4GfDMNIhaY47qd4a4P7VbpjP0skr2nzEpnws3BTOgR7HCfmVNcFFBERESuxpkD5iqdU5ZK56Yp0HWqXSqdlZHHeWp1jKXSceeNu1vTuWHZrnQKU0AREREpqgKVTg1LpdPD5mkzsvN4dm0MX2w/DkCnBtV5/e7W1PIs+5VOYQooIiIiVyonA76ZChELzXHdLpZKx9/mqfefSmHs4kj2nUrFyQkm9ghm/M3lp9IpTAFFRETkSpzZb6l0YjArnUeh6+N2qXRWRJiVTkaOWem8eXdrOpWzSqcwBRQREZHL2bUMvgq7UOnc9QE0uNnmadOzc3lmTSzLI8xKp3PD6swdWj4rncIUUERERC4lJwPWPw6Ri8xxvRvND/qzU6Uz5rNI9iekUsEJwno2Ymz3huW20ilMAUVEROSfnNkPX4yAhFjACbo+ZlY6FZxtnnrZ9mM8syaWjJw8anqaq3Q6NSjflU5hCigiIiKF7foCvgyDnDTwqGmeNWnQ3eZp07NzeXp1LCsizUqnS8MazB3ampqe7jbPXdYooIiIiPwlJwPWPwaRH5vjejeaq3Q8/Wyeep+l0jlgqXQe6dmIMap0LkkBRUREBOD0PnOVjrXSedysdWysdAzDYFnEcZ5ZE0NmTj61PN15c1goHepXt8tul1UKKCIiIjs/h68esVQ6tcxVOvW72TxtWlYuT6+OYWXUCQBuDDYrnRpVVOlcjgKKiIiUX9npZqUT9Yk5DroJ7vwQPH1tnnpvfApjPovgz9NpVHCCyb0b89+uDaigSueKKKCIiEj5dHofLBsBCbsBJ+g21bz5mh0qnS+2H+PZtbFk5uTj6+XOm3eH0l6VTpEooIiISPmzcyl8Nelvlc6HUL+rzdOmZeXy1OoYVlkqnZsa1WTukFZUV6VTZBWKsvH06dNxcnIq8OXnd+HKZsMwmD59OgEBAVSqVIlu3boRGxtbYI6srCzGjx9PjRo18PDwYMCAARw/ftw+r0ZEROTfZKfDmrGwapQZToK6wujNdgknf8Qnc9u8zayKOkEFJ3i0T2MW3n+9wslVKlJAAWjevDlxcXHWr+joaOtzs2fPZs6cOcybN49t27bh5+dHr169SElJsW4TFhbGqlWrWLp0KZs3byY1NZX+/fuTl5dnn1ckIiLyTxL+gA9uhqhPwakCdH8S7ltl8/UmhmGwdOtRbp/3CwdPp+HnVZGlD3dkbPeGut7EBkWueFxcXAqcNfmLYRi8/vrrPPnkk9x5550ALFq0CF9fXxYvXsyoUaNISkrio48+4pNPPqFnz54AfPrppwQGBrJp0yb69Olj48sRERH5BzuWwNeTICcdqvialU7QTTZPm5qVy5Orolmz4yQAXRvVZI4qHbso8hmU/fv3ExAQQFBQEHfffTcHDx4E4NChQ8THx9O7d2/rtu7u7nTt2pUtW7YAEBERQU5OToFtAgICCAkJsW7zT7KyskhOTi7wJSIiclnZabB6DKwebYaT+t3MSscO4WRPXDID3trMmh0nca7gxOO3NGGBKh27KdIZlPbt2/Pxxx/TqFEjTp06xQsvvECnTp2IjY0lPj4eAF/fgqfKfH19OXLkCADx8fG4ublRrVq1i7b56/v/yaxZs5gxY0ZRdlVERMq7hD/MVTqn/zArnW7T4MbJdlmls2TrMWZ8GUtWbj5+XhV5a3go19fzsdOOCxQxoPTt29f6/1u0aEHHjh1p0KABixYtokOHDgA4ORXs2wzDuOixwi63zbRp05g0aZJ1nJycTGBgYFF2XUREypMdi+HryX+rdD6CoBttnjY1K5cnVkazdqdZ6XRvXJPXhrTGx8PN5rmlIJuWGXt4eNCiRQv279/PwIEDAfMsib//hY+hTkhIsJ5V8fPzIzs7m8TExAJnURISEujUqdMlf467uzvu7jplJiIil5GdBusehR2fmeP63eHO96FKLZun3n0ymbGLIzl0Jg3nCk482qcxD99YXxfCFpMiX4Pyd1lZWezZswd/f3+CgoLw8/Nj48aN1uezs7MJDw+3ho+2bdvi6upaYJu4uDhiYmL+NaCIiIhcVsIec5XOjs/MSufmp+DelTaHE8Mw+Oz3Iwx85xcOnUnD37sinz/cgdG6K2yxKtIZlClTpnDbbbdRp04dEhISeOGFF0hOTmbEiBE4OTkRFhbGzJkzCQ4OJjg4mJkzZ1K5cmWGDx8OgLe3NyNHjmTy5MlUr14dHx8fpkyZQosWLayrekRERIos6jOz0snNgCp+MOgjqNfF5mlTMnOYtjKar3bFAXBzk1q8NrgV1VTpFLsiBZTjx48zbNgwzpw5Q82aNenQoQO//fYbdevWBeCxxx4jIyODMWPGkJiYSPv27dmwYQOenp7WOebOnYuLiwtDhgwhIyODHj16sHDhQpydbbtoSUREyqHsNDOY7Fxijut3hzs/gCo1bZ465kQS4xZHcvhsumWVTmMe6qJK51pxMgzDKOmdKKrk5GS8vb1JSkrCy8urpHdHRERKwqndsOx+OLP3wo3XukyCCjZdvYBhGHz6+1Ge/2o32bn5BHibq3Ta1tUqHVsV5f1bn8UjIiKli2GY15l8PcWsdDz9zVU69TrbPHVKZg5TV0bztaXS6dGkFq+q0ikRCigiIlJ6ZKWalc6upea4QQ9zlY5HDZunjjmRxNjFkRw5m46L5cZrD90YdNlbZUjxUEAREZHS4dRu88ZrZ/ZdWKXT+RH7VDq/HeH5r/aQnZfPdVUr8dbwUNrUqXb5b5Zio4AiIiKOzTAg6hNY99iFSmfQ/6Cu7benSM7MYeqKXayLNu9m3rOpWelUraxKp6QpoIiIiOPKSjU/5G/X5+bYjpVO9HGz0jl6zqx0pvZtwsguqnQchQKKiIg4plOx8MUIOLsfnJwtlU6YXSqdj389wotfX6h05g0PJVSVjkNRQBEREcdiGBD5Max/DHIzwTPAUul0tHnq5MwcHl++i/UxZqXTq5kvrw5qhXdlV5vnFvtSQBEREceRlQpfPQLRX5jjhr3gjvfAo7rNU+86fp5xi6M4ei4dV2cnpvZtyoOd66nScVAKKCIi4hjiY8wbr/1V6fR4GjpNtEuls2jLYV5ct4ecPIPa1Soxb3gbWgdWtctuS/FQQBERkZJlGBC5CNY/fqHSGbwA6nSweeqkDLPS+SbWrHR6N/PlFVU6pYICioiIlJysFEuls8wcB/eGge/apdLZeew845ZEcuxcBq7OTjxxa1Pu76RKp7RQQBERkZIRH22pdA5YKp1noNMEu1Q6C345zKz1Fyqdt4e3oZUqnVJFAUVERK4tw4CIhWalk5cFXteZq3TsUemk5/Do8p1s2H0KgFua+/HyoJZ4V1KlU9oooIiIyLWTlQJfhkHMcnMc3AfueBcq2/5JwTuOnWfc4kiOJ5qVzpO3NmWEKp1SSwFFRESujfho88Zr5/40K52e06HjOLtUOv/75TAvWSqdQB+z0mlZu6pddltKhgKKiIgUL8OAiAWwfqql0qltqXTa2zx1UnoOU5bvZKOl0ukb4sdLd6nSKQsUUEREpPhkJsOXEyF2pTm2Y6UTdTSRcYujOHE+AzfnCjzVvyn3dairSqeMUEAREZHiEbfTXKVz7iBUcIEez9qt0vlo8yFeWv8HufkGdXwq8/bwNrSo7W2f/RaHoIAiIiL2ZRiw/X/wzbQLlc7gBRB4g81Tn0/PZsqyXWzaY1Y6/Vr4M+uuFnhVVKVT1iigiIiI/RSudBr1hYHv2KXSiTyayPi/VTpP92/Kvap0yiwFFBERsY/ClU7PGdBxLNgYIAzD4MOfD/HyN2alU7e6WemEXKdKpyxTQBEREdsYBmz/yFLpZIN3IAxaAIHX2zx1Ylo2U5bt5Ls/EgDo19Kfl+5sgacqnTJPAUVERK5eZhKsnQC7V5tjO1Y6EUcSGb84kpNJmbi5VOCZ/s24p30dVTrlhAKKiIhcnZM7zEon8ZBZ6fR6DjqMsbnSyc83+ODng7zy7V5y8w3qVa/MPFU65Y4CioiIFI1hwLYP4dsnLJVOHXOVTu12Nk+dmJbN5GU7+d5S6fRv6c8sVTrlkgKKiIhcucwkWDsedq8xx437wcC3oVI1m6eOOHKOcYujiLNUOs/e1ozhN6jSKa8UUERE5MqcjLJUOoehgqul0vmv3Sqd2d/uJS/fIKiGB/OGh9I8QJVOeaaAIiIi/84wYOsHsOFJs9KpWgcGLYTabW2e+lxaNpO/2MEPe08DMKBVADPvbEEVd709lXf6DRARkUvLTII142DPWnPcpD/cPs8ulc72w+cYv+RCpTNjQHPuvj5QlY4ACigiInIpJyJh+QMXKp3ez0P70XapdN776SCvbjArnfo1PJg3vA3NArzss99SJiigiIhIQYYBW9+Hb5+E/Byz0hm8EK6zT6Uz6Ysd/GipdG5vHcCLd6jSkYvpN0JERC7IOA9rx8GeL81xk/5w+9tQqarNU287fI7xi6OIT87E3VLpDFWlI5eggCIiIqYTkeYqnfNHLJXOC9B+lF0qnXd/+pPXNuwzK52aHrw9vA1N/VXpyKUpoIiIlHeGAb+/BxueslQ6dc0br9mh0jmbmsWkL3YSvs+sdO4IvY4XBobgoUpHLkO/ISIi5VnGeVgzFv74yhw3vQ0GzLNLpbP10DnGL4nkVHIW7i4VeP72EAa3q61KR66IAoqISHl1PAKW3w/nj5qVTp8X4YaH7VLpzA//k9c27CXfgAY1PXj7njY08VOlI1dOAUVEpLwxDPj9Xdjw9N8qnYVwXRubpz6TmsUjn+/g5/1nALgz9DqeV6UjV0G/MSIi5UlGonnjNWulMwAGvGWXSue3g2eZsCSKhJQsKrpW4LnbQxjcVpWOXB0FFBGR8uLvlY6zG/SZCdc/ZJdK550fDzBn4z7yDWhYqwrv3NOGRr6e9tlvKZcUUEREyjrDgN/mw8ZnzEqnWj2z0gkItXnqiyqdNuYqncpuensR2+g3SESkLMtIhNVjYe/X5rjZ7WalU9H2Twr+9c+zTFx6odIxV+kE2jyvCCigiIiUXce3w7IHIMm+lU5evsHbPxzg9U2qdKT4KKCIiJQ1hgG/vg2bnoX8XKgWZKl0Wts89ekUs9LZfMCsdAa1rc1ztzdXpSN2p98oEZGyJP2ceeO1vevMcbOBMOBNu1Q6W/48w8SlOzidkkUlV2eeHxjCoLa1bZ5X5J8ooIiIlBXHtsHyByDpmFnp3DIL2o20S6Uz7/sDvPGdWekEWyqdYFU6Uowq2PLNs2bNwsnJibCwMOtjhmEwffp0AgICqFSpEt26dSM2NrbA92VlZTF+/Hhq1KiBh4cHAwYM4Pjx47bsiohI+WUYsGUeLLjFDCc+9eGhTXa53iQhJZP7PvqduZbrTQa3rc2acZ0VTqTYXXVA2bZtG++//z4tW7Ys8Pjs2bOZM2cO8+bNY9u2bfj5+dGrVy9SUlKs24SFhbFq1SqWLl3K5s2bSU1NpX///uTl5V39KxERKY/Sz8GSYbDhSfN6k+Z3wMPh4N/K5qm3HDjDrW9sZsufZ6nk6sxrg1vxyuBWut5EromrCiipqancc889fPDBB1SrVs36uGEYvP766zz55JPceeedhISEsGjRItLT01m8eDEASUlJfPTRR7z22mv07NmT0NBQPv30U6Kjo9m0aZN9XpWISHlwbCu8dxPsWw/O7tBvDgxaABVt+8ybvHyDuRv3cc9Hv3MmNYtGvlX4cnxn7tL1JnINXVVAGTt2LP369aNnz54FHj906BDx8fH07t3b+pi7uztdu3Zly5YtAERERJCTk1Ngm4CAAEJCQqzbFJaVlUVycnKBLxGRcis/H355Exb0/VulsxGut/16k4SUTO798Hfe+G4/hgFD2wWyZmwXGtZSpSPXVpHP0y1dupTIyEi2bdt20XPx8fEA+Pr6Fnjc19eXI0eOWLdxc3MrcOblr23++v7CZs2axYwZM4q6qyIiZU/6OVj9X9j3jTlufifc9obNZ00ANu8/Q9jnUZxJzaaymzMv3hHCHaE6ayIlo0gB5dixY0ycOJENGzZQsWLFS25X+IOhDMO47IdF/ds206ZNY9KkSdZxcnIygYG6W6GIlDPHtpo3Xks+blY6fV+Ctg/YZZXOG9/t563vzbMmjX09efueNjSsVcVOOy5SdEUKKBERESQkJNC2bVvrY3l5efz000/MmzePvXv3AuZZEn9/f+s2CQkJ1rMqfn5+ZGdnk5iYWOAsSkJCAp06dfrHn+vu7o67u3tRdlVEpOzIz4df34LvnjMvhPVpYN54zb/lZb/1chKSM5mwNIrfDp4D4O7rA3n2tuZUcnO2eW4RWxTpGpQePXoQHR3Njh07rF/t2rXjnnvuYceOHdSvXx8/Pz82btxo/Z7s7GzCw8Ot4aNt27a4uroW2CYuLo6YmJhLBhQRkXIr/RwsudvyQX+5EHIXjAq3Szj5ef9pbn3zZ347eI7Kbs68PrQ1L93VUuFEHEKRzqB4enoSEhJS4DEPDw+qV69ufTwsLIyZM2cSHBxMcHAwM2fOpHLlygwfPhwAb29vRo4cyeTJk6levTo+Pj5MmTKFFi1aXHTRrYhIuXb0d/PGa8knLJXOy9D2fpsrndy8fN74bj/zfjiAYUATP7PSaVBTlY44DrsvZn/sscfIyMhgzJgxJCYm0r59ezZs2ICn54UrwOfOnYuLiwtDhgwhIyODHj16sHDhQpydldpFRMjPhy1vmpWOkQfVG5qVjl8Lm6c+lZzJhCVR/H7IrHSG3VCHZ29rRkVX/f0Vx+JkGIZR0jtRVMnJyXh7e5OUlISXl+1XrouIOIy0s7B6NOzfYI5DBsFtr4O77ct8f9p3mkc+38HZtGw83JyZeWcLbm99nc3zilyporx/63aAIiKO4uhvsPxBs9JxqWhWOm1G2KXSeX3Tft7+0ax0mvp78fbwUOqr0hEHpoAiIlLS8vNhyxvw3fN/q3QWgV/I5b/3MuKTzEpn62Gz0hnevg7P9FelI45PAUVEpCSlnYVVo+CAZWVji8HQf65dKp1wS6VzzlLpzLqrJQNaBdg8r8i1oIAiIlJSjvxqVjopJy2Vzmxo8392qXTmbNzHOz/+CZiVzjv3tCGohoc99lrkmlBAERG51vLz4ZfX4fsXLJVOsGWVju2VTlxSBhOWRLHtcCIA93aow1P9VOlI6aOAIiJyLaWdsVQ6lk9vbznU/BRid9svWP1hbwKTPt9BYnoOVdxdeOmuFvRvqUpHSicFFBGRa+XIFlg+8kKlc+urEHqvXSqd1zbuY76l0mke4MXbw9tQT5WOlGIKKCIixS0/H36ZC9+/aFY6NRqZlY5vc5unPnnerHS2HzErnfs61OXJfk1V6Uipp4AiIlKc0s7Ayofhz+/MsT0rnT8SmPTFhUrn5bta0q+l/+W/UaQUUEARESkuh3+BFSMhJQ5cKsGtr9il0snJy+fVDXt5L/wgACHXeTFvmCodKVsUUERE7C0/Hza/Bj/MBCPfUuksAt9mNk998nwG45dEEWGpdEZ0rMsT/Zri7qJKR8oWBRQREXtKPQ2rHoY/vzfHrYaZF8PaodL5bs8pJi/byfn0HDzdXXh5UEtubaFKR8omBRQREXs5vNlcpZMab1Y6/V6D0HtsnjYnL59Xv93Lez+ZlU6L67yZNzyUutVV6UjZpYAiImKriyqdxjBkEdRqavPUJ85nMH5xJJFHzwNwf6d6TLu1iSodKfMUUEREbJF6Glb+Bw7+YI5bDYd+r4Kb7Wc3Nu02K52kjBw8K7rwyqCW3BKiSkfKBwUUEZGrdehnWPFQsVQ6s7/5gw9+PgRAy9rezBvWhjrVK9s8t0hpoYAiIlJU+Xnw82vw4yyz0qnZxLzxmh0qneOJ6YxfEkWUpdJ5oHM9pvZVpSPljwKKiEhRpCZYKp0fzXHre+HW2XapdDbuPsWUApVOK24J8bN5XpHSSAFFRORKHfrJUumcAtfKZqXTerjN02bnmpXOh5vNSqdVbW/mDW9DoI8qHSm/FFBERC4nPw9+ehXCX7JUOk0tlU4Tm6c+di6dcUui2HnsPAAPdg5iat8muLlUsHlukdJMAUVE5N+kJphnTQ6Fm+PQe6HvK+Bm+9mNDbHxTFm2k+TMXLwquvDK4Fb0aa5KRwQUUERELu1guHm9ibXSmQOth9k8bXZuPi+t/4P//WKpdAKrMm9YqCodkb9RQBERKSw/D356BX58CTDMSmfIIqjZ2Oapj51LZ9ziSHYeTwLgoS5BPHaLKh2RwhRQRET+LuUUrHzIvCAWIPQ+6DvbLpXONzHxPLp8JymZuXhXcuXVwa3o1czX5nlFyiIFFBGRvxz8EVb8B9ISwNUD+s+FVkNtnjY7N59Z6/ew4JfDALQOrMq84aHUrqZKR+RSFFBERPLzIHw2hL8MGFCrGQxeBDUb2Tz10bPpjFsSyS5LpfOfG4N4tI8qHZHLUUARkfIt5RSsGAmHfzbHbf4PbnnZTpVOHI8u32WtdF4b3IqeqnRErogCioiUX3/+YK7SSTtt10onKzePWev+YOGWwwCE1qnKvOFtuK5qJZvnFikvFFBEpPzJzzPrnPDZmJVOc/PGa3aqdMYujiT6hFnpjLqpPlP6NMbVWZWOSFEooIhI+ZISb954zVrpjIC+L4Or7Wc31kXH8fjyXaRk5VK1sitzhrTi5iaqdESuhgKKiJQff34PKx82Kx23KtD/dWg52OZps3LzmPn1Hhb9egSAtnWr8dawUAJU6YhcNQUUESn78vPMm6799ApggG+IWenUCLZ56iNn0xi7OJKYE8kAjOpanym9VemI2EoBRUTKtuQ4s9I5stkct70fbnnJLpXO17vimLrCrHSqVXZlzpDWdG9Sy+Z5RUQBRUTKsgPfmZVO+hmz0rntDWgxyOZpM3PyePHrPXzym1nptKtbjTdV6YjYlQKKiJQ9ebnw4yz4+TXMSqeFpdJpaPPUh8+YlU7sSbPS+W+3Bkzq1UiVjoidKaCISNmSHGfeeO3IL+a47QNwyyy7VDpf7jzJtJXRpP5V6QxtTffGqnREioMCioiUHcVY6bzw9W4+/e0oANfXMysdf29VOiLFRQFFREq/f6p0hiyC6g1snvrQmTTGfhbJ7jiz0hljqXRcVOmIFCsFFBEp3ZJPwvKRcHSLOW43EvrMBNeKNk+9dudJpq3YRVp2Hj4ebswZ0opuqnRErgkFFBEpvQ5sslQ6Z8HNEwa8ASF32TxtZk4ez321m8W/m5XODfV8eHNYKH7etoceEbkyCigiUvrk5cIPL8LmOebYrwUMtk+lc/B0KmMXR7EnLhknJxjbrSFhPYNV6YhcYwooIlK6JJ0wV+kc/dUcX/8Q9H7RLpXOmh0neGJlNGnZeVT3cGPu0Nbc1KimzfOKSNEpoIhI6bF/E6z6e6XzJoTcafO0mTl5zPhyN0u2mpVO+yCz0vH1UqUjUlIUUETE8eXlwg8vwOa55tivpXnjNTtUOn+eTmXsZ5H8EZ+CkxOM696QiT1U6YiUNAUUEXFsF1U6/4HeL9il0lkddYInVkWTbql0Xr+7NTcGq9IRcQRF+ifC/PnzadmyJV5eXnh5edGxY0fWr19vfd4wDKZPn05AQACVKlWiW7duxMbGFpgjKyuL8ePHU6NGDTw8PBgwYADHjx+3z6sRkbJl/0Z4t4sZTtw8zbMm/V61OZxk5uQxdcUuwj7fQXp2Hh3q+7Bu4o0KJyIOpEgBpXbt2rz00kts376d7du3c/PNN3P77bdbQ8js2bOZM2cO8+bNY9u2bfj5+dGrVy9SUlKsc4SFhbFq1SqWLl3K5s2bSU1NpX///uTl5dn3lYlI6ZWXAxufhc8GQcY58G8Fo8Kh+R02T30gIZWBb//C0m3HcHKCCT2C+eyhDrreRMTBOBmGYdgygY+PD6+88goPPvggAQEBhIWF8fjjjwPm2RJfX19efvllRo0aRVJSEjVr1uSTTz5h6NChAJw8eZLAwEDWrVtHnz59ruhnJicn4+3tTVJSEl5eXrbsvog4mqTjsPxBOPa7Ob7hYbPScXG3eepVUcd5clUM6dl51KjixutDQ+kSXMPmeUXkyhTl/fuqrwLLy8tj6dKlpKWl0bFjRw4dOkR8fDy9e/e2buPu7k7Xrl3ZssW8w2NERAQ5OTkFtgkICCAkJMS6zT/JysoiOTm5wJeIlEH7vjUrnWO/g7uXeW+TW1+xOZxkZOfx+PJdPPL5TtKz8+hYvzrrJtyocCLiwIp8kWx0dDQdO3YkMzOTKlWqsGrVKpo1a2YNGL6+vgW29/X15ciRIwDEx8fj5uZGtWrVLtomPj7+kj9z1qxZzJgxo6i7KiKlRV4OfP88/PKGOfZvDYMXgE99m6c+kJDC2M+i2HvKXKUz4eZgJvQIxrmCk81zi0jxKXJAady4MTt27OD8+fOsWLGCESNGEB4ebn3eyangf/SGYVz0WGGX22batGlMmjTJOk5OTiYwMLCouy4ijuiiSmcU9H7eLpXOiojjPLU6hoycPGpUceeNu1vTuaHOmoiUBkUOKG5ubjRs2BCAdu3asW3bNt544w3rdSfx8fH4+/tbt09ISLCeVfHz8yM7O5vExMQCZ1ESEhLo1KnTJX+mu7s77u62/7ESEQez9xtYPRoyEsHdG26fB80G2DxtRnYez6yJYVmEuUKwU4PqvH53a2p56kJYkdLC5jsRGYZBVlYWQUFB+Pn5sXHjRutz2dnZhIeHW8NH27ZtcXV1LbBNXFwcMTEx/xpQRKSMycuBDU/DkqFmOAkINVfp2CGc7D+Vwu1vb2ZZxHGcnOCRno34ZGR7hRORUqZIZ1CeeOIJ+vbtS2BgICkpKSxdupQff/yRb775BicnJ8LCwpg5cybBwcEEBwczc+ZMKleuzPDhwwHw9vZm5MiRTJ48merVq+Pj48OUKVNo0aIFPXv2LJYXKCIO5vwxs9I5vtUctx8NvZ6zS6WzPOI4T1sqnZqeZqXTqYEqHZHSqEgB5dSpU9x3333ExcXh7e1Ny5Yt+eabb+jVqxcAjz32GBkZGYwZM4bExETat2/Phg0b8PT0tM4xd+5cXFxcGDJkCBkZGfTo0YOFCxfi7Oxs31cmIo5n73pYNRoyz9u10knPzuWZNbEst1Q6XRrWYO7Q1tT0VDUsUlrZfB+UkqD7oIiUMnk58N0M2PKWOQ5oY67SqVbP5qn3nUph7GeR7E9IpYIThPVsxNjuDbVKR8QBFeX9W5/FIyLF6/xRS6WzzRx3GAM9Z4CLm81TL9t+jKfXxJCZk09NT3fevDuUjg2q2zyviJQ8BRQRKT6FK52B70DT/jZPm56dy1OrY1gZeQKAG4NrMGeIKh2RskQBRUTsLzfbrHR+nWeO7Vjp7I1PYeziSA5YKp1JvRoxpltDKqjSESlTFFBExL7OH4VlD8CJ7ebYTpWOYRgs236cZ9aalU4tT3feHBZKh/qqdETKIgUUEbGfP76G1f+FzCSo6A2326fSScvK5enVMayMulDpzB3amhpVVOmIlFUKKCJiu9xs2DQdfnvbHF/XFgYtgGp1bZ76j/hkxn4WyZ+n06jgBJN7N+a/XRuo0hEp4xRQRMQ2iUdg+QNwIsIcdxwHPZ61S6XzxfZjPLMmlqzcfHy9zFU67VXpiJQLCigicvUKVzoD50OTfjZPm5ZlrtJZZal0bmpUk7lDWlFdlY5IuaGAIiJFl5sNm56F394xx9e1g0H/s1ulM+azSA6eTsO5ghOTezdi9E2qdETKGwUUESmaxMPmKp2TkebYjpXO59uO8exas9Lx86rIW8NDub6ej+37LCKljgKKiFy5PV/C6rGQlQQVq1oqnVttnjY1K5cnV0WzZsdJALo1rsmcIa3x8bD9brMiUjopoIjI5eVmw8Zn4Pf55rj29WalU7WOzVPviTNX6Rw8Y1Y6U3o3ZtRN9VXpiJRzCigi8u8KVzqdxpuVjrOrTdMahsGSrceY8aVZ6fh7V+StYaG0U6UjIiigiMi/KVzp3PEuNO5r87SpWbk8sTKatTvNSqd745q8pkpHRP5GAUVELpabZal03jXHta83b7xWNdDmqXefTGbs4kgOWSqdx/o05j83qtIRkYIUUESkoHOHYNn9ELfDHHeaAD2esUuls3jrUWZ8uZvs3HwCvM1VOm3rqtIRkYspoIjIBbvXwJpxkJUMlarBHe9Boz42T5uSmcO0ldF8tSsOgJub1OK1wa2opkpHRC5BAUVEzEpnw9Ow9T1zHNjeXKXjXdvmqWNOJDFucSSHz6bjUsGJx25pzENdVOmIyL9TQBEp7wpXOp0nws1P26XS+fT3ozz/1d8rnTa0rVvN5l0WkbJPAUWkPItdDWvH/63SeR8a9bZ52pTMHKaujOZrS6XTs2ktXh3ciqqVVemIyJVRQBEpj3KzYMNTsPV9c2znSmfs4kiOWCqdqX2bMLJLEE5OqnRE5MopoIiUN+cOWiqdnea4cxjc/JR9Kp3fjvD8V3vIzsvnuqqVeGt4KG3qqNIRkaJTQBEpT2JXwdoJlkrHx7JKx/ZKJzkzh6krdrEuOh6Ank19eXVwS1U6InLVFFBEyoOcTNjwJGz70BwHdrBUOtfZPHX0cbPSOXpOlY6I2I8CikhZd/ZPs9KJ32WOu0yC7k+Cs23/+RuGwce/HuHFry9UOvOGhxKqSkdE7EABRaQsi1lpVjrZKWalc+cHENzT5mmTM3N4fPku1seYlU6vZr68OqgV3pVtu45FROQvCigiZVFOJnz7BGz/yBzX6Qh3fWSXSmfX8fOMWxzF0XPpuDo7Ma1vUx7oXE+VjojYlQKKSFlz9k9YNgLio82xHSudRVsO8+K6PeTkGdSuVom3h7ehVWBV2/dZRKQQBRSRsiRmBaydaFY6laubN16zQ6WTlGFWOt/EmpVOn+a+zB7UCu9KqnREpHgooIiUBTmZ8O002P4/c1ynEwz6CLwCbJ5657HzjFsSybFzGbg6O/HErU25v5MqHREpXgooIqXd2T/hixFwKhpwghsnQ7dpdql0FvxymFnrzUon0KcS84ap0hGRa0MBRaQ0i14OX06E7FSz0rnzfWhoh0onPYdHl+9kw+5TANzS3I+XB7VUpSMi14wCikhplJMB30yDiAXmuG5nuOtDu1Q6O46dZ9ziSI4nZuDmXIEn+zXl/zrWVaUjIteUAopIaXPmgHnjtb8qnZumQNepdql0/vfLYV6yVDp1fCrz9vA2tKjtbZfdFhEpCgUUkdKkQKVTw1Lp9LB52qT0HKYs38lGS6Vzaws/XrqrJV4VVemISMlQQBEpDXIy4JupELHQHNftYql0/G2eOupoIuMWR3HivFnpPNW/Kfd1UKUjIiVLAUXE0Z3Zb6l0YjArnUeh6+N2qXQ+2nyIl9b/QW6+Qd3qZqUTcp0qHREpeQooIo5s1zL4KuxCpXPXB9DgZpunPZ+ezZRlu9i0x6x0+rXwZ9ZdLVTpiIjDUEARcUQ5GbD+cYhcZI7r3Wh+0J8dKp3Io4mM/1ul8/Rtzbi3fR1VOiLiUBRQRBzN6X1mpZMQCzhB18fMSqeCs03TGobBhz8f4uVvzEqnXvXKzFOlIyIOSgFFxJHs+gK+DIOcNPCoaZ41adDd5mkT07KZsmwn3/2RAED/lv7MurMFnqp0RMRBKaCIOILsdFj/GER9Yo7r3Wiu0vH0s3nqiCOJjF8cycmkTNxcKvBM/2bco0pHRBycAopISTu9D5aNgITdmJXO42atY2Olk59v8MHPB3nl273k5hsE1fBg3vBQmgeo0hERx6eAIlKSdn4OXz1iqXRqmat06nezedrEtGwmL9vJ95ZK57ZWAcy8I0SVjoiUGgooIiWhcKUTdBPc+SF4+to8dcSRc4xbHEWcpdKZfltzht0QqEpHREqVCkXZeNasWVx//fV4enpSq1YtBg4cyN69ewtsYxgG06dPJyAggEqVKtGtWzdiY2MLbJOVlcX48eOpUaMGHh4eDBgwgOPHj9v+akRKg9N74cMelnDiBN2mwX2rbQ4n+fkG74b/yZD3fiMuKZP6NTxYPaYzw3W9iYiUQkUKKOHh4YwdO5bffvuNjRs3kpubS+/evUlLS7NuM3v2bObMmcO8efPYtm0bfn5+9OrVi5SUFOs2YWFhrFq1iqVLl7J582ZSU1Pp378/eXl59ntlIo5oxxJ4v5t5vYlHLfi/NdBtqs3Xm5xLy+bBRdt4af0f5OUbDGgVwNrxXWgW4GWf/RYRucacDMMwrvabT58+Ta1atQgPD+emm27CMAwCAgIICwvj8ccfB8yzJb6+vrz88suMGjWKpKQkatasySeffMLQoUMBOHnyJIGBgaxbt44+ffpc9ucmJyfj7e1NUlISXl76AyylQHY6rHsUdnxqjoO6mkuI7VDpbDt8jvGLo4hPzsTdpQLTBzTn7utV6YiI4ynK+3eRzqAUlpSUBICPjw8Ahw4dIj4+nt69e1u3cXd3p2vXrmzZsgWAiIgIcnJyCmwTEBBASEiIdZvCsrKySE5OLvAlUmok/AEf3GyGE6cK0P1JuG+VXSqdd348wN3v/0Z8sqXSGduZYTeo0hGR0u+qL5I1DINJkybRpUsXQkJCAIiPjwfA17fgH15fX1+OHDli3cbNzY1q1apdtM1f31/YrFmzmDFjxtXuqkjJ2bEEvp4EOelQxde8t0nQTTZPezY1i0lf7CR832kABrYO4IU7WlDFXde9i0jZcNV/zcaNG8euXbvYvHnzRc8V/tebYRiX/Rfdv20zbdo0Jk2aZB0nJycTGBh4FXstco1kp1kqnc/Mcf1uZqVTpZbNU289dI4JSy5UOs/d3pwh7VTpiEjZclUBZfz48axdu5affvqJ2rVrWx/38zPvehkfH4+//4UPNUtISLCeVfHz8yM7O5vExMQCZ1ESEhLo1KnTP/48d3d33N3dr2ZXRa69hD/MG6+d/sOsdLpNgxsn2+XGa/PD/2TOxn3k5Rs0qOnB2/e0oYmfrsMSkbKnSNegGIbBuHHjWLlyJd9//z1BQUEFng8KCsLPz4+NGzdaH8vOziY8PNwaPtq2bYurq2uBbeLi4oiJiblkQBEpNaI+gw+6m+Gkii/831q73BX2bGoW9y/cxivf7iUv3+CO0OtYO66LwomIlFlFOoMyduxYFi9ezJo1a/D09LReM+Lt7U2lSpVwcnIiLCyMmTNnEhwcTHBwMDNnzqRy5coMHz7cuu3IkSOZPHky1atXx8fHhylTptCiRQt69uxp/1coci1kp8HXU2DnYnNcvzvc+b5dKp3fD55lwtIoTiVnUdG1As8NCGFwu9qqdESkTCtSQJk/fz4A3bp1K/D4ggULuP/++wF47LHHyMjIYMyYMSQmJtK+fXs2bNiAp6endfu5c+fi4uLCkCFDyMjIoEePHixcuBBnZ9v+lSlSIhL2wBcj4MxeyyqdJ6DLZKhg0yI56yqdORv3kW9Ag5oevHNPWxr7eV7+m0VESjmb7oNSUnQfFHEYUZ/B15MhNwOq+MGgj6BeF5unPZOaxSOf7+Dn/WcAuLPNdTx/ewgeWqUjIqVYUd6/9ddO5Gpkp5nBZOcSc1y/u2WVTk2bp/7t4FkmLIkiIcWsdJ6/PYTB7bRqTUTKFwUUkaI6tRuW3f+3SudJ6DLJ5konL9/gnR8OMHeTWek0rFWFd+5pQyNfVToiUv4ooIhcKcOAqE/N+5vkZoCnP9z1EdTrbPPUp1PMSmfzAbPSGdS2Ns/d3pzKbvpPVETKJ/31E7kSWalmpbNrqTlu0MNcpeNRw+apt/x5holLd3A6JYtKrs48PzCEQW1rX/4bRUTKMAUUkcs5FWupdPaZlc7NT0HnR+xS6cz7/gBvfGdWOsGWSidYlY6IiAKKyCUZBkR9Yql0Ms1KZ9D/oK7tNxQ8nZJF2OdR/HLgLACD29ZmhiodEREr/TUU+SdZqeaH/O363Bzbs9I5cIaJn1+odF4YGMJdqnRERApQQBEp7FSseeO1s/vBydlS6YTZpdJ56/v9vPHdfgwDGvmalU7DWqp0REQKU0AR+YthQOTHsP4xS6UTYKl0Oto8dUJKJmFLd7DlT7PSGdoukOkDmlPJTXdPFhH5JwooIgBZKfDVIxC9zBw37AV3vAce1W2e+pcD5iqdM6lZVHZz5sU7QrgjVJWOiMi/UUARiY+BZSPg7AGz0unxNHSaaJdK543v9vPW92al09jXk7fvaUPDWlXstOMiImWXAoqUX4YBkYtg/eMXKp3BC6BOB5unTkjOZOLSHfx60Kx07r4+kGdvU6UjInKlFFCkfCpc6QT3hoHv2qXS2bz/DGGfR3EmNZvKbs7MvKMFA0Ovs3leEZHyRAFFyp/4aPPGa9ZK5xnoNME+lc6mfbz1wwEMA5r4mZVOg5qqdEREikoBRcoPw4CIhWalk5cFXteZq3TsUOmcSs5kwpIofj90DoBhN9Th2duaUdFVlY6IyNVQQJHyITMZvgqDmBXmOLgP3PEuVPaxeeqf9p3mkc93cDYtGw83Z2be2YLbW6vSERGxhQKKlH1xu8xK59yfZqXTczp0HGdzpZObl8/rm/bz9o9mpdPU34u3h4dSX5WOiIjNFFCk7DIMiFgA66daKp3alkqnvc1Tn0rOZPySKLZaKp3h7evwTH9VOiIi9qKAImVTZjJ8ORFiV5pjO1Y64ftOM+lvlc6su1oyoFWAzfOKiMgFCihS9sTttFQ6B6GCC/R41m6VztxN+3j7hz8Bs9J55542BNXwsMNOi4jI3ymgSNlhGLD9I/jmiQuVzuAFEHiDzVPHJ5mrdLYeNiudezvU4al+qnRERIqLAoqUDZnJ8OUEiF1ljhv1hYHv2KXS+XFvApO+2Mm5tGyquLvw0l0t6N9SlY6ISHFSQJHSL24nfDECEg+ZlU7PGdBxLDg52TRtbl4+r23cx/wfzUqneYAXbw9vQz1VOiIixU4BRUovw4BtH8K3T0BeNngHwqAFEHi9zVPHJWUwYUkU2w4nAnBfh7o82a+pKh0RkWtEAUVKp8wkWDsBdq82x41vhdvftkul88PeBCZ9voPE9ByquLvw8l0t6dfS3+Z5RUTkyimgSOlzcoe5SuevSqfXc9BhjM2VTk5ePq9t2Me74WalE3KdF/OGqdIRESkJCihSelxU6dQxV+nUbmfz1CfPZzB+SRQRR8xKZ0THujzRrynuLqp0RERKggKKlA6ZSbB2POxeY44b94OBb0OlajZP/f0fp5j0xU7Op+fg6e7Cy4NacmsLVToiIiVJAUUc38koS6VzGCq4Wiqd/9ql0nn1272899NBwKx03h7ehrrVVemIiJQ0BRRxXIYBWz+ADU+alU7VOjBoIdRua/PUJ85nMH5xJJFHzwNwf6d6TLu1iSodEREHoYAijikzCdaMgz1rzXGT/nD7PLtUOt/tOcXkZRcqndmDWtJXlY6IiENRQBHHcyISlj9wodLp/Ty0H22XSueVb/fyvqXSaVnbm3nD2lCnemU77LSIiNiTAoo4DsOAre/Dt09Cfo5Z6QxeCNfZXukcT0xn/JIoolTpiIiUCgoo4hgyzsPacbDnS3PcpL9547VKVW2eetNus9JJysjBs6ILrwxqyS0hqnRERByZAoqUvBMRsOwBOH/EUum8AO1H2VzpZOfmM/ubP/hw8yEAWtX2Zt7wNgT6qNIREXF0CihScgwDfn8XNjxtqXTqmjdes1OlM25xFDuOnQfgwc5BTO3bBDeXCjbPLSIixU8BRUpGRqK5SuePr8xx09tgwDy7VDobYuOZsmwnyZm5eFV04ZXBrejT3M/meUVE5NpRQJFr73gELL8fzh81K50+L8IND9ul0nn5mz/46K9KJ7Aq84aFqtIRESmFFFDk2vnHSmchXNfG5qmPnUtn3JIodloqnZFdgnj8FlU6IiKllQKKXBsXVToDYMBbdql0vo2N59G/VTqvDm5Fb1U6IiKlmgKKFL/jEeZn6SQdBWc36DMTrn/ILpXOrPV7WPDLYQBaB1blLVU6IiJlggKKFB/DgN/egY3PQH4uVKtnVjoBoTZPfexcOuMWR7LzeBIA/7kxiEf7qNIRESkrFFCkeGQkwuoxsHedOW52u1npVPS2eepvYuJ5dPlOUjJz8a7kyquDW9Grma/N84qIiONQQBH7O77dvPGanSudrNw8Zq37g4VbDgMQWsesdGpXU6UjIlLWKKCI/RgG/Po2bHrWUukEWSqd1jZPffRsOuOWRLLLUuk8fFN9Hu3TGFdnVToiImVRkf+6//TTT9x2220EBATg5OTE6tWrCzxvGAbTp08nICCASpUq0a1bN2JjYwtsk5WVxfjx46lRowYeHh4MGDCA48eP2/RCpISln4Olw2HDk2Y4aTYQRoXbJZx8ExNHv7d+ZtfxJKpWduWjEe144tamCiciImVYkf/Cp6Wl0apVK+bNm/ePz8+ePZs5c+Ywb948tm3bhp+fH7169SIlJcW6TVhYGKtWrWLp0qVs3ryZ1NRU+vfvT15e3tW/Eik5x7bBezeZ15s4u0G/18wzJzZeb5KVm8f0tbGM/jSSlMxc2tSpytcTbqRHU11vIiJS1jkZhmFc9Tc7ObFq1SoGDhwImGdPAgICCAsL4/HHHwfMsyW+vr68/PLLjBo1iqSkJGrWrMknn3zC0KFDATh58iSBgYGsW7eOPn36XPbnJicn4+3tTVJSEl5eXle7+2Irw4Bf58Gm6eZZE5/6ZjDxb2Xz1EfPpjN2cSTRJ8xKZ1TX+kzprUpHRKQ0K8r7t13/2h86dIj4+Hh69+5tfczd3Z2uXbuyZcsWACIiIsjJySmwTUBAACEhIdZtCsvKyiI5ObnAl5Sw9HOwZBhseMoMJ83vgIfD7RJO1kXH0e/Nn4k+YVY6/7u/HdP6qtIRESlP7HqRbHx8PAC+vgVPwfv6+nLkyBHrNm5ublSrVu2ibf76/sJmzZrFjBkz7LmrYotjW2H5g5B0DJzd4ZZZ0O5Bu6zSmfn1Hhb9av6utK1bjbeGhRJQtZI99lpEREqRYvknqVOhNyrDMC56rLB/22batGkkJSVZv44dO2a3fZUiyM+HX96EBX3NcOJTHx7aCNePtDmcHDmbxl3zt1jDyeiuDVj6cAeFExGRcsquZ1D8/MzPP4mPj8ff39/6eEJCgvWsip+fH9nZ2SQmJhY4i5KQkECnTp3+cV53d3fc3d3tuatSVOnnYPV/Yd835rj5nXDbG1DR9muAvt4Vx9QVu0jJyqVaZVfmDGlN9ya1bJ5XRERKL7ueQQkKCsLPz4+NGzdaH8vOziY8PNwaPtq2bYurq2uBbeLi4oiJiblkQJESdvR3ePdGM5w4u0P/uTDofzaHk8ycPJ5eHcPYxZGkZOXSrm411k28UeFERESKfgYlNTWVAwcOWMeHDh1ix44d+Pj4UKdOHcLCwpg5cybBwcEEBwczc+ZMKleuzPDhwwHw9vZm5MiRTJ48merVq+Pj48OUKVNo0aIFPXv2tN8rE9vl58Ovb8GmGWDkgU8DyyqdljZPffhMGmMXRxJ70rzgeUy3Bkzq1QgXXQgrIiJcRUDZvn073bt3t44nTZoEwIgRI1i4cCGPPfYYGRkZjBkzhsTERNq3b8+GDRvw9PS0fs/cuXNxcXFhyJAhZGRk0KNHDxYuXIizs7MdXpLYRfo5WDUa9n9rjkPuMisdd89//74r8OXOk0xbGU1qVi4+Hm7MGdKKbo111kRERC6w6T4oJUX3QSlmR3+H5Q9A8gmz0un7MrS93+YLYTNz8njh6918+ttRAG6o58Obw0Lx865oh50WERFHV5T3b30Wj1yQnw9b3oTvnjMrneoNzUrHr4XNUx86k8bYzyLZHWdWOmO7N+CRnqp0RETknymgiCntLKwaBQcsFy+HDILbXrdLpbN250mmrdhFWnYePh5uzB3amq6Nato8r4iIlF0KKAJHfjVvvJZyElwqmpVOmxF2qXSe+2o3i3+3VDpBPrx5tyodERG5PAWU8iw/H355Hb5/4W+VziLwC7F56oOnUxm7OIo9cck4OcG47g2Z2CNYlY6IiFwRBZTyKu0srHoYDmwyxy0Gm/c3sUOls2bHCZ5YGU1adh7VLZXOTap0RESkCBRQyqOLKp3Z0Ob/7FLpzPhyN0u2mpVO+yBzlY6vlyodEREpGgWU8uSiSifYskrH9krnz9OpjP0skj/iU3BygvHdGzJBlY6IiFwlBZTyIu2MZZWOpdJpORT6zQH3KjZPvTrqBE+siiY9O48aVdx4fWgoXYJr2DyviIiUXwoo5cGRLZZKJ86sdG59FULvtUulM31tLEu3mZ8u3aG+uUqnliodERGxkQJKWZafD5vnwA8vgpEPNRqZlY5vc5unPpCQyrjFFyqdCTcHM6FHMM4VbAs9IiIioIBSdqWeNlfp/Pm9ObZjpbMq6jhProqxVDruvHF3azo3VKUjIiL2o4BSFh3+BVaMtFQ6leDWV+xS6WRkm5XO59vNSqdTg+q8fndranmq0hEREftSQClL8vNh82vww8y/VTqLwLeZzVMfSEhh7GdR7D1lVjoTewQz/mZVOiIiUjwUUMqKwpVOq2HmxbB2qHRWRBznqdUxZOSYlc6bd7emkyodEREpRgooZcHhzbB8JKTGm5VOv9cg9B6bp83IzuOZNTEsizgOQOeG1Zk7VJWOiIgUPwWU0iw/D36eAz/+Vek0hiGLoFZTm6fefyqFsYsj2XcqlQpOENazEWO7N1SlIyIi14QCSmmVmgAr/wMHfzTHrYZDv1fBzcPmqZdHHOdpS6VT09NcpdOpgSodERG5dhRQSqNDP5urdFJP2bXSSc/O5Zk1sSy3VDpdGtZg7tDW1PR0t3luERGRolBAKU3y8+Dn1+DHWWalU7OJeeM1O1Q6+06lMPazSPYnmJXOIz0bMUaVjoiIlBAFlNKicKXT+l64dbZdKp1l24/x9JoYMnPyqeXpzpvDQulQv7rN84qIiFwtBZTS4NBPsOIhs9JxrWxWOq2H2zxtenYuT62OYWXkCQBuDDYrnRpVVOmIiEjJUkBxZPl58NOrEP6SpdJpaql0mtg89d54c5XOAUulM7l3Y/7btQEVVOmIiIgDUEBxVCmnYOVD5tkTMG9V3/cVcKts07SGYbBs+3GeWWtWOr5e7rx5dyjtVemIiIgDUUBxRAfDzUonLcFS6cyB1sNsnjYtK5enV8ewMsqsdG5qVJO5Q1pRXZWOiIg4GAUUR5KfBz+9Aj++BBhmpTNkEdRsbPPUf8QnM/azSP48naZKR0REHJ4CiqO4qNK5D/rOtkul88X2YzyzJpas3Hz8vCry5rBQbgjyscNOi4iIFA8FFEdw8EdY8R9LpeMB/edCq6E2T5uWlcuTq6JZveMkAF0b1WSOKh0RESkFFFBKUn4ehL8M4bMBA2o1g8GLoGYjm6feE5fM2MWRHDydhnMFJ6b0bsyom+qr0hERkVJBAaWkpMSbF8Ie/tkct/k/uOVlu1Q6S7cdY/raC5XOW8NDub6eKh0RESk9FFBKwp8/mHeFTTtt10on1VLprLFUOt0b1+S1Ia3x8XCzeW4REZFrSQHlWrqo0mlu3njNDpXO7pPJjFscycEzZqXzaJ/GPHyjKh0RESmdFFCulYsqnRHQ92VwrWTTtIZhsGTrMaZ/GUt2bj7+3hV5a1go7VTpiIhIKaaAci38+T2sfNisdNyqQP/XoeVgm6dNyczhiVUxfLnTrHRublKL1wa3opoqHRERKeUUUIpTXq75OTo/vQoY4BtiVjo1gm2eOvZkEuMWR3HIUuk8fktjHuqiSkdERMoGBZTikhxnVjpHNpvjtvfDLS/ZpdL57PejPPfVbrJz8wnwNlfptK2rSkdERMoOBZTicOA7s9JJP2NWOre9AS0G2TxtSmYO01ZG89WuOAB6NKnFq6p0RESkDFJAsae8XPhxFvz8Gmal08JS6TS0eeqYE0mMWxzJ4bPpuFRw4vFbmvDQjUE4OanSERGRskcBxV6S42DFSDjyizlu+wDcMssulc6nvx/l+S93k52Xz3VVK/HW8FDa1Klmh50WERFxTAoo9nBgE6wcZfdKJzkzh2krovk62qx0ejY1K52qlVXpiIhI2aaAYou8XPhxpqXSwax0hiyC6g1snjrmRBJjF0dyxFLpTO3bhJFdVOmIiEj5oIBytZJPwvKRcHSLOW43EvrMBNeKNk1rGAaf/HaEF77aY6105g0PJVSVjoiIlCMKKFdj/yZY9TCknwU3TxjwBoTcZfO0yZk5TF2xi3XR8QD0aubLq4Na4V3Z1ea5RUREShMFlKLIy4UfXoTNc8yxXwsYbHulYxgGUcfOE7Z0B0fPpePq7MTUvk15sHM9VToiIlIuKaBcqaQT5iqdo7+a4+sfgt4vFrnSycrNY/+pVHbHJbP7ZDK745LZE5dMSmYuALWrVWLe8Da0Dqxq5xcgIiJSeiigXIn9G80br2Wcs1Q6b0LInZf9tvPp2QWCyO6TyRxISCU337hoWzfnCvQJ8eOF20NU6YiISLlXogHlnXfe4ZVXXiEuLo7mzZvz+uuvc+ONN5bkLhWUlws/vACb55pjv5bmjdcKVTqGYXDsXAa745L+dlYkhRPnM/5xWu9KrjQP8KKZvxdN/b1oFuBFg5pVcHOpUMwvSEREpHQosYDy+eefExYWxjvvvEPnzp1577336Nu3L7t376ZOnToltVsXJB03V+kc+80cX/8f6P0CWU6u7D+RVOCsyJ64ZFKycv9xmjo+lWlmCSF//a+/d0VdWyIiIvIvnAzDuLhvuAbat29PmzZtmD9/vvWxpk2bMnDgQGbNmvWv35ucnIy3tzdJSUl4eXnZf+f2bYBVoyDjHLmuVfix8TOsy2vP7rh/r2ga+VUxQ4i/F80CvGni74lXRdU1IiIiULT37xI5g5KdnU1ERARTp04t8Hjv3r3ZsmXLRdtnZWWRlZVlHScnJxfLfqVnZBDz6WPccOJjAKLz6zEudQJHtvsBJ6zbVa18oaIxz4x4U7+mB67OqmhERETsoUQCypkzZ8jLy8PX17fA476+vsTHx1+0/axZs5gxY0ax75f7wQ3WcLIwtzczc+/Bv7o3ff3/FkYCvPDzUkUjIiJSnEr0ItnCb/KGYfzjG/+0adOYNGmSdZycnExgYKDd98e52QCi/e/ibM2ONG9zFxF+nniqohEREbnmSiSg1KhRA2dn54vOliQkJFx0VgXA3d0dd3f34t8xJydajPpf8f8cERER+VclctGEm5sbbdu2ZePGjQUe37hxI506dSqJXRIREREHUmIVz6RJk7jvvvto164dHTt25P333+fo0aOMHj26pHZJREREHESJBZShQ4dy9uxZnnvuOeLi4ggJCWHdunXUrVu3pHZJREREHESJ3QfFFsV+HxQRERGxu6K8f+vGHSIiIuJwFFBERETE4SigiIiIiMNRQBERERGHo4AiIiIiDkcBRURERByOAoqIiIg4HAUUERERcTgKKCIiIuJwSuxW97b46+a3ycnJJbwnIiIicqX+et++kpvYl8qAkpKSAkBgYGAJ74mIiIgUVUpKCt7e3v+6Tan8LJ78/HxOnjyJp6cnTk5Odp07OTmZwMBAjh07ps/5uQwdqyunY3XldKyunI7VldFxunLFfawMwyAlJYWAgAAqVPj3q0xK5RmUChUqULt27WL9GV5eXvpFvkI6VldOx+rK6VhdOR2rK6PjdOWK81hd7szJX3SRrIiIiDgcBRQRERFxOAoohbi7u/Pss8/i7u5e0rvi8HSsrpyO1ZXTsbpyOlZXRsfpyjnSsSqVF8mKiIhI2aYzKCIiIuJwFFBERETE4SigiIiIiMNRQBERERGHo4DyN++88w5BQUFUrFiRtm3b8vPPP5f0Ll1zP/30E7fddhsBAQE4OTmxevXqAs8bhsH06dMJCAigUqVKdOvWjdjY2ALbZGVlMX78eGrUqIGHhwcDBgzg+PHj1/BVFL9Zs2Zx/fXX4+npSa1atRg4cCB79+4tsI2OlWn+/Pm0bNnSeuOnjh07sn79euvzOk6XNmvWLJycnAgLC7M+puNlmj59Ok5OTgW+/Pz8rM/rOBV04sQJ7r33XqpXr07lypVp3bo1ERER1ucd8ngZYhiGYSxdutRwdXU1PvjgA2P37t3GxIkTDQ8PD+PIkSMlvWvX1Lp164wnn3zSWLFihQEYq1atKvD8Sy+9ZHh6ehorVqwwoqOjjaFDhxr+/v5GcnKydZvRo0cb1113nbFx40YjMjLS6N69u9GqVSsjNzf3Gr+a4tOnTx9jwYIFRkxMjLFjxw6jX79+Rp06dYzU1FTrNjpWprVr1xpff/21sXfvXmPv3r3GE088Ybi6uhoxMTGGYeg4XcrWrVuNevXqGS1btjQmTpxofVzHy/Tss88azZs3N+Li4qxfCQkJ1ud1nC44d+6cUbduXeP+++83fv/9d+PQoUPGpk2bjAMHDli3ccTjpYBiccMNNxijR48u8FiTJk2MqVOnltAelbzCASU/P9/w8/MzXnrpJetjmZmZhre3t/Huu+8ahmEY58+fN1xdXY2lS5datzlx4oRRoUIF45tvvrlm+36tJSQkGIARHh5uGIaO1eVUq1bN+PDDD3WcLiElJcUIDg42Nm7caHTt2tUaUHS8Lnj22WeNVq1a/eNzOk4FPf7440aXLl0u+byjHi9VPEB2djYRERH07t27wOO9e/dmy5YtJbRXjufQoUPEx8cXOE7u7u507drVepwiIiLIyckpsE1AQAAhISFl+lgmJSUB4OPjA+hYXUpeXh5Lly4lLS2Njh076jhdwtixY+nXrx89e/Ys8LiOV0H79+8nICCAoKAg7r77bg4ePAjoOBW2du1a2rVrx+DBg6lVqxahoaF88MEH1ucd9XgpoABnzpwhLy8PX1/fAo/7+voSHx9fQnvleP46Fv92nOLj43Fzc6NatWqX3KasMQyDSZMm0aVLF0JCQgAdq8Kio6OpUqUK7u7ujB49mlWrVtGsWTMdp3+wdOlSIiMjmTVr1kXP6Xhd0L59ez7++GO+/fZbPvjgA+Lj4+nUqRNnz57VcSrk4MGDzJ8/n+DgYL799ltGjx7NhAkT+PjjjwHH/b0qlZ9mXFycnJwKjA3DuOgxubrjVJaP5bhx49i1axebN2++6DkdK1Pjxo3ZsWMH58+fZ8WKFYwYMYLw8HDr8zpOpmPHjjFx4kQ2bNhAxYoVL7mdjhf07dvX+v9btGhBx44dadCgAYsWLaJDhw6AjtNf8vPzadeuHTNnzgQgNDSU2NhY5s+fz//93/9Zt3O046UzKECNGjVwdna+KAUmJCRclCjLs7+ukP+34+Tn50d2djaJiYmX3KYsGT9+PGvXruWHH36gdu3a1sd1rApyc3OjYcOGtGvXjlmzZtGqVSveeOMNHadCIiIiSEhIoG3btri4uODi4kJ4eDhvvvkmLi4u1ter43UxDw8PWrRowf79+/V7VYi/vz/NmjUr8FjTpk05evQo4Lh/rxRQMP94tm3blo0bNxZ4fOPGjXTq1KmE9srxBAUF4efnV+A4ZWdnEx4ebj1Obdu2xdXVtcA2cXFxxMTElKljaRgG48aNY+XKlXz//fcEBQUVeF7H6t8ZhkFWVpaOUyE9evQgOjqaHTt2WL/atWvHPffcw44dO6hfv76O1yVkZWWxZ88e/P399XtVSOfOnS+6DcK+ffuoW7cu4MB/r4rl0ttS6K9lxh999JGxe/duIywszPDw8DAOHz5c0rt2TaWkpBhRUVFGVFSUARhz5swxoqKirMutX3rpJcPb29tYuXKlER0dbQwbNuwfl6LVrl3b2LRpkxEZGWncfPPNZW7p3n//+1/D29vb+PHHHwssc0xPT7duo2NlmjZtmvHTTz8Zhw4dMnbt2mU88cQTRoUKFYwNGzYYhqHjdDl/X8VjGDpef5k8ebLx448/GgcPHjR+++03o3///oanp6f1b7aO0wVbt241XFxcjBdffNHYv3+/8dlnnxmVK1c2Pv30U+s2jni8FFD+5u233zbq1q1ruLm5GW3atLEuGS1PfvjhBwO46GvEiBGGYZjL0Z599lnDz8/PcHd3N2666SYjOjq6wBwZGRnGuHHjDB8fH6NSpUpG//79jaNHj5bAqyk+/3SMAGPBggXWbXSsTA8++KD1v6uaNWsaPXr0sIYTw9BxupzCAUXHy/TXfTpcXV2NgIAA48477zRiY2Otz+s4FfTll18aISEhhru7u9GkSRPj/fffL/C8Ix4vJ8MwjOI5NyMiIiJydXQNioiIiDgcBRQRERFxOAooIiIi4nAUUERERMThKKCIiIiIw1FAEREREYejgCIiIiIORwFFREREHI4CioiIiDgcBRQRERFxOAooIiIi4nAUUERERMTh/D9GjNySy3TPmQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(cumulative_desired, cumulative_desiredo)\n",
    "x = np.arange(0, 600)\n",
    "plt.plot(x, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "28d25528-bde4-4fb3-8b2b-e8507f989e14",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "602.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGdCAYAAAA1/PiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA85UlEQVR4nO3deXxU9b3/8feQZSAhGbJAwkBA0CBiAkJQFJWlbNoitbZFC0X7K+3FotgoCKXeWuxtQ6W3YntpuRevVyoudLG0tqISqiIICgRQArJJIAESwpJMEghZv78/kKNDgmSZ5Mzyej4e83icOeebk8+3BzPvfud7vsdhjDECAAAIAB3sLgAAAKCpCC4AACBgEFwAAEDAILgAAICAQXABAAABg+ACAAACBsEFAAAEDIILAAAIGOF2F9AS9fX1OnbsmGJiYuRwOOwuBwAANIExRuXl5XK73erQoWVjJwEZXI4dO6aUlBS7ywAAAC1QUFCgnj17tuhnAzK4xMTESDrf8djYWJurAQAATVFWVqaUlBTrc7wlAjK4XPh6KDY2luACAECAac00DybnAgCAgEFwAQAAAYPgAgAAAgbBBQAABAyCCwAACBgEFwAAEDAILgAAIGAQXAAAQMAguAAAgIBBcAEAAAGj2cHl6NGj+va3v62EhARFRUXpuuuuU05OjnXcGKMFCxbI7XarU6dOGjVqlHbt2uV1jqqqKs2aNUuJiYmKjo7WpEmTdOTIkdb3BgAABLVmBZeSkhLdfPPNioiI0Ouvv67du3fr17/+tbp06WK1WbRokZ566iktWbJEW7ZsUXJyssaNG6fy8nKrTWZmplatWqWVK1dqw4YNqqio0MSJE1VXV+ezjgEAgODjMMaYpjb+0Y9+pPfee0/r169v9LgxRm63W5mZmZo3b56k86MrSUlJevLJJzVjxgx5PB517dpVK1as0N133y1JOnbsmFJSUrR69WpNmDDhsnWUlZXJ5XLJ4/HwkEUAAHyovt7ofzccVE2d0QOjr/LpuX3x+d2sEZdXX31VQ4cO1Te/+U1169ZNgwcP1jPPPGMdz8vLU1FRkcaPH2/tczqdGjlypDZu3ChJysnJUU1NjVcbt9uttLQ0qw0AALDHiF+9razVe7Ri02G7S2lUs4LLwYMHtXTpUqWmpurNN9/U/fffr4ceekjPP/+8JKmoqEiSlJSU5PVzSUlJ1rGioiJFRkYqLi7ukm0uVlVVpbKyMq8XAADwnYMnKnTFj17TkZJKSVLJ2WqbK2pceHMa19fXa+jQocrKypIkDR48WLt27dLSpUt17733Wu0cDofXzxljGuy72Be1WbhwoZ544onmlAoAAJrAGKNXth3VnD9/6LV/z3/cZlNFX6xZIy7du3fXgAEDvPZdc801ys/PlyQlJydLUoORk+LiYmsUJjk5WdXV1SopKblkm4vNnz9fHo/HehUUFDSnbAAAcJHK6jr9JeeIBi5Y4xVaZo66Uvt/cftlBxzs0qzgcvPNN2vv3r1e+/bt26fevXtLkvr06aPk5GRlZ2dbx6urq7Vu3ToNHz5ckpSRkaGIiAivNoWFhcrNzbXaXMzpdCo2NtbrBQAAmu9cTZ3m/uVDXfP4G5rz5w9VXlVrHXvpe8M097b+igjz32XemvVV0cMPP6zhw4crKytLkydP1ubNm7Vs2TItW7ZM0vmviDIzM5WVlaXU1FSlpqYqKytLUVFRmjJliiTJ5XJp+vTpmj17thISEhQfH685c+YoPT1dY8eO9X0PAQAIcfX1Rv/cWajZf9qhmrqGNxP/x51p+saQnuoUGWZDdc3TrOBy/fXXa9WqVZo/f75+9rOfqU+fPnr66ac1depUq83cuXNVWVmpmTNnqqSkRMOGDdOaNWsUExNjtVm8eLHCw8M1efJkVVZWasyYMVq+fLnCwvz/fzAAAALFx4VlWvj6Hr2770SDY1ckROn57w5Tr4QoGypruWat4+IvWMcFAIDGlZ+r0e5jZXr0Lx8p//TZBsd//OX+mjw0RV2iItu9Nl98fjdrxAUAAPiv13cW6gcvbmuw/9bURM27rb8GdI9Vhw7+Oem2qQguAAAEgfp643V3UFKsUzdflajHJw6wZXSlrRBcAAAIAgdOVOhM9fln/v1z1i1K6+GyuaK24b/3OwEAgCbbdvj8+mg39U0I2tAiMeICAEBAq66t1+N/z9X6/SclSYN7dbG3oDZGcAEAwAb19UZna+pafZ639hRr5ZbPVpS/+arEVp/TnxFcAABoZzmHT+vrSzf59Jwj+nXVjBF9NfzKBJ+e198QXAAAaEc/+8du/d97eT49Z8eIDrp/ZF8NvzK4R1skggsAAO1i4ycndf+KHJWd++zZQEunDtHo/t1afe7wDg6F+/HzhXyJ4AIAQBuqrzd68OVtWr2zyGv/h4+PlysqwqaqAhfBBQCANlJbV69ZL2/X67mfhZbZ4/rpu7f0UbSTj+CW4H81AAB8zBijH6/aqZc3F3jt3/aTcYqPDp5VbO1AcAEAwIfOVtdqwONveu3rmxit1T+8VR0jwmyqKngQXAAAaKXisnPafOi0fvf2J/q4sMzr2L9mj9QVCdEKC/CHG/oLggsAAC1QerZa6/ad0B+3FGjjJ6caHE/v4dKrD94sh4PA4ksEFwAAmsgYo9dzi/Rfbx1oMLIiSVckROn29O4a1a+rMnrHEVraAMEFAIAvUFdvtGr7US3O3qejpZWNtnl0wtWaNMitlPiodq4u9BBcAABoRHVtvSb+13rtO17R6PFHJ1ytqcN6qUsUdwm1J4ILAACfk3vUo4Wvf6z3DjSct/JvI/pqxoi+io+O5GsgmxBcAAD4nGXvHvQKLSnxnfTcd67XlV07E1b8AMEFAIDP2V5QIun8V0GDU7rohj7xIfMcoEBAcAEA4FP7j5er4HSlHA5p2k29FduRZwn5G4ILACDgfXKiQp8UNz6Jtjne3X9CkuQM70Bo8VMEFwBAQHvpg3z9eNVOn57zgVFX+fR88B2CCwAgIO0tKtcz6w/qLzlHrH1DenVp9Xm7REVq8vUprT4P2gbBBQAQUOrrjWat3K7XPir02v9G5q3qnxxrU1VoLwQXAEDAOHiiQt9dvkWHTp219n19SE/dNaQHoSVEEFwAAH6vrt7o8b/n6sUP8r32b/33sUrs7LSpKtiB4AIA8GvF5ed0wy/+5bXvtmuTtfCudMVFs9x+qCG4AAD8UlVtnZ58fa/+7708r/2b5n9J3V2dbKoKdiO4AAD8ymsfFerlzfnacOCk1/6vDOyu39x9HavYhjiCCwDAVlW1dfrXx8X6n3cP6sOC0gbHO0WEac3DI5QSH9X+xcHvEFwAAO2uoqpWHxaU6ldv7tWORsKKJE2/pY++e0sf9ejC10L4DMEFANBuDhRXaOHqj/WvPcWNHr9jkFsPj03VFQnR6tCBJzGjIYILAKBNfXSkVK/kHNG+4xXadPCU17GuMU7dcEW8/m1EX6X3cBFWcFkEFwBAm/rJ33L14RGP176vDe6haTf11pBecTZVhUBFcAEAtJmz1bXKPVYmSfrBqCvV2Rmuu4b04HZmtBjBBQDQZj464lFdvVFybEfNu62/3eUgCBBcAAANlJ6tVkVVbavPs37/CUnSYB88tRmQCC4AgM8xxujFD/L173/L9el5mcsCXyG4AAAkSbuPlenLv13vta9jROtXqU3s7NRtacmtPg8gEVwAIOQZY/TInz7Uqu1Hvfb/c9YtSuvhsqkqoHEEFwAIYbuOefSV327w2vfEpGt1zw0pcoaH2VQVcGkEFwAIQcYYPfGP3Vq+8ZDX/u0/Gae46Eh7igKagOACACHovue26N19J6z3M0b01SPj+zHKAr/XrFlXCxYskMPh8HolJ3824coYowULFsjtdqtTp04aNWqUdu3a5XWOqqoqzZo1S4mJiYqOjtakSZN05MgR3/QGAHBJxhj96JWPdMWPXvMKLR/8eIzmf/kaQgsCQrOni1977bUqLCy0Xjt37rSOLVq0SE899ZSWLFmiLVu2KDk5WePGjVN5ebnVJjMzU6tWrdLKlSu1YcMGVVRUaOLEiaqrq/NNjwAAXs7V1Gnt7uPqM3+1Vm4psPZf645V3sIvKym2o43VAc3T7K+KwsPDvUZZLjDG6Omnn9Zjjz2mu+66S5L0hz/8QUlJSXrppZc0Y8YMeTwePfvss1qxYoXGjh0rSXrhhReUkpKitWvXasKECa3sDgDggoLTZzX7zx9qc97pBsfezByhq7p1lsPBQw0RWJo94rJ//3653W716dNH99xzjw4ePChJysvLU1FRkcaPH2+1dTqdGjlypDZu3ChJysnJUU1NjVcbt9uttLQ0q01jqqqqVFZW5vUCADS073i5Hv3zh5qxYqtuXfS2V2iJCHNoxoi++iTry7o6OUZhPIkZAahZIy7Dhg3T888/r379+un48eP6+c9/ruHDh2vXrl0qKiqSJCUlJXn9TFJSkg4fPixJKioqUmRkpOLi4hq0ufDzjVm4cKGeeOKJ5pQKACFje36J3tpTrDd3FWnf8YoGx7+Z0VPfyOipG/rEM8KCgNes4HL77bdb2+np6brpppt05ZVX6g9/+INuvPFGSWrwH4Ux5rL/oVyuzfz58/XII49Y78vKypSSktKc0gEgqKzff0KHT53Vk6/vUXkjzxT66nVupfdwaWS/rkpNirGhQqBttOp26OjoaKWnp2v//v268847JZ0fVenevbvVpri42BqFSU5OVnV1tUpKSrxGXYqLizV8+PBL/h6n0ymn09maUgEgKNTXG92z7H1tPtRw3sq3builrjFO/duIvursZLULBKdWPYSiqqpKH3/8sbp3764+ffooOTlZ2dnZ1vHq6mqtW7fOCiUZGRmKiIjwalNYWKjc3NwvDC4AgPP2FZd7hZbxA5L08Nh+2vvz27TwrnQ9Mq4foQVBrVn/uufMmaM77rhDvXr1UnFxsX7+85+rrKxM9913nxwOhzIzM5WVlaXU1FSlpqYqKytLUVFRmjJliiTJ5XJp+vTpmj17thISEhQfH685c+YoPT3dussIANBQTV29fvuv/co5XCJJiukYrq3/Ppa1VxBymhVcjhw5om9961s6efKkunbtqhtvvFHvv/++evfuLUmaO3euKisrNXPmTJWUlGjYsGFas2aNYmI++3518eLFCg8P1+TJk1VZWakxY8Zo+fLlCgvjPz4AuJQ3cov0X28dsN5/75a+hBaEJIcxxthdRHOVlZXJ5XLJ4/EoNjbW7nIAoM0teHWXlm88pBv7xuv2tO76ekZPvhJCwPHF5zf/6gGgjXxYUKqH/7RDkWGtmk4o6fxictL5Cbhfva5Hq88HBCqCCwC0gafW7NVvP/fVji84wzvopr4JPj0nEGgILgDgQ3uKyrT0nU/09x3HrH0PjUnVsD7xrT53r/godeO5QghxBBcA8AFjjB7/+y6teP+w1/4N80arZ1yUTVUBwYfgAgCt9MmJCt31+43yVNZY++4emqJpN/UmtAA+RnABgBaqqzf65esf65n1eV77188drZR4AgvQFgguANAC5edqlL5gjde+rwzsrt/eM5inLgNtiOACAM2wo6BUz27I0z8+POa1n1EWoH0QXADgMtbtO6F9ReX6zzV7VVVb73VsZL+uWv7/rv/CJ9wD8B2CCwBcpLauXhs/OaU9RWXKWr2n0TbfvbmPpt7YS1d27dzO1QGhjeACAJKqauu0p7Bcq7Yf1fKNhxpt87XBPdQ/OUbfvaWPInywGi6A5iO4AAhZdfVGhZ5KLXnrgFZuKWi0zQ1XxOvrGT30tcE9FRlOWAHsRnABEJK2Hjqtb/z3pgb7I8M7qFuMU1lfS9fwKxMUzsgK4FcILgBCjjFGP3hxm/Xe4ZCudcdq3m39dWtqVxsrA3A5BBcAIedoaaVOlFdJkp69b6jGXJNkc0UAmooxUAAhZ1t+qSRpYE8XoQUIMIy4AAgYWw+d1subC2SMadV5dheWSZIGp3TxQVUA2hPBBUBAyFr9sZa9e9Cn5xzWN8Gn5wPQ9gguAPza9vwS3bV0oz4/yDL9lj5Kju3YqvMmdI7Ubdcmt7I6AO2N4ALAL9XXGz32t1y9vDnfa/+Hj4+XKyrCpqoA2I3gAsDvfFhQqq/+7j2vff/x1Wt19/W9WAQOCHEEFwB+wxijea98pD9tPeK1f8tjY9U1xmlTVQD8CcEFgF+orze67Tfvat/xCmvfg6Ov0g/HpvJcIAAWggsAW9XXGz3ypx36245j1r6OER20Yd6XlNiZURYA3gguANpdbV29dhSU6pVtRxtMvr3+ijj9+f7hNlUGwN8RXAC0i9q6eh0+fVa/Wbtfr354rNE2qx+6Vf2SOrdzZQACCcEFQJuprzfaX1yhlzfna/nGQ422uTopRj+dNEDD+iQorIOjfQsEEHAILgB8bm9RuVZuyddz7x1qcCwyrIOuccdq1uir9KX+3dSBsAKgGQguAHzm+U2HtDh7n0rO1jQ4dud1bt01pKdG9OtqQ2UAggXBBYBP1NcbPf73XV77Jg1ya+qwXhrSO45bmgH4BMEFgE8cOPHZ+isvfW+Y0nu6FNORpfkB+BbBBYBPbM8vkSQN6xOv4Vcl2lwNgGDF2C0An9h2uFSSNKR3nL2FAAhqBBcAPrG94PyIy5BeBBcAbYfgAqDVys7VaH/x+Tkug3t1sbcYAEGN4AKg1T4sKJUxUq/4KJ4vBKBNMTkXCEHb80v0xq4iyfjmfDuPeiRJQxhtAdDGCC5AiNl/vFxf+/3GNjn30Cvi2+S8AHABwQUIIbV19Rq3+F3r/bQbe6tTZJhPzt0lKkLfyOjpk3MBwKUQXIAQMvG/Nljbv7wrXffc0MvGagCg+ZicC4SIUxVV2lNULknqnxxDaAEQkAguQAg4faZa0/+w1Xr/j1m32FgNALQcXxUBQa6mrl5D/iPbep/BAw8BBDD+egFB7o7PzWsZ3KuLFk++zr5iAKCVGHEBgtgft+Rb81p6xUdp1cybba4IAFqnVSMuCxculMPhUGZmprXPGKMFCxbI7XarU6dOGjVqlHbt2uX1c1VVVZo1a5YSExMVHR2tSZMm6ciRI60pBUAjFry629pe+8hIGysBAN9ocXDZsmWLli1bpoEDB3rtX7RokZ566iktWbJEW7ZsUXJyssaNG6fy8nKrTWZmplatWqWVK1dqw4YNqqio0MSJE1VXV9fyngCwGGN08y/fUmXN+f+mnvvO9YoM55thAIGvRX/JKioqNHXqVD3zzDOKi/vsSbDGGD399NN67LHHdNdddyktLU1/+MMfdPbsWb300kuSJI/Ho2effVa//vWvNXbsWA0ePFgvvPCCdu7cqbVr1/qmV0CIm/bsZh0trZQk9U6I0q2piTZXBAC+0aLg8sADD+grX/mKxo4d67U/Ly9PRUVFGj9+vLXP6XRq5MiR2rjx/BLjOTk5qqmp8WrjdruVlpZmtblYVVWVysrKvF4AGvenrQXacOCkJCm8g0PvzBmlcO4iAhAkmj05d+XKldq2bZu2bNnS4FhRUZEkKSkpyWt/UlKSDh8+bLWJjIz0Gqm50ObCz19s4cKFeuKJJ5pbKhByisvPae5fPrLef/wft8nhcNhYEQD4VrP+b1hBQYF++MMf6oUXXlDHjh0v2e7iP5TGmMv+8fyiNvPnz5fH47FeBQUFzSkbCAmb807rhl/8y3r/ZuYI1msBEHSaNeKSk5Oj4uJiZWRkWPvq6ur07rvvasmSJdq7d6+k86Mq3bt3t9oUFxdbozDJycmqrq5WSUmJ16hLcXGxhg8f3ujvdTqdcjqdzSkVCBnVtfWasWKr3t57wtp353VuXZ0cY2NVANA2mvV/x8aMGaOdO3dqx44d1mvo0KGaOnWqduzYob59+yo5OVnZ2Z+t0lldXa1169ZZoSQjI0MRERFebQoLC5Wbm3vJ4AKgcQeKK9Tv31/3Ci0Pj+2nn92ZZmNVANB2mjXiEhMTo7Q07z+I0dHRSkhIsPZnZmYqKytLqampSk1NVVZWlqKiojRlyhRJksvl0vTp0zV79mwlJCQoPj5ec+bMUXp6eoPJvgAaV1Vbp/mv7NRftx+19sVHR+rt2aPkioqwsTIAaFs+Xzl37ty5qqys1MyZM1VSUqJhw4ZpzZo1ion5bNh68eLFCg8P1+TJk1VZWakxY8Zo+fLlCgsL83U5QNDJOVyiry/1vgPvwdFXac6Eq22qCADaj8MYY+wuornKysrkcrnk8XgUGxtrdzlAuzDG6MerdurlzZ9NTo+KDNM7j45St5hLT5YHAH/hi89vnlUEBICdRzy6Y8kGr30/++q1uvemK+wpCABsQnAB/Fh1bb02HDih7y7fau2LcYYr+5GRSnYxygIg9BBcAD91+NQZjfzVO177Fn19oCZfn2JPQQDgB1idCvBTE//L+6uhX3wtjdACIOQx4gL4oZIz1So/VytJmjGyr+ZN6K8OHVi6HwAILoCfyT3q0V9yjkiS+naN1vzbr7G5IgDwHwQXwI9U1dbpW8veV3nV+dGWjF5xl/kJAAgtBBfAj+w6VqbyqlpFR4bp9vTumjn6KrtLAgC/QnAB/Mj2/FJJ0o19E/Sf3xxkbzEA4Ie4qwjwI9vySyRJQ3rzFREANIYRF8AHjpVWat2+E6pv5RM0NuedliQNTunig6oAIPgQXIBWqK2r16/W7NX/rDvos3OGdXBoIMEFABpFcAFaqPxcjW558m15KmusfUN6dVFiZ2erzju6fzd1dvKfJgA0hr+OQAv87/qD+vlrH3vty354hFKTYmyqCABCA8EFaIZjpZUaseht1dZ/NpdlcK8u+sv9wxXGyrYA0OYILkATGGP0zr4T+n/PbfHazygLALQvggtwGScrqjT052u99o0bkKRl0zLkcDDKAgDtieACXIIxRnknz+hLv17ntX/ajb310zsGEFoAwAYEF+Ai1bX1+uOWfC34x27VfW4uy11DeuiXdw1UZDjrNgKAXQgugKRzNXX6c84RZb32sSpr6hocnzGyL09pBgA/QHBByDpVUaW39hTrZ//crfJztQ2OOxzSz+9M07eu76UO3DEEAH6B4IKQUlx2TodOndXKzfn66/ajDY7HdAzXA6Ov0rdv7M0icADgh/jLjJBR5DmnGxf+q8H+K7tG62uDe+ir1/VQSnyUDZUBAJqK4IKQkb27yNru2zVazvAw/fae61iHBQACCMEFIWN7Qakk6RsZPfWf3xxkbzEAgBbhvk6EjO35pZKkrwzsbm8hAIAWI7ggJJw+U628k2ckSUNS4myuBgDQUgQXBL28k2f05d+sl3R+Iq4rKsLmigAALUVwQdD73dsHVFR2TpI0rG+CzdUAAFqD4IKgt+1wiaTzD0acO+Fqm6sBALQGwQVBreRMtQ5+Ordl0dcHqktUpM0VAQBag+CCoLbj01ug+yZGKy6a0AIAgY7ggqC2Lf/810SDe3EnEQAEA4ILgtqFtVsG9+piax0AAN8guCBo1dUb66uiIYy4AEBQILggaO0vLldFVa2iIsPUL6mz3eUAAHyA4IKgdeFrokE9uyg8jH/qABAMeMgi/E727uNau/t4q89jfU3Uu0urzwUA8A8EF/iV+nqj7z+/1afnvKEPq+UCQLAguMCv7C4ss7YfHttP4WGOVp0vObajRqQmtrYsAICfILjAr7y1p1iSFBUZph+OTbW5GgCAv2HGIvxGTV29Xnj/sCTpRh6GCABoBCMu8AvGGF37+JuqrquXJI3s19XmigAA/ogRF/iF3YVlVmi5OilGXxvSw+aKAAD+iBEX+IXsT29/7hjRQW8+PMLmagAA/qpZIy5Lly7VwIEDFRsbq9jYWN100016/fXXrePGGC1YsEBut1udOnXSqFGjtGvXLq9zVFVVadasWUpMTFR0dLQmTZqkI0eO+KY3CFgb9p+UxNwWAMAXa1Zw6dmzp375y19q69at2rp1q770pS/pq1/9qhVOFi1apKeeekpLlizRli1blJycrHHjxqm8vNw6R2ZmplatWqWVK1dqw4YNqqio0MSJE1VXV+fbniGgbD18/inOI1KZ2wIAuDSHMca05gTx8fH61a9+pe9+97tyu93KzMzUvHnzJJ0fXUlKStKTTz6pGTNmyOPxqGvXrlqxYoXuvvtuSdKxY8eUkpKi1atXa8KECU36nWVlZXK5XPJ4PIqNjW1N+fADnsoaDXpijSTp3UdHq1dClM0VAQDagi8+v1s8Obeurk4rV67UmTNndNNNNykvL09FRUUaP3681cbpdGrkyJHauHGjJCknJ0c1NTVebdxut9LS0qw2jamqqlJZWZnXC8Fj3b4T1nZKfCcbKwEA+LtmB5edO3eqc+fOcjqduv/++7Vq1SoNGDBARUVFkqSkpCSv9klJSdaxoqIiRUZGKi4u7pJtGrNw4UK5XC7rlZKS0tyy4cf+9fH5iblXJETJ4WjdSrkAgODW7OBy9dVXa8eOHXr//ff1gx/8QPfdd592795tHb/4g8cYc9kPo8u1mT9/vjwej/UqKChobtnwYxcm5o66upvNlQAA/F2zg0tkZKSuuuoqDR06VAsXLtSgQYP0m9/8RsnJyZLUYOSkuLjYGoVJTk5WdXW1SkpKLtmmMU6n07qT6cILwcEYo1NnqiVJ4wZc+t8AAACSDxagM8aoqqpKffr0UXJysrKzs61j1dXVWrdunYYPHy5JysjIUEREhFebwsJC5ebmWm0QWj4u/OyOs4zecV/QEgCAZi5A9+Mf/1i33367UlJSVF5erpUrV+qdd97RG2+8IYfDoczMTGVlZSk1NVWpqanKyspSVFSUpkyZIklyuVyaPn26Zs+erYSEBMXHx2vOnDlKT0/X2LFj26SD8F9nq2s195UPJUkRYQ51jAizuSIAgL9rVnA5fvy4pk2bpsLCQrlcLg0cOFBvvPGGxo0bJ0maO3euKisrNXPmTJWUlGjYsGFas2aNYmJirHMsXrxY4eHhmjx5siorKzVmzBgtX75cYWF8aIWa2X/6ULlHz98hxvotAICmaPU6LnZgHZfgcMWPXrO2188drZR41m8BgGBm6zouQGtUVn+2UvLbc0YRWgAATUJwgS0+v+jcFayUCwBoIoILbPH2nmJJUo8unVh0DgDQZAQX2GLjQZ4GDQBoPoIL2p0xRgWnKyVJo/tzNxEAoOkILmh3h06dtbZHs8w/AKAZCC5od2t2ffZYiGhns5YSAgCEOIIL2t2/Pp2YyxL/AIDmIrig3W3OOy1J+lJ/viYCADQPwQXt6jdr91vb43kaNACgmQguaDcFp89q8dp9kqTwDg5d1a2zzRUBAAINwQXt5rn3Dlnbbz48goXnAADNRnBBu/m/9/IkSdd0j9WVXRltAQA0H8EF7eJsda21PXtcPxsrAQAEMoIL2sWfthRY22Ou4W4iAEDLEFzQLlbvPL/oXGzHcOa2AABajOCCNneyokqbD51fu+WB0VfZXA0AIJARXNDm/vi5r4kmD02xsRIAQKAjuKBNfVhQql+9uVeSNHFgd8VFR9pcEQAgkBFc0GaMMfrq796z3v/biL42VgMACAYEF7SZZzfkWduPffkaDezZxb5iAABBgeCCNnPhYYodIzro+4y2AAB8gOCCNmGM0faCUknSC9OH2VsMACBoEFzQJo6UVOpEeZXCOziU1sNldzkAgCBBcEGbuDDacq07Vh0jwuwtBgAQNAguaBPb80skSYN7xdlcCQAgmBBc0Ca25ZdKkgb36mJrHQCA4EJwgc+dq6nT7mMeSdIQRlwAAD5EcIHP/eiVj1RTZ5TY2amecZ3sLgcAEEQILvCp0rPV+tuOY5KkG/vG8yRoAIBPEVzgUxfuJpKkX9yZbl8hAICgRHCBT23/dFLuXYN7yBUVYW8xAICgQ3CBT1m3QfdmUi4AwPcILvCZ+nqjHRdug07pYmstAIDgRHCBzxw4UaHyqlp1ighT/+QYu8sBAAQhggt85sLXRAN7uhQexj8tAIDv8ekCn9l2uFSSNIT5LQCANkJwgc9sL/h0Yi7zWwAAbYTgAp8oO1ej/cUVkhhxAQC0HYILfGJHfqmMkXrFRymxs9PucgAAQYrgAp/YztOgAQDtgOCCVnvvwEktXrtPEk+DBgC0LYILWu2F9w9b27emJtpYCQAg2BFc0CrGGG37dP2WJVMGq2/XzjZXBAAIZgQXtMoxzzkdL6tSeAeHxvRPsrscAECQa1ZwWbhwoa6//nrFxMSoW7duuvPOO7V3716vNsYYLViwQG63W506ddKoUaO0a9curzZVVVWaNWuWEhMTFR0drUmTJunIkSOt7w3a3YXVcq/pHqtOkWE2VwMACHbNCi7r1q3TAw88oPfff1/Z2dmqra3V+PHjdebMGavNokWL9NRTT2nJkiXasmWLkpOTNW7cOJWXl1ttMjMztWrVKq1cuVIbNmxQRUWFJk6cqLq6Ot/1DO3CWi2Xu4kAAO3AYYwxLf3hEydOqFu3blq3bp1GjBghY4zcbrcyMzM1b948SedHV5KSkvTkk09qxowZ8ng86tq1q1asWKG7775bknTs2DGlpKRo9erVmjBhwmV/b1lZmVwulzwej2JjY1taPnzgzt+9px0FpXr67ut05+AedpcDAPBjvvj8btUcF4/HI0mKj4+XJOXl5amoqEjjx4+32jidTo0cOVIbN26UJOXk5KimpsarjdvtVlpamtXmYlVVVSorK/N6wX5VtXXafez8teA2aABAe2hxcDHG6JFHHtEtt9yitLQ0SVJRUZEkKSnJe5JmUlKSdayoqEiRkZGKi4u7ZJuLLVy4UC6Xy3qlpKS0tGz4UO7RMlXX1Suxc6RS4jvZXQ4AIAS0OLg8+OCD+uijj/Tyyy83OOZwOLzeG2Ma7LvYF7WZP3++PB6P9SooKGhp2fChCxNzB/eKu+z1BQDAF1oUXGbNmqVXX31Vb7/9tnr27GntT05OlqQGIyfFxcXWKExycrKqq6tVUlJyyTYXczqdio2N9XrBfizzDwBob80KLsYYPfjgg/rrX/+qt956S3369PE63qdPHyUnJys7O9vaV11drXXr1mn48OGSpIyMDEVERHi1KSwsVG5urtUGgeHCwnPMbwEAtJfw5jR+4IEH9NJLL+nvf/+7YmJirJEVl8ulTp06yeFwKDMzU1lZWUpNTVVqaqqysrIUFRWlKVOmWG2nT5+u2bNnKyEhQfHx8ZozZ47S09M1duxY3/cQPldbV6/HVuWq0HNOYR0cGtjTZXdJAIAQ0azgsnTpUknSqFGjvPY/99xz+s53viNJmjt3riorKzVz5kyVlJRo2LBhWrNmjWJiYqz2ixcvVnh4uCZPnqzKykqNGTNGy5cvV1gYC5gFgrf2FOuPW8/PM0rv4VJUZLP+GQEA0GKtWsfFLqzjYq+Fqz/W/7x7UJL01uyRPJ8IANAktq/jgtB0YW7Lr74xkNACAGhXBBc0S01dvT46cn7hwcFMygUAtDOCC5rl48IyVdXWy9UpQn0To+0uBwAQYgguaJbPr93SoQOLzgEA2hfBBc1yYX7L4BS+JgIAtD+CC5rFWnSudxd7CwEAhCSCC5rsRHmVCk5XyuGQBqV0sbscAEAIIrigSYwx1kMVU7t1VmzHCJsrAgCEIpY8xWVV1dbp60s3KvdomSSeTQQAsA/BJYj97u0DeuH9w2rt2shFZees7fAODn05vXsrKwMAoGUILkHKGKNfvbnXp+ccPyBJT99zHc8mAgDYhk+gILXrWJm1/coPbpIzvHUPsAzr4FC/pBiFsXYLAMBGBJcgtfbj45KkqMgwZfSOt7kaAAB8g7uKgtR7B05Kkq6/gtACAAgeBJcgteXQ+VuXb01NtLkSAAB8h+AShMrP1VjbY69JsrESAAB8i+AShN47cMravoInOAMAggjBJQhdmJibEt/J5koAAPAtgksQWr//hCRpRGpXmysBAMC3CC5Bxhij42VVkqRxA5jfAgAILgSXIHK2ulbf/O9N1vthfRJsrAYAAN8juASRFZsOa+vh87dBX50Uo06RrVstFwAAf0NwCSLv7D1hba/43g02VgIAQNsguASRzYdOS5IenXC1usV0tLkaAAB8j+ASJM5W16qu3khiUi4AIHgRXILE5rzT1na/pBgbKwEAoO0QXIJE9u7zi851jXHaXAkAAG2H4BIkLkzMvfUqHqoIAAheBJcgYIzR0dJKSdJY5rcAAIIYwSUIfHLijLV9ayojLgCA4EVwCQJv7Tlubcd0jLCxEgAA2hbBJQis23d+fsv1V8TZXAkAAG2L4BLgDp86o/cOnJIk3crToAEAQY7gEsBq6+o18lfvWO/HX8vEXABAcCO4BLD1B05a2z8YdaX6J8faWA0AAG2P4BLA3vq4WJIU4wzXvNv621wNAABtj+ASwDZ8OuIy8mrmtgAAQgPBJYDlnTy/fsvoq7vZXAkAAO2D4BKgjpeds7ZH9ye4AABCA8ElQL21p9jajo+OtLESAADaD8ElQK399GnQ/ZNjbK4EAID2Q3AJUO/uP79a7sh+TMwFAIQOgksAqqmrV02dkcSicwCA0EJwCUDvHzxlbQ/s2cW+QgAAaGfNDi7vvvuu7rjjDrndbjkcDv3tb3/zOm6M0YIFC+R2u9WpUyeNGjVKu3bt8mpTVVWlWbNmKTExUdHR0Zo0aZKOHDnSqo6Eii2HTmvas5slnV94LiKM7AkACB3N/tQ7c+aMBg0apCVLljR6fNGiRXrqqae0ZMkSbdmyRcnJyRo3bpzKy8utNpmZmVq1apVWrlypDRs2qKKiQhMnTlRdXV3LexICzlbX6pv/vcl6/+2bettYDQAA7c9hjDEt/mGHQ6tWrdKdd94p6fxoi9vtVmZmpubNmyfp/OhKUlKSnnzySc2YMUMej0ddu3bVihUrdPfdd0uSjh07ppSUFK1evVoTJky47O8tKyuTy+WSx+NRbGzoPJ9n0pIN+uiIR5L0+MQBmjKslzpGhNlcFQAATeOLz2+ffs+Ql5enoqIijR8/3trndDo1cuRIbdy4UZKUk5OjmpoarzZut1tpaWlWGzT0cWGZFVquvyJO372lD6EFABBywn15sqKiIklSUpL3nS5JSUk6fPiw1SYyMlJxcXEN2lz4+YtVVVWpqqrKel9WVubLsgPCsncPWtsrpg+zsRIAAOzTJjM7HQ6H13tjTIN9F/uiNgsXLpTL5bJeKSkpPqs1ENTU1Vsr5f5wTCojLQCAkOXT4JKcnCxJDUZOiouLrVGY5ORkVVdXq6Sk5JJtLjZ//nx5PB7rVVBQ4Muy/VpNXb1SH3tdnsoaSdLd14dWaAMA4PN8Glz69Omj5ORkZWdnW/uqq6u1bt06DR8+XJKUkZGhiIgIrzaFhYXKzc212lzM6XQqNjbW6xUqJv/PZ3cR3XZtstxdOtlYDQAA9mr2HJeKigodOHDAep+Xl6cdO3YoPj5evXr1UmZmprKyspSamqrU1FRlZWUpKipKU6ZMkSS5XC5Nnz5ds2fPVkJCguLj4zVnzhylp6dr7NixvutZEMg7eUbb80slSde6Y/Xf0zLsLQgAAJs1O7hs3bpVo0ePtt4/8sgjkqT77rtPy5cv19y5c1VZWamZM2eqpKREw4YN05o1axQT89nDABcvXqzw8HBNnjxZlZWVGjNmjJYvX66wMOZufN5Xl2ywtv86s/HRKAAAQkmr1nGxSyis41JVW6f0n65RdV29Hp1wtR4YfZXdJQEA0Cp+t44LfGfXsTJV19VLkmaOutLmagAA8A8EFz9UU1evn/wtV5I09ppul72VHACAUEFw8UN/235Uu46dX2RvcK+4y7QGACB0EFz80PsHT1vbU4f1srESAAD8C8HFD20vOL843/99Z6i6REXaXA0AAP6D4OJnSs9W6+CJM5KkwSl8TQQAwOcRXPzMhQXn+iZGKy6a0RYAAD6P4OJntuef/5roul5d7C0EAAA/RHDxM9s+HXEZwt1EAAA0QHDxI3X1RjsKSiURXAAAaAzBxY8cKK5QRVWtoiLDdHVyzOV/AACAEENw8SPbPp3fMqhnF4V1YLVcAAAuRnDxIxcm5g7p3cXeQgAA8FMEFz9xrqZOf9t+TBLzWwAAuBSCix8wxujLv11vPQ36upQu9hYEAICfIrj4gUOnzlqr5U4e2lMJnZ02VwQAgH8iuPiBbYfPz23J6B2nRd8YZHM1AAD4L4KLH7jwUMXBfEUEAMAXIrj4gW2HSyVJQ3ozKRcAgC9CcLHZmapa7Skqk8TdRAAAXA7BxWYfHfGo3kjdXR2V7OpodzkAAPg1govNLqyWy2gLAACXR3Cx2YXVcgf36mJvIQAABACCi42MMdqeXypJGsyICwAAl0VwsVH+6bM6daZakWEdlNYj1u5yAADwewQXmxw+dUYPvrRdkjTAHStneJjNFQEA4P8ILjZ5aOUO7TzqkXR+xVwAAHB5BBcbVFTVaueRUknS2GuS9P1b+9pbEAAAASLc7gJC0UdHSlVvJLero/73vqF2lwMAQMBgxMUG1p1EfEUEAECzEFxs8Mq2I5J4qCIAAM1FcGln5edqdPDEGUlMygUAoLmY49JExhj9eNVOK3S01NnqOmv7OkZcAABoFoJLEz2/6bBe3lzgs/Pdc32KHA6Hz84HAEAoILg00U9f3WVt/27KkFadKyLMoeFXJba2JAAAQg7BpQkqqmqt7WfuHapxA5JsrAYAgNDF5Nwm2HjgpLVNaAEAwD4ElyZ4e2+xJKl3QpTNlQAAENoILpdRX2+0emeRJGn4lcxLAQDATsxx+QLVtfXq/5PXVW/Ovx/Tv5u9BQEAEOIYcfkCP/vnLiu0pHbrrFtSGXEBAMBOjLhcgjFGL7yfL0lKinXqzcwR6tCBdVcAALATIy6XsGr7UWt7xfRhhBYAAPwAweUSlm88ZG33S4qxrxAAAGAhuDTCGKOPjngkSY9PHGBzNQAA4AJbg8vvf/979enTRx07dlRGRobWr19vZzmWY55z1vZdQ3rYWAkAAPg824LLH//4R2VmZuqxxx7T9u3bdeutt+r2229Xfn6+XSVZsncVWdtdoiJtrAQAAHyebcHlqaee0vTp0/W9731P11xzjZ5++mmlpKRo6dKldpUkSXp33wkt+MduSdJ1KV1srQUAAHizJbhUV1crJydH48eP99o/fvx4bdy4sUH7qqoqlZWVeb3awonyKt37f5ut99/I6NkmvwcAALSMLcHl5MmTqqurU1KS9wMLk5KSVFRU1KD9woUL5XK5rFdKSkqb1FV6tvqz33lXuu6+vm1+DwAAaBlbJ+c6HN5roxhjGuyTpPnz58vj8VivgoKCNqknobNTD4y+Uiv/7UZ964ZeigjjpisAAPyJLSvnJiYmKiwsrMHoSnFxcYNRGElyOp1yOp1tXld8dKQendC/zX8PAABoGVuGFCIjI5WRkaHs7Gyv/dnZ2Ro+fLgdJQEAgABg27OKHnnkEU2bNk1Dhw7VTTfdpGXLlik/P1/333+/XSUBAAA/Z1twufvuu3Xq1Cn97Gc/U2FhodLS0rR69Wr17t3brpIAAICfcxhjjN1FNFdZWZlcLpc8Ho9iY2PtLgcAADSBLz6/uW0GAAAEDIILAAAIGAQXAAAQMAguAAAgYBBcAABAwCC4AACAgEFwAQAAAYPgAgAAAgbBBQAABAzblvxvjQuL/ZaVldlcCQAAaKoLn9utWbQ/IINLeXm5JCklJcXmSgAAQHOVl5fL5XK16GcD8llF9fX1OnbsmGJiYuRwOOwup9nKysqUkpKigoKCkHrWEv0OnX6HYp8l+h1K/Q7FPkut77cxRuXl5XK73erQoWWzVQJyxKVDhw7q2bOn3WW0WmxsbEj9g7+AfoeOUOyzRL9DSSj2WWpdv1s60nIBk3MBAEDAILgAAICAQXCxgdPp1E9/+lM5nU67S2lX9Dt0+h2KfZbodyj1OxT7LPlHvwNyci4AAAhNjLgAAICAQXABAAABg+ACAAACBsEFAAAEDIJLEx09elTf/va3lZCQoKioKF133XXKycmxjv/1r3/VhAkTlJiYKIfDoR07djQ4R1VVlWbNmqXExERFR0dr0qRJOnLkiFebkpISTZs2TS6XSy6XS9OmTVNpaalXm/z8fN1xxx2Kjo5WYmKiHnroIVVXV7dFt33S71GjRsnhcHi97rnnHq82/tTvL+pzTU2N5s2bp/T0dEVHR8vtduvee+/VsWPHvM4RbNe6qf0OtGt9uX5L0oIFC9S/f39FR0crLi5OY8eO1QcffOB1jmC73k3td6Bd78v1+fNmzJghh8Ohp59+2mt/MF7rz7tUv/3qWhtc1unTp03v3r3Nd77zHfPBBx+YvLw8s3btWnPgwAGrzfPPP2+eeOIJ88wzzxhJZvv27Q3Oc//995sePXqY7Oxss23bNjN69GgzaNAgU1tba7W57bbbTFpamtm4caPZuHGjSUtLMxMnTrSO19bWmrS0NDN69Gizbds2k52dbdxut3nwwQf9tt8jR4403//+901hYaH1Ki0t9WrjL/2+XJ9LS0vN2LFjzR//+EezZ88es2nTJjNs2DCTkZHhdZ5gu9ZN7XcgXeum9NsYY1588UWTnZ1tPvnkE5Obm2umT59uYmNjTXFxsdUm2K53U/sdSNe7KX2+YNWqVWbQoEHG7XabxYsXex0LxmvdlH7707UmuDTBvHnzzC233NKktnl5eY1+gJeWlpqIiAizcuVKa9/Ro0dNhw4dzBtvvGGMMWb37t1Gknn//fetNps2bTKSzJ49e4wxxqxevdp06NDBHD161Grz8ssvG6fTaTweT0u72Chf9NuY8//gf/jDH17yZ/2p383p8wWbN282kszhw4eNMcF/rS+4uN/GBNa1NqZl/fZ4PEaSWbt2rTEmdK73xf02JrCud1P7fOTIEdOjRw+Tm5trevfu7fUBHszX+ov6bYx/XWu+KmqCV199VUOHDtU3v/lNdevWTYMHD9YzzzzTrHPk5OSopqZG48ePt/a53W6lpaVp48aNkqRNmzbJ5XJp2LBhVpsbb7xRLpfLq01aWprcbrfVZsKECaqqqrrk0F9L+aLfF7z44otKTEzUtddeqzlz5lhP+Jb8q98t6bPH45HD4VCXLl0khc61vrjfFwTKtZaa3+/q6motW7ZMLpdLgwYNkhQa17uxfl8QKNe7KX2ur6/XtGnT9Oijj+raa69tcI5gvdaX6/cF/nKtCS5NcPDgQS1dulSpqal68803df/99+uhhx7S888/3+RzFBUVKTIyUnFxcV77k5KSVFRUZLXp1q1bg5/t1q2bV5ukpCSv43FxcYqMjLTa+Iov+i1JU6dO1csvv6x33nlHP/nJT/TKK6/orrvuso77U7+b2+dz587pRz/6kaZMmWI9cCwUrnVj/ZYC61pLTe/3P//5T3Xu3FkdO3bU4sWLlZ2drcTERKveYL3eX9RvKbCud1P6/OSTTyo8PFwPPfRQo+cI1mt9uX5L/nWtA/Lp0O2tvr5eQ4cOVVZWliRp8ODB2rVrl5YuXap77723Vec2xsjhcFjvP7/dmja+4Kt+f//737e209LSlJqaqqFDh2rbtm0aMmSIJP/pd3P6XFNTo3vuuUf19fX6/e9/f9lzB8u1/qJ+B9K1lpre79GjR2vHjh06efKknnnmGU2ePFkffPBBo3+oL1VvMPY7kK735fqck5Oj3/zmN9q2bVuzf28gX+um9tufrjUjLk3QvXt3DRgwwGvfNddco/z8/CafIzk5WdXV1SopKfHaX1xcbCXQ5ORkHT9+vMHPnjhxwqvNxcm0pKRENTU1DZJsa/mi340ZMmSIIiIitH//fkn+1e+m9rmmpkaTJ09WXl6esrOzvUYdgvlaf1G/G+PP11pqer+jo6N11VVX6cYbb9Szzz6r8PBwPfvss1a9wXq9v6jfjfHn6325Pq9fv17FxcXq1auXwsPDFR4ersOHD2v27Nm64oorrFqD7Vo3pd+NsfNaE1ya4Oabb9bevXu99u3bt0+9e/du8jkyMjIUERGh7Oxsa19hYaFyc3M1fPhwSdJNN90kj8ejzZs3W20++OADeTwerza5ubkqLCy02qxZs0ZOp1MZGRkt6t+l+KLfjdm1a5dqamrUvXt3Sf7V76b0+cKH9/79+7V27VolJCR4tQ/Wa325fjfGn6+11PJ/48YYVVVVSQre692Yz/e7Mf58vS/X52nTpumjjz7Sjh07rJfb7dajjz6qN998U1JwXuum9Lsxtl7rJk/jDWGbN2824eHh5he/+IXZv3+/efHFF01UVJR54YUXrDanTp0y27dvN6+99pqRZFauXGm2b99uCgsLrTb333+/6dmzp1m7dq3Ztm2b+dKXvtTobXQDBw40mzZtMps2bTLp6emN3k42ZswYs23bNrN27VrTs2fPNrmNzhf9PnDggHniiSfMli1bTF5ennnttddM//79zeDBg/2y35frc01NjZk0aZLp2bOn2bFjh9etgVVVVdZ5gu1aN6XfgXatm9LviooKM3/+fLNp0yZz6NAhk5OTY6ZPn26cTqfJzc21zhNs17sp/Q60692Uv2cXa+zummC71k3pt79da4JLE/3jH/8waWlpxul0mv79+5tly5Z5HX/uueeMpAavn/70p1abyspK8+CDD5r4+HjTqVMnM3HiRJOfn+91nlOnTpmpU6eamJgYExMTY6ZOnWpKSkq82hw+fNh85StfMZ06dTLx8fHmwQcfNOfOnfPLfufn55sRI0aY+Ph4ExkZaa688krz0EMPmVOnTvltv7+ozxdu+27s9fbbb1vtgu1aN6XfgXitL9fvyspK87Wvfc243W4TGRlpunfvbiZNmmQ2b97sdY5gu95N6XcgXu/L/T27WGPBJdiudWMu7re/XWuHMcY0fXwGAADAPsxxAQAAAYPgAgAAAgbBBQAABAyCCwAACBgEFwAAEDAILgAAIGAQXAAAQMAguAAAgIBBcAEAAAGD4AIAAAIGwQUAAAQMggsAAAgY/x8eK03ofCkX8wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ddf_name = 'COSMOS'\n",
    "mjds, night_mjd, cumulative_desired, cumulative_sched = optimize_ddf_times(\n",
    "            ddf_name,\n",
    "            ddfs[ddf_name][0],\n",
    "            ddf_grid,\n",
    "            season_unobs_frac=0.1,\n",
    "            sequence_time=60,\n",
    "        )\n",
    "plt.plot(night_mjd, cumulative_sched)\n",
    "print(cumulative_sched[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "57b509b9-2ff1-4f28-8cd7-2dd462b90de3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.404004105285132"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "first = np.where(night_season<0)\n",
    "good_first_season = night_season[first][np.where(vals[2][first])]\n",
    "(good_first_season.max() - good_first_season.min()) * 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "42fc8884-8b7c-45aa-a7c2-787912c22093",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "300.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGdCAYAAAA1/PiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAjElEQVR4nO3df3zP9f7/8fts9mYzY8a2N7MU+TXElB/5/dsJ6RdS5KSOjtBCpM456Xs6lE6olI5OEQmfin4cwpz8SCPMb0XU/N7Mj/1kv/f6/uH00ruGbba93j9u18vlfbm8n4/38/32ePbSdvd6v354GYZhCAAAwAVUsLoBAACAoiK4AAAAl0FwAQAALoPgAgAAXAbBBQAAuAyCCwAAcBkEFwAA4DIILgAAwGX4WN1ASRQUFOj06dMKCAiQl5eX1e0AAIAiMAxD6enpstvtqlChZPtOXDK4nD59WuHh4Va3AQAASuDEiROqU6dOid7rksElICBA0uWFV61a1eJuAABAUaSlpSk8PNz8PV4SLhlcfvl6qGrVqgQXAABczI0c5sHBuQAAwGUQXAAAgMsguAAAAJdBcAEAAC6D4AIAAFwGwQUAALgMggsAAHAZBBcAAOAyCC4AAMBlFCu4zJ07V82bNzevWNuuXTt99dVX5uuGYWjq1Kmy2+2qXLmyunTpogMHDjh8RnZ2tsaOHavg4GD5+/trwIABOnnyZOmsBgAAuLViBZc6dero5Zdf1o4dO7Rjxw5169ZNd999txlOZsyYoZkzZ2rOnDnavn27QkND1bNnT6Wnp5ufER0drRUrVmjp0qXavHmzMjIy1K9fP+Xn55fuygAAgNvxMgzDuJEPCAoK0quvvqpHH31Udrtd0dHRmjx5sqTLe1dCQkL0yiuvaNSoUUpNTVXNmjW1aNEiDR48WNKVOz2vWrVKvXv3LtKfmZaWpsDAQKWmpnKvIgAAXERp/P4u8U0W8/Pz9fHHH+vixYtq166d4uPjlZiYqF69eplzbDabOnfurNjYWI0aNUpxcXHKzc11mGO32xUZGanY2NirBpfs7GxlZ2eb47S0tJK2DQAACnEy+ZI+3Hpc2XmXvwEJrmLTk13rW9zV7xU7uOzbt0/t2rVTVlaWqlSpohUrVqhJkyaKjY2VJIWEhDjMDwkJ0bFjxyRJiYmJ8vX1VfXq1X83JzEx8ap/5vTp0/Xiiy8Wt1UAAFAE+QWG+r7+jdKz8szazTX93SO4NGzYULt371ZKSoo+/fRTPfLII9q4caP5+m9vVW0YxnVvX329OVOmTNH48ePNcVpamsLDw4vbOgAAKMT7m+PN0BIeVFkDWthV3c/X4q4KV+zg4uvrq/r1Lyew1q1ba/v27Xr99dfN41oSExMVFhZmzk9KSjL3woSGhionJ0fJyckOe12SkpLUvn37q/6ZNptNNputuK0CAIDrSMvK1T9W/WCO10Z3VmVfbws7urYbvo6LYRjKzs5WvXr1FBoaqpiYGPO1nJwcbdy40QwlUVFRqlixosOchIQE7d+//5rBBQAAlC7DMHQ+I1uPf7DDrH0x5k6nDi1SMfe4PPfcc+rbt6/Cw8OVnp6upUuXasOGDVq9erW8vLwUHR2tadOmqUGDBmrQoIGmTZsmPz8/DR06VJIUGBiokSNHasKECapRo4aCgoI0ceJENWvWTD169CiTBQIAgN977IMd+u/BJHN8Z/0aal6nmnUNFVGxgsuZM2c0bNgwJSQkKDAwUM2bN9fq1avVs2dPSdKkSZOUmZmp0aNHKzk5WW3atNHatWsVEBBgfsasWbPk4+OjQYMGKTMzU927d9eCBQvk7e3cCQ8AAHdwPiNb982N1dHzl8xa7WqV9eaDrSzsquhu+DouVuA6LgAAFN9nu04petluc1y/VhWtie4k7wrXPommtFh6HRcAAOA6/rnmkOasP2KOezUJ0VsPtSq30FJaCC4AALgxwzD0/Gf79dF3x83a60Nu013NwuTj7Xr3Wia4AADgpnYdT9acr484HIS7fHR7tapb/Rrvcm4EFwAA3ExefoE2/nhWI391qrMkrYnupIahAVd5l2sguAAA4EYysvP0j5U/aMm2K18N3R9VRyPa3+TyoUUiuAAA4DYysvPU5dX1OpeRY9b+2q+JRnaoZ2FXpYvgAgCAGzhx4ZK6z9yonLwCSZI9sJLmDW+tyNqBFndWugguAAC4uM2Hz+nh974zx90a1dJ7j7S+7k2OXRHBBQAAF1VQYGj04p1afSDRrE3q01Cju9S3sKuyRXABAMAF7T2ZoicWxel0apZZ++cDLXR/VB0Luyp7BBcAAFzMf/ae1piPdpnjKjYfrRzXQRE1/C3sqnwQXAAAcBGXcvK0aMsxTf/qoFm7t1Vt/eWuJgry97Wws/JDcAEAwAVk5uRrxPvbte3oBbM24/7muq9VHZe739CNILgAAODkDiWm6y+f7dP2o8lm7d/DW6tHkxALu7IGwQUAACdlGIbWHDijJz6Mc6h/M6mrwoP8LOrKWgQXAACcUEGBoYmf7NHynafMWof6wXrl/uaqXa2yhZ1Zi+ACAICTyckrUL83v9GPZzLM2qzBLXRPS/c+1bkoCC4AADiRDYeSNGL+dnMcFlhJq5/qpEC/ihZ25TwILgAAOIn1h5L0x1+Flg71g/XBo3d41FlD10NwAQDACZzLyHYILS8OaKqH20YQWn6D4AIAgBN4cN5W8/lf7mqs4e0i3PImiTeqgtUNAADg6f6z97QOJ10+ELd5nUCN7FCP0HIVBBcAACyUm1/gcN+hhY/eQWi5BoILAAAWen9zvPl83rAoVfPzjHsOlRTBBQAAC/1ygbkgf1/1ahpqcTfOj+ACAICFKvt6S5LGdatvcSeugbOKAAAoZXO+Pqy31v+kAsO47tzsvAJJUmTtwLJuyy0QXAAAKEUXs/P0z7U/Fus91f0qqkFIQBl15F4ILgAAlKJt8RfM51+O6aDq/te/VH8Nf5v5lRGujeACAEApiD93UQu+jTevx1K/VhU1q8PXP6WN4AIAwA1afzBJf1yw3aF2W3g1a5pxcwQXAABKKDe/QCt2ntKkT/eatab2qrr7Nrvuvq22hZ25L4ILAAAlkJ6Vq6lffK9Pd540a+O6N9CYrvXl68PVRsoKwQUAgGJKuZSjzq9uUGpmrll7a2gr3dU8zMKuPAPBBQCAYtj683kN+dWdnG+q4ad3hkWpUWhVC7vyHAQXAACKaPa6HzV73WFz3LtpiN55OIqbIpYjggsAANeRnpWrRxds1/ajyWZtztCW6tfcbmFXnongAgDANew6nqx73o51qG16pqvq1vCzqCPPRnABAOAq5n8brxe//N4cNwmrqkUj71CNKjYLu/JsBBcAAAox5+vDDvccerhtXU3p21j+Nn51Won/+gAA/MqRpAx9EHtUi7YeM2uvD7lN/ZvbVaECB+FajeACAIAkwzC0+0TK745n+fTP7RQVEWRRV/gtggsAwOMVFBh6Z9NPmrH6kFnrdGtNje5yC6HFyRBcAAAeLb/A0FNLd+k/exPM2ugut2hSn0YWdoWrIbgAADxWUnqW/vD6NzqXkWPWlv6prdreXMPCrnAtxboL1PTp03X77bcrICBAtWrV0sCBA3Xo0CGHOSNGjJCXl5fDo23btg5zsrOzNXbsWAUHB8vf318DBgzQyZMnBQBAeTlx4ZLu+Md/zdBSp3pl7Zvai9Di5IoVXDZu3Kgnn3xSW7duVUxMjPLy8tSrVy9dvHjRYV6fPn2UkJBgPlatWuXwenR0tFasWKGlS5dq8+bNysjIUL9+/ZSfn3/jKwIA4BoMw9BfPtunjjPWm7Xh7SK06ZmuCqhU0cLOUBTF+qpo9erVDuP58+erVq1aiouLU6dOncy6zWZTaGhooZ+Rmpqq9957T4sWLVKPHj0kSR9++KHCw8O1bt069e7du7hrAACgSH46m6FJn+xV3LErl+6/r1UdvTigKfcbchHF2uPyW6mpqZKkoCDHI643bNigWrVq6dZbb9Xjjz+upKQk87W4uDjl5uaqV69eZs1utysyMlKxsY6noP0iOztbaWlpDg8AAIpj8+Fz6v7aRofQsuxPbfXKfc0ILS6kxMHFMAyNHz9eHTp0UGRkpFnv27evFi9erK+//lqvvfaatm/frm7duik7O1uSlJiYKF9fX1WvXt3h80JCQpSYmFjonzV9+nQFBgaaj/Dw8JK2DQDwMDl5Bfok7qQefu87s9amXpC+mdRVbW6uIR/vG/o3PMpZic8qGjNmjPbu3avNmzc71AcPHmw+j4yMVOvWrRUREaGVK1fq3nvvvernGYZx1cQ7ZcoUjR8/3hynpaURXgAA15WbX6DJn+7Vil2nzNrEXrfq8U43y+bjbWFnKKkSBZexY8fqiy++0KZNm1SnTp1rzg0LC1NERIQOHz4sSQoNDVVOTo6Sk5Md9rokJSWpffv2hX6GzWaTzcYNrQAARZedl68pn+5zCC0z7muuQbfzD19XVqz9Y4ZhaMyYMVq+fLm+/vpr1atX77rvOX/+vE6cOKGwsDBJUlRUlCpWrKiYmBhzTkJCgvbv33/V4AIAQHE98M4WLf9VaFkT3YnQ4gaKtcflySef1EcffaTPP/9cAQEB5jEpgYGBqly5sjIyMjR16lTdd999CgsL09GjR/Xcc88pODhY99xzjzl35MiRmjBhgmrUqKGgoCBNnDhRzZo1M88yAgDgRny49Zj2nrx8AklAJR+tGH2n6teqYnFXKA3FCi5z586VJHXp0sWhPn/+fI0YMULe3t7at2+fFi5cqJSUFIWFhalr165atmyZAgICzPmzZs2Sj4+PBg0apMzMTHXv3l0LFiyQtzffNwIAbsyx8xf1l8/2m+Pdf+slb+7q7Da8DMMwrG6iuNLS0hQYGKjU1FRVrVrV6nYAAE7k6WW7zeNavnqqoxqH8XvCWZTG72/OAQMAuJW0zFxJUof6wYQWN8RNFgEAbiHlUo7e3vCT9vzv2JaH20ZY3BHKAsEFAOAWRszfrt0nUsxx3SA/65pBmeGrIgCAyzuYmGaGFl+fCnrjwZZqYudrInfEHhcAgNMqKDC071SqsnLzrzlvzvoj5vPYZ7spuAoXLXVXBBcAgNN6Zc1B/Wvjz0Wef0e9IEKLmyO4AACc0sHENIfQcktN/2vOr2Lz0Qv9m5R1W7AYwQUA4HReWX1Qczf8ZI43PtNFETWuHVzgGQguAACnkZtfoMc+2KGNP541ay8NjCS0wERwAQA4hR1HL2jE/O3KyM4za98+2021q1W2sCs4G4ILAMByR5LSdf87W8xxrQCb1kR3UnV/Xwu7gjMiuAAALJOelasPtx7XK6sPmrUR7W/S0z1vVWDlihZ2BmdFcAEAWOLwmXQ9vnCHjp6/ZNYGtLDrhf5N5OXF3ZxROIILAKDcHTidqrve2OxQm9ynkR7rWI/QgmsiuAAAyk1BgaHlu05p4sd7zNpdzcP017uaKDSwkoWdwVUQXAAA5cIwDP15cZzWHDhj1sZ1q6/xvRpa2BVcDcEFAFDmCgoM3f9OrHYeTzFr7zzcSn0iw6xrCi6J4AIAKFNf7UvQmCW7lF9gSJLqBvlpdXRH+fnyKwjFx98aAECZMAxDU5bv09LtJ8xas9qB+vzJO1WhAgfgomQILgCAUpeYmqU/L47Trl99NfTPB1ro7tvshBbcEIILAKBUbfzxrB55f5tDbdW4jmpir2pRR3AnBBcAQKn57w9nNPKDHeb4lpr+WjiyDfcbQqkhuAAAbljKpRwt3X5CL3915dL9E3reqhF33qSASly6H6WH4AIAuCGJqVnq9+ZmncvINmszB7XQva3qWNgV3BXBBQBQYrtPpGjgW9+aY39fb80e0lI9m4RY2BXcGcEFAFAiH8Qe1QtfHDDHXRvW1LzhrVXRu4KFXcHdEVwAAMWSm1+goe9u1fajyWbtpYGRerhthIVdwVMQXAAARXY2PVt9Zm/S+Ys5kiRf7wpaOa6DGoQEWNwZPAXBBQBQJIu/O6bnV+w3xzfX9NfnT97JWUMoVwQXAMB1fbzjhENoaR1RXR8+1kaVKnpb2BU8EcEFAHBNP53N0DOf7DXHbzzYUj0a1yK0wBIEFwDANb21/oj5/P9GtdMd9YIs7AaejnPWAABXlZCaqbUHzkiS+jUPI7TAcuxxAQAU6lxGttpN/9ocD2hht7Ab4DL2uAAACnXP21euiNuxQbDa1w+2sBvgMva4AAB+Z8G38TpxIVOS9HDbunppYDOLOwIuY48LAOB3Vh9INJ9P7d/Uwk4AR+xxAQAPMnPtIcX+dP66875PSJMkLXz0Dvlw7yE4EYILAHiIrT+f1xtfH7n+xP+p4CXVr1WlDDsCio/gAgAeICs3X0PmbTXH7zzc6rrviajhL3u1ymXZFlBsBBcAcHPHz1/S0/+32xw/94dG6hMZZl1DwA0guACAmzIMQwdOp6nfm5sd6iM73GxRR8CN44grAHBTa78/4xBa7rgpSOsndpF3BS8LuwJuDHtcAMDN5BcYOnHhkkYtijNrw9tF6MUBTeXlRWiBayO4AIAbSc/KVZ/Z3+hUSqZZe39Ea3VrFGJhV0DpKdZXRdOnT9ftt9+ugIAA1apVSwMHDtShQ4cc5hiGoalTp8put6ty5crq0qWLDhw44DAnOztbY8eOVXBwsPz9/TVgwACdPHnyxlcDAB7syz2n1WzqWofQEt2jAaEFbqVYwWXjxo168skntXXrVsXExCgvL0+9evXSxYsXzTkzZszQzJkzNWfOHG3fvl2hoaHq2bOn0tPTzTnR0dFasWKFli5dqs2bNysjI0P9+vVTfn5+6a0MADzIK6sPauySXeZ4QAu7Dv69j6J73GphV0Dp8zIMwyjpm8+ePatatWpp48aN6tSpkwzDkN1uV3R0tCZPnizp8t6VkJAQvfLKKxo1apRSU1NVs2ZNLVq0SIMHD5YknT59WuHh4Vq1apV69+593T83LS1NgYGBSk1NVdWqVUvaPgC4vFMpmZr6xQHFfH/GrL35YEvd1SxMFTgIF06mNH5/39BZRampqZKkoKAgSVJ8fLwSExPVq1cvc47NZlPnzp0VGxsrSYqLi1Nubq7DHLvdrsjISHPOb2VnZystLc3hAQCe7lBiuu58+WuH0LI6uqP6t7ATWuC2ShxcDMPQ+PHj1aFDB0VGRkqSEhMv35QrJMTx+9SQkBDztcTERPn6+qp69epXnfNb06dPV2BgoPkIDw8vadsA4PLyCwx9tuuUes/eZNZuDami/07orEah7IWGeytxcBkzZoz27t2rJUuW/O61355uZxjGdU/Bu9acKVOmKDU11XycOHGipG0DgEtLvpij55bvU/Sy3Watd9MQLR99p26pyX2F4P5KdDr02LFj9cUXX2jTpk2qU6eOWQ8NDZV0ea9KWNiVy0knJSWZe2FCQ0OVk5Oj5ORkh70uSUlJat++faF/ns1mk81mK0mrAOA2zmVkq/trG5WamWvWnureQE91b8BXQ/AYxdrjYhiGxowZo+XLl+vrr79WvXr1HF6vV6+eQkNDFRMTY9ZycnK0ceNGM5RERUWpYsWKDnMSEhK0f//+qwYXAPB08ecuqvVL68zQUjPAps+evFNP97yV0AKPUqw9Lk8++aQ++ugjff755woICDCPSQkMDFTlypXl5eWl6OhoTZs2TQ0aNFCDBg00bdo0+fn5aejQoebckSNHasKECapRo4aCgoI0ceJENWvWTD169Cj9FQKAi/v3Nz/rpZU/mOOeTUL0zsNRXLofHqlYwWXu3LmSpC5dujjU58+frxEjRkiSJk2apMzMTI0ePVrJyclq06aN1q5dq4CAAHP+rFmz5OPjo0GDBikzM1Pdu3fXggUL5O3tfWOrAQA3YhiGRi2K09pfnTX0dI9b9VSPBhZ2BVjrhq7jYhWu4wLA3f2QkKZRi+J0/MIlSVJFby/FPN1ZNwX7W9wZUHKl8fubexUBgJP5cs9ph6vgBlauqG+f7aYqNn5kA/xfAABOIis3X2+tP6I3vz5i1h6IqqO/9W9CaAH+h/8TAMAJHElK17glu/V9wpUrg796f3Pd07K2fLxv6CLngFshuACAxXYdT9Y9bzve8mT+H29X14a1LOoIcF4EFwCwSF5+gdZ+f0ajF+80ax0bBOulgZGKqMFBuEBhCC4AYIHc/AJFL9utlXsTzNoTnW/R5D4Nr3uLFMCTEVwAoJwZhqH+b27WwcR0s/bq/c31QGtuIAtcD8EFAMrZy18dNENL7WqVtWJ0e9WqWsnirgDXQHABgHL09oYj+temnyVJwVV89e2z3SzuCHAtnGMHAOVkz4kUzVh9yByvHNfRwm4A10RwAYBykJmTr7vf+tYcb3ymi0L4eggoNoILAJSDuRuuXA33uT804nRnoIQILgBQDtYcuHKH55EdbrawE8C1cXAuAJShggJDu0+mKCUzR5L0xoMt5V2B67QAJUVwAYASupSTp1PJmdec8+XeBL3x38PmuFntwLJuC3BrBBcAKIGcvAJ1f22jElKzivyeB++oq5tq+JVhV4D7I7gAQAkkpWeZoSXI3/eac/1t3po9+DZFRQSVR2uAWyO4AEAxfRJ3UtNW/SBJqhlg0/bne1jcEeA5CC4AUAx//8/3em9zvDm+Lbyadc0AHojgAgBF9O9vfnYILfOGRal74xALOwI8D8EFAK7ju5/Pa82BM3r/2yuhJebpTmoQEmBhV4BnIrgAwFXk5Rfoy72n9fSyPQ71r57qSGgBLEJwAYBCXLiYo5f+872W7zpl1h68o656NQlR47CqFnYGeDaCCwD8xtn0bHV5db0u5uSbtb8PjNSwthEWdgVAIrgAgIOj5y6qyz83mOPa1SrrX8OiFMkVbwGnQHABgP/Z+vN5DZm31Rz3aRqqtx9qpQrcWwhwGgQXAB7PMAw9umC71h86a9ae6d1QT3atb2FXAApDcAHg0fadTNVjC7frTFq2WXvzwZbq38JuYVcArobgAsBjrdh10uFU5+AqNq19utN17z0EwDoEFwAeJzMnX3PWH9Zb638yaw+1qavn/tBY/jZ+LALOjP9DAXiUExcuaeQH2/XjmQyz9vqQ29SvuV3eHIQLOD2CCwCPEXcsWffNjXWoLXz0DnW6taZFHQEoLoILALdXUGBo5b4EjV2yy6w1DquqecOiFB7kZ2FnAIqL4ALAreXkFWjaqh+0IPaoWXu4bV1N7d9UPt4VrGsMQIkQXAC4rbz8AvWZvUk/n7to1mYPvk0DW9a2sCsAN4LgAsAtbfzxrB55f5s5DguspIWP3sFdnQEXR3AB4Hbe/O9hvRbzoznu0biW/v3I7RZ2BKC0EFwAuI3kizl6dvlerTlwxqy9en9z3R9Vx8KuAJQmggsAt7DnRIrufutbh9rq6I5qFFrVoo4AlAWCCwCX90ncSU38+Mql+8ODKmvho21UL9jfwq4AlAWCCwCX9kHsUb3wxQFzfF+rOvpbvyYK9KtoYVcAygrBBYDL2nsyxSG0vNC/iUa0v0leXly6H3BXBBcALskwDA2Yc+WYliWPt1W7W2pY2BGA8sBlIwG4pBe//N58PrV/E0IL4CEILgBc0q4TKZKkGv6+eqT9TZb2AqD8FDu4bNq0Sf3795fdbpeXl5c+++wzh9dHjBghLy8vh0fbtm0d5mRnZ2vs2LEKDg6Wv7+/BgwYoJMnT97QQgB4BsMwNG7JLu35X3CZNzyKY1oAD1Ls4HLx4kW1aNFCc+bMueqcPn36KCEhwXysWrXK4fXo6GitWLFCS5cu1ebNm5WRkaF+/fopPz+/+CsA4FEWf3dcX+w5LUny9/VW/Zpcwh/wJMU+OLdv377q27fvNefYbDaFhoYW+lpqaqree+89LVq0SD169JAkffjhhwoPD9e6devUu3fv4rYEwEOcSsnUXz7bb463PtddAZU47RnwJGVyjMuGDRtUq1Yt3XrrrXr88ceVlJRkvhYXF6fc3Fz16tXLrNntdkVGRio2NrbQz8vOzlZaWprDA4Dn+cfKKwfkrhjdntACeKBSPx26b9++euCBBxQREaH4+Hj99a9/Vbdu3RQXFyebzabExET5+vqqevXqDu8LCQlRYmJioZ85ffp0vfjii6XdKgAnsfN4snYeS77uvFX7Lv+MuKNekFrWrX6d2QDcUakHl8GDB5vPIyMj1bp1a0VERGjlypW69957r/o+wzCueoDdlClTNH78eHOclpam8PDw0msagGVSM3N179uF7229mujuDcqoGwDOrswvQBcWFqaIiAgdPnxYkhQaGqqcnBwlJyc77HVJSkpS+/btC/0Mm80mm81W1q0CKGf5BYaGvfedOb77Nruud37QTcH+XLMF8GBlHlzOnz+vEydOKCwsTJIUFRWlihUrKiYmRoMGDZIkJSQkaP/+/ZoxY0ZZtwPASaRn5erpZXu092SqpMtf/7w+pKXFXQFwdsUOLhkZGTpy5Ig5jo+P1+7duxUUFKSgoCBNnTpV9913n8LCwnT06FE999xzCg4O1j333CNJCgwM1MiRIzVhwgTVqFFDQUFBmjhxopo1a2aeZQTAvS3bflzPLt8nw7hSm/tQK+saAuAyih1cduzYoa5du5rjX449eeSRRzR37lzt27dPCxcuVEpKisLCwtS1a1ctW7ZMAQFXrrUwa9Ys+fj4aNCgQcrMzFT37t21YMECeXt7l8KSADirggJDw9/fps1Hzpm1m2r46YuxHVSVM4QAFIGXYfz63zyuIS0tTYGBgUpNTVXVqlWtbgdAEew6nqxHF2xX8qVcSVJFby8t//OdalYn0OLOAJSX0vj9zd2hAZS5T+JOauLHe8yxPbCSVj/dib0sAIqN4AKgTC3cclR/+/yAOR7R/iZN6tNQfr78+AFQfPzkAFAmjp+/pLkbj2jJthNm7e2HWql301B5V+CmiABKhuACoNRtP3pBD7yzxaG2+LE2urN+sEUdAXAXBBcApSY3v0Ardp3SpE/2mrVGoQF6aWCkWt8UZGFnANwFwQVAqcjOy9efFsZp449nzdoj7SL0t/5N+WoIQKkhuAC4YTl5Beoxc6NOXMg0a28NbaW7modZ2BUAd0RwAXBD9p5M0YA535rj+rWq6NMn2ivQj1OdAZQ+gguAEpu59pDe+PrKLUD+0CxUbz8UZWFHANwdwQVAsZ1Nz9ZzK/Yp5vszZm1yn0Z6ovPNFnYFwBMQXAAUy87jybr37ViH2ldPdVSj0AB5eXEQLoCyRXABUCSGYej/dpzQ5E/3mbVbavrr/RG3K6KGv4WdAfAkBBcA13U+I1uz1x3Woq3HzNpDbepqct9G3G8IQLkiuAC4ptMpmfrDG98o5X93dZakv9/dVA+3jeCrIQDljuAC4Kp+Opuh7q9tNMdVbD6a/8fbdTtXwQVgEYILgEJ9EHtUL3xx5a7OfZqGavaQ21SporeFXQHwdAQXAL/z1NJd+nz3aXP85y63aHKfRhZ2BACXEVwAmI6dv6g/Ltiun89eNGufP3mnWoRXs64pAPgVggsASZcvKtf51Q3mOLiKrzZP7sZXQwCcSgWrGwDgHN7Z+JP5/N6WtfXNJEILAOfDHhcAijuWrPc2x0uSagXYNOP+5vLx5t81AJwPP5kAD5eVm6/75l65hP/kPo0ILQCcFj+dAA/34pdXTnn+U6ebNeA2u4XdAMC18VUR4MGOn7+kJdtOSJKC/H313B8aW9wRAFwbe1wAD/bq2kPm8+V/bm9hJwBQNAQXwINdys6TJHVpWFM3BXOHZwDOj6+KADdz4sIlDZm3Vecysq87Nye/QJL04B11y7otACgVBBfAjRiGoV6zNikzN7/I7/Hz9Vaz2oFl2BUAlB6CC+BGDp1JN0PL4NbhGtu9/nXfU93PV/42fhQAcA38tALcQG5+geZ8fURzN1y5+u30e5upQgUvC7sCgNJHcAFcXPy5i5r48R7FHUs2az2bhBBaALglggvgwr77+bwGz9vqUHv+D401sGVtizoCgLJFcAFcUH6BoY93nNCzy/eZtaiI6ppxf3PdUrOKhZ0BQNkiuAAuxjAMPbd8n5btOGHWHr2znv5yV2O+HgLg9ggugAsxDENjl+zSf/YmmLV3Ho5Sn8hQC7sCgPJDcAFcxM7jyXrk/W1Kz8oza99M6qrwID8LuwKA8kVwAVzAK6sPOpzqHFHDT2uiO6lSRW8LuwKA8kdwAZzYmbQsPb9in9b9kGTW/nJXY41of5N8vLnVGADPQ3ABnNT2oxf0wDtbHGoxT3dSg5AAizoCAOsRXAAnYxiGlmw7oedWXDnVuXFYVc0bFsXxLAA8HsEFcCJpWbl6+auD+ui742btj3fepAm9GqoK9xMCAIIL4CxOJl/SiPnbdSQpw6z9455IPdQmwsKuAMC5EFwAJ3AoMV29Z29yqC0f3V6t6la3qCMAcE4EF8Bi8zb9pGmrDprjLg1rauag2xTk72thVwDgnAgugIX+/GGcvtqfaI7HdW+g8T1vtbAjAHBuxb4QxKZNm9S/f3/Z7XZ5eXnps88+c3jdMAxNnTpVdrtdlStXVpcuXXTgwAGHOdnZ2Ro7dqyCg4Pl7++vAQMG6OTJkze0EMCVHEnKUOdX1zuElpXjOhBaAOA6ih1cLl68qBYtWmjOnDmFvj5jxgzNnDlTc+bM0fbt2xUaGqqePXsqPT3dnBMdHa0VK1Zo6dKl2rx5szIyMtSvXz/l5+eXfCWAi/jP3tPqMXOjjp2/JEkKrVpJh17qo6b2QIs7AwDn52UYhlHiN3t5acWKFRo4cKCky3tb7Ha7oqOjNXnyZEmX966EhITolVde0ahRo5SamqqaNWtq0aJFGjx4sCTp9OnTCg8P16pVq9S7d+/r/rlpaWkKDAxUamqqqlatWtL2gXKVl1+g6V8d1Hub483ag3eE64X+Tbl0PwCPUBq/v0v1muHx8fFKTExUr169zJrNZlPnzp0VGxsrSYqLi1Nubq7DHLvdrsjISHPOb2VnZystLc3hAbiS0ymZuu+dLQ6h5Z8PtNBLA5sRWgCgGEo1uCQmXv6+PiQkxKEeEhJivpaYmChfX19Vr179qnN+a/r06QoMDDQf4eHhpdk2UKa2/nxe7V/+WntOpJi1JY+31f1RdeRdwcu6xgDABZXJXdq8vBx/GBuG8bvab11rzpQpU5Sammo+Tpw4UWq9AmUlJ69AX+w5rSHztpq1m2r4aeMzXdTulhoWdgYArqtUT4cODQ2VdHmvSlhYmFlPSkoy98KEhoYqJydHycnJDntdkpKS1L59+0I/12azyWazlWarQJn762f7tWzHlZD9WId6mtCroSr78tUQAJRUqe5xqVevnkJDQxUTE2PWcnJytHHjRjOUREVFqWLFig5zEhIStH///qsGF8DVrD2Q6BBaZtzfXH/p14TQAgA3qNh7XDIyMnTkyBFzHB8fr927dysoKEh169ZVdHS0pk2bpgYNGqhBgwaaNm2a/Pz8NHToUElSYGCgRo4cqQkTJqhGjRoKCgrSxIkT1axZM/Xo0aP0VgZYJD0rV39aFGeOv5nUlbs6A0ApKXZw2bFjh7p27WqOx48fL0l65JFHtGDBAk2aNEmZmZkaPXq0kpOT1aZNG61du1YBAQHme2bNmiUfHx8NGjRImZmZ6t69uxYsWCBvb/41Ctc3Yv528/mCP95OaAGAUnRD13GxCtdxgbPaFn9Bg/61RZLUs0mI3h3e2uKOAMB5ON11XABP9/Gvjmt5a2grCzsBAPdEcAFKyZd7TuvjuMv33JrUp6F8ffjfCwBKGz9ZgVJwJCldY5fsMsfNa1ezrhkAcGOleh0XwB3tP5WqH8+kX3PO+P/bYz6fcX9ztecCcwBQJgguwDUkX8zRPW9/q9z8oh3D/ucut2hQa25JAQBlheACXEVWbr52n0hRbr4hm08Ftbn52ntR6tXw08ReDcupOwDwTAQXoBAXs/PU5Z8bdDY9W5JUN8hPCx+9w+KuAAAEF+A3tvx0Xg++e+XGiP6+3hrYsraFHQEAfkFwAX7l7//5Xu9tjjfHQ24P18v3NbewIwDArxFcAF2+v9CE/9ujtd+fMWv/uCdSD7WJsLArAMBvEVzg8b4/naZ7536rrNwCs7Z+YhfVC/a3sCsAQGEILvBon+8+paeW7jbHwVVs+nzMnapdrbJ1TQEArorgAo+Umpmr9zfH6/X/HjZrfZqGatq9zRTk72thZwCAayG4wOMkpWdp0DtbdPT8JbP2/B8a69EO9eRdwcvCzgAA10NwgUc5nZKp9i9/7VD74NE71PnWmhZ1BAAoDoILPIJhGPps9yk9vezKPYU6NgjWnAdbKdCvooWdAQCKg+ACj/DYBzv034NJ5vihNnX1j3uaWdgRAKAkCC5wa0npWRoyb6t+PnvRrH04so06NAi2sCsAQEkRXOC2Vuw66fDVUO1qlbV+Yhf5+lSwsCsAwI0guMDtGIahWesO641fnercvVEtvfVQK0ILALg4ggvcyokLl/Tilwe07ocrx7PMHNRC/ZrbCS0A4AYILnAbe06k6O63vnWoLftTW7W5uYZFHQEAShvBBS4vv8DQt0fOafj728xa/VpVNGdoSzUKrWphZwCA0kZwgUvLzMnXtFU/aNHWY2bt0TvraVKfhqpU0dvCzgAAZYHgApeVlZuvbq9tUEJqlln7y12N9VjHmy3sCgBQlggucEnHzl9Uz1mblJNXIEkKC6ykd4e3VmTtQIs7AwCUJYILXM5b64/o1TWHzHHXhjX1/ojb5eXFDRIBwN0RXOAy8gsMjVu6Syv3Jpi1Z/s20hOdb7GwKwBAeSK4wCXsO5mqJz6M06mUTLP21VMd1TiMs4YAwJMQXOD0jp2/qP5zNptjm08FbXimi8ICK1vYFQDACgQXOLW8/AL1mLnRHN99m11/uauJagbYLOwKAGAVgguc2sp9CcrNNyRJHRsEa+ag2+RdgYNwAcBTcfMWOK3c/AJNW/WDOX53eGtCCwB4OIILnFb00t06k5Yt6fLVcLkSLgCA4AKntPHHs1q578ppz0904Wq4AACOcYETysjO0yO/umHinr/1UqBfRQs7AgA4C/a4wOnMjvnRfP7u8NaEFgCAiT0uKDd7T6ZoxupDysrNv+a8HceSJUkRNfzUs0lIebQGAHARBBeUi6zcfA2Y822x3vPW0FZl1A0AwFURXFAunly803x+T8va6t302ntSwgIrc6dnAMDvEFxQplIu5ejtDT/pvweTJElVbD765wMtuB4LAKBECC4oM8fOX9Tgf21VYlqWWVsd3ZHQAgAoMYILysSPZ9LVa9Ymc+xTwUsLR96hOtX9LOwKAODqCC4oVYZhaPORcxr23pXrsPRoHKLXh9wmfxt/3QAAN4bfJChVw9/fpm8OnzPHU/o20qjOt1jYEQDAnZT6BeimTp0qLy8vh0doaKj5umEYmjp1qux2uypXrqwuXbrowIEDpd0GylliapY6zvjaIbS88WBLQgsAoFSVyZVzmzZtqoSEBPOxb98+87UZM2Zo5syZmjNnjrZv367Q0FD17NlT6enpZdEKysEncSfVdvp/deJCpiTpphp++vGlvhrQwm5xZwAAd1MmXxX5+Pg47GX5hWEYmj17tp5//nnde++9kqQPPvhAISEh+uijjzRq1KiyaAdlxDAM/XPtIb21/iez1rtpiN54sKV8fbibBACg9JXJb5fDhw/LbrerXr16GjJkiH7++WdJUnx8vBITE9WrVy9zrs1mU+fOnRUbG3vVz8vOzlZaWprDA9Y6dv6ihr+/zSG0vPFgS80Z2ko2H28LOwMAuLNSDy5t2rTRwoULtWbNGr377rtKTExU+/btdf78eSUmJkqSQkIcr5oaEhJivlaY6dOnKzAw0HyEh4eXdtsohl3Hk9X51Q0Ox7N88kQ7DWhhV0Vv9rQAAMpOqX9V1LdvX/N5s2bN1K5dO91yyy364IMP1LZtW0mSl5fjBcgMw/hd7demTJmi8ePHm+O0tDTCiwXy8gv0zeFz+uOC7WatSVhVvT7kNjUICbCwMwCApyjz06H9/f3VrFkzHT58WAMHDpQkJSYmKiwszJyTlJT0u70wv2az2WSz2cq6VVxDTl6B/vb5fi3dfsKsjep0syb0asjxLACAclPmv3Gys7P1ww8/KCwsTPXq1VNoaKhiYmLM13NycrRx40a1b9++rFtBCWXl5uu+ubEOoWVq/yaa8ofGhBYAQLkq9T0uEydOVP/+/VW3bl0lJSXppZdeUlpamh555BF5eXkpOjpa06ZNU4MGDdSgQQNNmzZNfn5+Gjp0aGm3glJwJCld/d/8Vpm5+ZIkm08FfTGmgxqG8tUQAKD8lXpwOXnypB588EGdO3dONWvWVNu2bbV161ZFRERIkiZNmqTMzEyNHj1aycnJatOmjdauXauAAH4ROpsl245ryvIr1+BpUSdQn/y5PQfgAgAs42UYhmF1E8WVlpamwMBApaamqmrVqla343Zy8ws05qOdWnPgjFn7W78merRDPQu7AgC4utL4/c29iuBg1/FkPb4wTucyss3a0j+1Vduba1jYFQAAlxFcYPpiz2mNW7LLHAfYfLTqqY4KD/KzsCsAAK4guEAZ2XlavPWYpn910KwNal1Hk/s0Uo0qnIYOAHAeBBcPdyYtS3+cv13fJ1y5jcKM+5vr/lZ1VKHC1S8KCACAFQguHuznsxnq9tpGh9p7j7RW98ZXvxggAABWIrh4qDUHEjVqUZw5bmqvqneHt5a9WmULuwIA4NoILh5oVsyPev2/h83xvS1r67VBLa55vygAAJwBwcXD7DuZ6hBa3nywpfq3sFvYEQAARUdw8SDx5y6q/5zN5jj22W58NQQAcClcu91DFBQY6v7aBnP8/+5uSmgBALgcgouHWPzdMRX87+YOQ24P17C2EdY2BABACRBcPEBWbr7++vkBc/z3gZEciAsAcEkEFw/wcdxJ8/m/h7fm7s4AAJfFbzAPsPNYsiSpgpfUowkXlwMAuC7OKnJRBQWG4s9fVMEvB65cQ2pmriRpeLubyrgrAADKFsHFRU36dK8++dVXQEXRIjywjLoBAKB8EFxcUE5egUNoCfL3ve57agXYdGf94LJsCwCAMkdwcUHHL1wyn++b2ksBlSpa2A0AAOWH4OJC1hxI1Nvrjyg9O0+SVLtaZUILAMCjEFxcxDsbf9LLXx10qDUMDbCoGwAArEFwcXI/n83Qm18f0Ypdp8zac39opEahVRUVUd3CzgAAKH8EFycWd+yC7pu7xaH24cg26tCAg2wBAJ6J4OKEcvMLtOHQWT2+cIdZa14nUH+/O1ItwqtZ1xgAABYjuDiZ7Lx8TVm+T8t3Xvlq6M9dbtGEnrfKh0v1AwA8HMHFiWTl5mvgW9/qYGK6Wfv7wEju5AwAwP8QXJzEDwlpuuftb5WVWyBJCqjkoxWj71T9WlUs7gwAAOdBcHECr645qLfW/2SOoyKq6/9GtZN3BS8LuwIAwPkQXCyUk1eg0Yt3at0PZ8waXw0BAHB1BBeL7Dh6QY8u2K60rDyztmFiF90U7G9hVwAAODeCiwVW7Dqpp5ftMcc1/H31VXRH1QqoZGFXAAA4P4JLOUrLytUH3x7VazE/mrVhbSM0vuetql6EOzwDAODpCC7lJCE1Uw/9+zv9fPaiWXvtgRa6t1VteXlxEC4AAEVBcCkHR5Iy1GPmRofa/D/erq4Na1nUEQAArongUoYMw9Dnu08retlus9aqbjW983CUalXleBYAAIqL4FJGDMPQuKW79eWe02ZtcOtwvXxfM74aAgCghAguZSAjO0/939ys+HNXjmd55+Eo9YkMtbArAABcH8GllH21L0F/XrzTHIcHVVbM051VqaK3hV0BAOAeCC6laPa6HzV73WFz3L1RLf1rWBR3dQYAoJQQXEpBYmqW/rHqB4fjWWbc31z3tarD/YYAAChFBJcblJ2Xr54zNyo9+8ql+z978k7dFl7NuqYAAHBTBJcbcOLCJT3zyR4ztNQN8tO7w1urYWiAxZ0BAOCeCC4llJtfoI4z1jvUVo7roIBKFS3qCAAA98dRoyU0fdVB83nLutW0bnxnQgsAAGWMPS4lkJSWpfe/jZck+ft6a/mf23NROQAAygF7XErgo23Hzedrnu5EaAEAoJxYGlzefvtt1atXT5UqVVJUVJS++eYbK9spspRLuZKkhiEBqlPdz+JuAADwHJYFl2XLlik6OlrPP/+8du3apY4dO6pv3746fvz49d9skRMXLmnou1v12e5TkqRBt4db3BEAAJ7FsuAyc+ZMjRw5Uo899pgaN26s2bNnKzw8XHPnzrWqpWsyDEO9Z29S7E/nHfa4AACA8mPJwbk5OTmKi4vTs88+61Dv1auXYmNjfzc/Oztb2dnZ5jgtLa1M+jqXka231h8p9LUDp9N0KSdfktTu5hp6/q7GamqvWiZ9AACAwlkSXM6dO6f8/HyFhIQ41ENCQpSYmPi7+dOnT9eLL75Y5n2lZeZq/rdHrztv/h9v56aJAABYwNLToX97No5hGIWeoTNlyhSNHz/eHKelpSk8vPSPL6nm56snu95y1de9vbzUv4Wd0AIAgEUsCS7BwcHy9vb+3d6VpKSk3+2FkSSbzSabzVbmfQX5++qZ3o3K/M8BAAAlY8nBub6+voqKilJMTIxDPSYmRu3bt7eiJQAA4AIs+6po/PjxGjZsmFq3bq127dpp3rx5On78uJ544gmrWgIAAE7OsuAyePBgnT9/Xv/v//0/JSQkKDIyUqtWrVJERIRVLQEAACfnZRiGYXUTxZWWlqbAwEClpqaqalVOSQYAwBWUxu9v7lUEAABcBsEFAAC4DIILAABwGQQXAADgMgguAADAZRBcAACAyyC4AAAAl0FwAQAALoPgAgAAXIZll/y/Eb9c7DctLc3iTgAAQFH98nv7Ri7a75LBJT09XZIUHh5ucScAAKC40tPTFRgYWKL3uuS9igoKCnT69GkFBATIy8vL6naKLS0tTeHh4Tpx4oRH3WuJdXvOuj1xzRLr9qR1e+KapRtft2EYSk9Pl91uV4UKJTtaxSX3uFSoUEF16tSxuo0bVrVqVY/6C/8L1u05PHHNEuv2JJ64ZunG1l3SPS2/4OBcAADgMgguAADAZRBcLGCz2fTCCy/IZrNZ3Uq5Yt2es25PXLPEuj1p3Z64Zsk51u2SB+cCAADPxB4XAADgMgguAADAZRBcAACAyyC4AAAAl0FwKaJTp07p4YcfVo0aNeTn56fbbrtNcXFx5uvLly9X7969FRwcLC8vL+3evft3n5Gdna2xY8cqODhY/v7+GjBggE6ePOkwJzk5WcOGDVNgYKACAwM1bNgwpaSkOMw5fvy4+vfvL39/fwUHB2vcuHHKyckpi2WXyrq7dOkiLy8vh8eQIUMc5jjTuq+15tzcXE2ePFnNmjWTv7+/7Ha7hg8frtOnTzt8hrtt66Ku29W29fXWLUlTp05Vo0aN5O/vr+rVq6tHjx767rvvHD7D3bZ3Udftatv7emv+tVGjRsnLy0uzZ892qLvjtv61q63bqba1geu6cOGCERERYYwYMcL47rvvjPj4eGPdunXGkSNHzDkLFy40XnzxRePdd981JBm7du363ec88cQTRu3atY2YmBhj586dRteuXY0WLVoYeXl55pw+ffoYkZGRRmxsrBEbG2tERkYa/fr1M1/Py8szIiMjja5duxo7d+40YmJiDLvdbowZM8Zp1925c2fj8ccfNxISEsxHSkqKwxxnWff11pySkmL06NHDWLZsmXHw4EFjy5YtRps2bYyoqCiHz3G3bV3UdbvSti7Kug3DMBYvXmzExMQYP/30k7F//35j5MiRRtWqVY2kpCRzjrtt76Ku25W2d1HW/IsVK1YYLVq0MOx2uzFr1iyH19xxWxdl3c60rQkuRTB58mSjQ4cORZobHx9f6C/wlJQUo2LFisbSpUvN2qlTp4wKFSoYq1evNgzDML7//ntDkrF161ZzzpYtWwxJxsGDBw3DMIxVq1YZFSpUME6dOmXOWbJkiWGz2YzU1NSSLrFQpbFuw7j8F/6pp5666nudad3FWfMvtm3bZkgyjh07ZhiG+2/rX/x23YbhWtvaMEq27tTUVEOSsW7dOsMwPGd7/3bdhuFa27uoaz558qRRu3ZtY//+/UZERITDL3B33tbXWrdhONe25quiIvjiiy/UunVrPfDAA6pVq5Zatmypd999t1ifERcXp9zcXPXq1cus2e12RUZGKjY2VpK0ZcsWBQYGqk2bNuactm3bKjAw0GFOZGSk7Ha7Oad3797Kzs6+6q6/kiqNdf9i8eLFCg4OVtOmTTVx4kTzDt+Sc627JGtOTU2Vl5eXqlWrJslztvVv1/0LV9nWUvHXnZOTo3nz5ikwMFAtWrSQ5Bnbu7B1/8JVtndR1lxQUKBhw4bpmWeeUdOmTX/3Ge66ra+37l84y7YmuBTBzz//rLlz56pBgwZas2aNnnjiCY0bN04LFy4s8mckJibK19dX1atXd6iHhIQoMTHRnFOrVq3fvbdWrVoOc0JCQhxer169unx9fc05paU01i1JDz30kJYsWaINGzbor3/9qz799FPde++95uvOtO7irjkrK0vPPvushg4dat5wzBO2dWHrllxrW0tFX/d//vMfValSRZUqVdKsWbMUExOj4OBgs1933d7XWrfkWtu7KGt+5ZVX5OPjo3HjxhX6Ge66ra+3bsm5trVL3h26vBUUFKh169aaNm2aJKlly5Y6cOCA5s6dq+HDh9/QZxuGIS8vL3P86+c3Mqc0lNa6H3/8cfN5ZGSkGjRooNatW2vnzp1q1aqVJOdZd3HWnJubqyFDhqigoEBvv/32dT/bXbb1tdbtSttaKvq6u3btqt27d+vcuXN69913NWjQIH333XeF/qC+Wr/uuG5X2t7XW3NcXJxef/117dy5s9h/ritv66Ku25m2NXtciiAsLExNmjRxqDVu3FjHjx8v8meEhoYqJydHycnJDvWkpCQzgYaGhurMmTO/e+/Zs2cd5vw2mSYnJys3N/d3SfZGlca6C9OqVStVrFhRhw8fluRc6y7qmnNzczVo0CDFx8crJibGYa+DO2/ra627MM68raWir9vf31/169dX27Zt9d5778nHx0fvvfee2a+7bu9rrbswzry9r7fmb775RklJSapbt658fHzk4+OjY8eOacKECbrpppvMXt1tWxdl3YWxclsTXIrgzjvv1KFDhxxqP/74oyIiIor8GVFRUapYsaJiYmLMWkJCgvbv36/27dtLktq1a6fU1FRt27bNnPPdd98pNTXVYc7+/fuVkJBgzlm7dq1sNpuioqJKtL6rKY11F+bAgQPKzc1VWFiYJOdad1HW/Msv78OHD2vdunWqUaOGw3x33dbXW3dhnHlbSyX/O24YhrKzsyW57/YuzK/XXRhn3t7XW/OwYcO0d+9e7d6923zY7XY988wzWrNmjST33NZFWXdhLN3WRT6M14Nt27bN8PHxMf7xj38Yhw8fNhYvXmz4+fkZH374oTnn/Pnzxq5du4yVK1cakoylS5cau3btMhISEsw5TzzxhFGnTh1j3bp1xs6dO41u3boVehpd8+bNjS1bthhbtmwxmjVrVujpZN27dzd27txprFu3zqhTp06ZnEZXGus+cuSI8eKLLxrbt2834uPjjZUrVxqNGjUyWrZs6ZTrvt6ac3NzjQEDBhh16tQxdu/e7XBqYHZ2tvk57rati7JuV9vWRVl3RkaGMWXKFGPLli3G0aNHjbi4OGPkyJGGzWYz9u/fb36Ou23voqzb1bZ3UX6e/VZhZ9e427YuyrqdbVsTXIroyy+/NCIjIw2bzWY0atTImDdvnsPr8+fPNyT97vHCCy+YczIzM40xY8YYQUFBRuXKlY1+/foZx48fd/ic8+fPGw899JAREBBgBAQEGA899JCRnJzsMOfYsWPGXXfdZVSuXNkICgoyxowZY2RlZTnluo8fP2506tTJCAoKMnx9fY1bbrnFGDdunHH+/HmnXfe11vzLad+FPdavX2/Oc7dtXZR1u+K2vt66MzMzjXvuucew2+2Gr6+vERYWZgwYMMDYtm2bw2e42/YuyrpdcXtf7+fZbxUWXNxtWxfmt+t2tm3tZRiGUfT9MwAAANbhGBcAAOAyCC4AAMBlEFwAAIDLILgAAACXQXABAAAug+ACAABcBsEFAAC4DIILAABwGQQXAADgMgguAADAZRBcAACAyyC4AAAAl/H/Ab+cQR1tJxpaAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ddf_name = 'ECDFS'\n",
    "mjds, night_mjd, cumulative_desired, cumulative_sched = optimize_ddf_times(\n",
    "            ddf_name,\n",
    "            ddfs[ddf_name][0],\n",
    "            ddf_grid,\n",
    "            season_unobs_frac=0.1,\n",
    "            sequence_time=60,\n",
    "        )\n",
    "plt.plot(night_mjd, cumulative_sched)\n",
    "print(cumulative_sched[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f813553b-e472-402d-b5e7-0554734fd61d",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}