fps_res.ipynb

{
 "cells": [
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "from fps_io import parse_fps_res, res2df",
   "execution_count": 1,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import marshaltools\nfrom marshaltools import ProgramList\nimport numpy as np\nimport pandas as pd\nimport sys",
   "execution_count": 2,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "from ztfquery import query",
   "execution_count": 3,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import seaborn as sns\nimport matplotlib.pyplot as plt\nsns.set_style('whitegrid')\nsns.set_context('notebook')",
   "execution_count": 4,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# `fps_second.txt` is the downloaded file I had from the forced photometry service \n# \n# Requested input R.A. = 153.1786446 degrees\n# Requested input Dec. = 39.373955 degrees\n# Requested JD start = 2458484.9111 days\n# Requested JD end = 2458584.9111 days\n# Authenticated program IDs = 2\ndata, names, rest = parse_fps_res('fps_second.txt')\ndf, requests = res2df(data, names, rest)",
   "execution_count": 5,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# requests is a series with ra, dec and jd start and end days\nrequests",
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 6,
     "data": {
      "text/plain": "ra         1.531786e+02\ndec        3.937396e+01\njdstart    2.458485e+06\njdend      2.458585e+06\ndtype: float64"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# df holds the obtained forced photometry results\ndf.head()",
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 7,
     "data": {
      "text/plain": "   index  field  ccdid  qid filter           pid  infobitssci  airmass  \\\n0      0   1710     15    2  ZTF_r  738456425815            0    1.002   \n1      1   1710     15    2  ZTF_r  739448415815            0    1.002   \n2      2   1710     15    2  ZTF_r  749493945815            1    1.074   \n3      3   1710     15    2  ZTF_r  755471655815            0    1.062   \n4      4    712      8    2  ZTF_i  784379673015            0    1.033   \n\n   moonalt  moonillf     ...      forcediffimfluxap  forcediffimfluxuncap  \\\n0 -71.1969  0.100470     ...              58.440480             69.554406   \n1 -68.2965  0.163419     ...             -51.131818             49.936020   \n2  27.6809  0.992139     ...            -227.959783            272.356399   \n3  48.7452 -0.653691     ...            1516.325832             85.909244   \n4  28.0442 -0.716601     ...             253.563648             97.903252   \n\n   forcediffimsnrap  aperturecorr  dnearestrefsrc  nearestrefmag  \\\n0          0.845871      1.048457        0.454230         16.460   \n1         -1.030843      1.037749        0.454230         16.460   \n2         -0.842628      1.061419        0.454230         16.460   \n3         17.769198      1.055294        0.454230         16.460   \n4          2.607384      1.100413        0.568666         16.148   \n\n   nearestrefmagunc  nearestrefchi  nearestrefsharp  procstatus  \n0             0.047          4.923            0.211           0  \n1             0.047          4.923            0.211           0  \n2             0.047          4.923            0.211           0  \n3             0.047          4.923            0.211           0  \n4             0.059          5.164            0.424           0  \n\n[5 rows x 53 columns]",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>index</th>\n      <th>field</th>\n      <th>ccdid</th>\n      <th>qid</th>\n      <th>filter</th>\n      <th>pid</th>\n      <th>infobitssci</th>\n      <th>airmass</th>\n      <th>moonalt</th>\n      <th>moonillf</th>\n      <th>...</th>\n      <th>forcediffimfluxap</th>\n      <th>forcediffimfluxuncap</th>\n      <th>forcediffimsnrap</th>\n      <th>aperturecorr</th>\n      <th>dnearestrefsrc</th>\n      <th>nearestrefmag</th>\n      <th>nearestrefmagunc</th>\n      <th>nearestrefchi</th>\n      <th>nearestrefsharp</th>\n      <th>procstatus</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>0</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>2</td>\n      <td>ZTF_r</td>\n      <td>738456425815</td>\n      <td>0</td>\n      <td>1.002</td>\n      <td>-71.1969</td>\n      <td>0.100470</td>\n      <td>...</td>\n      <td>58.440480</td>\n      <td>69.554406</td>\n      <td>0.845871</td>\n      <td>1.048457</td>\n      <td>0.454230</td>\n      <td>16.460</td>\n      <td>0.047</td>\n      <td>4.923</td>\n      <td>0.211</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>1</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>2</td>\n      <td>ZTF_r</td>\n      <td>739448415815</td>\n      <td>0</td>\n      <td>1.002</td>\n      <td>-68.2965</td>\n      <td>0.163419</td>\n      <td>...</td>\n      <td>-51.131818</td>\n      <td>49.936020</td>\n      <td>-1.030843</td>\n      <td>1.037749</td>\n      <td>0.454230</td>\n      <td>16.460</td>\n      <td>0.047</td>\n      <td>4.923</td>\n      <td>0.211</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>2</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>2</td>\n      <td>ZTF_r</td>\n      <td>749493945815</td>\n      <td>1</td>\n      <td>1.074</td>\n      <td>27.6809</td>\n      <td>0.992139</td>\n      <td>...</td>\n      <td>-227.959783</td>\n      <td>272.356399</td>\n      <td>-0.842628</td>\n      <td>1.061419</td>\n      <td>0.454230</td>\n      <td>16.460</td>\n      <td>0.047</td>\n      <td>4.923</td>\n      <td>0.211</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>3</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>2</td>\n      <td>ZTF_r</td>\n      <td>755471655815</td>\n      <td>0</td>\n      <td>1.062</td>\n      <td>48.7452</td>\n      <td>-0.653691</td>\n      <td>...</td>\n      <td>1516.325832</td>\n      <td>85.909244</td>\n      <td>17.769198</td>\n      <td>1.055294</td>\n      <td>0.454230</td>\n      <td>16.460</td>\n      <td>0.047</td>\n      <td>4.923</td>\n      <td>0.211</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>4</td>\n      <td>712</td>\n      <td>8</td>\n      <td>2</td>\n      <td>ZTF_i</td>\n      <td>784379673015</td>\n      <td>0</td>\n      <td>1.033</td>\n      <td>28.0442</td>\n      <td>-0.716601</td>\n      <td>...</td>\n      <td>253.563648</td>\n      <td>97.903252</td>\n      <td>2.607384</td>\n      <td>1.100413</td>\n      <td>0.568666</td>\n      <td>16.148</td>\n      <td>0.059</td>\n      <td>5.164</td>\n      <td>0.424</td>\n      <td>0</td>\n    </tr>\n  </tbody>\n</table>\n<p>5 rows × 53 columns</p>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df['SNR'] = df.forcediffimflux / df.forcediffimfluxunc",
   "execution_count": 8,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df['factor'] = df['SNR'] / df['forcediffimsnr']",
   "execution_count": 9,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.factor",
   "execution_count": 10,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 10,
     "data": {
      "text/plain": "0     0.89897\n1     0.89897\n2     0.89897\n3     0.89897\n4     0.89897\n5     0.89897\n6     0.89897\n7     0.89897\n8     0.89897\n9     0.89897\n10    0.89897\nName: factor, dtype: float64"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## What epochs do we have?"
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.groupby('filter').agg(dict(jd='count'))",
   "execution_count": 11,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 11,
     "data": {
      "text/plain": "        jd\nfilter    \nZTF_i    7\nZTF_r    4",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>jd</th>\n    </tr>\n    <tr>\n      <th>filter</th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>ZTF_i</th>\n      <td>7</td>\n    </tr>\n    <tr>\n      <th>ZTF_r</th>\n      <td>4</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df[['zpdiff', 'zpref', 'zpmaginpsci']]",
   "execution_count": 12,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 12,
     "data": {
      "text/plain": "     zpdiff   zpref  zpmaginpsci\n0   25.1376  26.275      25.1376\n1   26.2237  26.275      26.2237\n2   25.3055  26.275      25.3055\n3   26.1888  26.275      26.1888\n4   24.6872  25.660      24.6872\n5   24.4811  25.660      24.4811\n6   24.6440  25.660      24.6440\n7   24.3433  25.660      24.3433\n8   25.6178  25.660      25.6178\n9   25.6352  25.660      25.6352\n10  25.5981  25.660      25.5981",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>zpdiff</th>\n      <th>zpref</th>\n      <th>zpmaginpsci</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>25.1376</td>\n      <td>26.275</td>\n      <td>25.1376</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>26.2237</td>\n      <td>26.275</td>\n      <td>26.2237</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>25.3055</td>\n      <td>26.275</td>\n      <td>25.3055</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>26.1888</td>\n      <td>26.275</td>\n      <td>26.1888</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>24.6872</td>\n      <td>25.660</td>\n      <td>24.6872</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>24.4811</td>\n      <td>25.660</td>\n      <td>24.4811</td>\n    </tr>\n    <tr>\n      <th>6</th>\n      <td>24.6440</td>\n      <td>25.660</td>\n      <td>24.6440</td>\n    </tr>\n    <tr>\n      <th>7</th>\n      <td>24.3433</td>\n      <td>25.660</td>\n      <td>24.3433</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>25.6178</td>\n      <td>25.660</td>\n      <td>25.6178</td>\n    </tr>\n    <tr>\n      <th>9</th>\n      <td>25.6352</td>\n      <td>25.660</td>\n      <td>25.6352</td>\n    </tr>\n    <tr>\n      <th>10</th>\n      <td>25.5981</td>\n      <td>25.660</td>\n      <td>25.5981</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.columns",
   "execution_count": 13,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 13,
     "data": {
      "text/plain": "Index(['index', 'field', 'ccdid', 'qid', 'filter', 'pid', 'infobitssci',\n       'airmass', 'moonalt', 'moonillf', 'moonra', 'moondec', 'sciinpseeing',\n       'scibckgnd', 'scisigpix', 'scimaglim', 'zpmaginpsci', 'zpmaginpsciunc',\n       'zpmaginpscirms', 'clrcoeff', 'clrcoeffunc', 'ncalmatches', 'exptime',\n       'diffmaglim', 'zpdiff', 'programid', 'obsdate', 'jd', 'scifilename',\n       'diffilename', 'rfid', 'refmaglim', 'refbckgnd', 'refsigpix', 'zpref',\n       'refcreated', 'refjdstart', 'refjdend', 'reffilename',\n       'forcediffimflux', 'forcediffimfluxunc', 'forcediffimsnr',\n       'forcediffimchisq', 'forcediffimfluxap', 'forcediffimfluxuncap',\n       'forcediffimsnrap', 'aperturecorr', 'dnearestrefsrc', 'nearestrefmag',\n       'nearestrefmagunc', 'nearestrefchi', 'nearestrefsharp', 'procstatus',\n       'SNR', 'factor'],\n      dtype='object')"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# What we expect using ztfquery"
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "zquery = query.ZTFQuery()\nzquery.load_metadata(radec=[requests.ra, requests.dec], size=0.01,\n                     sql_query=\"obsjd>{0} and obsjd<{1}\".format(requests['jdstart'], requests['jdend']))",
   "execution_count": 14,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "expected =  zquery.metatable",
   "execution_count": 15,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "expected.groupby('filtercode').agg(dict(obsjd='count'))",
   "execution_count": 16,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 16,
     "data": {
      "text/plain": "            obsjd\nfiltercode       \nzg              4\nzi              9\nzr              4",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>obsjd</th>\n    </tr>\n    <tr>\n      <th>filtercode</th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>zg</th>\n      <td>4</td>\n    </tr>\n    <tr>\n      <th>zi</th>\n      <td>9</td>\n    </tr>\n    <tr>\n      <th>zr</th>\n      <td>4</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "So, we are missing the 'g' band epochs in the forced photometry. We are also missing 2 i band epochs. "
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "missing_g = expected.query('filtercode == \"zg\"')[['obsjd', 'field', 'ccdid', 'qid']]",
   "execution_count": 17,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "from astropy.time import Time",
   "execution_count": 18,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "missing_g.iloc[0]['obsjd']",
   "execution_count": 19,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 19,
     "data": {
      "text/plain": "2458508.0121065"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "missing_g['night'] = list(Time(t, format='jd').strftime('%Y-%m-%d') for t in missing_g['obsjd'].values)",
   "execution_count": 20,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "missing_g",
   "execution_count": 21,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 21,
     "data": {
      "text/plain": "           obsjd  field  ccdid  qid       night\n6   2.458508e+06   1710     15    3  2019-01-24\n8   2.458523e+06   1710     15    3  2019-02-08\n10  2.458507e+06   1710     15    3  2019-01-23\n12  2.458488e+06   1710     15    3  2019-01-04",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>obsjd</th>\n      <th>field</th>\n      <th>ccdid</th>\n      <th>qid</th>\n      <th>night</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>6</th>\n      <td>2.458508e+06</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>3</td>\n      <td>2019-01-24</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>2.458523e+06</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>3</td>\n      <td>2019-02-08</td>\n    </tr>\n    <tr>\n      <th>10</th>\n      <td>2.458507e+06</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>3</td>\n      <td>2019-01-23</td>\n    </tr>\n    <tr>\n      <th>12</th>\n      <td>2.458488e+06</td>\n      <td>1710</td>\n      <td>15</td>\n      <td>3</td>\n      <td>2019-01-04</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "y = expected.query('filtercode == \"zi\"')[['obsjd', 'field', 'ccdid', 'qid']].sort_values(by='obsjd')",
   "execution_count": 22,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "y.loc[[13, 5]].obsjd.values",
   "execution_count": 23,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 23,
     "data": {
      "text/plain": "array([2458538.8796759, 2458538.8801389])"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "x = df.query('filter == \"ZTF_i\"')[['jd', 'field', 'ccdid', 'qid']].sort_values(by='jd').rename(columns=dict(jd='obsjd'))",
   "execution_count": 24,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "x",
   "execution_count": 25,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 25,
     "data": {
      "text/plain": "           obsjd  field  ccdid  qid\n4   2.458539e+06    712      8    2\n5   2.458539e+06    711      5    3\n6   2.458544e+06    712      8    2\n7   2.458544e+06    711      5    3\n8   2.458548e+06    712      8    2\n9   2.458558e+06    712      8    2\n10  2.458562e+06    712      8    2",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>obsjd</th>\n      <th>field</th>\n      <th>ccdid</th>\n      <th>qid</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>4</th>\n      <td>2.458539e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>2.458539e+06</td>\n      <td>711</td>\n      <td>5</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>6</th>\n      <td>2.458544e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>7</th>\n      <td>2.458544e+06</td>\n      <td>711</td>\n      <td>5</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>2.458548e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>9</th>\n      <td>2.458558e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>10</th>\n      <td>2.458562e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>2</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "y.sort_values(by='obsjd')",
   "execution_count": 26,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 26,
     "data": {
      "text/plain": "           obsjd  field  ccdid  qid\n13  2.458539e+06    712      8    3\n5   2.458539e+06    711      5    4\n4   2.458544e+06    712      8    3\n3   2.458544e+06    711      5    4\n14  2.458548e+06    712      8    3\n9   2.458558e+06    712      8    3\n11  2.458562e+06    712      8    3\n0   2.458568e+06    712      8    3\n1   2.458568e+06    711      5    4",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>obsjd</th>\n      <th>field</th>\n      <th>ccdid</th>\n      <th>qid</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>13</th>\n      <td>2.458539e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>2.458539e+06</td>\n      <td>711</td>\n      <td>5</td>\n      <td>4</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>2.458544e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>2.458544e+06</td>\n      <td>711</td>\n      <td>5</td>\n      <td>4</td>\n    </tr>\n    <tr>\n      <th>14</th>\n      <td>2.458548e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>9</th>\n      <td>2.458558e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>11</th>\n      <td>2.458562e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>0</th>\n      <td>2.458568e+06</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>2.458568e+06</td>\n      <td>711</td>\n      <td>5</td>\n      <td>4</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "len(x) , len(y)",
   "execution_count": 27,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 27,
     "data": {
      "text/plain": "(7, 9)"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "expected.loc[[0, 1]]",
   "execution_count": 28,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 28,
     "data": {
      "text/plain": "   in_row_id       in_ra     in_dec          ra        dec  infobits  field  \\\n0          1  153.178645  39.373955  153.643493  39.073826         0    712   \n1          1  153.178645  39.373955  153.316472  39.081491         0    711   \n\n   ccdid  qid  rcid    ...            ra1       dec1         ra2       dec2  \\\n0      8    3    30    ...     154.175286  39.526503  153.056876  39.483411   \n1      5    4    19    ...     153.901112  39.493095  152.782302  39.532583   \n\n          ra3       dec3         ra4       dec4 ipac_pub_date  ipac_gid  \n0  153.119030  38.618998  154.223385  38.661275           NaN         2  \n1  152.738611  38.667051  153.843616  38.628239           NaN         2  \n\n[2 rows x 45 columns]",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>in_row_id</th>\n      <th>in_ra</th>\n      <th>in_dec</th>\n      <th>ra</th>\n      <th>dec</th>\n      <th>infobits</th>\n      <th>field</th>\n      <th>ccdid</th>\n      <th>qid</th>\n      <th>rcid</th>\n      <th>...</th>\n      <th>ra1</th>\n      <th>dec1</th>\n      <th>ra2</th>\n      <th>dec2</th>\n      <th>ra3</th>\n      <th>dec3</th>\n      <th>ra4</th>\n      <th>dec4</th>\n      <th>ipac_pub_date</th>\n      <th>ipac_gid</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>1</td>\n      <td>153.178645</td>\n      <td>39.373955</td>\n      <td>153.643493</td>\n      <td>39.073826</td>\n      <td>0</td>\n      <td>712</td>\n      <td>8</td>\n      <td>3</td>\n      <td>30</td>\n      <td>...</td>\n      <td>154.175286</td>\n      <td>39.526503</td>\n      <td>153.056876</td>\n      <td>39.483411</td>\n      <td>153.119030</td>\n      <td>38.618998</td>\n      <td>154.223385</td>\n      <td>38.661275</td>\n      <td>NaN</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>1</td>\n      <td>153.178645</td>\n      <td>39.373955</td>\n      <td>153.316472</td>\n      <td>39.081491</td>\n      <td>0</td>\n      <td>711</td>\n      <td>5</td>\n      <td>4</td>\n      <td>19</td>\n      <td>...</td>\n      <td>153.901112</td>\n      <td>39.493095</td>\n      <td>152.782302</td>\n      <td>39.532583</td>\n      <td>152.738611</td>\n      <td>38.667051</td>\n      <td>153.843616</td>\n      <td>38.628239</td>\n      <td>NaN</td>\n      <td>2</td>\n    </tr>\n  </tbody>\n</table>\n<p>2 rows × 45 columns</p>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Is it possible that these were marginal and did not actually fall on the CCD?"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "###  Marshalltools "
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "pl = ProgramList('Cosmology')\nt = pl.table.to_pandas().set_index('name')",
   "execution_count": 29,
   "outputs": [
    {
     "output_type": "stream",
     "text": "INFO:marshaltools.ProgramList:Initialized ProgramList for program Cosmology (ID 0)\nINFO:marshaltools.ProgramList:Loaded 834 saved sources for program Cosmology.\n",
     "name": "stderr"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "t.loc['ZTF18aaermez'].ra == requests.ra",
   "execution_count": 30,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 30,
     "data": {
      "text/plain": "True"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "t.loc['ZTF18aaermez'].dec == requests.dec",
   "execution_count": 31,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 31,
     "data": {
      "text/plain": "True"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Plots"
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "grouped = expected.groupby('filtercode')\nfig, ax = plt.subplots() #1, 3, figsize=(18,6))\nfor band in grouped.groups:\n    x = grouped.get_group(band)\n    ax.scatter(x['obsjd'], x['maglimit'],\n                label=band, marker='s')\n    ax.axvline(requests.jdstart, ls='--', color='k')\n    ax.axvline(requests.jdend, ls='--', color='k')\n    ax.invert_yaxis()\nax.legend(loc='lower right', framealpha=0.4)\nax.set_ylabel('mag_lim')\nax.set_xlabel('jd_obs')",
   "execution_count": 32,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 32,
     "data": {
      "text/plain": "Text(0.5, 0, 'jd_obs')"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEPCAYAAABsj5JaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcFNW99/HPzLAMDkGMIW4RvVzlJ49eXMaFhAgmwYUnBjUuxJhL1ECuu141aFyjESKaZ6ImN5oHFeVqXK5biCjRJ254CWo6EZ9B/KGioIKGqGTCNsPS94+qMeXQPdNdM9U13fN9v16+7K6uM/U79nh+c86pOqcqm80iIiJSrOq0AxARkfKkBCIiIrEogYiISCxKICIiEosSiIiIxKIEIiIisSiBiIhILEogIiISixKIiIjEogQiIiKxKIGIiEgsvdIOoCtlMpm+wAHACmBTyuGIiJSLGmAH4KX6+vrmQgtVVAIhSB5z0w5CRKRMHQw8X+jJlZZAVrS+2GuvvYouvNtuuwHwxhtvdF1EJdDY2BirvuVMde4ZVOfCdab9amlpYfHixRBpQwtRaQnkk2Grvn37Fl34iCOOiF02beUYc2epzj2D6lyYLmq/ihr6r7QE0inTp09POwQRkVjSaL90F1bEtGnTmDZtWtphiIgULY32Sz2QiEsvvRSAiy66KOVIRESKk0b7pR6IiIjEogQiIiKxKIGIiEgsmgMR6W6m7gQtq7c4vE9NP6h/P4WARHJTAok4//zz0w5BJGfyAKjZtK7EgUg5SaP9UgKJuO6669IOQUQkljTaL82BREyePJnJkyenHYaISNHSaL+UQCIaGhpoaGhIOwwRkaKl0X4pgYiISCxKICLdTZ/+OQ9vqulX4kBE2qdJdJHu5pL3ch5+OZOhvsShiLRHCaQHGXH3CNZsXLPF8bpedcw/aX4KEYlIOVMCiZgyZUraISQqV/Jo77iIlI802q9EE4iZDQDmAUe6+9tmdhhwPcH+u38CJrp7S5syA4G7gSHASuAEdy/J47dahVdEylUa7Vdik+hmdhDB3rpDI4dvA77l7nsBWwETchS9Bpjr7sOA6cCNScXY1qRJk5g0aVKpLici0mXSaL+SvAtrEnAmsDxyrAYYYGY1QC2Qa22GrxP0QADuAcaaWe8E4/zEjBkzmDFjRikuJSLSpdJovxJLIO4+0d3ntjl8BvAMQVL5HPBAjqI7Em7s7u4bgSZgUFJxiohIPCWbRDez7YFrgb2At4CG8J8z25xaleP95mKvl8lkio4xm83GLpu2QmKura5l/eb1OY9Xap0rjercM5RL+1XKu7AOBhrd/U0AM5sO3J/jvPeA7YF3zawX8Bngw2IvVl9f/B3zVVVVscumKZPJFBTzS/UvlSCa0ii0zpVEde4Z4ta5M+1Xc3MzjY2NRZcr5ZPojcCBZrZd+P4oIFeL9hj/mFwfTzChvqEE8YmISBFK1gNx90VmdjnwtJltBN4Avg9gZlcDy939FuBy4A4zWwisAk4qVYw333xzqS4lItKl0mi/Ek8g7r5r5PWdwJ05zrki8vojYFzSceWiW3hFpFyl0X5pMcWI8ePHM378+LTDEBEpWhrtl5YyiXjwwQfTDkFEJJY02i/1QEREJBYlEBERiUUJREREYlECERGRWDSJHvHQQw+lHYKISCxptF9KIBHjxqXy+ImIFGLqTtCyesvjffrn3Qa4J0mj/dIQVsTYsWMZO3Zs2mGISC65kkd7x3uYNNov9UAinnzyybRDEBGJJY32SwlEpLvJM1SzT00/qC/J7s4iBdEQlkh3k2dIpmZTrg08RdKjBCIiIrEogYhIeejTv7jjkjjNgUQ8//zzaYcgIvnoVt12pdF+KYFEjBgxIu0QRERiSaP90hBWxKhRoxg1alTaYUhPl2dIZlNNvxIHIuUkjfZLPZCIefPmpR2CSN6hmpczGepLHIqUjzTaL/VAREQkFiUQERGJJfEhLDMbAMwDjnT3t83sMOB6oAb4EzDR3VvalBkNPAS8Ex76s7ufknSsIiJSuER7IGZ2EPA8MDRy+DbgW+6+F7AVMCFH0f2Bn7r7PuE/Sh4iIt1M0j2QScCZwH9GjtUAA8ysBqgFcq3PcACwnZmdCLwNnOnu7+Q4r0stWbIk6UuIiCQijfarKpvNJn4RM3sbOCQcwjoauAdoAt4CRrt7c5vzbwGecPeHzOw04F/dfWRH18lkMruGP1NERIr3T/X19W8XenJJEwiwHngG+DpBQ98A9Hb3MzsovwrYxd3/1t550QRSX1/8DY+tZTKZTNFl05TJZGLVt5ypzj2D6ly4zrRfzc3NNDY2QpEJpNTPgRwMNLr7mwBmNh24P3qCmVUDPwSudfdNkY82Jh3cggULkr6EiEgi0mi/Sn0bbyNwoJltF74/CngpeoK7bwaOAY4FMLMJwAvuvqaUgYqISPtKmkDcfRFwOfC0mb1CcLfVhQBmdnU43wHwXeA8M1sInAJMLGWcIiLSsZIMYbn7rpHXdwJ35jjnisjrhcCXShGbiIjEoyfRRUQkFi2mGLFxY+Lz9CIiiUij/VICEakkP9q6nc/avQtepGgawooYNmwYw4YNSzsMEZGipdF+qQcS8frrr6cdgohILGm0X+qBiIhILEogIiISixKIiIjEojmQiOpq5VMRKU9ptF9KIBEtLS0dnyTSnelW3R4rjfZLf3KLiEgsSiARQ4YMYciQIWmHISJStDTaLw1hRSxbtiztEEREYkmj/VIPREREYlECERGRWJRAREQkFs2BRPTu3TvtEEREYkmj/VICiVi3bl3aIYiIxJJG+6UhLBERiSXRHoiZDQDmAUe6+9tmdjIwGdgEPAVc4O4b25QZCNwNDAFWAie4+/tJxtlqp512AuC9994rxeVERLpMGu1XYj0QMzsIeB4YGr434Brga+7+L0Bv4JwcRa8B5rr7MGA6cGNSMbb1wQcf8MEHH5TqciIiXSaN9ivJIaxJwJnA8vD9cOAP7r4ifP8ocHSOcl8n6IEA3AOMNTPNbouIdDOJDWG5+0SAoOMBwAKgwcx2JkgqxwHb5yi6I7Ai/BkbzawJGMQ/ElFBMplM0TFns9nYZdNWjjF3lurcM6jOhUmj/SrZXVjuvtjMLgZmAeuA+4EDc5xaleP95mKvV19fX3SMVVVVscumKZPJlF3MnaU69wyqc+E60341NzfT2NhYdLmSJRAzqwVedPd9w/fHA2/mOPU9gp7Ju2bWC/gM8GEpYqytrS3FZUREulwa7VcpnwOpA35vZnsCzcDZwC05znsMmABMBcYTTKhvKEWAq1evLsVlRES6XBrtV8meA3H3D4GrgPlAI/CUu/8awMyuNrPTwlMvB0aY2ULgDIKJeBER6WYS74G4+66R17cBt+U454rI64+AcUnHlcugQYMAWLlyZRqX73b2vHIOa5o3bXG8rm8NC686IoWIRCSfNNovLWUS8fHHH6cdQreSK3m0d1xE0pNG+6WlTEREJBYlEBERiaXgISwzGwF8NnrM3R/r8ohERKQsFJRAzOxe4GA+/TR4luCW24rRv3//tEMQEYkljfar0B7I/sAQd29OMpi0rVq1Ku0QupW6vjV578ISke4ljfar0ATyWnhuRScQ+TTdqisi7Sk0gfwKeMXM5gGfPBXu7qcmElVKBg4cCPTsnoie/RApT2m0X4UmkKuA35F77aqKoaVM9OyHSLlKo/0qNIFscvczEo1ERETKSqHPgWTM7OuJRiIiImWl0B7IocD3zawFaCHYoyPr7gMSi0xERLq1QhPIVxONopvYZptt0g5BRCSWNNqvdhOImX3V3Z8C8m1xtbTrQ0qPVuHVsx8i5SqN9qujHsiJwFMEmz+1lQUe6vKIJFW6VVdECtVuAnH3SeG/v1KacNLVuhSAbucVkXKTRvvV0RDWbwl6Gjm5eyobPyVl/fr1aYcgIhJLGu1XR0NYD3T0A8yszt3XdFE8IiJSJjoawrqzgJ8xF9iva8IREZFy0RVb2la196GZDQDmAUe6+9tmdjIwGdhEMEF/gbtvbFNmNMEE/TvhoT+7+yldEKuIiHSRrkggeedIzOwgYDowNHxvwDXAAe6+wsx+CZwDNLQpuj/wU3f/SRfEV7DtttuulJcTEekyabRfXZFA2jMJOBP4z/D9cOAP7r4ifP8ocDFbJpADgO3M7ETgbeBMd3+HhL333ntJX0JEJBFptF+JJhB3nwgQdDwAWAA0mNnOBLsbHgdsn6PoKuB+d3/IzE4D7gVGFnPtTCYTN+yy1NPqC6pzT6E6d19J90A+xd0Xm9nFwCxgHXA/cGCO806LvL7FzK41s63d/W+FXqu+Pt/D8/n169cPgHXr1hVdNk2ZTCZWfcuZ6twzqM6F60z71dzcTGNjY9HluiKBFPzUipnVAi+6+77h++Nps8eImVUDPwSudffomhqfmmhPwoYNGzo+SUSkG0qj/SoogZjZFW0OZYG1QKO7jyrienXA781sT4Ltcc8Gbome4O6bzewY4HXgfjObALygZ01EpEtM3Qlacvzd26c/XKJ50GIUuh/IvwD/BmwLDAROBY4Erjazywu9mLt/SLC74XygEXjK3X8NYGZXh/MdAN8FzjOzhcApwMRCryEi0q5cyaO945JXoUNY2wH17v4+gJlNAf4LOBjIAD9ur7C77xp5fRtwW45zroi8Xgh8qcDYur0Rd49gzcYtO1B1veqYf9L8FCIS6QbUEyh7hfZAtm1NHvBJT2Jbd28BKmbiYPDgwQwePLjLf26u5NHecZEeQT2BLpVU+9WeQnsgS8zsJ8D/JXjy/HvAm+GDgltuHlGmlixZknYIIiKxpNF+FdoDOQXYFfgz8CKwE8G8xH7AhYlEJiIi3VpBPRB3/yvB5lJt3dy14aSrT58+ALS0tKQciYgkpk///HMvZSyN9qvQ23i/SLDkSH+CIawa4J/cvbQDbgnbvHlz2iGISNIqdII+jfar0CGsWwlW1B0A3A00AQ8mFVSlqetVV9RxkR4h31/8Zd4T6EkKnUTPuvs0M/sc8BpwAvDH5MKqLLpVVySHCu0J9CSF9kD+Hv77TWAvd19HBd19JSIixSu0B/Kimd0HXA7MNrOhlGBtqlLbfffd0w5BRCSWNNqvQhPIecCI8PyfEyxj8p2kgkrLokWL0g5BRCSWNNqvQhPIT4GzCCbPWz0DfL6rAxIRkfJQaAI5FtgxXMKkYvXqFfzn2Lix4kbnRKTCpdF+FTqJvphgl0ARERGg8B7ITcCzZvY0kcUT3f3qRKISEZFur9AEcjHB/MfABGMREZEyUmgCqXP3LycaiYiIlJVCE4ib2XB3fyXRaFK29957px2CiEgsabRfhSaQwcAfzewtgr3MAXD34YlElZJMJpN2CCIisaTRfhWaQH6YaBTdxLJlywBKvquXiEhnpdF+FbofyLNxfriZXUmw8CLAbHefbGZjgAagH3Cfu1+Wo9xg4C6CBxUdOMndE9/ncsiQIYCeAxGpaBW6F3sa7Vehz4EULUwUhwH7AvsA9WZ2InA7cBQwDDjAzMbmKP5L4JfuvgfBqr+XJxWniPQw2ou9yySWQIAVwAXu3uLuG4BFwFDgdXd/y903EvQyjo8WMrPewCjggfDQHW3PERGR9BU6B1I0d1/Y+trMdicYyvo5QWJptQL4QpuinwOawgST7xwREUlZYgmklZntCcwGfkCwBPzQyMdVQNt9GKuBbJtjRe/VGOeOhGw2G7ts2sox5s5SnXuGrq5zfQmvFVe5tF+JJhAzG0mw9e157n6vmY0Gdoicsj2wvE2xvwBbm1mNu28Kz297Tofq69v7Nclt5MiRscumKZPJlF3MnaU69wyJ1Pm3+T/qDv9949a5M+1Xc3MzjY2NRZdLLIGY2c7AI8B4d38qPPxC8JHtBrwFfJtgUv0T7r7BzOYC44FfAxOAx5OKM+q5554rxWVEJE19+ue/C6uMpdF+JdkDuRCoBRrMrPXYLcDJBL2SWuAxwslyM7sVmOXus4AzgDvN7DJgGXBignF+Yv78YO/yESNGlOJyIpKGMr5Vtz1ptF9JTqKfC5yb5+Mtnrl394mR10uBQ5KJLL8vfzlY7kvPgbRvxN0jWLNxzRbH63rVMf+k+SlEJCJptF9J3sYrFSpX8mjvuIhUJiUQERGJRQlERERiUQIREZFYEn+QsJwceuihaYcgIhJLGu2XEkjE44+X5HGTslfXqy7vXVgiko402i8lkIhZs2YBMG7cuJQj6d50q24FqdClzXuiNNovJZCIb37zm4CeA5EeREubV4w02i9NoouISCxKICIiEosSiIiIxKIEIiIisWgSPeLYY49NOwSR0qrQpc17ojTaLyWQiPvuuy/tEERKS7fqVow02i8NYUVMnz6d6dOnpx2GiEjR0mi/1AOJOP300wGYNGlSypGIiBQnjfZLPRAREYlFCURERGJRAhERkVgSnQMxsyuBE8K3s919spmNARqAfsB97n5ZjnLfBa4FPoiUvTTJWEVEpDiJJZAwURwG7AtkgTlmdiIwDRgNvAPMNrOx7t52HeL9gfPd/Z6k4svllFNOKeXlRES6TBrtV5I9kBXABe7eAmBmi4ChwOvu/lZ47C7geKBtAjkA2N3MLgEWAGe7+8cJxgqgW3hFJKdsNsvatWtZv349mzdvTvRaffr0YeXKlUWXmzp1KkC7Zaurq6mtrWWrrbaiqqoqdoytEksg7r6w9bWZ7U4wlPVzgsTSagXwhRzFVwA/BeYBU4FfACclFWuradOmAXDRRRclfSkRKSNNTU1UVVUxcOBAqquru6TxzWft2rUMGjSo6HJz5swB4Igjjsj5eTabZfPmzaxZs4ampia23nrrTsUJJXgOxMz2BGYDPwA2EvRCWlUBW6Rzdz8mUv464M1ir5vJZIqO9ZJLLgFgzJgxRZdNW5z6ljvVuWfoDnXu3bs3AwYM4OOPEx8IAWDp0qVFl3nkkUcAGDZsWLvnZbNZmpqa2LBhQ6zYopKeRB8JPAic5+73mtloYIfIKdsDy9uU2Ro41d1/Fh6qIkg8Ramvry863ta/KuKUTVMmkym7mDtLde4ZukudV65cGatXEMfSpUvZZZddYpcvpGzb+jQ3N9PY2Fj0tRK7jdfMdgYeAb7t7veGh18IPrLdzKwG+DZbzn+sBiab2UHh+7OAh5OKU0RE4kmyB3IhUAs0mFnrsVuAkwl6JbXAY8ADAGZ2KzDL3WeZ2QnAzWbWD1gMTEgwThERiSHJSfRzgXPzfLx3jvMnRl7PBfZLKDQREekCWkwx4vzzz087BBGpEHteOYc1zZu2OF7Xt4aFV+W+U6pQM2fO5MEHHwRg/fr1vPPOO3zve99j9uzZHHXUUey///4899xzPPnkk526TkeUQCKuu+66tEMQkQqRK3m0d7wYEyZMYMKECWSzWc4++2yOP/54HnvsMa688kpGjx7NHXfcwaZNnb9OR7QWVsTkyZOZPHly2mGIiBTkxhtvpHfv3hx33HEsWbKEv/71r0DpdidUAoloaGigoaEh7TBERDo0Z84cnn76aaZOnUpNTQ0bNmxIfMiqLQ1hiYiUmUWLFjFt2jRmzpxJv379AOjXrx8ffvghAL/97W9LEocSiIhImbn++uvZuHEj55577idzHbvtthtvvvkmxxxzDGZGbW1t4nEogYiIJKCub03eu7A66/bbb9/i2OGHH85ee+3FjBkzeOKJJ1izZk2nr9MRJRARkQR09lbdYtXW1rJgwQKOPvpoBgwYwJQpUxK/phJIRCn+g4uIJOH0008H8q/GmwQlkAgt4y4i5aqUiaOVbuONmDRpEpMmTUo7DBGRos2cOZOZM2eW9JpKIBEzZsxgxowZaYchIlK0efPmMW/evJJeUwlERERi0RyIVJQkF7AT6c5ee+01dtxxx5JeUwlEKkqSC9iJdGd77LFHya+pBCIikoSpO0HL6i2P9+kPl7zXqR+dazn3vn37MnTo0E793GIpgUTcfPPNaYcgIpUiV/Jo73gRci3n/pvf/IYxY8Z0+mcXQwkkQrfwikg5aV3OfeLEiTz77LMMHz68pNfXXVgR48ePZ/z48WmHISLSoehy7gDLly/n0UcfLWkMifZAzOxK4ITw7Wx3nxwe7w3MAX7s7s/kKDcYuAv4PODASe7e+X5fB1rHFKV8JbmAnUh3kWs599WrV7N48eKSxpFYAjGzMcBhwL5AFphjZscArwK3A/u1U/yXwC/d/V4zuxy4HNA6I9Ih3aorPUGu5dxXrVpV8jiS7IGsAC5w9xYAM1sEDAa+CFwPnJerUNg7GQUcHR66A3gWJRARKSd9+ue/C6uTci3nftppp3X65xYrsQTi7gtbX5vZ7gRDWSPd/fXwWM4EAnwOaHL3jeH7FcAXkopTRCQRnbxVtxwkfheWme0JzAZ+0Jo8OlBNMOQVtbnY62YymWKLkM1mY5dNWznG3Fmqc8/QHercp08f1q5dW7LrLV26NNGyTU1NLFu2LPY1WiU9iT4SeBA4z93vLbDYX4CtzazG3TcBOwDLi712fX19sUV4+OGHY5dNUyaTKbuYO0t17hm6S51XrlzJoEGDSnKtpUuXsssuuxRdrnU/kELKtq1Pc3MzjY2NRV8zyUn0nYFHgPHu/lSh5dx9g5nNBcYDvwYmAI8nE+WnjRs3rhSXERHpcnvvvXfJr5lkD+RCoBZoMLPWY7e4+y25TjazW4FZ7j4LOAO408wuA5YBJyYY5yfGjh0LwOOPlyRfiYh0mZtuugmAc845p2TXTHIS/Vzg3HY+P6TN+4mR10uBQyixJ598stSXFBHpEq+++mrJr6mlTEpsxN0jWLNxzRbH63rVMf+k+SlEJCISj5YyKbFcyaO94yIi3ZUSiIiIxKIhLBGRBCQ5XJ1vP5CtttqKww8/nBtuuIFhw4Z16hqFUAKJeP7559MOQUQqRJLD1bn2A5kzZw7Dhw/niiuu6PTPL5QSSMSIESPSDkFEpGBt9wMZPXp0Sa+vOZCIUaNGMWrUqESvUderrqjjIiK5tN0PZNmyZTzyyCMljUE9kIh58+bl/WzPK+fk3WeimCXEdauuiHRWrv1A1q9fz/LlRa/61ClKIAXKlTzaOy4ikpSesB+IiEiPVderLu9dWJ1V8fuBiIj0ZD1huFqT6CIiEot6IBFLlixJOwQRkVimTJlS8msqgUQMHjw472d1fWvy3oUlIpWtqqqKbDZLVVVV2qHkte222xZ0XlfWQwkkonXns1xbaBZzq66IVJbevXvT1NRE//79qa6u7paJ5JprrgHgsssuy/l5Nptl8+bNrF69mt69e3fJNZVAIhYsWJB2CCLSDQ0YMIC1a9eyatUqNm/enOi1mpqaWLlyZdHlPvroI4B2y1ZXV1NbW8tWW20VO74oJRARkQ5UVVVRV1dHXV3yK0YsW7Ys1v7rrTsSNjQ0dHVIeekuLBERiUUJREREYqm0IaxPbolqbm4uuvAOO+wQu2zayjHmzlKdewbVuTCdab9aWlpaXxZ1W2lVNpst+mLdVSaT+TIwN+04RETK1MH19fUFb4xUaT2Ql4CDgRWAVjkUESlMDbADQRtasIrqgYiISOloEl1ERGJRAhERkViUQEREJBYlEBERiUUJREREYlECERGRWJRAREQklrJ9kNDMfgp8zt1PNrP9gF8BfYB3gO+4+yozGwjcDQwBVgInuPv7qQXdCWb2NPB5YEN46N+AfwYuA3oDN7j7f4TnjgEagH7Afe6ee4OAbi5PnY8ATgjfz3b3yeG5ZV/nXPV19xfCzz75fQ/f7wPcCgwAngNOc/eNJQ+6k/J8x58HrgTqgCfc/dzw3LL/jiFvnXcHfhi+f9zdLwzP7dbfc1n2QMzsa8B3I4duBK5w970BBy4Mj18DzHX3YcD08LyyY2ZVwFBgb3ffx933Ad4FpgBfBvYBvm9m/8vM+gG3A0cBw4ADzGxsSqHHlqfOnwEOA/YlqHO9mR1TCXXOVd9I8mj7+w5wF3CWuw8FqoBJJQ24C+T5jlcCtwBHA8OB/cxsbCV8x5C3zv8fuAkYDewNHBwmS+jm33PZJRAz+yxBwzk1criGIEMDbAWsC19/naAHAnAPMNbMumYrrtKy8N9PmNkCMzsLGAM85e4fufsa4AHgOOBA4HV3fyv8S+Uu4PhUou6cXHVeAVzg7i3uvgFYBAymMuqcq745f9/NbBegn7vPDw/dQfnVF3LX+RiC3sW74Xc8HniByviOIXedawja4jqC0YTewLpy+J7LcQjrV8ClwM6RY+cTfCE3AGuAg8LjOxI0Orj7RjNrAgYBy0sXbpfYBvg9cDbBL9czwH2EdQutIPifbMccx79Qkii7Vq46u7s/CWBmuxMMZY0E9qf867xFfc3Mge+z5e97JX/HNQT/L88i+OPgUeByKrvOTlDH14C1wLPAPGAE3bzOZZVAzGwi8I67/97MTg6P9QNuA8a4+4tmdj4wk6D30Xbj4iog2f0oE+DufwD+0PrezG4jGAu+JnJaa92qgWyO42UlT53/N/Ckme0JzAZ+4O6vm9mBlHmd89T3fmBG9Pc9VMnf8XCC3vUhwGpgFsHw3QYqt85nALsCuwB/I+hdXUiQRLp1ncsqgRB0Z3cws5eBzwL9Cf6jr3P3F8NzfgX8OHz9HrA98K6Z9SIYQ/+wtCF3npl9Gejr7r8PD1UBbxOsntlqe4Ke1bt5jpeVPHXeYGYjgQeB89z93vCzsq9znvr+ETgs+vtuZj8DbqDM6wt56/xX4P+5+8rwnIcJetb3Ubl1Phr4P+7+l/CcOwiSyv108zqXVQJx90NbX4d/kR0C/DvwmpmZuzvBJFvrksSPARMIxo/HE0yob6D8DASuNrMvEXR7vwt8B7jLzAYRDNsdSzDc8QpgZrYb8BbwbYLJx3KTq85nAY8A4939qci5L1D+dc5V39PCv1g/+X13938P3683s5Hu/t/AvwKPpxN2p+Sq8znAr8I7KP8OjCX4zivhO4bcdf4RcIyZ1REMYX0DeMndl3b377nsJtHbcvePgZOB+83sFeBU4JTw48uBEWa2kCCjn5lKkJ3k7o8SDNn8GcgAt4e/UJcCTwMvA7929xfdfT3Bf48HgVcJxlUfSCPuzsiTHR/fAAAEv0lEQVRVZ4KbBGqBBjN7OfzntEqoc57v+A/tFDkJ+JmZvUbQE78p+Si7Vp46PwtcBzxP8F0uJRjGK/vvGPLW+SqCm3wyBH8A9gauDYt06+9Z+4GIiEgsZd8DERGRdCiBiIhILEogIiISixKIiIjEUla38YqIlKvwGaafESz6+iFwqrsvzXPuZwjurvyeuz8THlsCNEVO+4a7v2NmPyC4+zQL3OruDR3EYQTPy20DvA98K7ybtWjqgYgUyMz2N7Mtbh01s0fbPCle7M99xsyO61RwUg7uBiaGCyjeTfu35P6CoIEHwMy2BVoiC23uEyaP3QgeUdgPOAA4NzyWU7iY4yzg2nDx2T8DF8etkHogIgVy9z8SPIsi0iEzOwQ4Odxyoi9wmbu/En78CsF6WLnKjSd4iPKVyOEDgCoz+2+C5ex/4u7/RdAJ6EPwfFRV+M+G8OdMAM4Lz8kQPAe3J7DG3eeEP3cqwcONsSiBiBQobBB+QbCk/J0EC/wtJdjboZDylwMnAhuBxQTLdLfuT3OMmV1MsJr03e4+JVx+5+cEC0ZuAJYAp7j76i6rlJSEuzcTrHGFmVUTPH3+SNvzzGwwQaP/VT791Hkt8DvgImA7YK6ZNbr7IjO7h2Bpo2qCIayl4Xpxk4Avuft6M/sJwfparwPvh2tw7UuwonXORFYIDWGJFO8/gPnuvifB0ht7dFTAzE4hWJbjAHcfDjQSLM/dagDB6qsjgO+Ee118kWC5nr3dvZ4ggQzvumpIEsI9al4m2AhqXLhiwozwsz4Ew1e9+PSWFK2J5TaCPyzWRT9z90fc/Rx3X+fubwMPEayTdgRQD+wU/nOAmZ0AfIVgk6r5YSxHEfye9iL4nbrZ3fcj+J1qd86kPeqBiBRvDOGmZe7+hpk91cH5ECSPGeHeLRBsbnZp2KBA8JfjRqApnGc5FLga2AS8YGa/Ax6MLBoq3ZS7Pww8HB3CAjCz/gTzDx8CR+VYl2+P8J/bgnludgNuNbNJBHuFvB8Oo8I/hqq+QfB7sTq8xj0EG1MtBu5393Mi1+5FkGxej/yce+jEkjDqgYgUL8untwooZIvRGj69NHc1wf/QrT9nU5vPNrj7KoId6i4MP7/PzM6IG7Sk7i7gDYLFQJvbfujur7r7zpGdCv9IMOn+NMFy71eaWbWZbQeMI1hTawFwpJnVWLBZ3liCxWSfIRgW/Xw4cX4zwdDYPGCQme0dXvYbBPMjsSiBiBRvDsHKx61j1l8psMyp4YqrEAx9PRdpSCaYWZWZbUOwUdYcMzuSYPOhee7+I4J9bg7oumpIktz9mUjvY1+CYaSRwJ/CYa3Hws9OM7OrO/hxtxBsKNVIsOHUReEtwLcSLC65kOCOqoXAne6+ALgKeCo8VkNw59U6gl0fp4eLzH4VuCBuHbWYokiBIpPoXwFmAP9MsBdJL4L/ae9op2zrxOmxBH+4vQGc7u7vmtkzBJOb9QR32Nzs7jeZWQ3BrZ5fI9hc6WNgUjgGLpI6JRCRApnZOOBSdz+ow5NFegBNoosUwMzOAi4DrmjnnLkEu17mcrC7/z2J2ETSoh6IiIjEokl0ERGJRQlERERiUQIREZFYlEBERCQWJRAREYlFCURERGL5HwEGLzSauXEWAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df['mag'] = -2.5 * np.log10(df.forcediffimflux) + df.zpdiff\ndf ['mag_err'] = 2.5 / np.log(10.) * df.forcediffimfluxunc / df.forcediffimflux",
   "execution_count": 33,
   "outputs": [
    {
     "output_type": "stream",
     "text": "/Users/rbiswas/soft/mypython3/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in log10\n  \"\"\"Entry point for launching an IPython kernel.\n",
     "name": "stderr"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "grouped = df.groupby('filter')\n#fig, ax = plt.subplots() #1, 3, figsize=(18,6))\nfor band in grouped.groups:\n    x = grouped.get_group(band)\n    ax.errorbar(x['jd'], x['mag'], yerr=x['mag_err'],\n                label=band, fmt='.')\n    ax.plot(x['jd'], x['mag'], 'ro', markersize='20', markerfacecolor='w', alpha=0.5)\n    ax.axvline(requests.jdstart, ls='--', color='k')\n    ax.axvline(requests.jdend, ls='--', color='k')\nax.legend(loc='lower right', framealpha=0.4)\nax.set_ylabel('mag')\nax.set_xlabel('jd_obs')\n# ax.invert_yaxis()\n#ax.set_ylim(bottom=35)",
   "execution_count": 34,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 34,
     "data": {
      "text/plain": "Text(0.5, 17.200000000000003, 'jd_obs')"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "lcs = pl.get_lightcurve('ZTF18aaermez').table.to_pandas()\ngrouped = lcs.query('magpsf < 98').groupby('filter')\nfor band in grouped.groups:\n    x = grouped.get_group(band)\n    ax.errorbar(x['jdobs'], x['magpsf'], yerr=x['sigmamagpsf'],\n                label=band, fmt='.')\n    ax.axvline(requests.jdstart, ls='--', color='k')\n    ax.axvline(requests.jdend, ls='--', color='k')\n    ax.legend(loc='best')",
   "execution_count": 35,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig",
   "execution_count": 36,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 36,
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEPCAYAAACp/QjLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNXd+PHPzGQmIQlgWMMStioHShRiQBaXWmvd0KLWVlusW8Hu4vbYPnV7tNb20VK7+Pxa64JiaUVFxYKiVVxRQEKgQONBDbuEJQ1L1klm5vfHmYEhmZnMmjuT+b5fL14kd2buPZcZznfuPed8vzafz4cQQggRzG51A4QQQqQfCQ5CCCE6kOAghBCiAwkOQgghOpDgIIQQogMJDkIIITqQ4CCEEKIDCQ5CCCE6yEn1AZRSvYAPgAuBLwL3Bz08BFiltb4w1e0QQggRvZQGB6XUZOBRYDSA1voV4BX/Y8XACuCmVLZBCCFE7FJ9W2k28CPg8xCPPQj8WWv9SYrbIIQQIkYpvXLQWs8CUEods10pdQJwJjAr2n1VVFTkApOA3YAnaY0UQojuzQEMAj4qLy9vifZFKR9zCON64P9praNuKCYwvJei9gghRHd3OvB+tE+2KjhcDJwT42t2A4wePRqXyxXTC0eOHAnAli1bYjyk9TZu3EhpaanVzehScs7ZQc45eon0YW63m82bN4O/D41WlwcHpVQ/oIfWOtaz9AC4XC5yc3NjeqHb7QaI+XXpIlPbnQg55+wg5xydJPVhMd2Ot+LKYRSwsysPuH///q48nBBCJJUVfViXBAet9Yign1cDU7riuAETJkwAYN26dV15WCGESAor+jCrxhy61Pr1661ughBCxM2KPkzSZwghhOhAgoMQQogOJDgIIYToQIKDEEKIDrJiQHrYsGFWNyGrLLxvNS1NbZzz3XEUj+ptdXOEyHhW9GFZceWwbds2tm3bZnUzskJN9UH276rncG0zix+qpKb6oNVNEiLjWdGHZUVwGD58OMOHD7e6GVlh1+Y68JmfPR6v+V0IkRAr+rCsuK20fft2q5uQNYaMLiLHacfj8eJw2BkyusjqJgmR8azow7IiOIiuUzyqNzNuKmPX5jqGjC6SMQchMpQEB5F07/xN09LUJlcNQmQwCQ4iqQID0vhg8UOVzLipTK4erLBnD1RWws6dsHcvuN3gcsGAATB0KJSVwcCBVrdSpDEJDiKpQg1IS3DoQnV1sGQJ7NtnAsDZZ0NxMeTmQksL1NTAli2wYAH07w8XXghFcoUnOsqK4DB+/Hirm5A1howuAhvgQwaku9qmTbB0KZx2GsycCfZ2kxHz8mDECPPnjDNg1Sp49FGYPh3GjbOixSJKVvRhWREcJFV31yke1Zt+QwplEVxX27QJli2Dq64yVwqdcThg2jQYNcpcRYAEiDRmRR+WFcGhX79+gBT96SqX33GK1U3ILnV15ooh2sAQrLjYXGXMnw+DB8stpjRlRR+WFYvgamtrqa2ttboZQqTGkiXmVlKsgSGguBhOPdXsR6QlK/qwrAgOQnRbe/aYwecpCRZXnDrVzGrasyc57RIZT4KDEJmsstLMSmo/+BzkwAsvcuCFFyPvx26Hk08GGZ8TfhIchMhkO3fCyJERn7L+sI3HP7dTsa2TPFcjR8KOHUlsnMhkWTEgLUTauH8IuOs7bJ7g6AHlNbHvb+/eiGMNFdvqmLtiJ2P3V3Pvvz/nrp9eQfnwMIPOxcVmf0KQJcHhrLPOsroJQhghAgOAw9MU5/7cZoFbGP9+YwX3vP8XnN422vQbbCgvofy7F4Z+sssFra3xtUOklBV9WFYEhzfffNPqJgiRmI0bYfFi2L4damtNJ+50Qn09vPEGjB9vUmO0c9L+z3B423Dgw+dp46T9n4U/httt9inSjhV9WFYEh169egFw6NAhi1siRIy2bIG5c6GpCc46Cy66yCxcy8+HxkZ46CFoa4MPPzQd+9Spx6xVOOGcM9j69ON4PW04XC6+cM4Z4Y9VUxMywAjrWdGHZUVwOHz4sNVNECJ2ixfDM8/AN74BM2aYVc1Ac1UVAHljx8KkSea20jnnwObN8OabUFoKY8YAkF9WRr9Z38VdXU2fa68hv6ws/PG2bIGSkpSfloidFX1YVgSH7m7Kgik0tDV02F6QU8DKmSstaJFI2OLF8NJLcN998IUvhH9eWZlJf3HGGTB2rMm0+u675jF/gBhw45zOj+fxwNq1cOWVSWi86A4kOHQDoQJDpO3CInv2wPYe0JgDLXngtYHdB7nNePLcOPbsMZ37li3miqGzwADm+f37myR606ZBnz4mUCxfbh6LNh3GypXmlpKk8RZ+KQ0OSqlewAfAhVrrrUqpc4AHAQewFpiltXansg1CpNqLc9cCcMktJ4d+QnAa7e/83awnaJdG+z+rV9M/kEb7vffMraTOAkPAhRea7KqjRpn99ukDX/yiGYe44ILOX19TAytWwOzZUZ6xyAYpWwSnlJoMvA+MDtr8OHCF1roUyAeuStXxhUgLmzaZjvsLX4Abb4Qvf9mkzM7LA5vtSBrt1tNOgzn+2z/79pk6DGE0uD00uD1HNxQVmbTbCxaYjh7MLSW3u/N1CzU15nXTp0vSPXGMVF45zAZ+BDwdtM0B9FJKOYA8IM7J3bG5+uqru+IwQhwrTBrtqg92AzB22qBjn+9wmKmqF18MBw+a3wsKOuy25WNN67830cPbenSAOZBue/58k0Rv6lRzJaF16BlIHo+5lbRihdRzyABW9GEpCw5a61kASqngzT8E3gYOAVuA51N1/GBPPvlkVxxGiKM6SaO9r76FjWt28I2J7WYHbd9upqsOHGi+1btcx6w9aKys5D/33YOtrY36F19g5FPzjg0QgwebW1grV5orlLo6M6PJ5TJXEoFKcGvXmqAxe7ZcMWQAK/qwLhuQVkoVA78GSjGB4bf+Pz+KZT8bN26M+djTpk0D4IMPPoj5temgoqIi4uN59jyavc0ht3f22nSVSe0+XG9WPQe3uf+yZfQoLaXV64XPPz/m+W63G5/Xi9vt4fOgx957/mOcB4cy2HMcA30+bPn5eD//nNY+fY48Z8XCVxna2oYdH21uN+8ufJVTjz/+6M5zcuDii7Ht3Yt99WqcTz7JgU2bsLW14cvJoaV3b5r69aPxxBNp69sXqqtT848SpUx6n5MlnnO2og/rytlKpwMbtdafASilHgWejXUnpaWl5EZIFxCK223GvMvLy2M9nOUqKio6bfdH5R8l52BpUpQ+mnNOJ9vfNgPS5eX+Aek9e6BHD3MFECJb6sGtu7G5W3C5HAwePBiADas+Y+Nbu/H1nsr6J6qZcVNvikf2haYm8nJzzfsAvJ0/nMvtOfi8bbTZc3g7fzjn9ezZsVE9e5orh4UL6f/QQyk570Rl2vucDPGecyJ9WEtLS1xfqrsyOGwE5iqlBmqt9wAzgCT1aiIhUpQ+Ie6mNlqa2qipPmjKokaRRru9/dsa8Xp8YLPjafOya3Od2VdhoUmR4b96OOnc0/nvnQc4af9n/KvfF7jq3NPD77SpyVxJCBGHLvvkaK2rlFJ3Am8ppdqAT4Hru+r4IgwpSp+QmuqD7N9VDz5Y/FAlM24qo3jnzoizjQ7ua6Rqcy26wMeo/oWUDy+i3/B87A4bvjYPjhw7Q0b7g2+PHiZ4+3178jAO7yhj+Q7FVdNG8+3Jw8I3rroa+vZN1qmKLJPy4KC1HhH081PAU6k+poiSFKVP2K7NdeAzP3s8/m/8EdJo11QfpOK17Xi9XoYBt2xfxdwfTGZQSQ8uueVkdj3yAkNGFJirBjg6kBzkshN6cdkJveg7PkJgANiwwQR1IeIgxX6yVWA2zcyZ8RelX7r0mG+12WjI6CKwmZ8dDv83/ghptHdtrsPn9WHHhh0obraxstrUBi4e1Zvy70yieM1rR1Nn22zg8x2zjwKXgwKXI3LDWltNttaLLkrk9EQWy4obkvfcc4/VTUg/ySxK/53vJLdtGaR4VG/6DSmkpamNc747znzjd7nMWE1e3rFPrq9nSH4DdrsPjxe8QI3Ty5TW/eTttEOvXiZpXo8e8I9/wKWXmsBgs8XesJdfNmskSkuTcp4pF6YIEq5C+Pmurm9PmrGiD8uK4HDXXXdZ3YT0EihKP3Nm2KeEXagVbOpUM58+kBMoS7l65ODqkXP0VtCAAWYQP3BLp6nJZExtbKS4pJizvz6cjz45THUezD2jhPI+OTTs2mVuA+Xnww9/CL/4hanRMGTIkZlKAXljx0Zu0GefwfPPw/33J/9kY7DtOyYBwvCn53f+5DBFkMJuzzJW9GFZcVvJ6XTilCImR0Uxm+bgvka2bthPTfXB8PuRovShDR1qZneBmRK8dq2ZbTRlCth3UzQsj6lnDee8M4ZQbq+CAxvILzgIkyeb59XVwXnnmcR7VVURK7118Nln5nVXXNFpbWmROazow7LiyqGtrc3qJqSXTmbT1FQfpPK1rXi9sH3DPmbcXH70W3F7I0eae9viqEAa7TFjTJAYP95MSQXweUK+xIbHBNuSEjNN2GYDrxduuw1mzYLLLjtSzyGk1lZzK+n5501gmDEjBScWhaC1MkWff47d5zVXMF28Vqa7saIPy4rgINrppCj9rtUb8Hm9gANPm4ddqzdQPOq00E+WovQdDRxobg+9+KKZCRYIDNEqLISTTjJXZGecYbK0LltmAvqJJ5rZYj16mNtV1dXmdtQbb5gxhvvvt+aKIcRaGW9hT7xOJ0y/IPRaGZHWJDhko06K0g9xbcDBCDz4cOBhiGsDECY4SFH68D75xCxgCxEcHDYv+XRMeXJEfb15vVJw552mhvQ//mEK+dTWmtKgOTlmHcOIEXDLLdYNPodbK7O20vwdZq1M3ogRkGUrpDOJBIdsFG42jV/xpHJmrJ3DruaxDOnxCcWT5obflxSl71jHYc8eU9/5uuvMN+VYpwsH0mhfdx28/rrZX2lpes48irBWxr19O+7qalwjRxxNDhi0Vmbg739v1smMG2dmJYWbrSQsIcEhG7WfTdNeySkUf+/3FG99D0Z8H0pOCb8vKUrfUWDA/8QTzbfo4DTaQH/nFjZ5h3d8Xag02vv3m9tL557bxScRhQiZZxsrK6l9/HFoa6P+3XcZNu+JY+tXFxfTdvnl5vWDB8t01TSUFcHhlVdesboJ6SUwmybS6tmSUyIHhQApSt9R8IB/+zTag51QUszHnl78a5+brw1sotx2AO/2Guy7V3dMo53OA/4R1so0rv7I3Pry+fC1ttK4+qNjgwPgHTBA1spEyYo+LCuCw/nnn291E9JLcFH6SDNgOiNF6UNrP+BfVGQ6vz17YN06Kipr+PnObbRiZ6HNx4IhTYzo24++V36l40yedB3w72StTP4pk/DlOKGtDbvTSf4pk0LvR9bKRMWKPiwrgoPD3wF6PKGnEWad9kXp45VpRem7KiV5uAH/gQPh3HNZ6foU9y6ND2i12VhZeiqDVH7oY6frgH8na2Xyy8pwf/cO7Fv+zfHXfK3DVcMRwWtl0vHWWZqwog/LiuDg9XqtbkL6aV+UPlYWFqW/dtm1AMw7bx5er5edO3fS0NAQ9vn2gwfp9e679KivxzdhAt6yMnNLw5+S3L53L/Zt27D93//RVFjIoTPOwNs7zLqOaJx3nlm8FibtxWn9vKgZg8CfGaNfYQv/+U8zBw+GWHDo85lOs6oqqkMXFBQwdOhQ7DGkC49LJ2tlAKr6DGdjTjGX9hlBxDlJ6XzrLE1Y0YdlRXAQIQQXpY93No1FRenr3fUcbj3Mur3rGMxgbDYbSqnQHeKmTWb652mn+Vcoh3hOSYmZUjljBsetWsWgd99NLCX5vn0mT1KE6cI1B5toafPSrzCXgtwc3G43rnZpMgAzq+zQIXOl1wmv18uuXbvYv38/A1I9SaCTtTIV2+r404ZaBrfaueVPJvNs+fAwn5V0vXWW5SQ4ZLOgovRbeiqaT5zI2NOGhH9+GhSlX7d3HbpO48PH7Ndn89jJjzH2+LHhA4MVKckDU4UjBAdXjgNXjoOC3E7+C7a0dMitFI7dbmfgwIFs27Yt9cGhk7UyK1d/zqWHXTgAT7P5PWxwSMGts5jyOomQJDhkO/9smh03z6e1ooai3SdQPMl/qykNi9Kv2bMGn7+AQqu3lZbWltA5ZyJMs+xUICX5/PlmplGs55qfbxaq9ewZX0bVAJ/PrJeIoWCP0+nsmlQLnayVKWmzsw2wYwN8lLRFuM0la2XSkgQHQU2dnQ15J4LPx6dv+5ix6z2KG3eYb3NOpwkIJSVmVpLFg88TB0488rPT7iTXkYstVAccNM1yX+M+APrnd35r5ohEUpI7neZPQ0PsqTOCNTQc3RfA7vXgC3Hv2WaHQePNj4kEo1h0slbmlMmD2fHebnxeyMmxc8rkweH3JWtl0lJWBAdfu2Ip4lhHq5nZ8Phg15jTKT5vhMWtCm3CgAm47C48Pg+3TboNlzfELZcoUpJHJZFplr17mwVsubnxfStubTUpNPr1O7otVGCItD2VOlkrUzyqN6NOH0zD/iZOmz4qfOJGkLUyUbCiD8uKlN0ispDVzNLUur3raPW24vF5eOCjB3B73B2fFGaaZasnxvvaiaQkz8kxAaK2Nvb76a2t5nW9e5v9pKOyMvPvHGFqZUHfPAaoosiBIbBWZsKEFDRSJCJNP3nJFRislCmtoYWsZpamOow5eFo6PilommWrp/XI8w+0HOC43ONwOmL4Jh/jNMtVq1bx5z//GafTyc6dOznrjDPIt9t544MPwGbjL48+yrJly1j04ks0NzWRl+ti7ty5DB06lFUrV3LfvffisNmYUFbGZ9u28fTTT0ff1q4UxVqZ/oVHB6zDFo9K0VoZT3093kOHaKysDL/GIoNY0YdlxZWDz+eTW0uduPyOU+jZJ48PX/zM6qZENHHgRGz+y5zAmEMHQdMsmz1HM5/68B3ze1SimGa55eAWthzccuT39evXc88997Bo0SIWLFxIn6FDeeHxx1HDhrF00SLeeP11/v70fF5duoQzv/QlFsyfT9uBA9x26608ePfdvLR4MTlhBnrTyoUXwvvvmzGDEMZOGxSxkqB9714z8y3J6bsbKytp+fhjWnftYvu119FYWZnU/VvBij4sK4KD6D4mDJiAKlIMKRzCo+c8issRYsyhk2mWMYljmuXo0aMZNGgQPXr0oKioiKmnngp9+zJ45EgOHTrE3DvuYOkzzzD33nt565//pPHgQbTW9O3fnzGTJ0NODpdddlly2p9KwWtlwgSIsGpqyFm4MCVrZRpXf2RmesGRvE4idllxW0l0L8997bkjP1fVhlg5HDTNMs+RR1NbEwA2bOQ5YvxGHsc0y/ZTawOpD3A42L13L5f/8IdceeWVnHHBBfQbMYKqqiooLMSbiVe3QWtljmSeDbHm5OC+RupqGinqn0vx3ipYsYI9ZWX0S8FamfxTJpkpxD4ftkh5nUREEhzEEe6mNlqa2qipPpjW4w6dCppm6XQ4j9yGinm8AZI+zXLDhg0MHz6ca665hubmZv7whz9QXFzMSP9VhdYapRT/+Mc/Or7YZg8/ldVK7TPPnnyyGavxr5Wp0bVUvrYNr9fH9nV7mFFaS/Hs2TRXV6ekOfllZeSOGYP30CEG/+bBbjHmYAUJDgIwdaP376oHHyx+qJIZN5VlboBoN80yx24+5jEHBkj6NMvTTjuNjz/+mAsuuACfz8ekSZP45JNPcDqdPPDAA/z0pz/FbrczcuRI8tqPO/jXMqSldplneeMNM1bT2squxuF4PSVgs+HBxuufl9DziS0MOzN1zXEUFuIoLJTAkICsCA4yGN25o2sdwOPxsmtzXeYGBwtTkk+ePJnJkycf+X358uVHfv7JT34S9nWHDx9m+fLl/O1vfyM/P5958+axZ8+e+NptJX/m2WBDqg/CgxXgA0eOndweWdHtJJUVfZi8SwIIWuvgS/+1Dp1qN82yKC/Oc4lymqXX58Xj89DY2ki+Mz+uQ9ntdo477jguu+wynE4nQ4YM4Ze//GVc+0o37adKp/uMOGGkNDgopXoBHwAXaq23KqWuAW4DPMBy4BatdcoTwcg6h85Fu9YhOF12WguVkjyWeg5RpiRvbG2kuc1Mj916aCsjeo2IO0Bcf/31XH/99XG9Nt25euTg6pGTuVejFrOiD0tZcFBKTQYeBUb7f1fAfcAkrfVupdT/A24AfpuqNgTIbaXodKv/wMHTLKdPh48+Mik1ysrMArni4iP1HI4kFlywwFxxTJpkkvZFMc2yofVoHQmfz0dDa0PcwUGIcLrbbaXZwI+AwBLPk4APtda7/b8vAX5GFwQHEZ1LbjnZ6iYk17hx8OmncM898I1vwLe/3XEMIi/PDFyPGGES9f3jH+b5V1wRVbruAmfBkZ9tNtsxv4ujgj9bgVlxh/ckMB7UCUnVnbiUBQet9SwAc8EAwHrgt0qpEuBz4DIg5hJkGzdujLtNFRUVcb/WSunU7r0H9tLobWThews5Pv/4lB0n2nPOyckJWwXO9u9/46yqwnfvvdgPH8b24YcwaBAcd5zJlupwmEHn+no4cAB278Y3YgTee+/FtmgRrYMH4/viFzttg8vuwuvz0t/VH5/bR4M7fFW6SCJVs4uF2+1Oq89MsMN72ti/05xn1VKAVfQcmF1Dn4m8N135vnbZu6K13qyU+hnwMtAEPAucEut+SktLyY1z9Wt5ecRihWmpoqIibdq9bu86dmzcgQ8fv9n2Gx4951EmDEh+wrRYzrmqqoqCghDf1uvq4K234Jprjo451NebW0jbtpl02B6PCRAFBaZy24QJR1NsX3MNzvnz4fjjO721lNNm/hv16dkn2lPsoKGhIfR5xMHlcjF+fHpOe61YthUw6xu8XuhpH0h5+Qgrm9SlEv3/HM9rW1pa4vpS3WXBQSmVB6zWWpf5f/8GINMWLDLu7mU0tHTMqFmQ62DTPeeFfE37pHdr9qxJSXBIiqB6DkcUFprOPhqJ1HPoRCAP08jeI5O630wQPCvObiezZ8V1c125tLIAeFMp1VMp5QJ+AizsigNL4r2OQgWGSNuhY9K74MI7aSVQz2HKlMT2M3WqmdWUiesN0lRgVlzPvnmMnV6QkskPL85dy4tz1yZ9v1bq1on3tNa1wD3ASmAjsFxr/beuOr5I3IQBE8hz5OGyu1J2SykpwtRziFki9RxEWK4eOfTsk5d1Yw2ZJuXvjtZ6RNDPjwOPp/qY7QVKJ8rVQ+IcdgcOHOkbGOCYeg4JS7Sew1lnkZ+fzxv+ffzlL3/h5Rdf5s1X38Tr9uJ0Opk7dy4DBw5k1apV3HfffTgcDiZMmMBnn32WvvUcRJeyog+TlN0iJkMLh9I7tzfr9qbxt+mgeg4Ji6Kew8jeI48ZPzimnsOCBfTp04cXXngBpRRLly7lg3c+4IE/PcCSJUs488wzWbBgAa2trdx22208+OCDvPTSS+SkawW4DOBuauPwf5qpqT5odVMymgQHEbV1e9eh6zS76ncx+/XZ6Rsg0q2ew9SpAAwePJhDhw7xs1/8jHf++Q5z587lrbfeorGxkU8//ZS+ffsyZswYgMyo55CGAgkkD9c2s/ihSgkQCZCvJ1ki3Oyk9gpywy9MypjZSkH1HBKWzHoOwO7du7npuzdx0Tcu4owzzqBfv35UVVXhcDiyJr1LYEFcKubsd6sEkhaT4JAlIgWGrb+eHtU+ArOVfPjSe7ZSUD2HhKWgnsOgoYO46IqL6JfTL7Z6DqJT3SqBpMUkOIioBUp0Hm49zK9P/3V6XjVAh3oOCUlyPYfJUyezftN6fnDFD/D5fEw5ZQpbPtsSXT2HLHD5Ix8CsPB7U+N6fbQJJEXnsiI4yCyl5Cl0FVLoKkzfwABpXc9hX+M+9jYeHeAekD+A/vn9u089hzTQrRJI+nW3xHtCWKNdPYe4RVnPIRbhEvV153oOXa3bJZC0SFYEB1nnkDxpX8chIFQ9h1hEWc8hVvnOfPJy8vD4PAwtHHpMeu/uXM9BJEbWOYiUCTcLKdLspIwWXM+hpia219bUHK0D0UnSvXjYbXacdqfUfRBpLSuuHARhk+l1a4F6DPPnmyR6U6dGTqnh8ZhbSStWmMAQRT0HIborCQ6iexs3DgYPNtlVV640uZJGjjS3mlwus44hUAlu7VozxjB7dkquGITIJBIcRPdXVGTSbu/ZY5LovfGGSYnR2moWuA0YYKarXnllUgefRZSCanvPeu9fOD1tUPtW6NreostIcBDZY+BAOPdcq1sh/ByHDsHTTx9T2/vQ2C/hdeXCiQM61va+8EK5outCWREcZJaSEGlm0yaGvPIKXHIJzJx5ZCzIu3+HeTy4tvcZZ5hpyY8+2vlY0P1DwF3fcburEH6+K+mn0VVknUMGmrJgCg1tHWv/FuQUsHLmSgtaJCKJpwJetObPn8+iRYsAaG5uZseOHcyfP5/7778fj8fDxIkTeffdd/nnP/+Z0HHSXmcd9KZNsGwZnpkzTWnWzjgcZr3KqFHmKgLCB4hQx420XYSVFVNZbTbbkXnCyRYqMETaLqwVTwW8aF111VUsXryYl156iRNOOIGbb76Z+++/nzlz5rB48WJKSkrweBI/TtqL1EHX1cHSpTBzJt5Yc1YVF5urjKVLzX6ySCr7sHCyIjgI0ZV+//vf43Q6ueyyy9i1axdf+tKXAPj6179uccvSQKja3rEIru0tUkqCgxBJtGzZMt566y3uv/9+HA6HjHcFq/dKbe8MIsFBiCSpqqrif//3f3n44Yfp0aMHPXv2pKSkhHfeeQdA0nDvc0ht7wwiA9JCJMmDDz5IW1sbc+bMOTK2cOedd/Lggw/yu9/9DqVUVqbhPqKl0CxAjGBrbQPV+xoY1b+Q8uERpq2Gq+3tKgw/GC5iIsEhQQU5BWFnK4n0U5DrCDtbKVFPPPFEh20PP/wwf/zjHxkwYACvv/46DQ1ZMFEhXAfd1jPiWEPFtjr+8m41bR4fb+m9LJg1JXyACFfbO4Onq6abrAgOqbzvK9NVM0tX55gaPHgw1113HTk5OfTq1Ss70nCH66D/538i1vZeWV1Lm8cUom1t87KyujZ8cIizBxjfAAAf7ElEQVSjtncmk3UOQnQzl156KZdeeqnVzUgPndT2njKqLzkOG20eH84cO1NG9Q2/rzhqe4vYRBUclFJ3tdvkAxqBjVrr15LeqiSTeg5CpIFOanuXDy/i+jNGUb2vgVmnj4o85pDk2t7pLp3rOZwIfA/oCxwHXAdcCNyrlLozRW0TQnQngdreEYzoW8BZYwZEDgyQ9NreoqNog8NAoFxrPUdrfTMwEXP1cDrwzVQ1TgjRjZSVmeyria4SD9T2jib1hohbtMGhr9b6SDktrXWtf5sbyJ5RISFE/IJreyciBbW9RUfRDkhXK6V+BfwFsAHfBT5TSk0Gwn4NUErdzdEri6Va69uUUmcDvwV6AAu11nfE3XohRGYJru0dj85qe3fTrKxWiPbK4VpgBFAJrAaGALOAk4FbQ73AHwTOAcqACUC5UupbwBPADGAsMEkpdX4C7Rcird1+++1s2LDB6makj6Da3vZQ6xQiiaa2t2RlTZqorhy01vuBb4V46E8RXrYbuMV/6wmlVBUwGvhEa73Fv+2vwDeAV2NpdKxklpKwSlasa4iVP92248knoaFBantHIW3XOSilpgI/Awoxt5UcwEit9bBwr9Fabwp6/QmY20t/xASNgN3A0FgavHHjxlie3i1UVFRY3YQuF+055+TkxLTqOP/3o7GF+BbpcxXSOGdz1PsJ5e9//zuLFy8GoKWlhZ07d1JSUsLPf/5zJk6c2Onrk7V62u12Z8RnxnHBBTS9/Tb5S5bgKyvDO3w4bY098TqdfF5djWPfPuzbtmGrrKSxsJD/TJuGp7kZIpxbeYTjpcu/Sbq0ozPRjjk8BswHLgP+DFwMLIrmhUqpccBS4L+ANszVQ4AN8EbbWIDS0lJyI6yyDCWT1zlUVFRQXh7pI9/9xHLOVVVVFBTEkKokzO0Fm7s+tv2EMGvWLGbNmoXP5+MnP/kJ3/zmN3nnnXfIy8vrdN8NDQ0JHz/A5XIxfvz4pOwrlSoqKhh2++1Ha3tv2EC/t/9FjtfD4NIhR2t733gjDBxI2G+iwSLkNkyH/0fx/n9OpA9raWmJ60t1tMHBp7X+X6VUP+BjzFXAms5epJQ6FRNEbtRaP6OU+hIwKOgpxcDnMbZZpJBUtktcoJ7DrFmzjmRkFREE1fZ+zPMhAF/+3lQrWySIPjgc9v/9GVCqtV6hlIo4WVkpVQK8BFyutV7u37zKPKSOB7YA38YMUIs0IZXtEhOo5/DMM89Y3ZTsJFlZkyba4LBaKbUQuBNYqpQajblFFMmtQB7wW6VUYNufgWswVxN5wCvA8zG2WYi0FKjnMH/+fHr06GF1c7KTTFdNmmiDw43AFP/z/4hJnXFlpBdorecAc8I8nP43RGM1b7r5+9ql1rZDWCZUPYfNmxMb5BbCKtEGh98APwYOBW17G8iIzFddXZhbpLEU3nYIVc9BiGSwog+LNjh8HRjsT5uRcbzemCZEie4sm247yGrhbsOKPiza4LAZOJDKhoj0IJXtupEMXC28UGYppY1og8MfgHeUUm8RlGhPa31vSlqVZJm8zqGryXRVIdKPFX1YtMHhZ5jxhuNS2BYhhBBpItrgUKC1Pi2lLRFCCJE2os3KqpVSJ6W0JUIIIdJGtFcOw4A1SqktQEtgo9ZaAoYQ6UpWC4sERBsc/julrUgxWecgspJMV+020nadg9Y6o7OHyToHEZDKxILz589n0SKTrLi5uZkdO3YwbNgwjj/+eD755BN+97vfMXbs2ISOIbKTFX1YtGMOQnQLqUwseNVVV7F48WJeeuklTjjhBG6++Wb69++PUorXXntNAoPIKFkRHGw2m9xaEl0mOGU3wEknydCcSIwVfVi0Yw5CiCiEStmdl5dnYYuEiI8EByGSRFJ2i+5EgoMQSSIpu0V3IsFBZJVUJhaUlN2iO8mK4GC3Z8W4u4iCJBYUmciKPiwrgkPgEl8IITKRFX1YVnylfvXVV3n11VetboYQQsTFij4sK64cLrjgAkDqOQghMpMVfVhWXDkIIYSIjQQHkZUuf+RDLn/kQ6ubIUTakuAghBCig6wYcxDdw7i7l9HQcuysjUe/NgjProOUDultUauE6J6yIjjk5GTFaXZ77QNDgDdNJhqsWbOGX/ziF8ds27x5M0OGDKGgoICDBw/S2NjIoEGDAHjggQeYN28eK1eupGfPnkfmsp955pncdNNNYY8zY8YMFi9enLoTEWnHij4sK3rN1tZWq5sg0szh5lYONbdRsa2O8uFFSdnnxIkTj+m0n3zySV588UWeffZZcnNzeeGFF1i9ejW//vWvj3ndDTfcwLnnnktBQXSrtCUwZB8r+rCsGHO49957uffee61uhkgTFdvq+LjmMDvrmpj52EoqttUl/Rhr1qzhT3/6E3/84x/Jzc1N6r6VUkndn0h/VvRhKb1yUErdDXzT/+tSrfVt/u1OYBnwC63126lsA8Ddd98NwF133ZXqQ4kMsLK6Fq//TlRrm5eV1bVJu3oAqK2t5eabb+a+++5j2LBhnT7/D3/4A/PmzTtyW2nBggUUFkqdZ3GUFX1YyoKDUups4BygDPABy5RSlwD/Bp4ATk7VsYWIZMqovtht4PWBM8fOlFF9k7Zvr9fLrbfeyvTp0/nqV78a1Wtiva0kRFdI5W2l3cAtWmu31roVqAKGAd8FHgRWpfDYohsqyHWE3G6PsUJW+fAixhT3ZGhRDxbMmpLUq4aHH34Yt9vNLbfckrR9CmGFlF05aK03BX5WSp2Aub10qtb6E/+2G+PZ78aNG+NuU0VFRdyv7czo+sMAbE7BMVLZ7nQV6pznf61/h205OQ5GHpdDQ0NsNaALXHYKXC7G9HPF/NpwVq5cybPPPstf//pXWlpaaGlpOebxlpYW2trajjleW1vbkefF0o5Iz3W73RnzmcmUdiZTIufclf9eKZ+tpJQaBywF/isQGBJRWloa9wBfeXl5oocP7189U3KMioqK1LY7DcVyzlVVVXHdjrHbzVVIMm/lPPXUU3i9XubMmXPM9iuuuIJvfetb5ObmkpOTc8wxc3JyjnyeY2lLpOe6XC7Gjx8fY+u7nny2YxfPa1taWuL6Up3qAelTgUXAjVrrZzp7fqoke7ZIOGv3ruXqp07ssL0gp0DqCKSZhd+bmvR9zp8/P+Ljl156KZdeeukx2wLTWmO5atBax944kdG6qg8LlsoB6RLgJeByrfXyVB0nGs3NzV1yHI8v9CKtUJXHhAinubmZyy+/PORjN9xwA1/5yle6uEXCal3VhwVL5ZXDrUAe8Nugedl/1lr/OYXHDOmaa64BzKIkIdJdXl6eLHQTx7CiD0vlgPQcYE6Ex89M1bHbe+qppwAJDkKIzGRFH5YVK6SFEELEJityK3UmVLZPMPPqN91zXvgX7tkDlZWwcye8e5CR9UXc2dPB5/ltfFrYzMYSJ7W95Z84Lc2bbv6+dqm17RAiTUnPRfhsn+G2U1cHS5bAvn1QVgZnnw2lDlpb91PnctG39hDDdu7lK+vWscm2l3fG90xh64UQIvkkOMRq0yZYuhROOw1mzgR/PhzqXAzKHcz3ymYefa7HwyWrVnHH+++b140bZ02bRZfoqpTdQnSFrAgOPXsm6Zv7pk2wbBlcdRUUF3f+fIcDpk2DUaNgwQKzTQJEemg5BM0HYcdqKDklKbvsqpTdIvskrQ+LQVYEh0OHDiW+k7o6c8UQbWAIVlxsrjLmz4fBg6Eoebl8RBx2rIY9G8Hnhae+Ble/nLQAERBI2f3cc88lfQHTlClTKC0tZd++fTz//PM4nc6k7l+kn6T0YTHKitlKX/nKVxJfOLRkibmVFGtgCCguhlNPNfsR1tr6ngkMAB63+T2J4knZfcUVVzBjxgxmzJhBfX19xOfX1dUxe/ZsFi9eLIEhSySlD4tRVlw5LF8eeYF2Qa4j7GwlwMxK2rfPfPtPxNSpsHKl2d/AgYntS8RvxOlgs5sA4XCZ35Okq1J2Z0LuJJE8nfVhqZAVwaEzEaergpmuWlZ2dPC5M5X+8YWydsHEboeTT4Z16+Dcc2NvqEiOklNgYKkZc/j6Y0m9pdRVKbvz8vJSun8hJDhEY+dOM101kv9Uw/7N0O+EyM8bORLeeCN5bRPxye1l/iQxMKxYsYLnnnuORYsWWVIQXohkkk9wNPbujTzWsGM1rPg9eNvgk3/C1B9Bn1Ghn1tcbPYnup1HHnkEj8fD7Nmzj9keSNktRCaR4BANtxsizTjZ+p4JDPjMAOf+zeGDg8sFra0paaawlqTsFt1JVgSHvn0TrBHsckFLC4S7zzvidLDngNdjBjj7jQ6/L7cbZIaJ9dI4bYak7BbtJdyHxSErgsP+/fsT28GAAVBTAyNGhH685BQ4dY65Ypj2E9gfoeBdTY3ZnxBhSMpu0V7CfVgcsmKdw4QJE5gwYUL8Oxg6FLZsifycPqNg9HmdD3Bu2QIlJfG3RQiRdRLuw+KQFcFh/fr1rF+/Pv4dlJWZ6ayeMIn4ouXxwNq10MVvshAisyXch8UhK4JDwgYOhP79YdWqxPazcqW5pSQL4Cx37bJruXbZtVY3Q4i0JcEhWhdeCO+/b8YM4lFTAytWmP0IIUSak+AQraIimD7dZFeNNUDU1JjXTZ8uSfeEEBkhK2YrJU0g3fb8+SaJ3tSpR1Nq1B4EvR0qHoNNK6CtDQZugeZmk0vpyislXbcQImNkRXCIJjNm1MaNM2m3lywxYwgnnADV1VCjTd6kk8+EySdAsxvqc6CpyQSQigrYvt3cVpKrB8vVu+s53HqYdXvXMWFA8iYIzJ07l9dee42ioiL69+/PWWed1WHhmxCxSmofFqWsCA7btm1L7g6LiuA734F33oHnnjOzmUrKTabP3bvh4Gfg8sG0y6Gw0Lzmq181A9qPPmpuL8lVhGXW7V2HrtP48DH79dk8es6jSQkQy5cvp6KigiVLltDU1MQll1zCWWedlYQWi2yX9D4sClkRHIYPHw4k+R940yZYswbuuMPkSzqSifV0qNxufg4EBpCqcGlkzZ41+PAB0OptZc2eNUkJDh988AHnn38+LpcLl8vF2Z0laxQiSinpwzqRFcFh+/btyd2hVIXLaBMHTsSGDR8+nHYnEwdOTMp+7XY7Xq83KfsSIljS+7AoyGyleEhVuIw2YcAEVJFiSOGQpN1SApg2bRqvv/46breb+vp63n77bWw2W1L2LURXy4orh6SSqnDdQqGrkEJXYVIHo88880wqKyu55JJL6N27NwMGDEh6/WghuopcOcQq1qpw4QRXhRPdQmVlJSNGjGDp0qU89dRTNDY2MmpUmNTtQqS5lF45KKXuBr7p/3Wp1vo2pdT1wA2AD1gDfE9r7U5lO5Iq1qpwkRLxSVU4y8w7b17S9zly5Egefvhh5s2bh8/n4+KLL2bMmDFJP44QXSFlwUEpdTZwDlCGCQTLlFI/BWYB5cBh4EngR8BDqWoHJLkYe6xV4a5+OXyAkKpw3cpxxx3H448/bnUzRDeU1D4sSqm8ctgN3BK4KlBKVQF5wA+11of82zYAKV/dsS6Zt25irQq39b3wwUGqwgkhopDUPixKKQsOWutNgZ+VUidgbi+dqrX+xL+tP/Bj4JpY9rtx48aY2xKYb/5GEm7hDDt4kNYtWzpUhevrNnfGDhWOoZ+/KpzP7qS2cAytn38eemfNzTgPHGB7RUXEY1Z08nh3FO055+TkxFRiM50l6zzcbnfGfGYypZ3JFM85J7MPi1bKZysppcYBS4H/CgoMQ4BXgce11m/Hsr/S0tKYZ4AcOHAAgPLy8pheF1JlpRlMHjz42O17XAD0L5sO/zFV4WzTfkL/SGMOW7fCuHH0j9CuioqK5LQ7g8RyzlVVVRQUFKS4RanX0NCQtPNwuVyW3IaIlXy2o5dIH9bS0hLXl+qUzlZSSp0KvAn8TGv9lH/bGOAD4Cmt9S9SefyUkKpw3cK271zFtu9cZXUzhEhbKQsOSqkS4CXg21rrZ/zbegKvA3doreem6tgpJVXhRAQbNmzg9ttvt7oZQiQslbeVbsUMQP9WKRXYthAYCNyilLrFv+1lrfVdKWxHcgVXhZs2Lf79SFW4bunEE0/kxBNPtLoZQiQslQPSc4A5IR76VaqO2WUuvNBkVx01Kr4UGoGqcLNnJ79tIiqe+nq8hw7RWFlJfllZ0va7atUqHn74YZ5++umk7VMIK2RF+oykp00Orgo3c2ZsAUKqwlmusbKSlo8/Bp+P7ddex7B5TyQ1QAiRbFakfs+K4PDmm28mf6ftq8LleSOn1PB4zK2kFSuknoPFGld/BD6TstvX2krj6o8kOIi0lpI+rBNZERx69eoFwKFDh5K74+CqcCv+BmUToGgrtLjBmWNKhNbUmFlJa9eaMYbZs+WKwWL5p0wCmw18PmxOp/ldiDSWsj4sgqwIDocPH07dzgNV4Wqfgk/fAccX4N0nwQt8cYcJCCUlpoa0DD6nhfyyMnLHjMF76BCDf/OgXDWItJfSPiyMrAgOXaJ3rvlz7SxwvGi2Xftza9skwnIUFuIoLJTAIEQYEhyESKLJkyczefJkq5shRMKknoMQQogO5MpBZKXhT8+3uglCpLWsCA5XX3211U0QQoi4WdGHZUVwePLJJ61ughBCxM2KPiwrxhzy8vLIa1d/QQghMoUVfVhWXDm0tLRY3QSRZl6cuxaAS2452eKWCNE5K/qwrLhyEEIIEZusuHLoctcutboFQgiRELlyEEII0YEEB5GV3E1tHP5PMzXVB5O631WrVnHZZZdx6aWX8tOf/jSp+xaiK2XFbaV77rnH6iaINFJTfZD9u+rBB4sfqmTGTWUUj+qdtP1v3bqVt956i549eyZtnyK7WdGHZUVwuOuuzKlCKlJv1+Y6MOUc8Hi87Npcl9TgMHLkSAkMIqms6MOy4raS0+nE6XRa3QyRJoaMLgKb+dnhsJvfk0jW1Ihks6IPy4orh7a2NqubINJI8aje9BtSSEtTG+d8d1xSrxqESAUr+rCsCA5CtOfqkYOrR44EBiHCkOAgRBJJPQfRXUhwSBZZ+JZRJG2GEJFlxYC0EEKI2GTFlcMrr7xidROEECJuVvRhWREczj//fKubIFLI5/Nhs9msbkZa8Pl8VjdBpIAVfVhKg4NS6m7gm/5fl2qtb1NK/QD4MWam+VLgNq11Sj/RDocDAI/Hk8rDCAvk5eVRW1tL3759sz5A+Hw+amtrZZ1FN2RFH5ay4KCUOhs4ByjDrEddppS6CfghMAFoBt4Fvgq8nqp2AHi93lTuXlho6NCh7Ny5k3379lndlIS43W5cLlfC+8nLy2Po0KFJaJFIJ1b0Yam8ctgN3KK1dgMopaoAL/BFrXWrUqov0Bs4kMI2iG7O6XQycuRIq5uRsIqKCsaPH291M4Q4ImXBQWu9KfCzUuoEzO2lU/2BYTbwG2A1sC5VbRBCCBEfW6oHsJRS4zBjC3drrZ8K2p4DzAN2aK1/3tl+KioqRgBb4mnDxIkTAVizZk08LxdCCEslqQ8bWV5evjXaJ6d6QPpUYBFwo9b6GaVUCTBMa71Ca92mlHoG+EGUu3MAjB49OuZ7s4MGDQKgtLQ0ptelg40bN2ZkuxMh55wd5Jyjl0gf5na72bx5M/j70Gil7MrBHwjWApdrrZf7t5UCSzAD0geBJwCttf51Z/urqKg4DXgvJY0VQoju7/Ty8vL3o31yKq8cbgXygN8qpQLb/gz8CvgAaMN09nOj3N9HwOmYgW6ZkyqEENFxAIMwfWjUUj7mIIQQIvNIbiUhhBAdSHAQQgjRgQQHIYQQHUhwEEII0YEEByGEEB1IcBBCCNGBBAchhBAdpGWxH6XUb4B+WutrlFInA48ALmAHcKXW+oBS6jhgATAK2Ad8U2tdY1mjE6CUegsYALT6N30P+AJwB+AEfqe1/j//c88Gfgv0ABZqre/o+hYnJsz5nke72h/+52b8+ULoc9Zar/I/duTz7v99AvAY0AuT1v77Wuu2Lm90gsK8zwOAu4EC4HWt9Rz/c7vt+wycAPy3//dXtda3+p+b1u9z2l05KKW+AlwdtOn3wF1a6/GAxqy8BrgPeE9rPRZ41P+8jKOUsgGjgfFa6wla6wnATuCXwGmYVCPXK6W+qJTqgUk5MgMYC0xSSmVUmbsw59uTo7U/JgDlSqlLusP5QuhzDgoM7T/vAH8Ffqy1Ho0pijW7SxucBGHe532YLAkXAycBJyulzu/O7zOwAfgD8CVgPHC6PxBCmr/PaRUclFJ9MJ3i/UGbHZjICpAPNPl/no65cgD4O3C+UsrZFe1MskBukdeVUuuVUj8GzgaWa63/o7VuAJ4HLgNOAT7RWm/xf8P4K/ANS1odv1Dne6T2h9a6FagChtE9zhdCn3PIz7tSajjQQ2u90r/pSbrPOV+CuSrY6X+fLwdW0b3fZwemny3A3AVwAk2Z8D6n222lR4DbgZKgbTdj/rF/BzQAk/3bB2M6FfwZXg8B/YHPu665SVEEvAn8BPPBeRtYiP/c/HZj/gMNDrE908p+hTpfrbX+Jxxb+wOYSOafL4Q4Z6WUBq6n4+e9O7zHEPp9dmD+L7+MCf5LgDvp3uesMef4MdAIvIPJLTeFND/ntAkOSqlZmNoObyqlrvFv6wE8DpyttV6tlLoZmI+5amhfMNiGqTSXUbTWHwIfBn5XSj2Oufd6X9DTAudmx5Rcbb89Y4Q53wuAfwbV/vgvrfUnSqlTyPDzhbDn/CwwL/jz7pfx7zGEPeeTMFfFZwL1wMuYW2qtdN9z/iEwAhiOyUT9V8yt8Q9I83NOm+CAucQcpJRaB/QBCjH/oE1a69X+5zwC/ML/8y6gGNjpLxzUE6jt2iYnTil1GpCrtX7Tv8kGbMVkUQwoxlwR7QyzPWOEOd/W9rU//I9l/PlC2HNeA5wT/HlXSj0E/I7ue877gTe01vv8z3kRc0W8kO57zhcDc7XWe/3PeRITMJ4lzc85bYKD1vqrgZ/936TOBG4CPlZKKa21xgxYBdLOvgJchblfezlmcLqVzHMccK9SahrmUvRq4Ergr0qp/phbaV/H3IL4F6CUUsdjquJ9GzOQl0lCne+PgZcIqv3ht4rMP18Ifc7f93/TPPJ511rf5P+9WSl1qtZ6BfAd4FVrmp2QUOd8A/CIf6bhYeB8zPvend/n/wEuUUoVYG4rXQR8pLXelu7vc1oNSLenta4DrgGeVUr9C7gOuNb/8J3AFKXUJkwk/pEljUyQ1noJ5lZKJVABPOH/sNwOvIWpsf03rfVqrXUz5t9jEfBvzH3M561od7xCnS9msD1Q+2Od/8/3u8P5Qtj3+MMIL5kJPKSU+hhzBf2H1LcyucKc8zvAA8D7mPdzG+bWWnd+n+/BTJipwHy5cwKB4mZp/T5LPQchhBAdpPWVgxBCCGtIcBBCCNGBBAchhBAdSHAQQgjRQdpMZRVCiEzlX6fzECZBaC1wndZ6W5jn9sTMQvyu1vpt/7Zq4FDQ0y7SWu9QSv0XZpamD3hMa/3bTtqhMOvBioAa4Ar/rM+YyZWDyHpKqYlKqQ5TJ5VSS9qtXo51v28rpS5LqHEiUywAZvmT7S0g8rTUhzGdNwBKqb6AOygp4wR/YDgeM03/ZGASMMe/LSR/4r+XgV/7E5VWAj+L94TkykFkPa31GsxaCyGiopQ6E7jGX1YgF7hDa/0v/8P/wuRXCvW6yzELAP8VtHkSYFNKrcCkLP+V1vo5zJd3F2YNkM3/p9W/n6uAG/3PqcCs8xoHNGitl/n3ez9mYV5cJDiIrOf/j/4wJm34U5hEcNswefmjef2dwLeANmAzJg1zoLbIJUqpn2EyCi/QWv/Sn+7lj5jkgq1ANXCt1ro+aScluozWugWTMwmllB2zKvql9s9TSg3DdOhncexq6DzgNeCnwEDgPaXURq11lVLq75h0OnbMbaVt/hxks4FpWutmpdSvMPmaPgFq/DmdyjDZjUMGqWjIbSUhjvo/YKXWehwm1cOYzl6glLoWkwZiktb6JGAjJv1yQC9MBs4pwJX+OgVTMelhxmutyzHB4aTknYZIFX+dkXWYIj1f86/mn+d/zIW5pZTDsWUHAkHjccwXh6bgx7TWL2mtb9BaN2mttwIvYPJunQeUA0P8fyYppb4JfBlTQGilvy0zMJ/VHMzn6k9a65Mxn6uIYxSRyJWDEEedjb+YlNb6U6XU8k6eDyYwzPPX3QBTdOp2f0cB5tteG3DIP67xVeBewAOsUkq9BiwKSi4p0pjW+kXgxeDbSgBKqULM/f5aYEaIPG9j/H8eN2PGHA88ppSajan1UOO/vQlHbx9dhPls1PuP8XdM0aDNwLNa6xuCjp2DCSSfBO3n7ySQhkSuHIQ4ysexqeCjKdno4NjUy3bMf9TAfjztHmvVWh/AVAW71f/4QqXUD+NttEgLfwU+xSSPbGn/oNb631rrkqAKcWswA9hvYVJ6362UsiulBgJfw+RoWg9cqJRyKFPI7HxM4tG3MbcrB/gHof+EuV31AdBfKTXef9iLMOMRcZHgIMRRyzDZbwP3h78c5Wuu82fdBHM76t2gDuIqpZRNKVWEKWK0TCl1IaYozAda6//B1CiZlLzTEKmmtX476KqhDHNr51Rgrf9W0yv+x76vlLq3k939GVPsZyOmGNBP/dNgH8MkItyEmXm0CXhKa70euAdY7t/mwMxQasJU23vUn5D0LOCWeM9REu+JrBc0IP1lYB7wBUwtiRzMf8YnI7w2MAD5dcyXrU+BH2itdyql3sYMEpZjZqH8SWv9B6WUAzPV8SuYojd1wGz//WYh0oIEB5H1lFJfA27XWk/u9MlCZAkZkBZZTZki8HcAd0V4znuYSoOhnK61PpyKtglhJblyEEII0YEMSAshhOhAgoMQQogOJDgIIYToQIKDEEKIDiQ4CCGE6ECCgxBCiA7+P0XHfG/rktzTAAAAAElFTkSuQmCC\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, axyz = plt.subplots(1, 3, figsize=(18, 6), sharey=True, squeeze=True)\n\nfor (ax, band) in zip(axyz, 'gri'):\n    e = expected.query('filtercode == \"z{}\"'.format(band))\n    f = df.query('filter == \"ZTF_{}\"'.format(band))\n    m = lcs.query('filter == \"{}\" and magpsf < 98'.format(band))\n    ax.scatter(e.obsjd, e.maglimit, marker='s', color='k', label='mag lim')\n    ax.errorbar(f.jd.values, f.mag.values, yerr=f.mag_err, fmt='o', color='b', ms=15, markerfacecolor='w', alpha=0.4,\n                label='FPS')\n    ax.errorbar(m.jdobs, m.magpsf, m.sigmamagpsf, fmt='d', color='r', label='Marshall')\n    ax.set_title('{} band'.format(band))\nax.invert_yaxis()\naxyz[-1].legend(loc='lower left')\nfor ax in axyz:\n    ax.set_xlabel('JD')\naxyz[0].set_ylabel('mag')\nfig.savefig('fps_marshall_ztfquery_ZTF18aaermez.pdf')",
   "execution_count": 37,
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 1296x432 with 3 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCQAAAGECAYAAADusvFAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X983XV99//HyY82JCSIQANNC1QqH9oEaw2IyGVvoC3W1k4LCNuUDYrj8hJRdAzl2i79uunKnHUbuO8Y+1IvRaYTxCFGwSJiJ9dFhayDNu3eIJZr/cHSjl2YNmmbpOd8/3gnbdr0R5Ke3+dxv93OjZzP+XzO53VC+snJ87zfr3cqk8kgSZIkSZKUT1WFLkCSJEmSJFUeAwlJkiRJkpR3BhKSJEmSJCnvDCQkSZIkSVLeGUhIkiRJkqS8M5CQJEmSJEl5ZyChspUkyaVJkqzP4fPfmiTJ/8zV80tSKUuS5OUkSS7I0XOfmiSJ65ZLEpAkyQVJkjx4mO2+F1bRqyl0AZIkSZKkiQkhPAtcVeg6pIkwkFBJSpLkM8ANwE5gNfD+EMLZh9n1xKHEeCbwGnBjCOGFJEnOBf4aaATOAP4FuCaEsCdJkj3AHcDlQ499KYTwN0mS1AJ3AguA7UA38OscvkxJKlpJklwK/BXQC5wIXBhC2HvIbjclSTIHmAysCCGsTJKkCvgL4G3Ea3AK+HAI4amhT9p6gPOB6cDzwO+EEHYlSXIF8EWgD3gm169PkkrF0PX4qyGEtsM87HthFTWnbKjkJEnybuA64EKgnXghPZLpwFdCCG8G/h64b2j77wFfDyG8jXiBngEsHnpsMvAfIYS3E9Pmv0iSpA74KHAuMJt4IT4ziy9LkkpRG/BbIYQ3HSaMANgdQngL8Zq5PEmSVuAiYCpwcQhhNvB14DMjjmkHFgKzgLOBDyRJ0gysBK4MIbQD/ydXL0iSyozvhVXUDCRUihYBD4QQXgshZIjp7pE8H0L4X0Nf/0/ggiRJTgI+DexIkuQ24G+Ib45PHHHcw0P//WfiRbkBmA/8fQihP4TQC9yfrRckSSVqcwjhaOHA3wKEELYBPwbeFUL438AfAf81SZIvE9/sjrz+PhpC2BtCGADWAa8H/guwLoSwYeTzSpKOyffCKmoGEipFg8QhvsP2HWXfQx/LAAPAt4AbiZ+y/QXxYjvyOXcDDAUejHhs5D6D46paksrPrmM8PvIaXAUMJEmyGOgY2vYwcDeHuf4OyeD1V5KOh++FVdQMJFSKOoArh9JdiL0kjtRtfU6SJG8e+vq/Aj8PIfQB7wb+OITwD0OPXQRUH+O8PwJ+J0mSuqFha9dM+BVIUmW4DiBJkjOJn6z9hDjM95EQwt8AzwLv59jX39VA61A/iv3PK0k6Jt8Lq6jZ1FIlJ4TwRJIkfwf87yRJ+oAuYpOzw9kIfC5JkjcQm+/87tD2/w58L0mSXmIznp8R588dzd8O7bMeeBV48bheiCSVv7okSf4ZmATcPNRI7W7gW0mSrCO+D/kxMWQ+4ockIYQdSZL8NnB/kiT9xGu2JOnYfC+sopbKZFzGW6VlaF37t4cQ7hy6/yngohCCKa0kSZIklQhHSKgUvQB8OkmSG4lTNf6NOAdOkiRJklQicj5CIkmSJuB/Ae8lLhHzpyMebgHWhBDem9MiJEmSJElSUcnpCIkkSS4C/o64Xi0hhB8CPxx67HTgKeCTuaxBkiRJkiQVn1yvsvF7wE3AtsM89ufA3SEEm6FIkiRJklRhcjpCIoTwYYAkSQ7aniTJG4FLgQ+P9bk6OzsnAxcCrzB6PV1JKpRq4Azgmfb29r2FLibXvBZLKlIVcy32OiypSE3oOlyoppY3Av9vCGE8vzAuBP4pR/VI0vF6B/DzQheRB16LJRWzSrgWex2WVMzGdR0uVCDxfuDycR7zCsC5557LpEmTsl8RsH79etra2nLy3IXiayoNvqbid6TX09/fzwsvvABD16gKkPNr8USV289ctvn9OTq/P0dX7N+fCrsW+574CEq1bijd2q07/4q19oleh/MeSCRJcipwQghh0zgP3QcwadIkJk+enP3ChuTyuQvF11QafE3F7xivp1KGzeblWjxRxVhTMfH7c3R+f46uRL4/lXAt9j3xUZRq3VC6tVt3/hV57eO6Due6qeXhvAHYUoDzSpIkSZKkIpGXERIhhLNHfP0L4G35OK8kSZIkSSpOhRghIUmSJEmSKpyBhCRJkiRJyjsDCUmSJEmSlHcGEpIkSZIkKe8MJCRJkiRJUt4ZSEiSJEmSpLwzkJAkSZIkSXlnICFJkiRJkvLOQEKSJEmSJOWdgYQkScWqqwva2uJ/JUmSyoyBhCRJxai3F9rbYxixeHG8L0mSVEYMJCRJKkbLlkF/f/y6uxtuuKGw9UiSJGWZgYQkScVm5Uro6IBMJt7fswceeSRulyRJKhMGEpIkFZvbbx89RaOvL26XJEkqEwYSkiQVm+XLoaHh4G319XDHHYWpR5IkKQcMJCRJKjbLlsVGlnV18X5dHSxZAtdfX9i6JEmSsshAQpKkYrRyJUyZAqkUNDfDvfcWuiJJkqSsMpCQJKkYNTTAnXfCpEnwV381egqHJElSiTOQkCSpGPX2wjXXwN698IlPjG5yKUmSVOIMJCRJKkbLlkF/f/y6uxtuuKGw9UiSJGWZgYQkScVm5Uro6IBMJt7fswceeSRulyRJKhMGEpIkFZvbbx89RaOvL26XJEkqEwYSkiQVm+XLRzexrK+HO+4oTD2SJEk5YCAhSVKxWbYMFi+OS34C1NXBkiVw/fWFrUuSJCmLDCQkSSpGK1fGJT8Bmpvh3nsLW48kSVKWGUhIklSMGhqgsxNaW2ODy0OncEiSJJW4mkIXIEmSjqC1FdavL3QVkiRJOWEgIUmSJElj0N0N27fDwADU1sKUKXFWnaSJMZCQJEmSpCNIpyEE6OqCfftg+vQYRuzZAxs2QHV1HNCWJFDlhHhpXAwkJEmSJOkwBgZg1aoYSlx8MUydemABJIBMBrZuhbVrYdMmWLAghhWSxsYMT5IkSZIOkU7HMKK+HhYtgpaWoTDiy1+ON+L9adPiSs319QfCC0ljYyAhSZIkSYcIIYYL8+YdMhXj3/89BhJdXfs3VVXF/fbti8dJGhsDCUmSJEk6RFcXzJ17SBjR2wv33hs7Wy5eHO8PqaqK+4/IKSQdg4GEJEmSJI3Q3R1HO0ydesgDy5bBrl2xeUR3N9xww0EPt7TA4GB8SNKxGUhIkiRJ0gjbt8fVNEY2sGTlSujoiIkDxGU2Hnkkbh+SSsXjduzIb71SqTKQkCRJkspUU1MTqVRq1K2pqanQpRW1gYHDrJZx++0HTdEAoK8vbh9h0iTo789tfVK5MJCQJEmSytTOnTvHtV1RbW0MJQ6yfDk0NBy8rb4e7rjjoE39/TGUkHRsBhKSJEmSNMKUKbB5c2wVsd+yZbGRZU1NvF9XB0uWwPXX798lk4nHnXZafuuVSpWBhCRJkiSN0NwM1dWwbdshD6xcCSeeeGCne+896OGtW+Poiubm/NQplToDCUmSJEnF59JL461AWlth7VpIp0dsbGiIK2s0N8cGlyOmcKTTcf/Zs/Nfq1SqDCQkSZIk6RBJAlVVsHr1IaHE6afDrbfGxGJIOh33q6mJx0kam5pCFyBJkiQpNxobGw/bwLKxsbEA1ZSWqipYsABWrYqDIebOhZYWSN166/59Mpk4TWPt2hhGzJ8fj5M0NgYSkiRJUpnq6ekpdAklrbYWFi6EEODpp2FwEKZPP7C05+bNcZ/Zsw+MqJA0dgYSkqSK1NTUdMRPDX0DXwKG55U/+WQhq5CUI+vXQ9jxDvakJ1H33fjHfltbYWqpqoJZs+Ktuxt27IhhxEknwcyZNrCUjoeBhCSpIh0ujDjadklSbg0OwqOPwlNPxVEHrUs/QEMd7B2E73wH7r8fLrkkjlioKdBfMc3NBhBSNuX0n3KSJE3A/wLeG0J4OUmSy4E/B6qBfwY+HELoz2UNkiRJkopbXx985StxRc0PfSiORhj4m6ep/dKfUvXDDtIfaGXDBvjxj+Ff/gU+9Smory901ZKOV85mOSVJchHwc+DcEZvvBX4zhNAG1AO/k6vzS5IkjVLgZQQljTY4GMOI6dPhYx+Li1dU7e5l0uf/B6nN/waLF1O1u5e2Nvj4x+N+X/lKPE5Sactl25XfA24Cto3YVg00JUlSDdQBu3N4fkmSJElF7tFH48iID35wxFSMZcvg16+RymRi44YbbgDi4x/8YBwd8eijhatZUnbkLJAIIXw4hPBPh2z+KPAkMaQ4FXgwV+eXJEmSFD3ySLwVo6eeistr7g8jVq6Ejg5S/UMzu/fsicWvXAnE/S6/PB4nqbSlMplMTk+QJMnLwKXAHmIYsRjYBHwFqA0h3DSW5+ns7Dx76DhJKkYz2tvbXy50EblWTtfiefPm0dfXN2p7fX09q1evLkBFGo9zb7wRgBfuuSerx5137bXUbdrEv3796+w555zjK1KFUPbX4olch199tYZVq04mnU7x1rf2cPLJg5xySnHMd/jlL+v4xS/O5DOf6d+/ZObJ551H1X/8x6h906eeyv/913+NX6fhT/90Em97278xc+aefJYs6ejGdR3OZ3/adwDrQwgvASRJ8nfAd8b7JG1tbUyePDnbtQHQ2dlJe3t7Tp67UHxNpcHXVPyO9Hr27t3L+vXrC1BRYeXyWjxR4/2Z6+3tzWE1xafc/k3S2Agw/td0hOM6OztpP+88+NWvYO9eWm+7Dbq6oKEhK+WWumL/+anEa/GxrsPpNIQQf4z37Yt9GWprY/+FzZuhpyduSxL2BwHD8vn/++WX4a1vhdNOG7Hxz/4sNosYeZ2ur6fqS1/ilFNO2b/prW+FSZNaGS612H9Oj6ZUa7fu/CvW2id6Hc5lD4lDrQfemiTJ8EI57wOeyeP5JUlSuejthWeeiX9tZcuyZTA8RHzEnHWp1AwMxP4KL70EF18MV18Nc+dCWxu8/e3x/tveBr/8ZdxvYKBwte7ZA6NylWXLYPFiMpMmxft1dbBkCVx//UG7nXBCXJ1DUunKWyARQtgI/A/gp0mSPA9cANyar/NLkqQy0dsL69bFv0QWLz74U9QJOuX734eODhieynrInHWpVKTTsGpVbPq4aBG0tEAqdfA+qRRMmxb/+dTXx/3T6cLUW1cHe/ce5oGVK8k0vY5MKgXNzXDvvaN22b3bpT+lUpfzQCKEcHYI4eWhr78eQpgdQnhTCOGKEMLoyWGSJElHczwjGY4wsqLlrrtGBxt9fXD77cdZrJRfIcRwYd680VMxDlVVFffbty8eVwhJEv85jgpEGhrY/Yd/QnramTEsPGT6VDoNGzfCG9+Yv1olZV8+p2xIkiQdn6Hu+xMayXCUkRVbb755dL+I+nq4444sFi/lXldXnJ5xrDBiWFVV3D+bs5/Go60tThnZuHH0Y+kzz6bvL++JzS4OsWFDDFLa2vJQpKScMZCQJEml4/bbJz6S4SgjK179jd+IIcXw2PYjzFmXill3d/wjferU8R3X0gKDg/H4QrjkkjhtZPCQhT9qa+PtUIOD8OMfx/4YkkqbgYQkSSody5dPbCTDWEZWrFwZlyA4ypx1qZht337gR3g8Uql43I4duanrWBYuhF274P77R4cShxocjPvt3h2Pk1TaDCQkSVLpGOq+P+6RDGMZWdHQAD/8Icyefdg561KxGxg4/IiCsZg06cAAonyrqYFPfSouR3rnnbB+fewRUTc53iDeX78+Pr5lC3zyk/E4SaXNf8aSJGVRd3f8lHL4D4MpU+KH7cqilSvh4Ydja/6xjmRYvhw+/vGDQ4nDjaxobY1/9UglqLY2Dv6ZiP5+OOmk7NYzHvX18JnPxGVI778/vpZZ+5o54YQUu1+NPSb27YvTNBYuNIyQyoX/lCVJOk7pdOxQ39UV3zBPn37gD4MNG6C6Ov6dmyRjbzSno2hogPPPj9/csY5kWLYMHnsMHnggTtuwR4TK0JQp8Z9FJjO+aRuZTBydMHNm7mobi5oaeO974239enjxNx+gd99k6v/0TVx9tQ0spXJkICFJ0nEYGIjN2NLp+Mnd1KmQ+oNb4ZvfhJ/8hMzFrWzdCmvXwqZNsGDBxIdUa4SGBrjwwsN23z+iiYyskEpIc3MMQLdti40qx2rr1nhdKqbRXG1t0Lb+C4UuQxqlqamJnTt3jtre2NhIT09PASoqbX5OI0nSBKXTMYyor4dFi+IfAKm+3viH7vbtsHgxqb5epk2LbQ/q6w+EFyqA4ZEV9fX2iFDZam2NAeih15mXX463Q6XTcf/Zs/NRnTQBl14ab0XicGHE0bbr6AwkJEmaoBDim/l580ZMxVi2LLaLz2QOWlqyqirut29fPE4FMpGRFVIJGZ4atnr1scPPdDruV1MTj5OkfHPKhiRJE9TVFadp7A8jhpeWHF63buTSksuWUVUFc+fC00/DrFkFK1tSGauqilPDVq2Kl6Pvfe+zVFf/O6nUbwDQ0fF9MhnYt+90li79Y2pqYP58+9tIKgwvPZIkTUB3dxztMHXqiI1jWFqypSXmFd3d+alTUuWprY0rUcycCQ8/3M2aNWezeXMT3d1NbN58JmvWnM3DD/87M2fCu99tXxtJheMICUmSJmD79riaxkGd7MewtGQqFY/bsaO4GshJKi9VVXEk1n/+5z38539OYd26k4Fq4NvADmA7s2bdU9giJVU8R0hIkjQBAwOH+VRx2TJYvJhMzVDef4SlJSdNgv7+/NQpSVEXsAEYx3qgkkZpbGwc13YdnSMkJEmagNra2CJilJUryTz2Y/j1a6SOsLRkfz+cdFLuayxrTz5Z6AqkopZOxwa6r3/9jbS0nMm0aTOHwtBT2bIlxdat/8bGjQeaYEpFq7cXNmyIjZuKoCGxS3tml4GEJEkTMGVKfH+UyRwybaOhgRfu/TlnffoaTnj4H0YtLZnJwObNcW63JOXCwMCBJYYXLjyJ+vr/Qyo1nIL+Gy0t0Nv7On75S9i0KTbBtI+EilJvL6xbB3v3xvWzu7pcsrnMmIdKkjQBzc1QXQ3bto1+7LwrWznhl+sP+0nO1q3xjb/9IwrkyScdXaGylk7HMKK+HhYtglNO2XNwaEoMUU89dQ+LF8f9hsMLqegsW3ZgjuOIpbRVPhwhIUnSBLW2wtq1cMYZYxvynE7H/WfPzn1tkipTCPFaM29evC7deeedADz3XHx8zpyD9583Ly4PGoLLEavIrFxJ+gcdVGUy8f6ePaS//whVQ0tpqzw4QkKSpAkannu9evWxP11Mp+N+NTXxOEnKha4umDt37H0hqqri/l1dua1LGqt0GjZuhD2fup2qvoOX0q7a3ceeT93Oxo2O6ikXBhKSJE1QVVWce93XFz9h3LIl9ogYKZOJ2zs6YPdumD/fBnIqMpdeGm8qed3dsG8fTJ06vuNaWmBwMB4vFdLAADz6KLz0Euz+7HIyh/Zhqq+n77N38Mtfxv0GBgpUqLLGKRuSJB2H2lpYuDAOd3766fimfvr0A0t7bt4c95k92272knJr+/Z4/Tm0Z8SxpFLxuB077G+jwhnZ/2TePKh67zJY8xiZhx4iNTgIdXWklizh9Z+6nsVDow5XrYq/g/3dWroMJCRJOk5VVXHu9axZ8RPGHTsOLO05c6Zv8CXlx8DAxFfLGA5RpUI5tP8JcMSltKuq7H9SLsySpFJ24YVx6SMnfkpFo7kZ2trgLW+J/zWMUNHr7YVnnvF3SRmorZ34EPb+/hhKSIVy2P4nDQ3s/uANpE9rjunDiCkc9j8pDwYSUqkaXpe5ry+uy9zbe+xjJEkayd8lZWXKlDhN7NBeNseSycTjTjstN3VJx3K0/ifpKafTd9Oth11K2/4npc9AQipVrsssSTpe/i4pK83NUF0N27aN77itW+PoCkd0qVCy0f9EWZTHZscGElIpWrkyDlsbsS4zjzwSt0vKP1cpUCnyd0lZam2FtWvHviRiOh33nz07t3VJR2P/k8plICGVottvHz2stq8vbpeUf87BVynyd0lZGl7NZ/XqY4cS6aGVCmpq4nFSodj/pHIZSEilaPnyg5r6AHGNpDvuKEw9UiVzDr5Klb9LylJVFSxYEC9JHR2wZcvonhKZTNze0QG7d8P8+S6bqMKy/0nlctlPqRQtWwaPPQbf+U68P3kyLFkC119//M/9utfF/7722vE/l1QJDjcH/9vfLmxN0lgM/y554IH4rr6uLnu/S1RQtbWwcGFcDvHpp2PTv92746fIu3bFP+Bqa+M0jeERFVIhjex/0tJy8GO/ev+tAMw5zHH2Pyl9Xn6kUnXXXQe+HhyEO+/MzvPu2wc7dzr0XBoL5+Cr1K1ceWCsc3Mz3HtvYetR1lRVwaxZcNVVcNllcfDWjh1w0knx/pVXxscNI1Qs7H9SmbwESaXq5psPfF1dDR//+PE/Z29vvKXTDj2XxsI5+Cp1DQ1w/vlxqkZHx+gpHCoLzc1w1llwzjnQ1lb8nyanUqkj3lS+7H9SmQwkpFI0/KnssP7+7Hwqu2zZgU96Xf5NOjbn4KscNDTAhRfGjyclqUDsf1KZ7CEhlaKjfSq7bNnEnvPQkGPk0POJPqdU7pyDr3IzvHztk08WsgpJFcr+J0Witxc2bIhTuHMcVhtISKVo+fI4RWNkKHG8n8rmIuTIJd80q1isXAkPPwx79zoHX6XJ62hFWLKk0BVIYzPc/2TWrDhgd8eOA0t7zpxZ/FOOSt7w6mF798Yp3F1dOZ3OZ6YklaJly+IFYnguZTY+lV2+PPaiGKm62qHn0rE4B1+SpJxobo59T97yltLof1IWDrd6WA4ZSEilKl+d0ce7ILRUiZ55Js7Bv+mmQlciHZ/e3vjz7EpLklR5CrB6mIGEVKqy/ans7bfHJT9H2rfP1QIkqVIMD9Pt63OlJUmqRAVYPcxAQiplzzwTLxrZaDbjagHS8fGTZZW6PA/TlY4kk8kc8SYphwrw94CBhKRouC/FMFcLkMbOT5ZV6gowTFeSVGRy0afuGAwkJB2wciVMnhy/drUAaez8ZFml7kjDdD/ykQOrGkmSyl+++tQNMZCQdICrBUjj5yfLKgdHGqY7Y0Zh6pEkFUae/x4wkJB0sIaGuFpANvpSSJWgAA2gpKw70jDdM84obF2SpPzL498DBhKSSpMNBFUsbAircpHnYbqSJBlISDrYk0/GWzGzgaCKSQEaQEk54bQ9SVKeGUhIKj02EFSx8ZNllQun7UmS8shAQuXp0kvtCl6ubCCoYuQny5IkSeOW00AiSZKmJEnWJ0ly9tD965Ik2ZAkybokSf4qSZKaXJ5fUhmygaCKlZ8sqxwN9eupe+mlQlciSSpDOQskkiS5CPg5cO7Q/QT4AvCuEML5QC3w8VydXxXOhoflywaCKlal0H9FOpaRP8cj+vXMvOUW+/VIUqXI43uaXI6Q+D3gJmDb0P03Af87hPDK0P0fAO/P4flVqWx4WN5sIChJ+TGiX0/tf/6n/XokSVmXs0AihPDhEMI/jdj0HPC2JEmmJ0lSDVwFnJ6r86uC2fCw/NlAUJJy65B+PVV799qvR5KUdanMcGO4HEmS5GXg0hDCy0mSfBC4FdgNfAf4cAihbSzP09nZeTawKUdlqkyc8v3vM/1LX6J6z5792/bV1bH5ttt49Td+o4CVKdvOu/Za6jZt4l+//nX2nHNOocsBmNHe3v5yoYvINa/FUmV404IF1P7f/ztq+8DJJ/P8qlUFqGjMyv5a7HVYUpEb13U4b00lkySpA34RQpg7dP8DwLg7JLW1tTF58uRslwdAZ2cn7e3tOXnuQqm417RoUVx1YYTqPXs4++67Ofvzn89DdRNTcf+fsmHKFJgyhdarr87dOUY40uvZu3cv69evz0sNxSSX1+KJKsd/R9nk9+fo/P4c4stfho9//OBpj/X11K5YUZTfp0q8FvueeLRSrRtKt3brzr9irX2i1+F8LvvZAPwkSZLGJEkmATcD/5DH86sS2PBQx8PlYiUpOqRfT3ryZPv1SJKyLm+BRAjhVeDzwNPAeuCJEMLf5+v8qhA2PJQkKTtG9OsZeP3r7dcjScq6nE/ZCCGcPeLrewF/mwmApqYmdu7cOWp7Y2MjPT09E3/ilSvh4Ydh714bHkqSNFENDXD++bBhA7/8y7+k9dARiLkyPFLNZXQlqezlc8qGdJDDhRFH2z5mDQ1QWwtVVbFDeL7eQEmSVG4aGuDCC4ulebAkqcwYSKg8VVdDYyO0tha6EkmSNB69vfDMM9DVVehKJEk5ZiCh8rRvH+zc6ZsZSZJKSW8vrFsHfX2xJ9TIVT4kSWXHQELlp7cXBgYgnfbNTDl78knnF0tSuVm2DPr749fd3XDDDYWtR5KUUwYSKj++mZEkqfSsXBl7P2Uy8f6ePfDII3G7JKksGUioYBobG8e1fUx8MyNJUmm6/fbRoxr7+uJ2SVJZMpBQwfT09JDJZEbdjmvJT9/M6HjYSE2SDpbP6XHLl49eGau+Hu64Iz/nlyTlnYGE8qKpqYlUKjXq1tTUlN0T+WZGE2UjNUkqrGXL4vU3lYr36+pgyRK4/vrC1iVJyhkDCeXFzp07x7V9wnwzo4my94gkFd7KlTBpUvy6uRnuvbew9UiScspAQuXHNzMaL3uPSFJxaGiA88+Poxs7OkaPepSkIvHII/Gm42MgofLjmxmNl71HJKl4NDTAhRdCa2uhK5Ek5ZiBhMqTb2Y0HvYekSRJkvLOQEKS7D0iScUjnyt7SJIKykBCedHY2Diu7cfNNzMaL3uPSJIkSXllIKG86OnpIZPJjLr19PQUujQpsveIJEmSlFcGEpI0zN4jkiRJUt4YSEiSJEmSpLwzkJAkSZIkSXlnICFJkiRJkvLOQEKSJEmSJOWdgYQkSZIkSco7AwlJkiRJkpR3BhKSJEmSJCnvagpdgCQVjSefLHQFkiRJUsVwhIQkSZIkSRWuqamJVCo16tbU1JSzcxpISJIkSZJU4Xbu3Dmu7dlgICFJyolCpOySJKk8+D6iMthDQpKUE4VI2SVJUnko5vcR3d3w8suwbx+sWwdTpkBzc6GrKk0GEpIkSZIkHUXc8TVnAAAgAElEQVQ6DSFAV1cMIhobobYWdu6EDRuguhpaWyFJoMp5CGNmICFJkiRJ0hEMDMCqVTGUuPhimDoVnn8+PjZnDmQysHUrrF0LmzbBggUxrNCxmd1IkiRJknQY6XQMI+rrYdEiaGmBVOrgfVIpmDYNFi+O+w2HF6WmsbFxXNuzwUBCkiRJkqTDCCGGC/PmHXsqRlVV3G/fvnhcqenp6SGTyYy69fT05OycBhKSpJwoRMouSZLKQ7G8j+jqgrlzx94Xoqoq7t/Vldu6yoU9JCRJOZHLNF2SJJW3Yngf0d0dRztMnTq+41paYHAwHu/qG0fnCAlJkiRJkg6xfTtMnz66Z8SxpFLxuB07clNXOTGQkCRJkiTpEAMDE18tY9Ik6O/Pbj3lyEBCkiRJkqRD1NbGUGIi+vtjKKGjM5CQJEmSJOkQU6bA5s2QyYzvuEwmHnfaabmpq5wYSEiSJEmSdIjmZqiuhm3bxnfc1q1xdIUNLY/NQEKSJEmSpMNobYW1ayGdHtv+6XTcf/bs3NZVLgwkJEmSJEk6jCSBqipYvfrYoUQ6HferqYnH6dhqxrJTkiSfPWRTBugD1ocQHst6VZIkSZIkFVhVFSxYAKtWQUcHzJ0LLS0H75PJxGkaa9fGMGL+/Hicjm1MgQRwPvB24EFgH7AUeBm4OkmSt4YQ/iQ35UmSJEmSVDi1tbBwIYQATz8Ng4Owe3dcRWPXrtjAsrY2TtMYHlGhsRlrINEMtIcQ/h0gSZIvAg8A7wA6AQMJSZIkSVJZqqqCWbPirbsbHngAfv1rmDMHZs60geVEjTW7OWU4jAAIIbw6tK0fmODKrJIkSZIklZbmZjjrLDjnHGhrK58woqmpiVQqNerW1NSUs3OOdYTEr5IkWQ7cA6SAG4CXkiS5iDiF47CSJPkccPXQ3Y4Qwm1JkswHvgKcAPxDCOGPJly9JEmSJEk6bjt37hzX9mwY6wiJ64GzgbXAL4AW4MPAW4BbD3fAUPBwOTAXeDPQniTJbwErgfcBs4ALkyR5z3HUL0mSJEmSStCYRkiEEP4D+K3DPPQ3RznsFeD3h6Z1kCTJRuBc4MUQwqahbd8EPgD8aDxFS5IkSZKk0pbKZDLH3ClJkouBzwAnEqdsVAMzQghnjuUkSZK8EXgKuAtIQggfGto+H7gthHD5sZ6js7PzbGDTWM4nSQUwo729/eVCF5FrXoslFbmyvxZ7HZaKw+rVJwEwb96vC1xJ9lxwwQVHfOzZZ58d69OM6zo81h4S/x/wDeAq4G7g/cB3x3JgkiStQAfwB8AgcZTEsBSQHmuxAG1tbUyePHk8h4xZZ2cn7e3tOXnuQvE1lQZfU/E70uvZu3cv69evL0BFhZXLa/FEldvPXLb5/Tk6vz9HV+zfn0q8FvueeLRSrRtKt/ZKrnvbtvjffL/8Qn3Pj3XOiV6Hx9pDIhNC+DPgSeBfiY0qjzmqIUmSS4CfAJ8JIXwd2AKcMWKX04Ft4ylYkiRJkiRlV2Nj47i2Z8NYA4nhtpovAW0hhN0cZXUNgCRJpgP/CPx2COHbQ5vXxIeSmUmSVAO/jf0jJEmSJEkqqJ6eHjKZzKhbT09Pzs451ikbv0iS5B+A/wF0JElyLnH6xdHcCtQBX0mSZHjb3cB1xOkedcAPgQfHWbMkSZIkSSpxYw0kbgHeNrT/XcB7gQ8d7YAQwieATxzh4TljLVCSpFLR1NR02LW6Gxsbc/rpgiRJUikaayDxZeBjwMh3U08CU7JdkCRJpepwYcTRtkuSJFWysQYSVwJTQwiv5rIYSVL5eeSR+N8lSwpbhyRJkorLWJtavgC8lstCJEmSJElS5RjrCIk7gZ8lSfJTYGB4Ywjhj3NSlSRJkiRJKmtjDSQ+Q+wf8boc1iJJkiRJkirEWAOJhhDCf8lpJZIklbjGxsYjrrIhSZKkg421h0RIkuRNOa1EkqQS19PTQyaTGXVzyU9JKn1NTU2kUqlRt6ampkKXpgJYssSG3dkw1hESZwLPJkmyCdg7vDGEYEghSZIkqey5tLOUfWMNJG7PaRWSJEmSJKmijCmQCCH8LNeFSJIkSZKkyjHWHhKSJEmSJElZYyAhSZIkSZLyzkBCkiRJko7hSEs4u7SzNHFjbWopSZIkSRXLJZyl7HOEhCRJkiRJyjsDCUmSJEmSlHcGEpIkSZIkKe8MJCRJkiRJUt4ZSEiSJEmSpLwzkJAkSZIkSXlnICFJkiRJkvLOQEKSJEmSJOWdgYQkSZIkSco7AwlJkiRJkpR3BhKSJEmSJCnvDCQkSZIkSVLeGUhIkiRJkqS8M5CQJEmSJEl5ZyAhSZIkSZLyzkBCkiRJkiTlnYGEJEmSJEnKOwMJSZIkSZKUdwYSkiRJkiQp7wwkJEmSJElS3hlISJIkSZKkvDOQkCRJkiRJeWcgIUmSJEmS8s5AQpIkSZIk5Z2BhFSCmpqaSKVSo25NTU2FLk2SJEmSxsRAQipBO3fuHNd2SZIkSSo2BhKSJEmSJCnvDCQkSZIkSVLeGUhIkiRJkqS8M5CQJEmSJEl5V5PLJ0+S5HPA1UN3O0IItw1trwUeBf4khPBkLmuQylFjY+NhG1g2NjYWoBpJkiRJGr+cBRJJkswHLgfmAhng0SRJlgIbgJXAW3J1bqnc9fT0FLoESZIkSTouuZyy8Qrw+yGE/hDCALAROBO4AfhzYE0Ozy1JkiRJkopYzkZIhBC6hr9OkuSNxKkbl4QQXhzadkuuzi1JkiRJkopbKpPJ5PQESZK0Ah3A50IIXx+x/Ung/xlrD4nOzs6zgU05KFGSsmFGe3v7y4UuItcmci1evfokAObN+3UOKpKkg5T9tdj3xJKK3Liuw7luankJ8F3glhDCt7PxnG1tbUyePDkbTzVKZ2cn7e3tOXnuQvE1lQZfU/E70uvZu3cv69evL0BFhTWea/G2bfG/uf5xKLefuWzz+3N0fn+Orti/P5V4LfY98WilWjeUbu3WnX/FWvtEr8O5bGo5HfhH4JoQwhO5Oo8kSZIkSSo9uRwhcStQB3wlSZLhbXeHEO7O4TklSZIkSVIJyGVTy08AnzjK45fm6tySJEmSJKm45XLZT0mS9mtqaiKVSo26NTU1Fbo0SZIkFYCBhCQpL3bu3Dmu7ZIkSSpvBhKSJEmSJCnvDCQkSZIkSVLeGUhIkiRJkqS8M5CQJEmSJEl5V7aBhN3cJam4NDY2jmu7JEmSyltNoQvIFbu5S1Jx6enpKXQJkiRJKiJlO0JCkiRJkiQVLwMJSZIkSZKUdwYSkiRJkiQp7wwkJEmSJElS3pVtIGE3d0kqvO5uePlleOklWLcu3pckSZKgjFfZsJu7JBVGOg0hQFcX7NsHjY1QWws7d8KGDVBdDa2tkCRQVbaxuCRJko6lbAMJSVL+DQzAqlUxlLj4Ypg6FZ5/Pj42Zw5kMrB1K6xdC5s2wYIFMayQJElS5fGzKUlSVqTTMYyor4dFi6ClBVKpg/dJpWDaNFi8OO43HF5IkiSp8hhISJKyIoQYLsybd/SpGM89F/tJzJsXp3SEkL8aJUmSVDwMJCRJWdHVBXPnjr0vRFVV3L+rK7d1SZIkqTgZSEiSjlt3dxztMHXq+I5raYHBQVffkCRJqkQGEpKk47Z9O0yfPrpnxLGkUvG4HTtyU5ckSZKKl4GEJOm4DQxMfLWMSZOgvz+79UiSJKn4GUhIko5bbW0MJSaivz+GEpIkSaosBhKSpOM2ZQps3gyZzPiOy2Ticaedlpu6JEmSVLwMJCRJx625GaqrYdu28R23dWscXdHcnJu6JEmSVLwMJCRJWdHaCmvXQjo9tv3T6bj/7Nm5rUuSJEnFyUBCkpQVSQJVVbB69bFDiXQ67ldTE4+TJElS5TGQkCRlRVUVLFgAfX3Q0QFbtozuKZHJQHc3/PSnsHs3zJ8fj5MkSVLlqSl0AZKk8lFbCwsXQgjw9NMwOBiDh0mTYNeu2MBywwY4/XR497sNIyRJkiqZgYQkKauqqmDWrHjr7oYHHoBf/xrmzIGZM6Gh4cB+kiRJqlwGEpKknGluhrPOil+3tRW2FkmSJBUXP5+SJEmSJEl5ZyAhSZIkSZLyzkBCkiRJkiTlnYGEJEmSJEnKOwMJSZIkSZKUdwYSkiRJkiQp7wwkJEmSJElS3hlISJIkSZKkvDOQkCRJkiRJeVdT6AIkSZWjuxtefhn27YN162DKFGhuLnRVkiRJKgQDCUlSTqXTsHUrPPhgDCIaG6G2FnbuhA0boLoaWlshSaDKcXuSJEkVw0BCkpQzAwPw3HMxhLj4Ypg6FZ5/Pj42Zw5kMjGsWLsWNm2CBQtiWCFJkqTy52dRkqScSKdh1So44wy49FJoaYFU6uB9UimYNg0WL4b6+rh/Ol2QciVJkpRnBhKSpJwIIYYLF1xw7KkYVVUwb16c0hFCfuqTJElSYRlISJJyoqsL5s4de1+Iqqq4f1dXbuuSJElScchpD4kkST4HXD10tyOEcFuSJDcCHwcywLPAfw0h9OeyDklSfnV3x9EOU6fCf/zH2I9raYHBwXi8q29IkiSVt5yNkEiSZD5wOTAXeDPQniTJp4E/AN4OvGno/DflqgZJUmFs3w7Tp4/uGXEsqVQ8bseO3NQlSZKk4pHLERKvAL8/PPohSZKNQB3w0RBCz9C2dcCZOaxBklQAAwMTXy1j0iTod9ycJElS2ctZIBFC2D8LOEmSNxKnblwSQnhxaNtpwMeA63JVgySpMGprYc+eiR3b3w8nnZTdeiRJklR8UplMJqcnSJKkFegAPhdC+PrQthbgR8ADIYQ/GcvzdHZ2ng1sylWdknScZrS3t79c6CJybazX4ldfrWHDhtN4//sH+OUvTwDgjW/cDcCLLx58f6RMBr73vUm0tm7nlFMGs1e4pEpR9tdi3xNLKnLjug7nuqnlJcB3gVtCCN8e2nYe8BhwZwhhxXifs62tjcmTJ2e30CGdnZ20t7fn5LkLxddUGnxNxe9Ir2fv3r2sX7++ABUV1liuxT09caTEaafF+2edFf/72msH3x9pyxY480y4/PIzxl1Tuf3MZZvfn6Pz+3N0xf79qcRrse+JRyvVuqF0a7fu/CvW2id6Hc5ZIJEkyXTgH4FrQghPDG1rBH4M/GEI4b5cnVuSVHitrbB2bVw5YyxLf6bTcf/Zs3NfmySVo4GBAbZs2cKeic6ZG1JTU8PGjRuzVFX+HG/ddXV1TJs2jdqJNkGSNG65HCFxK7GJ5VeSJBne9g9AM/D7SZL8/tC274cQPpvDOiRJBZAksGkTPPssXHDB0fdNp2H1aqipicdJksZvy5YtNDY2cvbZZ5Ma7zJHI/T29tLQ0JDFyvLjeOrOZDK8+uqrbNmyhRkzZmS5MklHksumlp8APnGYh5bn6pySpOJRVQULFsDy5fDTn8apGy0tB++TycDWrXFkRE0NzJ8/ttEUkqTR9uzZc9xhRKVKpVKccsop7HDdaSmvctpDQpJU2WprYc6cGDo8/TQMDsLu3XFpz127YPPmuM/s2XFkhGGEJB0fw4iJ83sn5Z9v/SRJOVVVBdOnw1VXwWWXQW8v7NgRl/a87DK48kqYNcswQpI0Ptdeey1r1qxh3bp1/OEf/mGhy5E0AY6QkCTlTXPzgdU12toKW4skqTycf/75nH/++fT29ha6FEnjZCAhSZIkKavWrFnD3XffTW1tLVu2bOGd73wn9fX1PP744wDcc889nHrqqXzzm9/k4YcfZvfu3dTW1rJixQre8IY3sGbNGr7whS9QXV3Nm9/8Zl566SXuu+/wi/StWbOGr371q9x9991ce+21zJ49m87OTvbu3cutt97KN77xDV566SWuu+46rrvuujx+FyQdiwNkJUmSpArU1NREKpUadTv99NOz8vzPPfccn//85/nud7/L/fffz+tf/3oeeughkiSho6ODXbt28fjjj3Pffffxgx/8gEsvvZT777+fgYEBbrvtNv78z/+cf/zHf6SmZnyfoWYyGR588EHe/e5384UvfIGvfvWr3H///fz1X/91Vl6XpOwxkJAkSZIq0M6dOw+7fdeuXVl5/nPPPZczzjiDE044gZNPPpmLL74YgKlTp9LT08OJJ57IihUr6OjoYMWKFfz0pz+lr6+PF154gVNOOYXzzjsPgKuuumpc5503b97+88yZM4cTTjiBlpYWenp6svK6JGWPgYQkSZKkrKutrT3ofnV19UH3X3nlFa655hp27tzJvHnzWLp0KZlMhurqatLpdFbOO97RFZLyy0BCkiRJUt6tW7eOs846i+uuu47zzz+fxx9/nH379vGGN7yBnp4eQggAPPLIIwWuVFKuGBlKkiRJyrtLLrmEb33rWyxatIhMJsOFF17Iiy++yKRJk/jSl77Epz/9aaqqqpgxYwZ1dXWFLldSDhhISJIkScqqiy66iIsuumj//SeeeGL/1zfffPP+r7/2ta+NOjadTvPEE0/w93//99TX1/O1r32N7u7uUfuNXHXjoosuore396BtV1xxBVdcccX++8MjLiQVDwMJSZIkqQI1NjYetrHliSeeWIBqDqiqquJ1r3sdV111FbW1tbS0tPDFL36xoDVJyg0DCUmSJKkCHWnVid7e3jxXMtqNN97IjTfeWOgyJOWYTS0lSZIkSVLeOUJCkiRJqlDd3bB9OwwMQG0tTJkCBZ6xIamCGEhIkiRJFSSdhhCgqwv27YPp02MYsWcPbNgAg4M1vOUtkCRQ5XhqSTlkICFJkiRViIEBWLUqhhIXXwxTp0IqdeDxTAZ+9asYTGzaBAsWxLBCknLBzFOSJEmqAOl0DCPq62HRImhpiWHEc8/FG8T706bB4sVxv+HwQpJywRESkiRJUgUIIYYL8+YdeypGVVXcr6MjHjdr1vjPt2XLFhYuXMg555xz0PbW1lZ+9rOfceqppwKwZ88eFi5cyCc/+UkAHn30Ue655x4GBwfJZDK8733v48Mf/vD4C5BU9AwkJEmSpArQ1RWnaYy1L0RVFcydC08/PbFAAmDKlCk8/PDDB2276667+M3f/E1uvvlmAPr6+li0aBEXXHAB5557Ln/2Z3/GQw89xMknn0xvby/XXnstM2bM4F3vetfEipBUtAwkJEmSpDLX3R0bWE6dOr7jWlpgcDAe39ycm9rq6+t505vexIsvvshpp53GwMAAe/bsAaChoYE77riDyZMn5+bkkgrKQEKSJEkqc9u3x9U0RjawHItUKh63Y8fEAont27fzvve9b//9JUuWjNpn69at/PM//zO/+7u/y3nnnce73vUu5s+fz6xZs7joootYsmQJZ5111vhPLqnoGUhIkiRJZW5gYOKrZUyaBP39Ezv2SFM2vv3tb/P444+TTqeprq7mIx/5CO3t7QB8/vOf56Mf/Sg///nP+fnPf87VV1/Nl7/8ZS6//PKJFSGpaBlISJIkSWWuthaGZkGMW38/nHRSdusZ2UNipCeffHJ/T4krr7ySK6+8ku985zs8+OCDBhJSGXLZT0lSXi1ZEm+SpPyZMgU2b4ZMZnzHZTLxuNNOy01dh6qrq2PFihVs2bJl6PwZNm7cyKyJdtWUVNQcISFJkiSVueZmqK6Gbdtio8qx2ro1jq7IVUPLQ73tbW/jYx/7GB/5yEcYGBgA4B3veAc33XRTfgqQlFcGEpIkSVIFaG2FtWvhjDPGtvRnOh33nz17YuebNm0aTzzxxKjth5uqMdLSpUtZunTpxE4qqaQ4ZUOSJEmqAEkSg4jVq2PYcDTpdNyvpiYeJ0m5YCAhSZIkVYCqKliwAPr6oKMDtmyJPSLmzIk3iPe3bImP794N8+ePbTSFJE2EUzYkSZKkClFbCwsXQgjw9NMwOAjTpx9Y2nPz5jg6Yu7cAyMqJClXDCQkSZKkClJVBbNmxVt3N+zYcWBpz5kz4cQTB2lomFzoMiVVAAMJSZIkqUI1N49eQaO3tzC1SKo8DsKSJEmSKtmll8abJOWZgYQkSZKkrNuyZQtJkvDZz372oO0bN24kSRIeeuih43r+a6+9ljVr1ox5/7vuuou77roLgMSlQ6SiYCAhSZIkVbLeXnjmGejqyvpTv+51r+Of/umf2Ldv3/5tP/zhD3n961+f9XNJKj0GEpIkSVKl6u2FdeviWqCLF2e9gURDQwOzZs3imWee2b/tqaee4u1vfzsA3/zmN/nABz7Ae9/7XpYuXcqvfvUrAN75zndyyy238O53v5tt27Zx4403csUVV3DFFVfwk5/8ZP9zPfjggyxdupR3vetd/OxnPwPghRde4Nprr+XKK6/ksssu41vf+lZWX5Ok7DGQkCRJkirVsmVxiQ2IS27ccEPWT/Ge97yHxx57DIDnn3+eJEmora1l165dPP7449x333384Ac/4NJLL+X+++/ff9y8efN47LHHWLNmDS0tLTz00EN88Ytf5Nlnn92/T2NjI9/73vf4oz/6I/7u7/4OgAceeICPfvSjfPe73+Ub3/gGX/rSl7L+miRlh4GEJEmSVIlWroSODshk4v09e+CRR6i5776snuad73wnq1evJp1O86Mf/Yj3vOc9AJx44omsWLGCjo4OVqxYwU9/+lP6+vr2HzdnzhwA5s6dy+OPP85HP/pR1q1bx0033bR/n/nz5wMwc+ZMXnvtNQA+85nPsHfvXv72b/+Wv/zLvzzoOSUVFwMJSVJOLVkSb5KkInP77aOnaPT1UXtIE8rj1dDQwHnnnUdnZydPP/30/ukar7zyCtdccw07d+5k3rx5LF26lMxwOAJMnjwZgLPPPpsf/ehHLFmyhGeffZarrrqKdDoNQHV1NQCpVGr/sbfccgurVq3inHPO4ZZbbsnqa5GUXTWFLkCSJElSASxfDh//+MGhRH09A3/8x0zO8qne8573sGLFCtra2qipqRk6VT1nnXUW1113HXv27OHOO+/k9NNPH3XsN7/5TTZv3sztt9/OvHnzuOyyy9i1a9cRz/XUU0/xox/9iObm5v1TQEY21ZRUPBwhIUmSJFWiZctiI8tUKt6vq4MlSxi89tqsn+qyyy5j48aNLFq0aP+22tpa0uk0ixYtYunSpcyYMYMtW7aMOvb9738/mzZtYsmSJXzwgx/kD/7gD2hqajriuW6++WZ++7d/m0WLFtHZ2UlLS8thn1dS4TlCQpIkSapUK1fCww/D3r3Q3Az33pu1p542bRpPPPEEEKdtPPfcc/sfu+OOOwD40Ic+dNhjh4+D2GvinnvuGbXPfSN6XUybNo2Ojg4Arr/+eq6//vpR+9988837vw4hjOelSMoRR0hIkiRJlaqhAc4/H+rrY4PLhoZCVySpgjhCQpIkSapkzzxT6AokVShHSEiSJEmSpLwzkJAkSZIkSXmX0ykbSZJ8Drh66G5HCOG2JEn+G/AxIAV0ALeFEDJHeg5JkiRJY5PJZEgNr5qhcclk/JNEyrecjZBIkmQ+cDkwF3gz0J4kySeBTwFvBc4H3g4syFUNkiRJUqWoq6vj1Vdf9Q/rCchkMrz66qvU1dUVuhSpouRyhMQrwO+HEPoBkiTZCKSB2SGEgSRJTgFOAl7LYQ2SJElSRZg2bRpbtmxhx44dx/U8/f39TJo0KUtV5c/x1l1XV8e0adOyWJGkY8lZIBFC6Br+OkmSNxKnblwyFEb8HvBl4BfAv+SqBkmSJKlS1NbWMmPGjON+ns7OTubMmZOFivKrVOuWKlkq10O6kiRpJfaK+FwI4esjttcAXwM2hxD++7Gep7Oz82xgU67qlKTjNKO9vf3lQheRa16LJRW5sr8Wex2WVOTGdR3OdVPLS4DvAreEEL6dJMl04MwQwlMhhMEkSb4N/LfxPGdbWxuTJ0/ORbl0dnbS3t6ek+cuFF9TafA1Fb8jvZ69e/eyfv36AlRUWLm8Fk9Uuf3MZZvfn6Pz+/P/t3fvQXdV5R3Hv2/CJRRUoIVyabBjIb8iI0RD0JGLVCkjFqRpLbG2hhhuGYHGVqJOBWOqAtYaYrEDNCQpM8EoF8lQsFQlpEVAlEASJfgIlKSGy5ShViFASODtH2u95ORknze3N/vy7t9nJpNz9llnv8/aZ6/nrLP22nsPru7bp4252H3izTU1bmhu7I67fHWNfXvz8E4bkMiDD4uAiRGxOC9+E3C9pLHAr4APAT/YylWOhHRu2M60bt26nbr+KrhOzeA61V9RfTpy0shSg6lOKbl4ew23fW6oefsMzttncHXePi3Lxe4TD6KpcUNzY3fc5atj7Nubh3faKRuSvgZMAR7vWHw10A9MAzYAd5NmT6zf0vqWLl16XC5vZlZHx48bN25rB1gby7nYzGpu2Odi52Ezq7ltysM786KW00gDD0Wu2Y5V/hg4nnT3jle3Ny4zsyE2EjiQlKPawLnYzOqoTbnYedjM6mi78vBOv6ilmZmZmZmZmVm3EVUHYGZmZmZmZmbt4wEJMzMzMzMzMyudByTMzMzMzMzMrHQekDAzMzMzMzOz0nlAwszMzMzMzMxK5wEJMzMzMzMzMyudByTMzMzMzMzMrHS7VB1AFST9A/BbETFZ0inAl/NLPwHOi4gXJO0NXA+8BXgWOCMinqkm4t4k3QXsD6zPi84Dfg+4GNgVmB0R/5TLngTMAvYAvhURF5cf8ZYV1Ski7pf0RuBe4NSIWJXLNrZOwPuBM/Lz2yPiU7lsk+v0R8CHgH5gbkTMymUbW6eIuD+/9nreyM/HAtcCbwT+E5gaERtKD9oAkPQRCvJeW0maQVd+aUo7LFNXf8BtOpN0GjAD2BP4bkRM8/4z/DW1PXTFPQOYAvwyvzynjt8H29J/r5MecV/evWyg71QXTc5pPWKfDxwHrM3FZkbELVXFuKNaN0NC0vuAM/PjvYHrgA9HxJHAcuDSXPSLwN0RcTgwB/haBeEOSlIfMAY4KiLGRsRYYA3wJdJOOhY4V9JbJe0BzANOBw4HxufBmFopqlMejHgn8IP82kDZxtYJeMS3RgoAAAqSSURBVANwMvB20uc0TtKEhtdpFPBe4EjgaOBCJY2tU8dgxOt5o8MC4IKIGAP0AeeUGrC9TtLBFOS9aqOqTu5kdeeXP6cB7bBMBe3abRqQ9BbgauCPSfn8HXlf8f4zjDW1PRTEfTSpXz/wPV7HH/Vb3X+vMMzN9Ij7R93LajgY0dicNkjsRwMndGzzxg5GQMsGJCTtS2rsA4MOhwGrI2Jlfn4b6QOHdKT3+vx4IXCKpF3LinUrKf//XUnLJV0AnAQsjoj/jYi1wE2kI9bHAI9GxBN5hHsB8GeVRD24ojpB+iI8H3iqo2yT6/Q08MmIeCUi1gOPAIfQ4DpFxH8Af5Dj3p80A2stDa4TFOYNJL0Z2CMifpgX/Qv1rFNb9Mp7bVWUX8bQjHZYiu527Ta9iQmko4Vr8v4zEXgR7z/DVlPbQ9H3M+mH2t9KWiHp65JGVRPdoLal/14nRXH36rfXSZNzWlHsy0m/Gebl/XympEb/pm/bKRvXAJ8FRufnjwKjJR0VEctJ01sPyK8dROrUEREbJP0a2I9NfxBXbR/gTuBC0vSuJcC3yHFnT5N+EB5UsPx3Soly22xWJ0kREWcDSOos29g6ARER3wOQdBhp3zuW9EXayDrlz+l7kmYCFwE3Ak8Cx9PgOgHnsmnegObse21R9HkcU1EslYuIhwced+SXK/E+26m7P+A2vdGhwCuSbiV1em8DHsbbZzhranvYJG5JewEPAdOBx0gDKZfkMnWyLf33OimKu2d/sKogCzQ5pxXFPh9YDHwc+FVedhZpRn8jtWZAQtLZwC8i4k5JkwEi4v8kTQL+OY8szQFeyW/p61pFH/BaWfFujYi4D7hv4LmkuaRzob7YUWwg7hGk8/q7l9dKjzp9AChKbI2vk6QjgNuB6RHxqKRjaHidImKGpC8D/0qa2bKW5tbpBmB+Z97IGrHvtYg/jwKd+QXYQMcpb7R4GxX1B/A+1GkX4ATgROAF4FbgJbx9hqWmtoce/foXSH2RgTJfJU3Lr9WAxDb232ujR9yHRMSkrmW9+u1VaXJOK4r9sYiYMFBA0pXAJDwg0QgTgQMlLQP2BfaSdAXpYjfvBJA0Hng8l3+SNFtijaRdSOf8P1d+2L1JOg7YPSLuzIv6gFXAgR3FDiDN6ljTY3mt9KjT+h7FG10nSccCNwOfiIhv5teaXKeDJY2NiGUR8aKkb5POd7uR5tbpAeDkgrwxmwbUqUXWkGbiDGj959GdXyS9B++zAzbrD5A6pt4+yTPA9yPiWQBJt5CmMr/aUabN22e4aWp7KOrXzydd/21eLjNYH7Iy29h/r40ecR8l6X1b2W+vSpNzWlHsZ0p6PiJuzmXquM23SWsGJCLiDwce55HUE4FPAv+tdMHEp4C/IU2ZAvgOabTpUlLSuzufu1MnewN/J+ndpGlSZwJ/CSyQtB/pyPSfkqacrwAk6VDgCeAjpFHjuimq09QeZe+nuXW6AFgETIyIxR1lm1ynG4A5+Qurn3ShoHk0u05T8xGB1/NGRPx1fv6ypGMj4h7go8C/VRO2Ad8HPl+Q91pJ0mg2zy9NaYc7XVF/ICI+JumnbtNAmv57ndKFv58HTiGdz/4Z7z/DT1PbQ49+/aeAR5TuBLGKdO2xOl7sb1v673VSFPdNwFe2st9elSbntKLYFwGzJS0mzZo4l3SThsZq9AUwdlREvEa6Xc0dQJBuEfSV/PIlwLskPUw6R+f8SoIcRETcRpqO+xCwFJiXvzg+C9wFLAO+ERE/ioiXgcmkI2YrgZ+RGmOt9KjTfT3KNrZOpAsVjQJmSVqW/01tcp0i4ktdy+6NiG82vE6F+172F8AVkn5GOqL0jzs/SisSEU9SkPeqjapSF9GVX0htcDI1b4cVc5sG8hXy/550Z6uVwGrgKrz/tE3j2kM+inwe6ZTRIB05/mqlQRXYlv57dVFurkfcXyhYNljfqXRNzmk9Yr8SuAy4Jy9bFhELKwtyCPT19/dvuZSZmZmZmZmZ2RBq9QwJMzMzMzMzM6uGByTMzMzMzMzMrHQekDAzMzMzMzOz0nlAwszMzMzMzMxK15rbfppZs0k6FrgC2A14DpgSEat7lH0D6SrVZ0XEkrzsv4BfdxQ7LSJ+IWk6MIV0q9JrI2LWFuIQcA2wD+n+0B+OiF/uSN3MzOqk7HwraT5wHOl2hwAzI6Ln7Rqdh82sDdqSiz0gYa0iqR8YD9wP/CQvHkG6j+/siLihqthsi64HPhgRKyRNId2C7PQeZb9OSo4ASPpN4JWIGNtZKN9/+uPAW0n7wUpJt0bEY0UrldQH3ApMi4g7JF0OfAb49I5Vzaw9nIcboex8ezRwQkQ8vaXAnIfNhoZzcSO0Ihd7QMLa6qXOBirpzcCdkl6NiJsrjMsySScCkyNisqTdgYsjYkV+eQVwYY/3TQSez2UGjAf6JN0D7AFcFhE3khLxbsAo0v3K+4D1eT2TgE/kMkuB84EjgLURcUde76XA3kNSYbP2cR6uiSrzraTfAA4B5kk6GLiFdFTuNedhs1I4F9dEW3OxryFhBuTpT58Dplcdi20uItZFxAIASSOAzwOLustJOoSUMLs/x1HAvwMnAX8CzJJ0eET8HFgIrAJWAzdHxGpJRwDnAO/OX9L/A1wEHAo8I2mupAeBq0hHEsxsBzkP10PZ+RY4AFhMmj78LuB44CznYbNqOBfXQ5tysWdImG20HHhb1UG0naQJwAxgL2BfScuAhyLiY5J2A64j5a5Lu943ApgLXBARL6XT2pKIWMTGJL5K0reBk/NRgHHAwaQB2u9IOgPYHzgM+GFez27Ag8ATwImk6WwPSPoCMAuYPNTbwaylnIdLVId8m6eFT+hY95XAJGB3nIfNquJcXKK252IPSJht1A+8WHUQbZcvnnNL57Q1AEl7kc5Vew44PSLWd7319/O/uTlhHgpcK+kcYE/gmYh4IJcdODXjNNLI8Av5bywE3gP8HLghIv6q42/vQkrgj3asZyFw05BuALN2cx4uUR3yraRHgDEdU8MHyo/EedisKs7FJWp7LvYpGzasSXq7pLflx315cXdjHjCejRf1sfpZADwGTIyIdd0vRsTKiBgdEWPzlLIHgLMj4i7gd4EZkkZI+m3gg8DtpCMAp0oaKWlX4BTgx8ASYIKk/fN+cxVpOty9wH6Sjsp/9jTSuXRm1oPzcCOVmW/7gNmS9snLzyWdu7wE52GzIeNc3EityMWeIWHD3anAGOCjwDuAZ9n09jcASBoDXAJMKzU66ynfsmgJpC9R0lWFVwIP5lHgpyLiA5KmAgdFxOcGWd3VwJHAT0kDsZ/O14q4ljSy/DCwgZSor4uIfkkzSefSjSDdRunyiHg5T6ubI2lPYA1p3zKz3pyHa64G+fYy4B5gV9KRu4U5Fudhs6HjXFxzbc3Fff39/dvzPrNGkLQP8A1gNKlxTSddgfZxNo78vga8DFyRrz5rZmZDxHnYzKx6zsVWVx6QMDMzMzMzM7PS+RoSZmZmZmZmZlY6D0iYmZmZmZmZWek8IGFmZmZmZmZmpfOAhJmZmZmZmZmVzgMSZmZmZmZmZlY6D0iYmZmZmZmZWek8IGFmZmZmZmZmpfOAhJmZmZmZmZmV7v8B72cg51vSsV0AAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python",
   "version": "3.7.1",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}