KPF-CC-Scheduler_BJ.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0a1da87e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'1.3.3'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "from datetime import datetime, timedelta\n",
    "import astropy\n",
    "import calendar\n",
    "from astropy.time import Time,TimezoneInfo,TimeDelta\n",
    "from astropy import units as u\n",
    "from astropy import coordinates as co\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pylab as pl\n",
    "from astroplan import Observer, FixedTarget, is_observable, is_always_observable, months_observable, observability_table\n",
    "from astroplan import (AltitudeConstraint, AirmassConstraint,\n",
    "                       AtNightConstraint)\n",
    "\n",
    "\n",
    "np.random.seed(2)\n",
    "pd.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3ec431b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "nnights = 365\n",
    "max_obs = 60\n",
    "target = FixedTarget.from_name('HD155060')\n",
    "t = Time('2021-02-2') + arange(nnights) * u.day # UT date\n",
    "df = pd.DataFrame(dict(utdate=t.to_value('iso','date')))\n",
    "keck = Observer.at_site(\"Keck\", timezone=\"US/Hawaii\")\n",
    "sunset = keck.sun_set_time(t,horizon=-12*u.deg,which='next')\n",
    "sunrise =  keck.sun_rise_time(t,horizon=-12*u.deg, which='next')\n",
    "moonrise = keck.moon_rise_time(t,horizon=0*u.deg, which='next')\n",
    "nightlength = (sunrise - sunset)\n",
    "constraints = [AltitudeConstraint(20*u.deg, 85*u.deg)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "33b8f3ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "savecols = 'utdate quarter start stop start_jd stop_jd kpfcc dark'.split()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "61192c88",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1012/1460 quarters remain after dropping dark all, early, late\n",
      "362/1460 quarters are KPF-CC\n",
      "272/365 nights KPF-CC onsky\n"
     ]
    }
   ],
   "source": [
    "df2 = []\n",
    "for i in range(1,5):\n",
    "    df = pd.DataFrame(dict(utdate=t.to_value('iso','date')))\n",
    "    df = pd.DataFrame(dict(utdate=Time(t.to_datetime(timezone=keck.timezone)).to_value('iso','date')))\n",
    "    df['quarter'] = i\n",
    "    start = sunset + (i-1) / 4.0 * nightlength\n",
    "    stop = sunset +  i / 4.0 * nightlength\n",
    "    df['nightlength'] = nightlength.value\n",
    "    df['start_jd'] =  start.jd\n",
    "    df['stop_jd'] = stop.jd\n",
    "    df['start'] = pd.Series(start.to_value('iso','date_hm')).str.slice(start=11,stop=16)\n",
    "    df['stop'] = pd.Series(stop.to_value('iso','date_hm')).str.slice(start=11,stop=16)\n",
    "    df['month'] = start.ymdhms['month']\n",
    "    df['moon_alt_start'] = keck.moon_altaz(start).alt.deg\n",
    "    df['moon_alt_stop'] = keck.moon_altaz(stop).alt.deg\n",
    "    df2.append(df)\n",
    "\n",
    "df = pd.concat(df2)\n",
    "df = df.sort_values(by=['start_jd'])\n",
    "df = df.reset_index(drop=True)\n",
    "df['kpfcc'] = False\n",
    "df['dark'] = (df.moon_alt_start < 0) & (df.moon_alt_stop < 0)\n",
    "\n",
    "dark_early = df[df.quarter.isin([1,2])].groupby('utdate').filter(lambda x: x['dark'].sum()==2)\n",
    "dark_late = df[df.quarter.isin([3,4])].groupby('utdate').filter(lambda x: x['dark'].sum()==2)\n",
    "dark_all = df.groupby('utdate').filter(lambda x: x['dark'].sum()==4)\n",
    "\n",
    "\n",
    "dfrem = df.copy()\n",
    "rs = 0 \n",
    "nstart_quarters = len(df)\n",
    "nkpfcc_quarters = 365 \n",
    "ndark_quarters = 365 # remove this many quarters for darktime\n",
    "\n",
    "\n",
    "dfrem = dfrem.drop(dark_all.index,errors='ignore')\n",
    "dfrem = dfrem.drop(dark_early.index,errors='ignore')\n",
    "dfrem = dfrem.drop(dark_late.index,errors='ignore')\n",
    "\n",
    "print(\"{}/{} quarters remain after dropping dark all, early, late\".format(len(dfrem),nstart_quarters))\n",
    "\n",
    "# schedule 1/8 of nights are first half\n",
    "\n",
    "for i in range(4):\n",
    "    if i==0:\n",
    "        # schedule 1/8 of nights as Q1+Q2\n",
    "        nnights_select = int(0.125 * nkpfcc_quarters)\n",
    "        dates = dfrem.utdate.drop_duplicates().sample(nnights_select,random_state=rs)\n",
    "        df.loc[df.utdate.isin(dates) & (df.quarter.isin([1,2])),'kpfcc'] = True\n",
    "\n",
    "    if i==1:\n",
    "        # schedule 1/8 of nights as Q3+Q4\n",
    "        nnights_select = int(0.125 * nkpfcc_quarters)\n",
    "        dates = dfrem.utdate.drop_duplicates().sample(nnights_select,random_state=rs)\n",
    "        df.loc[df.utdate.isin(dates) & df.quarter.isin([3,4]),'kpfcc'] = True\n",
    "\n",
    "    if i==2:\n",
    "        # schedule 1/5 of nights as Q1 \n",
    "        nnights_select = int(0.25 * nkpfcc_quarters)\n",
    "        dates = dfrem.utdate.drop_duplicates().sample(nnights_select,random_state=rs)\n",
    "        df.loc[df.utdate.isin(dates) & df.quarter.isin([1]),'kpfcc'] = True\n",
    "\n",
    "    if i==3:\n",
    "        # schedule 1/4 of nights as Q4\n",
    "        nnights_select = int(0.25 * nkpfcc_quarters)\n",
    "        dates = dfrem.utdate.drop_duplicates().sample(nnights_select,random_state=rs)\n",
    "        df.loc[df.utdate.isin(dates) & df.quarter.isin([4]),'kpfcc'] = True\n",
    "\n",
    "    dfrem = dfrem[~dfrem.utdate.isin(dates)]\n",
    "    \n",
    "\n",
    "print(\"{}/{} quarters are KPF-CC\".format(df.kpfcc.sum(),nstart_quarters))\n",
    "print(\"{}/{} nights KPF-CC onsky\".format((df.groupby('utdate')['kpfcc'].sum()>0).sum(),nnights))\n",
    "\n",
    "# df[savecols].to_csv('kpfcc_schedule_v2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7af581db",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/bfulton/anaconda3/lib/python3.8/site-packages/xarray/plot/plot.py:1451: FutureWarning: Conversion of the second argument of issubdtype from `str` to `str` is deprecated. In future, it will be treated as `np.str_ == np.dtype(str).type`.\n",
      "  and not np.issubdtype(x.dtype, str)\n",
      "/Users/bfulton/anaconda3/lib/python3.8/site-packages/xarray/plot/plot.py:1466: FutureWarning: Conversion of the second argument of issubdtype from `str` to `str` is deprecated. In future, it will be treated as `np.str_ == np.dtype(str).type`.\n",
      "  and not np.issubdtype(y.dtype, str)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f93786f63d0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAARwCAYAAAB+ejU9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABFCUlEQVR4nO39f5TleVkf+r6fdFDsAAGDmGYGg4nTngA2IAS513MTAocfmiDmqPciJ4AeXUSuJOg1hx96icTgChqPSVAjdy56AY/ReAR0iJBahMBRPGAYYGgcJ3aPRGHsWRLUCKQjZrqf+0ftMZWm9nTVTNX+1Lc+r9datbr2d3/3/n6+fil4fJ7P5/lUdwcAgPn8idEDAABgDIEgAMCkBIIAAJMSCAIATEogCAAwKYEgAMCkBIIAABtUVT9eVR+rql9d835V1auq6taqOltVX3pYYxEIAgBs1muTPO0u3v+KJNetfp6X5EcPayACQQCADeruX0zye3dxyjOSvL63vSfJ/avq1GGMRSAIAHC0XJPkozte37Y6duD+5GF86XHwWfXZfe/8qdHDAICN+GR+/+Pd/Xmbut5T/+qf6t/9vUubutxGve/sp29O8oc7Dl3f3dfv4ytql2OHsiewQHCNe+dP5cvqSaOHAQAb8a/7Z39rk9f73d+7lH+79QWbvOTGnDh1/g+7+7H34CtuS/KQHa+vTXLhno1qd0rDAABHyw1JnrNaPfz4JH/Q3bcfxoVkBAEANqiqfirJE5I8sKpuS/LdSe6VJN396iRvSfKVSW5NcjHJNx7WWASCAMDGdZLLuTx6GEN099df5f1O8q2bGIvSMADApASCAACTEggCAExKIAgAMCmLRQCAATqXes7FIkeJjCAAwKRkBAE4EFsXPjh6CNwDJ06NHgEjyAgCAExKRhAA2LjthtI9ehjTEwgCR4ryIsDmKA0DAExKIAgAMCmlYeBIeeqDHzl6CNxNyvrs1+XoIziajCAAwKQEggAAk1IaBo4U5UWAzREIAgAb1+lcan0ER1MaBgCYlIwgcKRYNQyjnB89AAaQEQQAmJRAEABgUkrDABwIK76X7cSpzV/zciwWGU1GEABgUgJBAIBJKQ0DcCCs+F46q4ZnJBAEADauk1wyR3A4pWEAgEkJBAEAJiUQBACYlDmCAMAQ+giOJxAE4EBoKL1sIxpKM57SMADApASCAACTUhoG4EBoKL10m20o3UkutTmCo8kIAgBMSiAIADAppWEADoRVw8tm1fCcZAQBACYlIwgADHF59AAQCAJwMKwaXrrNrhrmaFAaBgCYlEAQAGBSSsNrnD5zMVtbVsAtkfIUwNHX6VyKhtKjyQgCAExKIAgAMCml4TXOnT2pxAiwDxpKL5uG0nMSCAIAm9fJJVMEh1MaBgCYlIwgcKQoLwJsjowgAMCkZAQBgI3r2Gv4KJARBACYlEAQAGBSAkEAgEmZIwgcKRq5wyjnRw+AAQSCAMAAlUup0YOYntIwAMCkZAQBOBCagS+bvYbnJCMIADApGUEAYOM6yeUePQoEgmucPnMxW1vKHLBpVg0vl2e3dFYNz0hpGABgUgJBAIBJKQ2vce7sSWUOgH2wanjZRqwa1kdwPBlBAIBJCQQBACYlEAQAmJQ5ggDAxnXMETwKZAQBACYlI7iGhtLLZbU3jOFvb+k0lJ6RjCAAwKQEggAAk1IaBgCGuNwWi4wmIwgAMCmBIADApASCAACTMkcQANg4DaWPBhlBAIBJCQQBACYlEAQAmJQ5gmucO3vSdkkAcEg6lUvyUcN5AgAAkxIIAgBMSiAIADApcwQBgCHsNTzeoWUEq+ohVfWOqrqlqm6uqheujn9uVb2tqs6v/n3A6viTq+p9VfWh1b9PXB0/WVW/UFX/bvU9r7yLaz5m9flbq+pVVVWr43+5qt5fVXdU1dce1j0DACzJYWYE70jyHd39/qq6b5L3VdXbknxDkrd39yur6iVJXpLkxUk+nuTp3X2hqh6RZCvJNavv+oHufkdVfVaSt1fVV3T3W3e55o8meV6S9yR5S5KnJXlrko+srvt3D+leAaa3deGDo4fAPXDi1OgRMMKhZQS7+/bufv/q908muSXbgd0zkrxuddrrknz16pwPdPeF1fGbk9y7qj67uy929ztW5/xRkvcnufbK61XVqST36+53d3cnef2O7/7N7j6b5PJh3CsAwBJtZLFIVT00yaOT/EqSz+/u25PtYDHJg3b5yNck+UB3f/qK77l/kqcnefsun7kmyW07Xt+W/5pR3Os4n1dVN1bVjf8ln776BwAAFuzQF4tU1X2SvCHJt3X3J1bT9u7q/Icn+b4kT7ni+J9M8lNJXtXdH97to7sc6/2MtbuvT3J9kjz2kffurS1lDtg0jdyXy7NbuvMbvVonubTr/3SzSYeaEayqe2U7CPzJ7n7j6vDvrMq4d5ZzP7bj/GuTvCnJc7r7N674uuuTnO/uf7I690RV3bT6+Z5sZwB3loyvTXIhAADs6tAygqsVuz+W5Jbu/sEdb92Q5LlJXrn69+dX598/yS8keWl3//IV3/WKJH86yTffeay7LyV51BXnfbKqHp/tEvRzkvzQgd4UAMAxcpil4S9P8uwkH6qqm1bHvjPbAeDPVNU3ZXs179et3ntBki9K8rKqetnq2FOSfFaS70ry75K8f1Va/uHufs0u13x+ktcm+ZxsrxZ+a5JU1V/KdqbxAUmeXlV/v7sffmB3ypGiPAUAe3NogWB3vyu7z9tLkiftcv4rkrxizfl7mkTQ3TcmecQux9+bXVYaAwCjVC61Dc5G8wQAACYlEAQAmJRAEABgUofeRxAA4Eqd5LJ81HACwTXOnT1p9SkMYL9aGMNew3MSigMATEogCAAwKaXhNU6fuRh7DcPmmZIBo2x2r+HEXsNHgYwgAMCkBIIAAJNSGubYUVoEgL2REQQAmJSMIACwcd2VSy0fNZonAAAwKYEgAMCkBIIAAJMyR3ANew0DwOG6rKH0cDKCAACTEggCAExKaXgNew0vl5I+AOyNQBAA2LhOcklhcjhPAABgUgJBAIBJCQQBACZljiAAMIC9ho8CgeAaGkoDAMedUBwAYFICQQCASSkNA3Agti5owr9kJ06NHgEjCAQBgI3rJJcVJofzBAAAJiUjCBwpyosAmyMjCAAwKRlBAGCIS12jhzA9GUEAgEkJBAEAJiUQBACYlDmCwJFij28Y5fxGr9apXJKPGs4TAACYlEAQAGBSSsNrnD5zMVtbGtsukdIiAOyNQBAAGOJyK0yO5gkAAExKIAgAMCmBIADApMwRBAA2rhN9BI8AgeAa586etPoUYB+2Lui0sGQnTo0eASMIxQEAJiUQBACYlNIwAAfCdJql2+xewxwNAkEAYOM6lUtdo4cxPaVhAIBJyQiuYa/h5VKeAoC9kREEAJiUjCAAMMRl+ajhBIJraCgNY2hKDGNoKD0noTgAwKQEggAAk1IaXsOq4eVS0l82zw9G2WxD6e7kUstHjeYJAABMSiAIADAppWGOHatOl01pGGBzBIIAwACVy7HX8GhKwwAAk5IRXENDaYD9MS1j2TSUnpOMIADApASCAACTUhpeQ0Pp5VLShzH87S3dhhtKR0Ppo8ATAACYlEAQAGBSSsNrWDUMsD9WDS+bVcNzEggCAENcUpgczhMAAJiUjCAAB8J0mqXb7KphjgYZQQCASckIAgAb16lc7ho9jOnJCAIATEogCAAwKYEgAMCkzBEEAIbQR3A8TwAAYFICQQCASSkNr3H6zMVsbdk3c4k0tQWAvZERBACYlIwgALBxneRyy0eNJhBc49zZk0qMAPuwdcF0miU7cWr0CBhBKA4AMCmBIADApJSGgSNFeRFmUbmUGj2I6ckIAgBMSiAIADAppeE1NJReLqu9l83zg1HOjx4AAwgEAYCN00fwaPAEAAAmJSO4hobSAPtjxfeyaSg9JxlBAIBJyQgCAEPoIzieQHANq4ZhDFMylsuzWzqrhmekNAwAMCmBIADApASCAACTMkcQANi47tJQ+gjwBAAAJiUjuIaG0gD7o6H0smkoPScZQQCASckIAgBDXDJHcDhPAABgUgJBAIBJCQQBACZljiAAB0KnhaXb7F7DneRyaqPX5DPJCAIATEogCAAwKaXhNU6fuZitLc1RYdOUFwE2RyAIAAxQ+ggeAZ4AAMCkBIIAABtUVU+rql+vqlur6iW7vP+nq+rNVfXBqrq5qr7xsMYiEAQA2JCqOpHkR5J8RZKHJfn6qnrYFad9a5Jf6+5HJnlCkv+1qj7rMMYjEAQA2JzHJbm1uz/c3X+U5KeTPOOKczrJfauqktwnye8lueMwBmOxyBrnzp60ehFgH7Yu6LSwZCdObfZ6neRyH9uG0g+sqht3vL6+u69f/X5Nko/ueO+2JF92xed/OMkNSS4kuW+S/0d3Xz6MgQoEAQAO1se7+7Fr3tst+u0rXj81yU1JnpjkLyR5W1X9Und/4uCGuE1pGABgc25L8pAdr6/NduZvp29M8sbedmuSf5/kvzuMwcgIAkeK8iJwzL03yXVV9YVJfjvJM5M864pzPpLkSUl+qao+P8kXJ/nwYQxGIAgADHFpwsJkd99RVS9IspXkRJIf7+6bq+pbVu+/Osk/SPLaqvpQtkvJL+7ujx/GeASCd0FmYpks8lk+z3CZti580LNbtPOjBzCN7n5LkrdccezVO36/kOQpmxjLfKH4Hp0+c3H0EGBKAonl8uxgeQSCAACTUhpeQx9BGMOUDBhj830E6zj3EVwMGUEAgEkJBAEAJqU0DBwppmTAKFYNz0ggCAAMcVlhcjhPAABgUjKCABwIK76XbdOrhjkaZAQBACYlEAQAmJTSMACwcd3JJQ2lh5MRBACYlEAQAGBSSsNrnD5zMVtbVsAtkYbEMIa/vaXTUHpGAkEAYIjL5ggOpzQMADApGcE1zp09qcwBA2hKDGNoKD0nGUEAgEnJCAIAG9epXG75qNEEgsCRYkrGcinrw/IIxQEAJiUQBACYlNIwcKQoL8I8LkUfwdFkBAEAJiUQBACYlNIwcKRYNQyj2Gt4RjKCAACTkhEEADauk1xui0VGEwiucfrMxWxtWb0Im6Y0DLA5SsMAAJMSCAIATEppeI1zZ08qUcEAGkrDGCdObfqKlcstHzWaJwAAMCmBIADApJSGgSPFlIzlUtaH5REIAgBDXI4+gqMpDQMATEpGEIADoay/dPYanpGMIADApGQEAYCN604u2Wt4OBlBAIBJCQQBACYlEAQAmJQ5gmucPnMxW1uaoy6RlYsAsDcCQQBgiMutMDmaJwAAMCkZQY4d+50um9I+wObICAIATOrQAsGqekhVvaOqbqmqm6vqhavjn1tVb6uq86t/H7A6/uSqel9VfWj17xN3fNf3VtVHq+pTV7nmrudV1V+uqvdX1R1V9bWHcb8AwN51Kpf7eP4syWGWhu9I8h3d/f6qum+S91XV25J8Q5K3d/crq+olSV6S5MVJPp7k6d19oaoekWQryTWr73pzkh/O1TdCXHfeR1bX/bt7Hfy5syeVqAD2wbSMZTtxavQIGOHQAsHuvj3J7avfP1lVt2Q7sHtGkiesTntdkncmeXF3f2DHx29Ocu+q+uzu/nR3vydJqu46yl53Xnf/5ur45XtyTwAAx8lG5ghW1UOTPDrJryT5/FWQeGew+KBdPvI1ST7Q3Z/exPjuVFXPq6obq+rG/5KNXhoAYOMOfdVwVd0nyRuSfFt3f+JqWb2qeniS70vylMMe25W6+/ok1yfJYx9579ZQepmU9GEMf3tLd7XZVwfvcpY1n+44OtSMYFXdK9tB4E929xtXh3+nqk6t3j+V5GM7zr82yZuSPKe7f+Mq332iqm5a/XzP4dwBAMDxdWgZwdpO/f1Yklu6+wd3vHVDkucmeeXq359fnX//JL+Q5KXd/ctX+/7uvpTkUQc7agCAeRxmafjLkzw7yYeq6qbVse/MdgD4M1X1Tdlezft1q/dekOSLkrysql62OvaU7v5YVX1/kmclOVlVtyV5TXe//MoLrjuvqv5StjOND0jy9Kr6+9398LsavFXDMIaVpzCGVcNzOsxVw+9K1hb/n7TL+a9I8oo13/WiJC/awzV3Pa+735vk2qt9HgDYjE4W13PvOLKzCADApOw1DBwppmQsl7I+LI+MIADApASCAACTUhpe4/SZi9FQepmUFmEMf3tLN6ChdMtHjeYJAABMSiAIADAppeE1NJSGMaw8hTE0lJ6TQBAA2LwuDaWPAKVhAIBJyQiuYdUwjGFKBoyy+VXDjCcjCAAwKRlBAGDjOsnlmCM4mkBwDauGYQyrhmEMq4bnpDQMADApgSAAwKSUhoEjxZSM5VLWZ7/0ERxPRhAAYFICQQCASQkEAQAmJRAEAJiUxSIAwMZ1LBY5CgSCa9hreLmsOoUx/O0tnb2GZ6Q0DAAwKYEgAMCklIbXsNcwwP5oKL1sI/YaNkdwPBlBAIBJCQQBACalNLyGVcMwhikZy+XZLZ1VwzMSCAIAG9cpcwSPAKVhAIBJyQhy7ChPAcDeyAgCAExKRhAAGOJyzBEcTSC4hobSAPujofSyjWgozXhKwwAAkxIIAgBMSml4DQ2ll0tJH8bwt7d0GkrPSCAIAGxeR0PpI0BpGABgUjKCa1g1DAAcdzKCAACTkhEEADauY47gUSAQBOBAaCi9bBpKz0lpGABgUgJBAIBJKQ0DAEOYIziejCAAwKQEggAAkxIIAgBMyhxBAGDjOmWO4BEgIwgAMCkZQQAOhP3Zl+786AEwgIwgAMCkBIIAAJNSGgYAhmiLRYaTEQQAmJRAEABgUkrDAByIrQsfHD0E7oETp0aPgBEEggDAEJdjjuBoSsMAAJOSEVzj9JmL2dpS5lgiTW1hDH97S6eh9IxkBAEAJiUjCABsXHdyWR/B4QSCHDtWLi6b8iLA5igNAwBMSiAIADAppeE1zp09qUQFAIfIXsPjyQgCAExKIAgAMCmBIADApASCAACTslgEABigNJQ+AmQEAQAmJRAEAJiUQBAAYFLmCK5x+szFbG3ZsxY2TSN3mIeG0uPJCAIATEogCAAwKaVhjh2lRQDYG4EgALBxnegjeAQIBDl2ti5Y5LNkMroAm2OOIADApASCAACTUhrm2FFaBFiATrpHDwIZQQCASQkEAQAmpTS8xrmzJ5UYAfbBiv1lO3Fq9AgYQUYQAGBSMoIAwBCXo6H0aALBNU6fuZitLWWOJVLShzH87S3d+dEDYAClYQCASQkEAQAmpTQMAGxcJ+k2R3A0GUEAgEkJBAEAJqU0zLGjqe2yWXkKsDkCQQBggMplcwSHUxoGAJiUQBAAYFICQQCASZkjCAAM0T16BAgE1zh39qTViwD7YMX+sp04NXoEjKA0DAAwKYEgAMCklIbXOH3mYra2lDlg00zJWC7PbunOjx4AAwgEAYAhWkPp4ZSGAQAmJSPIsaM8BQB7IyMIADApGUEAYOO6zRE8CmQEAQAmJRAEAJiUQBAAYFLmCAIAQ1w2R3A4GUEAgEkJBAEAJqU0vMa5syc1JgbYh60L9mdfshOnRo+AEQSCAMAQ3aNHgNIwAMCkZATXOH3mYra2lDmWSEkfxvC3t3TnRw+AAWQEAQAmJRAEAJiU0vAaVg0D7I9Vw8s2YtVwayg9nIwgAMCkBIIAAJMSCAIATMocQQBg4zpljuARICMIADApgSAAwKQEggAAkzJHEAAYokcPABlBAIBZCQQBACYlEAQAmJQ5ggDA5rW9ho8CGUEAgEkJBAEAJqU0DMCBeOqDHzl6CNwj50cPgAFkBAEAJiUjCACMoaP0cDKCAACTEggCAExKIAgAMClzBAE4EFsXPjh6CNwDJ05t/poaSo8nIwgAMCmBIADApASCAACTMkcQABii9REcTkYQAGBSAkEAgA2qqqdV1a9X1a1V9ZI15zyhqm6qqpur6v84rLEoDQMAbEhVnUjyI0menOS2JO+tqhu6+9d2nHP/JP8sydO6+yNV9aDDGo9AEADYuM60fQQfl+TW7v5wklTVTyd5RpJf23HOs5K8sbs/kiTd/bHDGoxAcI3TZy5ma0tz1CV66oMfOXoIMCV/e0t3fvQAjpMHVtWNO15f393Xr36/JslHd7x3W5Ivu+Lzp5Pcq6remeS+Sf5pd7/+MAYqEAQAOFgf7+7HrnlvtzToleun/2SSxyR5UpLPSfLuqnpPd587wDH+8YUAANiM25I8ZMfra5Nc2OWcj3f3f0ryn6rqF5M8MolAcFPOnT2pzAEAHLT3Jrmuqr4wyW8neWa25wTu9PNJfriq/mSSz8p26fgfH8ZgBIIAwOZ1kgkXi3T3HVX1giRbSU4k+fHuvrmqvmX1/qu7+5aq+ldJzia5nOQ13f2rhzEegSAAwAZ191uSvOWKY6++4vU/SvKPDnssGkoDAExKIAgAMCmlYQBgiL6yaQobJxBcQ0Pp5bLaGwD2RmkYAGBSAkEAgEkpDa+hoTTA/mxdMJ1myU6cGnBRcwSHkxEEAJiUQBAAYFJKwwAcCNNplu786AEwgEAQABig0hPuNXzUKA0DAExKRnANDaWXS3kKAPbm0DKCVfWQqnpHVd1SVTdX1QtXxz+3qt5WVedX/z5gdfzJVfW+qvrQ6t8n7viu762qj1bVp65yzcesPn9rVb2qqmp1/M9V1dur6mxVvbOqrj2s+wYAWIrDLA3fkeQ7uvsvJnl8km+tqocleUmSt3f3dUnevnqdJB9P8vTu/pIkz03yEzu+681JHreHa/5okucluW7187TV8R9I8vruPpPke5L8w3tyYwAAx8GhlYa7+/Ykt69+/2RV3ZLkmiTPSPKE1WmvS/LOJC/u7g/s+PjNSe5dVZ/d3Z/u7vckySrBt6uqOpXkft397tXr1yf56iRvTfKwJN++OvUdSX7uHt8gR5amtsumtA8T0VB6uI0sFqmqhyZ5dJJfSfL5qyDxzmDxQbt85GuSfKC7P72Py1yT5LYdr29bHUuSD66+M0n+RpL7VtWf2cd3AwAcO4ceCFbVfZK8Icm3dfcn9nD+w5N8X5K/td9L7XLszv9f4+8m+StV9YEkfyXJb2e7dH3ltZ9XVTdW1Y3/4Xcv7fPyAADLcqirhqvqXtkOAn+yu9+4Ovw7VXWqu29flXM/tuP8a5O8Kclzuvs3rvLdJ5K8b/XyhmzPD9y5COTaJBeSpLsvJPkfV5+7T5Kv6e4/uPI7u/v6JNcnyf3qc1uJCgA4zg4tEFyt2P2xJLd09w/ueOuGbC8GeeXq359fnX//JL+Q5KXd/ctX+/7uvpTkUVdc85NV9fhsl6Cfk+SHVscfmOT3uvtykpcm+fF7cm8AwD3U0VD6CDjM0vCXJ3l2kidW1U2rn6/MdgD45Ko6n+TJq9dJ8oIkX5TkZTvOf1CSVNX3V9VtSU5W1W1V9fI113x+ktckuTXJb2R7oUiyvTjl16vqXJLPT/K9B3yvAACLc5irht+V3eftJcmTdjn/FUlesea7XpTkRXu45o1JHrHL8Z9N8rNX+zwAwExsMQcAMClbzAEAY+gjOJxAEDhSNASHMU6cGj0CRlAaBgCYlEAQAGBSSsMAwCD6CI4mIwgAMCmBIADApASCAACTEggCAEzKYhEAYAwNpYcTCAJHylMf/MjRQ4BJnR89AAZQGgYAmJRAEABgUkrDABwI+0Qv25C9hs0RHE5GEABgUgJBAIBJKQ0DcCCs+F46q4ZnJBAEADavk3SNHsX0lIYBACYlIwgcKVaewhhDVg0znIwgAMCkZAQBgCFaH8HhBIJrnD5zMVtbSlRLZOXisnl+MIpVwzNSGgYAmJRAEABgUkrDa5w7e1KJCmAfrPheNnsNz0lGEABgUgJBAIBJKQ0DcCBMp1k6q4ZnJCMIADApGUEAYIyu0SOYnkAQOFKsPIUx7DU8J6VhAIBJCQQBACalNLyGvYaXy8rFZfP8YJTNrxouDaWHkxEEAJiUQBAAYFICQQCASZkjCABsXq9+GEpGkGNHHzoA2BuBIMeOVacAsDcCQQCASZkjuMa5sydllmAApX0YY/NbzJW9ho8AGUEAgEkJBAEAJqU0vIYt5mAMUzJglM1vMcd4MoIAAJOSEQQAxtBQejgZQQCASQkEAQAmJRAEAJiUOYJraCgNAIfMHMHhZAQBACYlEAQAmJTS8BoaSi+Xkj4A7I1AEAAYwxzB4ZSGAQAmJSO4hlXDAPuzdcF0miU7cWr0CBhBRhAAYFIyggDA5nWSrtGjmJ5AEIADYTrN0p0fPQAGUBoGAJiUQBAAYFICQQCASZkjCAAMURpKDycjCAAwqT0HglX1OVX1xYc5GAAANmdPgWBVPT3JTUn+1er1o6rqhkMcFwAAh2yvcwRfnuRxSd6ZJN19U1U99HCGBABMwRzB4fZaGr6ju//gUEcCAMBG7TUQ/NWqelaSE1V1XVX9UJL/8xDHBQDAHlXV26rq/jteP6Cqtq72ub0Ggn87ycOTfDrJP0/yB0leeDfGCQDAwXtgd//HO1909+8nedDVPrTXOYJ/rbu/K8l33Xmgqr4uyf++z0ECcExtXfjg6CFwD5w4NXoE3EOXq+oLuvsjSVJVfy57mIW514zgS/d4DACAzfuuJO+qqp+oqp9I8ovZQ6x2lxnBqvqKJF+Z5JqqetWOt+6X5I57MFgAAA5Id/+rqvrSJI9PUkm+vbs/frXPXa00fCHJjUm+Ksn7dhz/ZJJvv5tjXYTTZy5ma0uZY4me+uBHjh4CTMnf3tKdHz0A7oGq+htJ/k13/8vV6/tX1Vd398/d1efuMhDs7g8m+WBV/fPu/i8HNloAYHr2Gj5Q393db7rzRXf/x6r67iQ/d1cf2utikYdW1T9M8rAk995xkT9/NwYKAMDB2m3dx1XjvL0Ggv+/JN+d5B8n+atJvjHb9edj69zZk8ocAMBS3FhVP5jkR7K9Wvhv57+d1rerva4a/pzufnuS6u7f6u6XJ3ni3R0pAAAH6m8n+aMk/yLJzyT5z0m+9Wof2mtG8A+r6k8kOV9VL0jy29lDk0IAAA5PVf1Edz87yTd390v2+/m9BoLfluRkkr+T5B9kOxv43P1ebEmsGl4uJX2AhehjPctsUx6zah79P1fV63PF1L3u/r27+vCeAsHufu/q109le34gAADjvTrJv0ry57M9J3BnINir42tdraH0m3MX25N091fteZgAAByo7n5VkldV1Y929/P3+/mrZQR/YPXv/5jkzyb531avvz7Jb+73Ykti1TCMYb9aGMNew8vW3c9f7Szy32c7ifeu7v7A1T53tYbS/0eSVNU/6O6/vOOtN1fVL96TAQMAE+vcRc2R/aqqlyX5vyd54+rQa6vqf+/uV9zV5/a6WOTzqurPd/eHVxf7wiSfd7dHCwDAQXpWkkd39x8mSVW9Msn7kxxIIPjtSd5ZVR9evX5okr9198YJsJ4pGTCKvYYX7jezvfvbH65ef3aS37jah/YaCL4jyXVJ/rvV63+3z8EBAHB4Pp3k5qp6W7aL7k9O8q6qelWSdPff2e1Dew0E393dX5rkj2dxV9X7k3zpPRoyADAvcwQP0ruSvGnH63cmuV+ST9zVh67WPubPJrkmycnVSpQ73S/bDaYBIIkV30tn1fDiPSvJc7v7Q0lSVV+f5G9295fd1YeulhF8apJvSPLg/NdWMknyySQvvdtDBQDgIH1tkp+tqv8p2y1knpPkKVf70NXax7wuyeuq6kPZTjHu7Fb9JflvU5AAAAzQ3R+uqmcm+bkkH03ylO7+z1f73F7nCL52x+/3TvLXk9yyzzECXJXyIsyjzBG8x1bJup3/l/zcJCeS/EpVpbvP3NXn97rX8P96xUV/IMkN+xwrAAAH66/fkw/vNSN4pZO5yibGAAAcru7+rXvy+T0FglekHU9ke1eR77knFz7qTp+5mK0tJSrYNA2lYRQNpWe014zgzrTjHUl+p7vvOITxAACwIXudI3iP0o4AAJ/BYpHh7u4cwWPv3NmTSlQAwLH2J0YPAACAMQSCAACTUhoG4EBoBr5sQ/YaNkdwOBlBAIBJCQQBACalNLyGhtLLZbU3jOFvb+k0lJ6RQBAA2Ljq7R/GUhoGAJiUQBAAYFICQQCASZkjCACM0TV6BNOTEQQAmJRAEABgUgJBAIBJCQQBACZlsQgAMIaG0sPJCAIATEogCAAwKYEgAMCkzBEEAIYocwSHEwiuce7syTz1wY8cPQyYztaFD44eAkzpxKnRI2AEpWEAgEkJBAEAJqU0vMbpMxeztaVEtURK+svm+cEo5zd/SXMEh5MRBACYlEAQAGBSSsNrWDUMABx3AkEAYPNaH8GjQGkYAGBSMoJrWDUMY5iSAbA5MoIAAJMSCAIATEppmGNHaRFgISwWGU5GEABgUgJBAIBJKQ1z7GxdsNp7yZT2ATZHIAgAjGGO4HBKwwAAk5IRXMNewzCG0j6MceLU6BEwgowgAMCkZAQBgCHKHMHhBIJr2GsYxjAlA0Y5P3oADKA0DAAwKYEgAMCklIbXsGoYADjuZAQBACYlEAQAmJTSMAAHQjPwZdNQek4yggAAk5IRBADG0FB6OIEgAAdCp4Wl01B6RkrDAACTEggCAExKaRgA2LxOyhzB4WQEAQAmJRAEAJiU0vAap89czNaW5qhLZOUiAOyNQBAAGMMcweGUhgEAJiUjuMa5syeVGGEA+9XCGPYanpOMIADApGQEAYAxzBEcTiAIHCmmZCyXsj4sj9IwAMCkBIIAAJNSGl5DQ2kYQ2l4uTy7pTs/egAMIBAEADaukpTFIsMpDQMATEpGcA0NpQH2x6rhZdNQek4yggAAk5IRBADGMEdwOIHgGlYNwximZCyXZ7d0Vg3PSGkYAGBSAkEAgEkpDa9h1TDA/lg1vGwbXzXc+ggeBYeWEayqh1TVO6rqlqq6uapeuDr+uVX1tqo6v/r3AavjT66q91XVh1b/PnHHd31vVX20qj51lWs+ZvX5W6vqVVVVq+P/uKpuWv2cq6r/eFj3DQCwFIdZGr4jyXd0919M8vgk31pVD0vykiRv7+7rkrx99TpJPp7k6d39JUmem+QndnzXm5M8bg/X/NEkz0ty3ernaUnS3d/e3Y/q7kcl+aEkb7yH9wYAsHiHVhru7tuT3L76/ZNVdUuSa5I8I8kTVqe9Lsk7k7y4uz+w4+M3J7l3VX12d3+6u9+TJKsE366q6lSS+3X3u1evX5/kq5O89YpTvz7Jd19t/FYNL5eSPozhb2/prBqe0UbmCFbVQ5M8OsmvJPn8VZCY7r69qh60y0e+JskHuvvT+7jMNUlu2/H6ttWxneP4c0m+MMm/WTPO52U7o5gvuMb0SQA4VOYIDnfoq4ar6j5J3pDk27r7E3s4/+FJvi/J39rvpXY5duV/xJ6Z5Ge7+9JuX9Dd13f3Y7v7sZ/3Z07s8/IAAMtyqIFgVd0r20HgT3b3nfPyfmdVxr2znPuxHedfm+RNSZ7T3b9xle8+sWMByPdkOwN47Y5Trk1y4YqPPTPJT92TewIAOC4Orf65WrH7Y0lu6e4f3PHWDdleDPLK1b8/vzr//kl+IclLu/uXr/b9q6zeo6645ier6vHZLkE/J9sLQ+5874uTPCDJu+/2TQEAHCOHORHuy5M8O8mHquqm1bHvzHYA+DNV9U1JPpLk61bvvSDJFyV5WVW9bHXsKd39sar6/iTPSnKyqm5L8prufvku13x+ktcm+ZxsLxLZuVDk65P8dHfvaUaCPoIA+6OP4LJtvI8gR8Jhrhp+V3aft5ckT9rl/FckecWa73pRkhft4Zo3JnnEmvdefrXPAwAbZLHIcLaYAwCYlB4pABwI02mWTh/BGckIAgBsUFU9rap+fbUl7kvu4ry/VFWXquprD2ssMoIAwBA14RzBqjqR5EeSPDnbre/eW1U3dPev7XLe9yXZOszxyAgCAGzO45Lc2t0f7u4/SvLT2d5+90p/O9u9mD+2y3sHRiAIAHCwHlhVN+74ed6O965J8tEdr3fbEveaJH8jyasPe6BKwwAAB+vj3f3YNe/tZUvcf5Lkxd19aXt/jsMjEFzj9JmL2drSHHWJrFwEWIgJ5whmOwP4kB2vd9sS97FJfnoVBD4wyVdW1R3d/XMHPRiBIADA5rw3yXVV9YVJfjvJM7O9e9of6+4vvPP3qnptkn95GEFgIhAEANiY7r6jql6Q7dXAJ5L8eHffXFXfsnr/0OcF7iQQXMNewwD7Y6/hZbPX8OZ091uSvOWKY7sGgN39DYc5FoEgALB5nVnnCB4p2scAAExKRnANq4aXS0kfxvC3t3T2Gp6RjCAAwKQEggAAk1IaXsOqYRjDylMYY8Sq4bJYZDgZQQCASQkEAQAmpTS8hlXDy6Wkv2yeH4xi1fCMBIIAwBjmCA6nNAwAMCkZQY4dq06XTWkYYHNkBAEAJiUjCAAMoY/geALBNTSUBgCOO6VhAIBJCQQBACalNLyGhtIwhikZMBFzBIeTEQQAmJRAEABgUkrDa1g1DLA/mrkv24lTo0fACDKCAACTkhEEADavY7HIESAQBOBAmE6zdOdHD4ABlIYBACYlEAQAmJTSMHCkWHkKY2x61XCtfhhLRhAAYFICQQCASSkNA0eKlafLpawPyyMQBADG0EdwOKVhAIBJCQQBACYlEAQAmJQ5ggDAEGWO4HACwTVOn7mYrS0r4JbIqlMYw9/e0tlreEZKwwAAkxIIAgBMSml4jXNnTypzAOyDhtLLtum9hjkaBIIAwBgWiwynNAwAMCkZQQAOhOk0S2fV8IxkBAEAJiUjCACMYY7gcDKCAACTEggCAExKIAgAMClzBAGAzeukzBEcTkYQAGBSAkEAgEkJBAEAJmWOIAAwhjmCw8kIAgBMSkYQgAOxdeGDo4fAPXDi1OgRMIKMIADApASCAACTUhoG4EA89cGPHD0E7pHzG7+ihtLjyQgCAExKIAgAMCmlYQAOhFXDy2bV8JwEggDAGOYIDqc0DAAwKRnBNU6fuZitLWUO2DQrT5fLs1u6za8aZjwZQQCASckIAgBD6CM4nkBwjXNnTypzAOyDVcPLZtXwnJSGAQAmJRAEAJiU0vAaVg0vl5I+jOFvb+k2vGq4o4/gESAjCAAwKYEgAMCklIbXsGoYxrDyFMawanhOMoIAAJOSEQQAxrBYZDgZQQCASQkEAQAmJRAEAJiUOYJraCi9XFZ7L5vnB6NstqF0JSlzBIeTEQQAmJRAEABgUgJBAIBJmSMIAIxhjuBwMoIAAJOSEVzDXsMA+2Of6GWz1/CcZAQBACYlIwgADFFtkuBoMoIAAJMSCAIATEogCAAwKXME17DX8HJZ7Q1j+Ntbus3uNczRIBAEADavo6H0EaA0DAAwKRnBNTSUBtgfDaWXTUPpOckIAgBMSkYQABiizBEcTkYQAGBSAkEAgEkJBAEAJmWOIHCkWHkKEzFHcDgZQQCASQkEAQAmpTQMHCkauS+Xsj4sj0AQABhCH8HxlIYBACYlIwjAgVDWX7rzowfAADKCAACTEggCAExKaXiN02cuZmvLCrglUp4CWAiLRYaTEQQAmJRAEABgUkrDa5w7e1KJEQA41gSCAMDmtYbSR4HSMADApGQEgSPFfrUwxolTo0fACDKCAACTkhEEAMYwR3A4geAaGkovl9Xey+b5wSj2Gp6R0jAAwKQEggAAk1IaBgA2rqKP4FEgIwgAMCmBIADApJSG17DXMMD+aAa+bBpKz0lGEABgUjKCa+gjuFwyuTCGv72lG9BHsK0WGU1GEABgUgJBAIBJCQQBACZljiAAMISG0uPJCAIATEpGcA19BAGA405GEABgUjKCAMDm9eqHoWQEAQAmJRAEAJiUQBAAYFLmCAIAQ9Tl0SNARhAAYFICQQCASQkEAQAmJRAEAJiUxSIAwBgaSg8nIwgAMCmBIADApASCAACTMkcQABiizBEcTkYQAGBSAkEAgEkpDa9x+szFbG19cPQwuBue+uBHjh4CACyCQBAA2LxO0iYJjqY0DAAwKYEgAMCkBIIAAJMyRxAAGEIfwfEEgmucO3vS6lOAfdi6oNPCkp04NXoEjKA0DAAwKYEgAMCklIaBI0V5ESZijuBwMoIAAJMSCAIATEogCAAwKYEgAMCkDi0QrKqHVNU7quqWqrq5ql64Ov65VfW2qjq/+vcBq+NPrqr3VdWHVv8+ccd3fW9VfbSqPnWVaz5m9flbq+pVVVWr499QVf+hqm5a/XzzYd03AHB1le2G0sfxZ0kOc9XwHUm+o7vfX1X3TfK+qnpbkm9I8vbufmVVvSTJS5K8OMnHkzy9uy9U1SOSbCW5ZvVdb07yw0nOX+WaP5rkeUnek+QtSZ6W5K2r9/5Fd7/gwO4OOBQaucMoV/ufWI6jQ8sIdvft3f3+1e+fTHJLtgO7ZyR53eq01yX56tU5H+juC6vjNye5d1V99uq993T37Xd1vao6leR+3f3u7u4kr7/zuwEA+EwbmSNYVQ9N8ugkv5Lk8+8M6lb/PmiXj3xNkg9096f3cZlrkty24/Vt+a8ZxST5mqo6W1U/W1UPWTPO51XVjVV143/Jfi4NALA8h95Quqruk+QNSb6tuz+xmrZ3V+c/PMn3JXnKfi+1y7E7K/VvTvJT3f3pqvqWbGcin/gZJ3dfn+T6JHnsI+/dW1sa2y6R0iLAAnRv/zDUoWYEq+pe2Q4Cf7K737g6/DurMu6d5dyP7Tj/2iRvSvKc7v6Nq3z3iR2LP74n2xnAa3eccm2SC0nS3b+7I7v4/03ymHt+dwAAy3aYq4YryY8luaW7f3DHWzckee7q9+cm+fnV+fdP8gtJXtrdv3y17+/uS939qNXP31uVmT9ZVY9fXfs5O7771I6PflW25ysCAEztMEvDX57k2Uk+VFU3rY59Z5JXJvmZqvqmJB9J8nWr916Q5IuSvKyqXrY69pTu/lhVfX+SZyU5WVW3JXlNd798l2s+P8lrk3xOtlcL37li+O9U1VdleyXz72V75fJdOnf2pBIjAHCsHVog2N3vyu7z9pLkSbuc/4okr1jzXS9K8qI9XPPGJI/Y5fhLk7z0ap8HADZnaT33jiM7iwAATOrQVw0D7MfWBav1YYQTp65+DsePjCAAwKRkBAGAMcwRHE4gCBwpVusvl7I+LI/SMADApASCAACTUhoG4EAo6y/d+dEDYACBIAAwhIbS4ykNAwBMSiAIADApgSAAwKTMEQQANq+TXDZJcDSBIHCkaEoMY9hreE5KwwAAkxIIAgBMSmkYOFI0JV4uZX32zRTB4WQEAQAmJRAEAJiU0vAap89czNaWMgdsmtLwcnl2S2ev4RkJBAGAIew1PJ7SMADApGQE1zh39qQyB8A+WDW8bBpKz0lGEABgUgJBAIBJKQ2vYdUwjGFKxnJ5dks3YNVwWy0ymowgAMCkBIIAAJNSGubYUZ4CgL0RCAIAQ2goPZ7SMADApGQEOXY0tV02pX2AzZERBACYlIwgALB5vfphKBlBAIBJCQQBACYlEAQAmJQ5ghw7Vp0CHH2VpCbda7iqnpbknyY5keQ13f3KK97/n5K8ePXyU0me392H0hJDRhAAYEOq6kSSH0nyFUkeluTrq+phV5z275P8le4+k+QfJLn+sMYjEAQA2JzHJbm1uz/c3X+U5KeTPGPnCd39f3b3769evifJtYc1GKXhNc6dPanECADcHQ+sqht3vL6+u+/M6l2T5KM73rstyZfdxXd9U5K3HvD4/phAEADgYH28ux+75r3a5diukyWr6q9mOxD87w9qYFcSCAIAY1wePYAhbkvykB2vr01y4cqTqupMktck+Yru/t3DGoxAcI3TZy5ma8uetbBppmQAx9x7k1xXVV+Y5LeTPDPJs3aeUFVfkOSNSZ7d3ecOczACQQCADenuO6rqBUm2st0+5se7++aq+pbV+69O8veS/Jkk/6yqkuSOuyg13yMCQQCADerutyR5yxXHXr3j929O8s2bGItAcA2rhgHgcM3aUPoo0UcQAGBSAkEAgEkJBAEAJmWOIACweZ01bZTZJBlBAIBJCQQBACYlEAQAmJQ5ggDAAJ3oIzicjCAAwKQEggAAkxIIAgBMyhzBNU6fuZitrQ+OHgZMxx7fAJsjEAQAhihrRYZTGgYAmJSM4Brnzp5UogLYh60LptMs2YlTo0fACDKCAACTkhEEAMbQUHo4gSAAB8J0mqU7P3oADKA0DAAwKYEgAMCklIbX0FAaxlBehEl0UpdHDwIZQQCASQkEAQAmJRAEAJiUOYIAwBj6CA4nIwgAMCkZQY4dq04BYG9kBAEAJiUQBACYlNLwGufOnlRiBNiHrQua8C/ZiVMDLmqtyHAyggAAkxIIAgBMSmkYOFKUFwE2RyAIAAxRGkoPpzQMADApGUHgSLFaH0Y5P3oADCAjCAAwKRlBAGAMcwSHEwiucfrMxWxtWb24REqLALA3SsMAAJMSCAIATEppeA17DQPsj2bgy7bxvYY7yeUNX5PPICMIADApGUEADoxKypLpIzgjGUEADoQgEJZHIAgAMCmlYQAOjAUjy7XpxSKVTmkoPZyMIAAHQhAIyyMQBACYlNLwGraYWy4T1mEMf3tLZ9XwjASCAMAY5ggOpzQMADApGcE1bDEHABx3MoIAAJOSEQQAxjBHcDiBIAAHQh/BZdt0Q2mOBqVhAIBJCQQBACalNLyGhtLLZbU3jOFvb+k23FC6k1ze7CX5TDKCAACTEggCAExKaXgNDaUB9seq4WWzanhOMoIAAJOSEQQAhigNpYcTCAJwIEynWboNrxrmSFAaBgCYlEAQAGBSSsMAHAirhpdtyKphcwSHkxEEAJiUQBAAYFJKw2vYa3i5rFyEMfztLZ1VwzMSCAIAA7Q5gkeA0jAAwKRkBNew1zDA/lg1vGz2Gp6TjCAAwKRkBAGAzeuYI3gECATXsGp4uZT0YQx/e0tn1fCMlIYBACYlEAQAmJTS8BpWDQPsj1XDy2bV8JwEggDAGJdHDwClYQCASckIrmHV8HIp6cMY/vaWzqrhGckIAgBMSkYQABiiNJQeTiC4hlXDAPtj1fCyWTU8J6VhAIBJCQQBACalNAzAgTCdZukGrBo2R3A4GUEAgEkJBAEAJiUQBACYlDmCAMDmdZLL5giOJiMIADApgSAAwKQEggAAkxIIAgBMymIRAGCA1lD6CJARBACYlEAQAGBSAkEAgEmZIwgAjGGO4HAyggAAk5IRXOP0mYvZ2vrg6GFwNzz1wY8cPQQAWAQZQQCASckIAgBjmCM4nIwgAMCkBIIAAJMSCAIATMocwTXOnT1p9SnAPmxd0GlhyU6c2vAFO8llcwRHkxEEAJiUQBAAYFJKw8CRorwIsDkyggAAk5IRBAAG6KQvjx7E9ASCa9hreLms9l42zw9GOT96AAygNAwAMCmBIADApJSGAYAxWkPp0WQEAQAmJRAEAJiU0jDHjobEy2bVMMDmCAQBgM3rJJfNERxNaRgAYFIygmucO3tSiQpgH0zLWLYTp0aPgBEOLSNYVQ+pqndU1S1VdXNVvXB1/HOr6m1VdX717wNWx59cVe+rqg+t/n3iju/63qr6aFV96irXfMzq87dW1auqqlbHv2V1/KaqeldVPeyw7hsAYCkOszR8R5Lv6O6/mOTxSb51FYC9JMnbu/u6JG9fvU6Sjyd5end/SZLnJvmJHd/15iSP28M1fzTJ85Jct/p52ur4P+/uL+nuRyX5/iQ/eE9uDAA4AN3H82dBDq003N23J7l99fsnq+qWJNckeUaSJ6xOe12SdyZ5cXd/YMfHb05y76r67O7+dHe/J0lWCb5dVdWpJPfr7nevXr8+yVcneWt3f2LHqX8q21NU75K9hpdLSR/G8Le3dPYantFG5ghW1UOTPDrJryT5/FWQmO6+vaoetMtHvibJB7r70/u4zDVJbtvx+rbVsTvH8K1J/l9JPivJEwMAMLlDXzVcVfdJ8oYk33ZFZm7d+Q9P8n1J/tZ+L7XLsT/O/HX3j3T3X0jy4iT/7zXXfl5V3VhVN/6H3720z8sDACzLoQaCVXWvbAeBP9ndb1wd/p1VGffOcu7Hdpx/bZI3JXlOd//GVb77xGrxx01V9T3ZzgBeu+OUa5Nc2OWjP53tkvFn6O7ru/ux3f3Yz/szJ/Z0jwAAS3VopeHVit0fS3JLd+9cnHFDtheDvHL178+vzr9/kl9I8tLu/uWrfX93X0ryqCuu+cmqeny2S9DPSfJDq+PXdfedkx/+WkyEAIDxFraw4jg6zIzglyd5dpIn7sjcfWW2A8AnV9X5JE9evU6SFyT5oiQv23H+g5Kkqr6/qm5LcrKqbquql6+55vOTvCbJrUl+I8lb7/zuVQubm7I9T/C5B32zAABLc5irht+V3eftJcmTdjn/FUlesea7XpTkRXu45o1JHrHL8Rde7bMcH5raLpuVpwCbY4s5AIBJ2WIOABhgec2XjyMZQQCASQkEAQAmJRAEAJiUOYJrnDt70upFgH2wYn/ZTpza8AU7yeXLG74oV5IRBACYlEAQAGBSSsPAkaK8CLA5AkEAYAx9BIdTGgYAmJSMIHCkWK0Po5wfPQAGkBEEAJiUQBAAYFJKwwAcCCu+l23jDaUTi0WOABlBAIBJCQQBACalNAwcKcqLAJsjEAQABujksjmCoykNAwBMSkZwjdNnLmZrS4kKNk1DaRhFQ+kZyQgCAExKRhAA2LxOui+PHsX0BIIcO0qLALA3SsMAAJMSCAIATEppGAAYQx/B4WQEAQAmJRAEAJiU0jDHjr1ql82qb4DNkREEAJiUjCAAMEZbLDKaQHCNc2dPKlEB7INpGct24tToETCC0jAAwKQEggAAk1IaBo4U5UWYRHdy+fLoUUxPRhAAYFICQQCASQkEAQAmZY4gADCGPoLDyQgCAExKRnCN02cuZmvL6sUl0gh82Tw/GOX86AEwgIwgAMCkZAQBgCFaH8HhBIJr2GsYADjulIYBACYlEAQAmJRAEABgUuYIAgADtIbSR4CMIADApASCAACTEggCAEzKHME1bDEHY+jfCZPoJJfNERxNRhAAYFICQQCASQkEAQAmZY4gADBGXx49gunJCAIATEpGcI1zZ09avQgDbF2wWh9GOHFq9AgYQUYQAGBSMoIAwMZ1ktZHcDiB4BoaSi+Xkv6yeX4wyvnRA2AApWEAgEkJBAEAJiUQBACYlDmCAMDmdWsofQTICAIATEpGcA0NpQH2RzPwZdNQek4yggAAk5IRBACG0FB6PBlBAIBJCQQBACYlEAQAmJQ5ggAcCJ0Wlm7AXsOT9hGsqqcl+adJTiR5TXe/8or3a/X+Vya5mOQbuvv9hzEWGUEAgA2pqhNJfiTJVyR5WJKvr6qHXXHaVyS5bvXzvCQ/eljjEQgCAGzO45Lc2t0f7u4/SvLTSZ5xxTnPSPL63vaeJPevqkPp9CgQBADYnGuSfHTH69tWx/Z7zoEwR3CNT+b3P/Wv+2d/ffQ4DtEDk3x89CAO0XG+v+N8b4n7Wzr3t1xfvMmLfTK/v/Wv+2cfuMlrbtC9q+rGHa+v7+7rV7/XLudf2VBxL+ccCIHger/e3Y8dPYjDUlU3ur9lOs73lri/pXN/y3VF4HLouvtpm7zeEXJbkofseH1tkgt345wDoTQMALA5701yXVV9YVV9VpJnJrnhinNuSPKc2vb4JH/Q3bcfxmBkBAEANqS776iqFyTZynb7mB/v7pur6ltW7786yVuy3Trm1my3j/nGwxqPQHC9669+yqK5v+U6zveWuL+lc3/LdZzv7Ujp7rdkO9jbeezVO37vJN+6ibHU9rUAAJiNOYIAAJOaPhCsqqdV1a9X1a1V9ZJd3n9CVf1BVd20+vl7I8Z5d1TVj1fVx6rqV9e8X1X1qtW9n62qL930GO+JPdzfkp/dQ6rqHVV1S1XdXFUv3OWcxT6/Pd7fkp/fvavq31bVB1f39/d3OWfJz28v97fY55ds7/5QVR+oqn+5y3uLfXZ3usr9LfrZsT9TzxHcsc3Lk7O9VPu9VXVDd//aFaf+Unf/9Y0P8J57bZIfTvL6Ne/v3MLmy7K9hc2XbWRkB+O1uev7S5b77O5I8h3d/f6qum+S91XV2674z+aSn99e7i9Z7vP7dJIndvenqupeSd5VVW9d7RBwpyU/v73cX7Lc55ckL0xyS5L77fLekp/dne7q/pJlPzv2YfaM4F62eVms7v7FJL93F6dsbAubw7CH+1us7r79zg3Gu/uT2f4v7Cu7yi/2+e3x/hZr9Uw+tXp5r9XPlROyl/z89nJ/i1VV1yb5a0les+aUxT67ZE/3x0RmDwT3uoXL/2VVAnlrVT18M0PbiI1tYTPQ4p9dVT00yaOT/MoVbx2L53cX95cs+PmtSm83JflYkrd197F6fnu4v2S5z++fJHlRkstr3l/0s8vV7y9Z7rNjn2YPBPeyhcv7k/y57n5kkh9K8nOHPagN2tgWNoMs/tlV1X2SvCHJt3X3J658e5ePLOr5XeX+Fv38uvtSdz8q2zsCPK6qHnHFKYt+fnu4v0U+v6r660k+1t3vu6vTdjm2iGe3x/tb5LPj7pk9ELzqFi7d/Yk7SyCrvj/3qqrjsjfixrawGWHpz2419+oNSX6yu9+4yymLfn5Xu7+lP787dfd/TPLOJFdup7Xo53endfe34Of35Um+qqp+M9vThZ5YVf/bFecs+dld9f4W/Oy4G2YPBK+6zUtV/dmqqtXvj8v2/81+d+MjPRwb28JmhCU/u9W4fyzJLd39g2tOW+zz28v9Lfz5fV5V3X/1++ck+R+S/LsrTlvy87vq/S31+XX3S7v72u5+aLb/N+HfdPffvOK0xT67vdzfUp8dd8/Uq4b3uM3L1yZ5flXdkeQ/J3lmL6QLd1X9VJInJHlgVd2W5LuzPal741vYHIY93N9in122/7/2Zyf50GoeVpJ8Z5IvSI7F89vL/S35+Z1K8rpVZ4I/keRnuvtf1qAtpA7BXu5vyc/vMxyjZ7er4/zsuGt2FgEAmNTspWEAgGkJBAEAJiUQBACYlEAQAGBSAkEAgEkJBIEjraruX1X/z9HjADiOBILAkbXqU3f/JPsKBFeNfv33G8BV+C9K4MBU1XdV1a9X1b+uqp+qqr9bVe+sqseu3n/gamurVNVDq+qXqur9q5//6+r4E6rqHVX1z5N8KMkrk/yFqrqpqv7R6pz/pareW1Vnq+rv7/i+W6rqn2V7r9SHfOYIAdhp6p1FgINTVY/J9pZVj872f7e8P8ldbWz/sSRP7u4/rKrrkvxUkseu3ntckkd097+vqoeufn/U6jpPSXLd6pxKckNV/eUkH0nyxUm+sbuVkgH2QCAIHJT/W5I3dffFJKmqG65y/r2S/HBVPSrJpSSnd7z3b7v736/53FNWPx9Yvb5PtgPDjyT5re5+z90bPsB8BILAQdptz8o78l+nodx7x/FvT/I7SR65ev8Pd7z3n+7iGpXkH3b3/+e/ObidObyrzwFwBXMEgYPyi0n+RlV9TlXdN8nTV8d/M8ljVr9/7Y7z/3SS27v7cpJnJzmx5ns/meS+O15vJfmfq+o+SVJV11TVgw7mFgDmIiMIHIjufn9V/YskNyX5rSS/tHrrB5L8TFU9O8m/2fGRf5bkDVX1dUnekTXZvO7+3ar65ar61SRv7e7/par+YpJ3V1WSfCrJ38x2eRmAfaju3So5APdMVb08yae6+wdGjwWA3SkNAwBMSkYQAGBSMoIAAJMSCAIATEogCAAwKYEgAMCkBIIAAJMSCAIATOr/D13J+5Yam6JlAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x1440 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "da = df.groupby(['utdate','quarter'])['kpfcc'].first().to_xarray()\n",
    "da.coords['utdate'] = Time(list(array(da.utdate))).datetime\n",
    "da.sortby(da.utdate,ascending=False).plot.pcolormesh(figsize=(10,20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f5415c65",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2459249.106688459 1\n",
      "2459250.059417007 1\n",
      "2459253.09286689 2\n",
      "2459260.125596199 2\n",
      "2459261.0190041875 2\n",
      "2459262.111134746 1\n",
      "2459266.950185211 2\n",
      "2459268.1099689654 1\n",
      "2459270.084943889 2\n",
      "2459271.0654437374 1\n",
      "2459276.1148567637 1\n",
      "2459280.026814646 2\n",
      "2459281.1220216947 1\n",
      "2459289.059139398 1\n",
      "2459291.137326309 1\n",
      "2459293.0679493668 1\n",
      "2459295.0775366384 1\n",
      "2459297.070549879 1\n",
      "2459299.0796948383 1\n",
      "2459302.0698245405 1\n",
      "2459304.1239861357 1\n",
      "2459307.077082991 2\n",
      "2459308.127139249 1\n",
      "2459313.1389600346 1\n",
      "2459313.8816464585 2\n",
      "2459320.8061936377 1\n",
      "2459321.8260970474 2\n",
      "2459324.8235150143 2\n",
      "2459326.1029022182 1\n",
      "2459326.7806885336 2\n",
      "2459328.063117215 1\n",
      "2459328.987803917 2\n",
      "2459332.06854561 1\n",
      "2459333.0030463687 2\n",
      "2459334.769926875 1\n",
      "2459336.964664036 2\n",
      "2459337.9560520076 2\n",
      "2459339.084497761 1\n",
      "2459340.7990041072 1\n",
      "2459347.7741028676 1\n",
      "2459348.7849066523 1\n",
      "2459350.0649792324 1\n",
      "2459350.771438949 1\n",
      "2459351.7908531507 1\n",
      "2459352.9619807624 2\n",
      "2459354.7755374103 1\n",
      "2459357.1155770943 1\n",
      "2459358.094146289 2\n",
      "2459359.086323545 1\n",
      "2459360.0443681367 1\n",
      "2459360.780567459 1\n",
      "2459362.0846117944 1\n",
      "2459363.0861548753 1\n",
      "2459363.7582882107 1\n",
      "2459364.825284251 1\n",
      "2459366.0777728837 2\n",
      "2459367.0116591067 2\n",
      "2459367.8135418436 1\n",
      "2459369.098510151 1\n",
      "2459369.8718242906 2\n"
     ]
    }
   ],
   "source": [
    "timestamps = []\n",
    "dates_observed = []\n",
    "f = open('timestamps_C.csv', 'w')\n",
    "date_groups = df[df['kpfcc']].groupby('utdate', as_index=False)\n",
    "for grp in date_groups:\n",
    "#     if np.random.uniform() > 0.50:  # 33% chance of observing on a given CC night\n",
    "#         continue\n",
    "    num = len(grp[1])\n",
    "    block = grp[1].sample(1)\n",
    "    date = block['utdate'].values[0]\n",
    "    if date not in dates_observed:\n",
    "        obs = Time(np.random.uniform(block.start_jd, block.stop_jd), format='jd')[0]\n",
    "        if keck.target_is_up(obs, target):\n",
    "            print(obs, num)\n",
    "            timestamps.append(obs)\n",
    "            dates_observed.append(date)\n",
    "            print(obs, file=f)\n",
    "        if len(timestamps) >= max_obs:\n",
    "            break\n",
    "\n",
    "\n",
    "f.close()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "423557b9",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "timestamps = []\n",
    "f = open('timestamps_A.csv', 'w')\n",
    "\n",
    "date_groups = df[df['kpfcc']].query('3 <= month <= 10')\\\n",
    "        .sort_values(by='utdate')\\\n",
    "            .groupby('utdate', as_index=False)\n",
    "random_dates = date_groups.sum().sample(10)['utdate'].values\n",
    "\n",
    "i = 0\n",
    "while len(timestamps) < max_obs:\n",
    "    for grp in date_groups:\n",
    "        ind = grp[0]\n",
    "        if ind not in random_dates:\n",
    "            continue\n",
    "        num = len(grp[1])\n",
    "        block = grp[1].sample(1)\n",
    "        date = block['utdate'].values[0]\n",
    "        start = \"{:s} {:s}\".format(date, block.start.values[0])\n",
    "        end = \"{:s} {:s}\".format(date, block.stop.values[0])\n",
    "        range = Time([start, end])\n",
    "        can_obs = is_observable(constraints, keck, target, range)[0]\n",
    "        tab = observability_table(constraints, keck, [target], range, time_grid_resolution=60*u.second)\n",
    "        obs = Time(np.random.uniform(block.start_jd, block.stop_jd), format='jd')[0] + 365*i\n",
    "        timestamps.append(obs)\n",
    "        timestamps.append(obs+600/86400)\n",
    "        timestamps.append(obs+1200/86400)\n",
    "        print(obs, file=f)\n",
    "        print(obs+600/86400, file=f)\n",
    "        print(obs+1200/86400, file=f)\n",
    "        if len(timestamps) >= max_obs:\n",
    "            break\n",
    "\n",
    "        \n",
    "#         frac = tab['fraction of time observable']\n",
    "#         hours_observable = ((range[1] - range[0]) * frac) * 24\n",
    "#         hours_observable = hours_observable[0]\n",
    "#         print(start, end, hours_observable)\n",
    "#         if can_obs and hours_observable >= 1:\n",
    "#             # print(obs, num, can_obs, hours_observable)\n",
    "#             timestamps.append(obs)\n",
    "#             print(obs, file=f)\n",
    "#         if len(timestamps) >= max_obs/3:\n",
    "#             break\n",
    "    i+=1\n",
    "\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "41725250",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-13-d911b7d5a3c9>:21: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
      "  ax.xaxis.set_ticklabels([x-2450000 for x in ax.get_xticks()])\n",
      "<ipython-input-13-d911b7d5a3c9>:21: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
      "  ax.xaxis.set_ticklabels([x-2450000 for x in ax.get_xticks()])\n",
      "<ipython-input-13-d911b7d5a3c9>:21: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
      "  ax.xaxis.set_ticklabels([x-2450000 for x in ax.get_xticks()])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcjUlEQVR4nO3df3xU9Z3v8fcXpUS7ARTdakEJRb0CCUmWqPy4hUSuLWKlkKsRpV1K10cfIm0FHkWhXr0U3cqlXrDqxXZVULeWIO6qVVi3agmtugqhDT8UFShhiSh1sQSoRfnxuX/MmTCZzCSZ/Jj5nsnr+XjkkZlzzvec7yczk3fOj5yvMzMBAOCbbpnuAAAAiRBQAAAvEVAAAC8RUAAALxFQAAAvnZrKwmeddZbl5eV1Ulc611/+8hd9/vOfz3Q3Ok0215fNtUnUF3bU134bN278LzM7O356SgGVl5en6urqjutVGlVVVam0tDTT3eg02VxfNtcmUV/YUV/7Oed2J5rOIT4AgJcIKACAlwgoAICXCCgAgJcIKACAlwgoAICXCKgkRi18pcPbN7fO+Hmxz9vbl/j1duT6km0jHdK1HbQsnZ+XlpYN0+clne/hMH5eUvo/qK7k/QNHOrx9c+uMnxf7vL19aW0fwrSNdG4HLUvn56WlZcP0eUnneziMnxf2oAAAXnKpDFhYUlJi2XwniVELX0n4V0bf3jl6be7YFreRrH0ifXvnSEr9r5pkfWmuvub61draWtLen11zYmvrzO1kSljvRNDa1yJZfal8Xro56UQbxlZty/si1c9Lqq9fOt/DHbGtNN1JYqOZlTSZYWat/ho2bJiF1dq1a1Navv9tL7Rre4naN7fO+Hmxz1vTl9bW1/+2F9pdW2u20ZGS1dbZdaRLqu9NHzX3WrSmvlQ+Ly0t25Hvi9Z8Xtrz+qXzPdzWbaXj/Smp2hJkDof4AABe4iKJJKKH4DqyfXPrjJ8X+7y9fWltH8K0jXRuBy1L5+elpWXD9HlJ53s4jJ8XzkFliWyuL5trk6gv7Kiv/ZKdg+IQHwDASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEsEFADASwQUAMBLBBQAwEvhCKi19yR+3JFtl49vOj/6vKXpy8c3Xkfs80TrjW1zz3lN1xO7neg6YuetvafpV+z2E82TpCX5TfsT218A8Eg4AmrdwsSPO7Lt7teazo8+b2n67tcaryP2eaL1xrb59GDT9cRuJ7qO2HnrFjb9it1+onmSVL+naX9i+wsAHglHQAEAuhxnZq1euKSkxKqrqzuxOzHW3tPy3tKYuVLZvLa3rf0dew/x+o+Spq3JdC8aqaqqUmlpaaa70WmoL9yor/2ccxvNrCR++qmdutX2KJt3Mnzm95Lm1zd9nELbqqoqlVZ9Pa5tTLjN7xV8r288bX59023GT4/93rBMguex4tsk206jNvWJp0uqKn0uUl8iidrFT2vpZwoAacYhPgCAl/zdg4o1Zm7ixx3Ztv8oKe/LidvGt4uf3n9U43XU/q7x9Pj1xs7r0VMafnPjabHbia4r0fZbU0OsXudJvc9vvEyybQBAhvl7DqqDcZw4vLK5Non6wo762i/ZOSgO8QEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALxEQAEAvERAAQC8REABALx0aqY7AKDrOXr0qOrq6nTkyJFMd6XdevXqpW3btmW6G52mI+vLyclRv3791L1791YtT0ABSLu6ujrl5uYqLy9PzrlMd6ddDh06pNzc3Ex3o9N0VH1mpv3796uurk4DBgxoVRsO8QFIuyNHjqhPnz6hDye0nnNOffr0SWmvmYACkBGEU9eT6mtOQAEAvERAAQC8REABCI0lL73XIes5cOCAli5dKknau3evrrnmmg5ZbyI1NTVas2ZNp62/rZYsWaKcnBzV19dnuitJEVAAQuOnr2zvkPXEBtQXv/hFPf300x2y3kR8DagVK1bokksu0TPPPJPpriRFQAHocubOnaudO3eqqKhI1157rfLz8yVJjz32mCZOnKirr75aAwYM0IMPPqjFixeruLhYw4cP18cffyxJ2rlzp8aNG6dhw4bpq1/9qt555x1J0qpVq5Sfn6/CwkKNHj1an332me68806tXLlSRUVFWrlypdavX6+RI0equLhYI0eO1LvvvpvStktLSzVz5kyNHDlS+fn5Wr9+vSRp3bp1KioqUlFRkYqLi3Xo0KGk9e/cuVOHDx/W3XffrRUrViRdbtGiRRo+fLgKCws1d+5cSdL999+vwYMHa+jQoZo8ebJOnDihvLw8HThwoKHdBRdcoH379rXx1YlhZq3+GjZsmIXV2rVrM92FTpXN9WVzbWZds76333671e0X//pd63/bC02+Fv/63Tb3adeuXTZkyJAmj5cvX24DBw60gwcP2p/+9Cfr2bOnPfTQQ2ZmNnPmTFuyZImZmV1++eX23nvvmZnZK6+8YmVlZWZmlp+fb3V1dWZm9uc//7lhnTNmzGjYdn19vR09etTMzF566SUrLy9PadtjxoyxG2+80czM1q1b19D3r33ta/bqq6+amdmhQ4catpHIXXfdZQsWLLDjx49b//79bd++fU2WWbNmjY0YMcI+/PBDMzPbv3+/mZmde+65duTIkUY1fv/737dly5aZmdkbb7xhY8eOTbrtRK+9pGpLkDn8oy4Ar8264iLNuuIiSVLe3NWqXXhVp26vrKxMubm5ys3NVa9evXT11VdLkgoKCrR582YdPnxYr7/+uq699lpJ0okTJ3T06FFJ0qhRo/Stb31LFRUVKi8vT7j++vp6TZ06Vdu3b5dzrqFta7Yddf3110uSRo8erYMHD+rAgQMaNWqUZs+erSlTpqi8vFz9+vVLWmNlZaWeeeYZdevWTeXl5Vq1apVmzJjRaJmXX35Z06ZN0+mnny5JOvPMMyVJQ4cO1ZQpUzRx4kRNnDhRknTddddpwYIFmjZtmiorK3Xddde17ofdAg7xAUCMHj16NDzu1q1bw/Nu3brp2LFjOnHihHr37q2amhrV1NTotddea7gV0M9+9jPdfffd2rNnj4qKirR///4m67/jjjtUVlamrVu36vnnn2/0j6stbTsq/v+JnHOaO3euHnnkEf31r3/V8OHDGw47xtu8ebO2b9+uK664Qnl5eaqsrEx4mM/MEv7f0urVqzVjxgxt3LhRw4YN07FjxzRixAjt2LFDH330kZ599tmk4ZwqAgpAaNwy9sIOWU9ubm6z52ia07NnTw0YMECrVq2SFPlFvmnTJkmRczuXXXaZFixYoLPOOkt79uxpsq36+nr17dtXUuS8U1usXLlSkvTqq6+qV69e6tWrl3bu3KmCggLddtttKikpSRpQK1as0Pz581VbW6va2lrt3btX77//vnbv3t1oua985StatmyZPvnkE0nSxx9/rBMnTmjPnj0qKyvTokWLdODAAR0+fFjOOU2aNEmzZ8/WoEGD1KdPnzbVFY+AAhAa0UN97dWnTx+NGjVK+fn5mjNnTsrtn3zyST366KMqLCzUpZdequeee06SNGfOHBUUFCg/P1+jR49WYWGhysrK9PbbbzdcJHHrrbdq3rx5GjVqlI4fP96m/p9xxhkaOXKkbrrpJj366KOSpPvuu6/hAo3TTjtNV155ZcK2lZWVmjRpUqNpkyZNUmVlpaqrq3XjjTdKksaNG6cJEyZozJgxKioq0r333qvjx4/rG9/4hgoKClRcXKxZs2apd+/ekiKH+X7xi1902OE9SVwkkS2yub5srs2sa9aXykUSvjt48GBatzdmzBjbsGFD2rbX0fWlcpEEe1AAAC9xFR8AhEhVVVWrltuyZYu++c1vNprWo0cPvfnmm53Qq85BQAFAFiooKFBNTU2mu9EuHOIDAHiJgAIAeImAAgB4iYACEB5r78l0D5BGBBSA8Fi3sENW05XHg6qtrdVpp52moqIiFRYWNrqjum8IKABdTlcfD2rgwIGqqanRpk2bNHXqVP34xz/OdJcSIqAA+G3tPdL8XpEv6eTjdhzu6+rjQcU6ePCgzjjjjITznnjiCY0YMUKFhYUN/1MVX6MkXXbZZXrrrbca2pWWlmrjxo0pvSYJJbq9RLIvbnXkr2yuL5trM+ua9bX5Vkf/u2f7OhPoyuNB7dq1y3JycqywsNC+9KUv2TnnnGO7d+9ustzWrVvtoosusl27dpnZyfGgEtW4ePFiu/POO83MbO/evXbhhRcm+9FzqyMAaKvomExnn312kzGZamtrG40HVVRUpJkzZ+qDDz6QdHI8qIcffjjpjWDr6+sb9tpmzZrVaM+jpW1HNTce1P33368DBw7o1FOT34cheohv586duu+++/Sd73ynyTK/+c1vdM011zTcmTw6HlSiGisqKhru7v7UU081jJXVXgQUgPAYM7fTN5Ht40HFmzBhgn772982mW5JxoNKVGPfvn3Vp08fbd68WStXrtTkyZNbte2WEFAAwqNsXoespiuPBxXv1Vdf1cCBA5tMHzt2rJ566qmGkI09/xZfoyRNnjxZixYtUn19vQoKCtpUVzwCCkCX05XHg5LUcIFIYWGhfvjDH+qRRx6RFLnkfvz48ZKkIUOG6Pbbb9f48eNVWFio2bNnJ61Rkq655hpVVlaqoqKiTTUllOjEVLIvLpLwVzbXl821mXXN+hgPqu0YDwoAgAxjuA0ACBHGgwIAhBrjQQEA0EkIKACAlwgoAICXCCgAobG0Zmmmu4A0IqAAhMZDmx7qkPV05fGgJOm9997T+PHjdcEFF2jQoEGqqKjQvn37Mt2tJggoAF1OVx4P6siRI7rqqqs0ffp07dixQ9u2bdP06dP10UcfZbprTRBQALy2tGapCh4vUMHjkfu7RR+353BfVx4P6pe//KVGjBjRcKd0KXIX9ejPIBbjQaVJV7ydTLbI5trMumZ9bb3VUf5j+e3sTURXHg9q1qxZdt9997X4M/JhPKiM/KNu9C+fm4tu1tKapbq56OaGedNenKZLzrmkYV6s6HLRNtNenKbl45YnXHe8wRrcaLuxbWOnxz+P9mfDhxu0fNzyhufx4muIitYSnXbJOZfo2R3PauIFE7Xhww2NlgGQedExmXJzc5uMybR58+ZG40FJ0okTJ3T06FFJJ8dKqqioUHl5ecL119fXa+rUqdq+fbuccw1tW7PtqObGg5oyZYrKy8vVr1+/dv0cWhoPKrbGiooKXXHFFfrRj37UoeNBZSSgoic6by66WQ9teqjRL+fqfdWq3lfdMC9WdLlom+p91UnXHe+B/g802m5s29jp8c+j/YnvX7z4GmIfx24v+j22n9FlADRveuH0Tt9GKuNBSdKhQ4eUm5srKTJW0ptvvqnVq1erqKgo4Z0couNBPfPMM6qtrVVpaWmrtx2VbDyoq666SmvWrNHw4cP18ssv6+KLL26y/SFDhmjdunUt/hysmfGg4muMHw/q5z//eYvrbw3OQQEIjY76Q64rjwd1ww036PXXX9fq1asbpr344ovasmVLo+V8GA8qbXtQS2uWNtm7iT3pGa+lac21TeR7u7/XYttUn3fEvPhlphdOZ28K6GSx40ENGjQo5fZPPvmkpk+frrvvvluffvqpbrjhBhUWFmrOnDnavn27zExjx45VYWGhzj//fC1cuFBFRUWaN2+ebr31Vk2dOlWLFy/W5Zdf3qb+R8eDOnjwoJYtWyYpMh7U2rVrdcopp2jw4MFJx4M67bTT9MILL2jmzJmaOXOmunfvrqFDh+qnP/2p9u7dqxtvvFFr1qxpNB5U9+7dVVxcrMceeyxhjVJkPKhbbrlFd9xxR5tqSsRFzk+1TklJiVVXNz28laroL+MtU7eo4PECbZm6Jem8WNHlom3i28a2j/dA/wcaQiq+bew2m+tDoj7F9y1RH5prl6h9W1RVVTU6VJBNsrk2qWvWt23btjYFg49iD/GlQ2lpqe69916VlJSkZXsdXV+i1945t9HMmhTEIT4AgJcycpFE7InO+JOeJV8oabhKLtkJ0ej0ki80/Qsi6UnUA43nxbaNbxO/XPQqvvj+JRO77uiy0WnJruIDgNboSuNBZeQQXyZ0xcMo2SKba5O6Zn3btm3TxRdfnPAqsbBJ9yG+dOvI+sxM77zzDof4APgrJydH+/fvVyp/ICPczEz79+9XTk5Oq9swoi6AtOvXr5/q6uq8vP9bqo4cOZLSL92w6cj6cnJyUvoHYgIKQNp1795dAwYMyHQ3OkRVVZWKi4sz3Y1Ok8n6OMQHAPASAQUA8BIBBQDwUkqXmTvnPpK0u/O606nOkvRfme5EJ8rm+rK5Non6wo762q+/mZ0dPzGlgAoz51x1ouvss0U215fNtUnUF3bU13k4xAcA8BIBBQDwUlcKqH/KdAc6WTbXl821SdQXdtTXSbrMOSgAQLh0pT0oAECIEFAAAC+FJqCcc7c457Y6595yzs0Mpt3lnNvsnKtxzv3aOffFmOXnOed2OOfedc59NWb6MOfclmDe/S7J/f6TtU9nfTHzfuCcM+fcWS31L0z1OefmO+feD16/Gufc+DDWl+y1c859L9j+W865RWGsLdheotduZczrVuucq8my+oqcc28E9VU75y7NsvoKnXP/EfT3eedcTy/rMzPvvyTlS9oq6XRFbnD7sqQLJfWMWeb7kn4WPB4saZOkHpIGSNop6ZRg3npJIyQ5Sf8m6coE20vaPp31BfPOk/TvivyD9FnZVJ+k+ZJ+kEr/fKuvmdrKgsc9guX+Nmy1tfTejFnm/0q6M5vqk/TraP8kjZdUlWX1bZA0Jljm25Lu8rG+sOxBDZL0hpl9YmbHJK2TNMnMDsYs83lJ0Ss+vi6p0sw+NbNdknZIutQ5d64iofYfFvlpPiFpYoLtJWzfKZVFJKwvmLdE0q06WVvS/oW0vkTCVF+y2qZLWmhmn0qSmf0phLU1V58kKfgrukLSiub6F8L6TFJ0r6KXpL3N9S+E9f03Sb8NlnlJ0v9srn+Zqi8sAbVV0mjnXB/n3OmK/EVzniQ55/7RObdH0hRJdwbL95W0J6Z9XTCtb/A4fnq8ZO07S8L6nHMTJL1vZpta2b9Q1RfM+66LHKZd5pw7o4X++VhfstoukvRl59ybzrl1zrlLWuibj7VJzb92kvRlSfvMbHsL/QtbfTMl/ST43XKvpHkt9C9s9W2VNCFY5lqdfE29qi8UAWVm2yT9H0WS/kVFdiGPBfNuN7PzJD0p6btBk0THRq2Z6fFau1yHaKa+23UydFvTv7DV95CkgZKKJH2gyKGi5vrnXX3N1HaqpDMkDZc0R9JTwd5GaGqTmv/sBa7Xyb2n5voXtvqmS5oV/G6ZJenRFvoXtvq+LWmGc26jpFxJn7XQv4zUF4qAkiQze9TM/s7MRkv6WNL2uEV+qZO7qXVq/FdeP0V20euCx/HT4yVr32kS1FeryDHcTc652qAPv3fOndNM/8JU33Yz22dmx83shKSHdfJQQKjqS/LerJP0rxaxXtIJRW66GarapOSfPefcqZLKJa1sRf/CVt9USf8aLLJKIX1vSkk/e++Y2VfMbJgif2DsbKF/mamvrSev0v2lkyeZz5f0jiJ/nV4YM/97kp4OHg9R4xN1f9TJE30bFPmrNnqib3yCbSVtn8764ubX6uRFEllRn6RzY+bPUuTYdejqS1LbTZIWBNMvUuSwhwtbbc29NyWNk7Sutf0LU32StkkqDaaPlbQxy+qLTuumyPmkb/tYX6f9UDrhh/w7SW8HxY8Npv2LIsdSN0t6XlLfmOVvV+SvgncVc7WJpJKgzU5JD+rk3TQmKPiF0lz7dNYXN79WQUBlS32S/lnSluD1+5UaB1Zo6ktS2+ck/SLo6+8lXR7G2pp7b0p6TNJNCZYPfX2S/rukjcG0NyUNy7L6bpH0XvC1MNpX3+rjVkcAAC+F5hwUAKBrIaAAAF4ioAAAXiKgAABeIqAAAF4ioJC1nHOHnXN5zrm/Ouf+4Jzb5pxb75yb2oZ1zXbOvR3clukV51z/uPk9XeTO7A/GTHvMObcr5q7fRcF0F9wNekewvr+LaTMuuAv0Dufc3JjpZzrnXnLObQ++nxEzL613xwbShYBCV7DTzIrNbJCkyZJmOeempbiOP0gqMbOhkp6WtChu/l2K3Igz3hwzKwq+aoJpVypyR+kLJX1HkVs+yTl3iqT/F8wfLOl659zgoM1cSa+Y2YWSXgmeK5g/WZF/kBwnaWmwHiD0CCh0KWb2R0mzFRmeJZV2a83sk+DpG4q57YtzbpikLygyRENrfF3SExbxhqTewd2iL5W0w8z+aGafSaoMlo22eTx4/LhO3kk63XfHBtKGgEJX9HtJF7ej/T8ocqsXOee6KXKT2zlJlv3H4DDeEudcj2Bac3eMTnYn6C+Y2QeSFHz/2xbWBYQeAYWuKNEdl1vX0LlvKHLLl58Ek26WtMbM9iRYfJ4iQXiJpDMl3dbM9lO5Y3SjLrWhDRAKp2a6A0AGFCtyM9BGnHPLg3l7zWx8gvn/Q5H7jI2xYCBCRUYY/bJz7mZJfyPpc865w2Y2N7rHI+nTYN0/CJ4nu+Pz55JMl6R9zrlzzeyD4HBgdADEtN8dG0gX9qDQpTjn8hQZgO6B+HlmNi24mCFROBVL+rmkCXZydFyZ2RQzO9/M8hQJoCfMLHoBw7nBd6fIOaOtQbNfSfr74Gq+4ZLqgzDbIOlC59wA59znFLn44VcxbaJXH06V9FzM9MnOuR7OuQGKXHixPvWfDOAf9qCQlYKxiqJ7OQOdc3+QlCPpkKQHzGx5iqv8iSJ7SKsieaP/NLMJzTfRk865sxU5DFejyBAckrRGkZFNd0j6RNI0STKzY86570r6d0mnSFpmZm8FbRYqMujhP0j6T0VGQZWZveWce0qRu1UfkzTDzI6nWBvgJe5mjqzknCuU9LCZcUUbEFIc4kPWcc7dpMgoof8r030B0HbsQQEAvMQeFADASwQUAMBLBBQAwEsEFADASwQUAMBL/x+nG/4jz/lGrQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def jd2date(t):\n",
    "    return Time(t, format='jd').datetime\n",
    "\n",
    "def date2jd(t):\n",
    "    try:\n",
    "        time = Time(t, format='datetime').jd\n",
    "    except:\n",
    "        time = Time(2455565, format='jd')\n",
    "        pass\n",
    "    return time\n",
    "\n",
    "i = 0\n",
    "for fname in ['timestamps_A.csv', 'timestamps_B.csv', 'timestamps_C.csv']:\n",
    "# for fname in ['timestamps_A.csv']:\n",
    "    times = pd.read_csv(fname, header=None)\n",
    "    dts = [jd2date(t[0]) for t in times.values]\n",
    "    times = [date2jd(t) for t in dts]\n",
    "    pl.plot(times, np.zeros_like(times)+i, '+', label=fname)\n",
    "    ax = pl.gca()\n",
    "    ax.set_xlabel('JD - 2450000')\n",
    "    ax.xaxis.set_ticklabels([x-2450000 for x in ax.get_xticks()])\n",
    "    ax.yaxis.set_ticks([])\n",
    "    ax.set_xlim(Time('2021-01-01', format='iso').jd, Time('2022-12-31', format='iso').jd)\n",
    "    i -= 1\n",
    "\n",
    "\n",
    "# axyrs = ax.twiny()\n",
    "# xl = np.array(list(ax.get_xlim()))\n",
    "# pl.locator_params(axis='x', nbins=24)\n",
    "# start = 1\n",
    "# end = 24\n",
    "# axyrs.set_xlim(start, end)\n",
    "# months = [calendar.month_abbr[int(x)] for x in axyrs.get_xticks()] * 2\n",
    "# axyrs.xaxis.set_ticklabels(months)\n",
    "# axyrs.set_xlabel('Month', fontweight='bold')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.grid(True)\n",
    "legend = ax.legend()\n",
    "pl.savefig('observing_schedule.pdf')\n",
    "pl.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cc243762",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ True]\n",
      "[False]\n",
      "target name ever observable always observable fraction of time observable\n",
      "----------- --------------- ----------------- ---------------------------\n",
      "   HD155060            True             False                         0.5\n"
     ]
    }
   ],
   "source": [
    "start = \"{:s} {:s}\".format(block.utdate.values[0], block.start.values[0])\n",
    "end = \"{:s} {:s}\".format(block.utdate.values[0], block.stop.values[0])\n",
    "time_range = Time([start, end])\n",
    "print(is_observable(constraints, keck, target, time_range))\n",
    "print(is_always_observable(constraints, keck, target, time_range))\n",
    "print(observability_table(constraints, keck, [target], time_range))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "738bf398",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-34-960176e4d6cc>:1: MatplotlibDeprecationWarning: Passing the minor parameter of get_ticklocs() positionally is deprecated since Matplotlib 3.3; the parameter will become keyword-only two minor releases later.\n",
      "  axyrs.get_xticks([])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.06, -0.05, -0.04, -0.03, -0.02, -0.01,  0.  ,  0.01,  0.02,\n",
       "        0.03,  0.04,  0.05,  0.06])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "axyrs.get_xticks([])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "c563245e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9, 2, 1, 9, 3, 5, 6, 0, 3, 5])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(10, size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8cf4610",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "836c0756",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "f03a80e608be7842e601eaf57f5b09fad8736b450ed14f8b12320f64d023db2c"
  },
  "kernelspec": {
   "display_name": "Python 3.8.5 64-bit ('base': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}