notebook_retrieve_field_star_radec_coords.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a4532f31",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "#sys.path.append('/Users/hhsieh/anaconda3/envs/astroconda/lib/python3.6/site-packages')\n",
    "#sys.path.append('/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages')\n",
    "#print(sys.version)\n",
    "import datetime,ssl\n",
    "import math\n",
    "import numpy as np\n",
    "import subprocess\n",
    "from os import path\n",
    "import astropy.units as u\n",
    "import astropy.coordinates as coord\n",
    "from astropy.time import Time, TimeDelta\n",
    "from astropy.table import Table\n",
    "from astropy.coordinates import SkyCoord,Angle\n",
    "from astropy.units import Quantity\n",
    "from astroquery.gaia import Gaia\n",
    "from astroquery.jplhorizons import Horizons\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.figure as fig\n",
    "from matplotlib.offsetbox import TextArea, DrawingArea, OffsetImage, AnnotationBbox\n",
    "import matplotlib.image as mpimg\n",
    "from matplotlib.ticker import NullFormatter, MaxNLocator, ScalarFormatter, MultipleLocator, FormatStrFormatter\n",
    "from urllib.parse import urlencode\n",
    "from urllib.request import urlretrieve\n",
    "import IPython.display\n",
    "from PIL import Image\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5be2f795",
   "metadata": {},
   "outputs": [],
   "source": [
    "def output_write_log_entry(log_filepath,log_entry):\n",
    "    print('{:s} - {:s}'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),log_entry))\n",
    "    with open(log_filepath,'a') as of:\n",
    "        of.write('{:s} - {:s}\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),log_entry))\n",
    "    return None\n",
    "\n",
    "def get_horizons_obj_desig(object_name):\n",
    "    horizons_obj_desig = object_name\n",
    "    if object_name == '66P':  horizons_obj_desig = '90000688;'\n",
    "    if object_name == '16P':  horizons_obj_desig = '90000273'\n",
    "    if object_name == '43P':  horizons_obj_desig = '90000507'\n",
    "    if object_name == '104P': horizons_obj_desig = '90000954'\n",
    "    if object_name == '213P': horizons_obj_desig = '90001150'\n",
    "    if object_name == 'C/2019 Q4': horizons_obj_desig = '90004424'\n",
    "    return(horizons_obj_desig)\n",
    "\n",
    "def remove_extra_characters(obj_name):\n",
    "    name_length = len(obj_name)\n",
    "    obj_name_ed = ''\n",
    "    for idx in range(0,name_length):\n",
    "        if obj_name[idx] != ' ' and obj_name[idx] != '/' and obj_name[idx] != '-':\n",
    "            obj_name_ed = '{:s}{:s}'.format(obj_name_ed,obj_name[idx])\n",
    "    return obj_name_ed\n",
    "\n",
    "def observatory_code_lookup(o_observatory_id):\n",
    "    # Specify user-readable observatory codes\n",
    "    if   o_observatory_id == '568': o_observatory_code = 'MKO'  # Mauna Kea\n",
    "    elif o_observatory_id == '807': o_observatory_code = 'CTIO' # Cerro Tololo\n",
    "    elif o_observatory_id == '809': o_observatory_code = 'LSO'  # La Silla\n",
    "    elif o_observatory_id == '309': o_observatory_code = 'VLT'  # VLT, Paranal\n",
    "    elif o_observatory_id == 'VLT': o_observatory_code = 'VLT'  # Paranal\n",
    "    elif o_observatory_id == 'I33': o_observatory_code = 'SOAR' # SOAR, Cerro Pachon\n",
    "    elif o_observatory_id == 'I11': o_observatory_code = 'GemS' # Gemini-South, Cerro Pachon\n",
    "    elif o_observatory_id == '304': o_observatory_code = 'LCO'  # Las Campanas\n",
    "    elif o_observatory_id == 'G37': o_observatory_code = 'DCT'  # Discovery Channel Telescope\n",
    "    elif o_observatory_id == 'F51': o_observatory_code = 'PS1'  # Pan-STARRS1\n",
    "    else: o_observatory_code = o_observatory_id\n",
    "\n",
    "    return o_observatory_code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3f522a39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def query_horizons(object_id,obs_code,obs_datetime,time_step,log_filepath):\n",
    "    # Create and execute JPL Horizons query for a given object, observatory code, date, and time step\n",
    "    output_write_log_entry(log_filepath,'Executing Horizons query...')\n",
    "    date_start = Time('{:s}'.format(obs_datetime),scale='utc',format='iso')\n",
    "    dt = TimeDelta(1,format='jd')\n",
    "    date_end = date_start+dt\n",
    "    obj = Horizons(id=object_id, location=obs_code,epochs={'start':date_start.iso,'stop':date_end.iso,'step':time_step})\n",
    "    query_successful = True\n",
    "    max_retries = 50\n",
    "    ephems = ''\n",
    "        \n",
    "    # try comet mode\n",
    "    output_write_log_entry(log_filepath,'Trying Horizons in comet mode...')\n",
    "    for _ in range(max_retries):\n",
    "        try:\n",
    "            if (not os.environ.get('PYTHONHTTPSVERIFY','') and getattr(ssl,'_create_unverified_context',None)):\n",
    "                ssl._create_default_https_context = ssl._create_unverified_context\n",
    "            ephems = obj.ephemerides(closest_apparition=True,no_fragments=True)\n",
    "            print(obj.uri) # print Horizons query URL\n",
    "        except Exception:\n",
    "            print('.',end='')\n",
    "        else:\n",
    "            print('')\n",
    "            #print(ephems)\n",
    "            break\n",
    "    else:\n",
    "        output_write_log_entry(log_filepath,'Horizons query (comet mode) failed. Maximum retries reached.')\n",
    "        query_successful = False\n",
    "        \n",
    "    # try asteroid mode\n",
    "    if not query_successful:\n",
    "        output_write_log_entry(log_filepath,'Trying Horizons in asteroid mode...')\n",
    "        for _ in range(max_retries):\n",
    "            try:\n",
    "                if (not os.environ.get('PYTHONHTTPSVERIFY','') and getattr(ssl,'_create_unverified_context',None)):\n",
    "                    ssl._create_default_https_context = ssl._create_unverified_context\n",
    "                ephems = obj.ephemerides()\n",
    "            except Exception:\n",
    "                print('.',end='')\n",
    "            else:\n",
    "                print('')\n",
    "                #print(ephems)\n",
    "                break\n",
    "        else:\n",
    "            output_write_log_entry(log_filepath,'Horizons query (asteroid mode) failed. Maximum retries reached.')\n",
    "            output_write_log_entry(log_filepath,obj.uri)\n",
    "            query_successful = False\n",
    "            \n",
    "    return ephems\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "01b96f2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_nightly_ephems(object_id,obs_code,obs_datetime,time_step_min,log_filepath):\n",
    "    # Sends a query to JPL Horizons for a given object, observatory code, date, and time step,\n",
    "    #   and outputs date/time, RA, and Dec to an astropy Table\n",
    "    output_write_log_entry(log_filepath,'Retrieving nightly ephemeris table...')\n",
    "    time_step = '{:d}m'.format(time_step_min)\n",
    "    ephems = query_horizons(object_id,obs_code,obs_datetime,time_step,log_filepath)\n",
    "    num_lines_all = len(ephems)\n",
    "    num_lines     = 0\n",
    "    \n",
    "    # Count ephemeris entries that meet night-time (or astronomical twilight) and airmass requirements\n",
    "    for idx in range(0,num_lines_all):\n",
    "        if (ephems[idx]['solar_presence'] == '' or ephems[idx]['solar_presence'] == 'A') and ephems[idx]['airmass'] < 2.7:\n",
    "            num_lines += 1\n",
    "    \n",
    "    # Initialize ephemeris arrays\n",
    "    eph_datetime_str = [0 for idx in range(num_lines)]\n",
    "    eph_ra           = [0 for idx in range(num_lines)]\n",
    "    eph_dec          = [0 for idx in range(num_lines)]\n",
    "\n",
    "    # Save ephemeris entries that meet night-time (or astronomical twilight) and airmass requirements\n",
    "    idx2 = 0\n",
    "    for idx1 in range(0,num_lines_all):\n",
    "        #print(ephems[idx1]['solar_presence'],ephems[idx1]['airmass'])\n",
    "        if (ephems[idx1]['solar_presence'] == '' or ephems[idx1]['solar_presence'] == 'A') and ephems[idx1]['airmass'] < 2.7:\n",
    "            eph_datetime_str[idx2] = ephems[idx1]['datetime_str']\n",
    "            eph_ra[idx2]           = ephems[idx1]['RA']\n",
    "            eph_dec[idx2]          = ephems[idx1]['DEC']\n",
    "            idx2 += 1\n",
    "\n",
    "    # Create ephemeris table for entries that meet night-time (or astronomical twilight) and airmass requirements\n",
    "    #print(eph_datetime_str,eph_ra,eph_dec)\n",
    "    ephem_table,distance_traveled_arcsec,start_time,end_time = '',0,0,0\n",
    "    if eph_datetime_str != []:\n",
    "        ephem_table = Table()\n",
    "        ephem_table = Table([eph_datetime_str,eph_ra,eph_dec],names=('datetime','ra','dec'),dtype=('S17','f8','f8'))\n",
    "    \n",
    "        ephem_table['clear'] = True\n",
    "        ephem_table['nearby_star_idx'] = -1\n",
    "        ephem_table['nearby_star_mag'] = 30.0\n",
    "        ephem_table['nearby_star_dist'] = 9999.0\n",
    "\n",
    "        c1 = SkyCoord((ephem_table[0]['ra']/180*math.pi)*u.radian,(ephem_table[0]['dec']/180*math.pi)*u.radian,frame='icrs')\n",
    "        c2 = SkyCoord((ephem_table[len(ephem_table)-1]['ra']/180*math.pi)*u.radian,(ephem_table[len(ephem_table)-1]['dec']/180*math.pi)*u.radian,frame='icrs')\n",
    "        sep = c1.separation(c2)\n",
    "        distance_traveled_arcsec = sep.arcsecond\n",
    "\n",
    "        # Confirm start and end times and positions of object when it is observable\n",
    "        output_write_log_entry(log_filepath,'Start time and position: {:s} {:f} {:f}'.format(ephem_table[0]['datetime'],ephem_table[0]['ra'],ephem_table[0]['dec']))\n",
    "        output_write_log_entry(log_filepath,'End time and position:   {:s} {:f} {:f}'.format(ephem_table[len(ephem_table)-1]['datetime'],ephem_table[len(ephem_table)-1]['ra'],ephem_table[len(ephem_table)-1]['dec']))\n",
    "        output_write_log_entry(log_filepath,'Distance traveled:       {:.1f} arcsec'.format(distance_traveled_arcsec))\n",
    "        \n",
    "        start_time = ephem_table[0]['datetime']\n",
    "        end_time   = ephem_table[len(ephem_table)-1]['datetime']\n",
    "    \n",
    "    return ephem_table,distance_traveled_arcsec,start_time,end_time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a79e237a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_star_field_gaia(output_path,ephem_table,fov_arcmin,log_filepath):\n",
    "    output_write_log_entry(log_filepath,'Retrieving Gaia sources...')\n",
    "    # Compute/retrieve parameters needed for field star retrieval\n",
    "    num_lines  = len(ephem_table)                  # number of ephemeris entries\n",
    "    center_idx = int(num_lines/2)                  # approximate midpoint of observing window\n",
    "    start_ra   = ephem_table[0]['ra']              # RA at start of observing window\n",
    "    start_dec  = ephem_table[0]['dec']             # Dec at start of observing window\n",
    "    mid_ra     = ephem_table[center_idx]['ra']     # RA at approximate midpoint of observing window\n",
    "    mid_dec    = ephem_table[center_idx]['dec']    # Dec at approximate midpoint of observing window\n",
    "    end_ra     = ephem_table[num_lines-1]['ra']    # RA at end of observing window\n",
    "    end_dec    = ephem_table[num_lines-1]['dec']   # Dec at end of observing window\n",
    "    field_radius_deg = (fov_arcmin / 60 / 2) * (2**0.5)\n",
    "\n",
    "    query_successful = True\n",
    "    max_retries = 50\n",
    "    for _ in range(max_retries):\n",
    "        try:\n",
    "            if (not os.environ.get('PYTHONHTTPSVERIFY','') and getattr(ssl,'_create_unverified_context',None)):\n",
    "                ssl._create_default_https_context = ssl._create_unverified_context\n",
    "            adql_query = \"SELECT * FROM gaiadr1.gaia_source WHERE CONTAINS(POINT('ICRS',gaiadr1.gaia_source.ra,gaiadr1.gaia_source.dec),CIRCLE('ICRS',{:f},{:f},{:f}))=1;\".format(mid_ra,mid_dec,field_radius_deg)\n",
    "            job = Gaia.launch_job_async(adql_query,dump_to_file=False)\n",
    "            gaia_results_table = job.get_results()\n",
    "        except Exception:\n",
    "            print('.',end='')\n",
    "        else:\n",
    "            print('')\n",
    "            break\n",
    "    else:\n",
    "        output_write_log_entry(log_filepath,'Gaia source search failed. Maximum retries reached.')\n",
    "        query_successful = False\n",
    "\n",
    "    if query_successful:\n",
    "        output_write_log_entry(log_filepath,'{:d} Gaia sources found.'.format(len(gaia_results_table)))\n",
    "    \n",
    "    return gaia_results_table,query_successful,mid_ra,mid_dec\n",
    "\n",
    "def retrieve_star_field_gaia_radec(output_path,ra_coord,dec_coord,fov_arcmin,log_filepath):\n",
    "    output_write_log_entry(log_filepath,'Retrieving Gaia sources...')\n",
    "    # Compute/retrieve parameters needed for field star retrieval\n",
    "    mid_ra     = ra_coord\n",
    "    mid_dec    = dec_coord\n",
    "    field_radius_deg = (fov_arcmin / 60 / 2) * (2**0.5)\n",
    "\n",
    "    query_successful = True\n",
    "    max_retries = 50\n",
    "    for _ in range(max_retries):\n",
    "        try:\n",
    "            if (not os.environ.get('PYTHONHTTPSVERIFY','') and getattr(ssl,'_create_unverified_context',None)):\n",
    "                ssl._create_default_https_context = ssl._create_unverified_context\n",
    "            adql_query = \"SELECT * FROM gaiadr1.gaia_source WHERE CONTAINS(POINT('ICRS',gaiadr1.gaia_source.ra,gaiadr1.gaia_source.dec),CIRCLE('ICRS',{:f},{:f},{:f}))=1;\".format(mid_ra,mid_dec,field_radius_deg)\n",
    "            job = Gaia.launch_job_async(adql_query,dump_to_file=False)\n",
    "            gaia_results_table = job.get_results()\n",
    "        except Exception:\n",
    "            print('.',end='')\n",
    "        else:\n",
    "            print('')\n",
    "            break\n",
    "    else:\n",
    "        output_write_log_entry(log_filepath,'Gaia source search failed. Maximum retries reached.')\n",
    "        query_successful = False\n",
    "\n",
    "        \n",
    "    if query_successful:\n",
    "        output_write_log_entry(log_filepath,'{:d} Gaia sources found.'.format(len(gaia_results_table)))\n",
    "    \n",
    "    return gaia_results_table,query_successful\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "83f32a00",
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_refcat_sources(ra_deg,dec_deg,radius_arcmin,source_list_filepath):\n",
    "    # Retrieve refcat sources\n",
    "    refcat_table = 0\n",
    "    source_list_filtered_filepath = source_list_filepath[:-4] + '_filtered.txt'\n",
    "    print('Retrieving refcat sources within {:.1f} arcmin of RA={:.6f},Dec={:.6f}'.format(radius_arcmin,ra_deg,dec_deg))\n",
    "    with open(source_list_filepath,'w') as outf:\n",
    "        outf.write(' objRA       objDec      gmag   gerr  rmag   rerr  imag   ierr  zmag   zerr\\n')\n",
    "    with open(source_list_filepath,'a') as outf:\n",
    "        outf.write(' objRA       objDec      gmag   gerr  rmag   rerr  imag   ierr  zmag   zerr\\n')\n",
    "        radius_deg = radius_arcmin / 60\n",
    "        cmd_ra  = '{:f}'.format(ra_deg)\n",
    "        cmd_dec = '{:f}'.format(dec_deg)\n",
    "        cmd_radius = '{:f}'.format(radius_deg)\n",
    "        cmd_refcat,cmd_dir = '/users/hhsieh/astro/tools/refcat','-dir'\n",
    "        if os.path.exists('/volumes/Fantom12a/refcat_catalog/'):\n",
    "            directories = '/volumes/Fantom12a/refcat_catalog/00_m_16_binary,/volumes/Fantom12a/refcat_catalog/16_m_17_binary,/volumes/Fantom12a/refcat_catalog/17_m_18_binary,/volumes/Fantom12a/refcat_catalog/18_m_19_binary,/volumes/Fantom12a/refcat_catalog/19_m_20_binary'\n",
    "        elif os.path.exists('/volumes/RavPower_1TB_A/refcat_catalog/'):\n",
    "            directories = '/volumes/RavPower_1TB_A/refcat_catalog/00_m_16_binary,/volumes/RavPower_1TB_A/refcat_catalog/16_m_17_binary,/volumes/RavPower_1TB_A/refcat_catalog/17_m_18_binary,/volumes/RavPower_1TB_A/refcat_catalog/18_m_19_binary,/volumes/RavPower_1TB_A/refcat_catalog/19_m_20_binary'\n",
    "            \n",
    "        var_fields = 'RA,Dec,g,dg,r,dr,i,di,z,dz'\n",
    "        cmd_rad,cmd_var,cmd_nohdr = '-rad','-var','-nohdr'\n",
    "        cmd = [cmd_refcat,cmd_ra,cmd_dec,cmd_dir,directories,cmd_rad,cmd_radius,cmd_var,var_fields,cmd_nohdr]\n",
    "        print(cmd)\n",
    "        subprocess.call(cmd,stdout=outf)\n",
    "    print('Refcat catalog sources written to {:s}.'.format(source_list_filepath))\n",
    "            \n",
    "    # parse local file into astropy.table object\n",
    "    objRA  = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(0))\n",
    "    objDec = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(1))\n",
    "    gmag   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(2))\n",
    "    gerr   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(3))\n",
    "    rmag   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(4))\n",
    "    rerr   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(5))\n",
    "    imag   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(6))\n",
    "    ierr   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(7))\n",
    "    zmag   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(8))\n",
    "    zerr   = np.genfromtxt(source_list_filepath,skip_header=1,usecols=(9))\n",
    "    refcat_table = Table([objRA,objDec,gmag,gerr,rmag,rerr,imag,ierr,zmag,zerr], \\\n",
    "        names=('right_ascension','declination','gp1_mag','gp1_err','rp1_mag','rp1_err','ip1_mag','ip1_err','zp1_mag','zp1_err'), \\\n",
    "        dtype=('f8','f8','f8','f8','f8','f8','f8','f8','f8','f8'))\n",
    "    num_refcat_sources = len(refcat_table)\n",
    "    if num_refcat_sources > 0:\n",
    "        idx = 0\n",
    "        # Remove sources that are likely to be saturated or have dummy values (i.e., -999)\n",
    "        while idx < len(refcat_table):\n",
    "            if refcat_table[idx]['gp1_mag'] < 10 or refcat_table[idx]['rp1_mag'] < 10 or refcat_table[idx]['ip1_mag'] < 10 or refcat_table[idx]['zp1_mag'] < 10:\n",
    "                refcat_table.remove_row(idx)\n",
    "            else:\n",
    "                idx += 1\n",
    "        num_refcat_sources = len(refcat_table)\n",
    "    if num_refcat_sources > 0:\n",
    "        print('{:d} refcat catalog sources retrieved.'.format(num_refcat_sources))\n",
    "        with open(source_list_filtered_filepath,'w') as outf:\n",
    "            outf.write('starID   objRA       objDec      gmag   gerr  rmag   rerr  imag   ierr  zmag   zerr\\n')\n",
    "            for idx in range(num_refcat_sources-1):\n",
    "                outf.write('{:>6d}   {:11.7f} {:11.7f} {:6.3f} {:5.3f} {:6.3f} {:5.3f} {:6.3f} {:5.3f} {:6.3f} {:5.3f}\\n'.format(idx+1, \\\n",
    "                             refcat_table[idx]['right_ascension'],refcat_table[idx]['declination'],refcat_table[idx]['gp1_mag'],refcat_table[idx]['gp1_err'], \\\n",
    "                             refcat_table[idx]['rp1_mag'],refcat_table[idx]['rp1_err'],refcat_table[idx]['ip1_mag'],refcat_table[idx]['ip1_err'],\\\n",
    "                             refcat_table[idx]['zp1_mag'],refcat_table[idx]['zp1_err']))\n",
    "    else:\n",
    "        print('No refcat calibration sources found.')\n",
    "                \n",
    "    return refcat_table\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f07d244e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_finder_chart_gaia(output_path,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,gaia_results_table,log_filepath):\n",
    "    # Generates finder chart showing star field (from source catalog), object path\n",
    "    #   during the nightly observing window, and close encounters with background stars\n",
    "    output_write_log_entry(log_filepath,'Generating Gaia finder chart for {:s} for {:s}, {:s}...'.format(object_id,obs_code,obs_datetime))\n",
    "    \n",
    "    majorFormatter0 = FormatStrFormatter('%d')\n",
    "    majorFormatter1 = FormatStrFormatter('%.1f')\n",
    "    majorFormatter2 = FormatStrFormatter('%.2f')\n",
    "    majorFormatter3 = FormatStrFormatter('%.3f')\n",
    "    majorFormatter4 = FormatStrFormatter('%.4f')\n",
    "    \n",
    "    finder_chart_name = 'finder_chart_' + remove_extra_characters(object_id) + '_' + obs_code + '_'+ obs_datetime[0:4] + obs_datetime[5:7] + obs_datetime[8:10]\n",
    "\n",
    "    ra_toplot  = coord.Angle((ephem_table['ra']/180*math.pi)*u.rad)\n",
    "    dec_toplot = coord.Angle((ephem_table['dec']/180*math.pi)*u.rad)\n",
    "\n",
    "    stars_ra_toplot  = coord.Angle((gaia_results_table['ra']/180*math.pi)*u.rad)\n",
    "    stars_dec_toplot = coord.Angle((gaia_results_table['dec']/180*math.pi)*u.rad)\n",
    "\n",
    "    # initialize and populate star size array\n",
    "    star_size = [0 for idx in range(len(gaia_results_table))]\n",
    "    for idx in range(0,len(gaia_results_table)):\n",
    "        star_size[idx] = 55 * (10**(0.15*(21 - gaia_results_table[idx]['phot_g_mean_mag'])))\n",
    "    \n",
    "    num_lines  = len(ephem_table)\n",
    "    center_idx = int(num_lines/2)\n",
    "    mid_ra     = ephem_table[center_idx]['ra']\n",
    "    mid_dec    = ephem_table[center_idx]['dec']\n",
    "\n",
    "    a_mid_ra   = Angle(mid_ra,u.deg)\n",
    "    a_mid_dec  = Angle(mid_dec,u.deg)\n",
    "\n",
    "    x_min  = mid_ra-(0.8*fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    x_max  = mid_ra+(0.8*fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    y_min  = mid_dec-(0.8*fov_arcmin/2/60)\n",
    "    y_max  = mid_dec+(0.8*fov_arcmin/2/60)\n",
    "    xlim   = [x_min,x_max]\n",
    "    ylim   = [y_min,y_max]\n",
    "\n",
    "    major_xticks = np.arange(int(xlim[0]-1),int(xlim[1]+1),0.05)\n",
    "    minor_xticks = np.arange(int(xlim[0]-1),int(xlim[1]+1),0.01)\n",
    "    major_yticks = np.arange(int(ylim[0]-1),int(ylim[1]+1),0.01)\n",
    "    minor_yticks = np.arange(int(ylim[0]-1),int(ylim[1]+1),0.002)\n",
    "    \n",
    "    fig = plt.figure(figsize=(15,15))\n",
    "    rect_plot = [0.10,0.10,0.60,0.60]\n",
    "    axPlot  = plt.axes(rect_plot)\n",
    "    axPlot.plot(ra_toplot.deg,dec_toplot.deg,ls='--',zorder=4)\n",
    "    axPlot.scatter(stars_ra_toplot.deg,stars_dec_toplot.deg,marker='.',color='#000000',s=star_size,rasterized=True,zorder=2)\n",
    "    window_start = axPlot.scatter(ra_toplot[0].deg,dec_toplot[0].deg,marker='+',s=150,lw=3,color='#00a822',zorder=7)\n",
    "    window_end   = axPlot.scatter(ra_toplot[len(ra_toplot)-1].deg,dec_toplot[len(dec_toplot)-1].deg,marker='x',s=80,lw=3,color='#cc0e00',zorder=7)\n",
    "    axPlot.xaxis.set_major_formatter(majorFormatter2)\n",
    "    axPlot.yaxis.set_major_formatter(majorFormatter2)\n",
    "    axPlot.tick_params(axis='both',which='both',direction='in',labelsize=10)\n",
    "    axPlot.set_xlabel('Right ascension (deg)',fontsize=14,labelpad=10)\n",
    "    axPlot.set_ylabel('Declination (deg)',fontsize=14,labelpad=10)\n",
    "    axPlot.set_xlim(xlim)\n",
    "    axPlot.set_ylim(ylim)\n",
    "    axPlot.set_xticks(major_xticks)\n",
    "    axPlot.set_xticks(minor_xticks,minor=True)\n",
    "    axPlot.set_yticks(major_yticks)\n",
    "    axPlot.set_yticks(minor_yticks,minor=True)\n",
    "    axPlot.set_xlim(xlim)\n",
    "    axPlot.set_ylim(ylim)\n",
    "    axPlot.invert_xaxis()\n",
    "\n",
    "    for idx in range(len(stars_ra_toplot)):\n",
    "        if (stars_ra_toplot.deg[idx] > x_min) and (stars_ra_toplot.deg[idx] < x_max) and (stars_dec_toplot.deg[idx] > y_min) and (stars_dec_toplot.deg[idx] < y_max):\n",
    "            axPlot.text(stars_ra_toplot.deg[idx]-0.002,stars_dec_toplot.deg[idx]-0.0007,'{:d}'.format(idx+1),fontsize=10,color='#0000ff',zorder=3)\n",
    "\n",
    "    print(start_time,end_time)\n",
    "    \n",
    "    props1 = dict(facecolor='white',alpha=1.0,boxstyle='round,pad=0.50',edgecolor='black')\n",
    "    props2 = dict(facecolor='white',alpha=1.0,boxstyle='round,pad=0.50',edgecolor='white')\n",
    "    axPlot.text(0.03,0.035,'{:s} ({:s}, UT {:s} - {:s})\\nField center: {:s} {:s}\\nField size: {:.2f} arcmin x {:.2f} arcmin'.format(object_id,observatory_code_lookup(obs_code),start_time,end_time[-5:],a_mid_ra.to_string(unit=u.hour,sep=':'),a_mid_dec.to_string(unit=u.degree,sep=':'),fov_arcmin,fov_arcmin),fontsize=12,transform=axPlot.transAxes,bbox=props1)\n",
    "\n",
    "    axPlot.grid(True,zorder=1)\n",
    "\n",
    "    axPlot.legend([window_start,window_end],['Observing window start','Observing window end'],ncol=1,scatterpoints=1,handletextpad=0.2,columnspacing=0.,fontsize=13,edgecolor='#000000',framealpha=1,labelspacing=0.4,loc='upper right',bbox_to_anchor=(0.99,0.99))\n",
    "\n",
    "    plt.draw()\n",
    "    plot_filename = finder_chart_name\n",
    "    plot_filepath = output_path + plot_filename + '_gaia.pdf'\n",
    "    plt.savefig(plot_filepath,format='pdf',transparent=True,dpi=400,bbox_inches='tight')\n",
    "    plt.clf()\n",
    "    plt.cla()\n",
    "    plt.close()\n",
    "    output_write_log_entry(log_filepath,'>>> Finder chart written to {:s}.'.format(plot_filepath))\n",
    "\n",
    "    return plot_filepath\n",
    "\n",
    "def generate_finder_chart_gaia_radec(output_path,ra_coord,dec_coord,fov_arcmin,gaia_results_table,log_filepath):\n",
    "    # Generates finder chart showing star field (from source catalog), object path\n",
    "    #   during the nightly observing window, and close encounters with background stars\n",
    "    output_write_log_entry(log_filepath,'Generating Gaia finder chart for ra={:.6f}, dec={:.6f}...'.format(ra_coord,dec_coord))\n",
    "    \n",
    "    majorFormatter0 = FormatStrFormatter('%d')\n",
    "    majorFormatter1 = FormatStrFormatter('%.1f')\n",
    "    majorFormatter2 = FormatStrFormatter('%.2f')\n",
    "    majorFormatter3 = FormatStrFormatter('%.3f')\n",
    "    majorFormatter4 = FormatStrFormatter('%.4f')\n",
    "    \n",
    "    finder_chart_name = 'finder_chart_ra{:.6f}_dec{:.6f}'.format(ra_coord,dec_coord)\n",
    "\n",
    "    ra_toplot  = coord.Angle((ra_coord/180*math.pi)*u.rad)\n",
    "    dec_toplot = coord.Angle((dec_coord/180*math.pi)*u.rad)\n",
    "\n",
    "    stars_ra_toplot  = coord.Angle((gaia_results_table['ra']/180*math.pi)*u.rad)\n",
    "    stars_dec_toplot = coord.Angle((gaia_results_table['dec']/180*math.pi)*u.rad)\n",
    "\n",
    "    # initialize and populate star size array\n",
    "    star_size = [0 for idx in range(len(gaia_results_table))]\n",
    "    for idx in range(0,len(gaia_results_table)):\n",
    "        star_size[idx] = 55 * (10**(0.15*(21 - gaia_results_table[idx]['phot_g_mean_mag'])))\n",
    "    \n",
    "    a_mid_ra   = Angle(ra_coord,u.deg)\n",
    "    a_mid_dec  = Angle(dec_coord,u.deg)\n",
    "\n",
    "    x_min  = mid_ra-(0.8*fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    x_max  = mid_ra+(0.8*fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    y_min  = mid_dec-(0.8*fov_arcmin/2/60)\n",
    "    y_max  = mid_dec+(0.8*fov_arcmin/2/60)\n",
    "    xlim   = [x_min,x_max]\n",
    "    ylim   = [y_min,y_max]\n",
    "\n",
    "    major_xticks = np.arange(int(xlim[0]-1),int(xlim[1]+1),0.05)\n",
    "    minor_xticks = np.arange(int(xlim[0]-1),int(xlim[1]+1),0.01)\n",
    "    major_yticks = np.arange(int(ylim[0]-1),int(ylim[1]+1),0.01)\n",
    "    minor_yticks = np.arange(int(ylim[0]-1),int(ylim[1]+1),0.002)\n",
    "    \n",
    "    fig = plt.figure(figsize=(15,15))\n",
    "    rect_plot = [0.10,0.10,0.60,0.60]\n",
    "    axPlot  = plt.axes(rect_plot)\n",
    "    axPlot.plot(ra_toplot.deg,dec_toplot.deg,ls='--',zorder=4)\n",
    "    axPlot.scatter(stars_ra_toplot.deg,stars_dec_toplot.deg,marker='.',color='#000000',s=star_size,rasterized=True,zorder=2)\n",
    "    window_start = axPlot.scatter(ra_toplot.deg,dec_toplot.deg,marker='+',s=150,lw=3,color='#00a822',zorder=7)\n",
    "    axPlot.xaxis.set_major_formatter(majorFormatter2)\n",
    "    axPlot.yaxis.set_major_formatter(majorFormatter2)\n",
    "    axPlot.tick_params(axis='both',which='both',direction='in',labelsize=10)\n",
    "    axPlot.set_xlabel('Right ascension (deg)',fontsize=14,labelpad=10)\n",
    "    axPlot.set_ylabel('Declination (deg)',fontsize=14,labelpad=10)\n",
    "    axPlot.set_xlim(xlim)\n",
    "    axPlot.set_ylim(ylim)\n",
    "    axPlot.set_xticks(major_xticks)\n",
    "    axPlot.set_xticks(minor_xticks,minor=True)\n",
    "    axPlot.set_yticks(major_yticks)\n",
    "    axPlot.set_yticks(minor_yticks,minor=True)\n",
    "    axPlot.set_xlim(xlim)\n",
    "    axPlot.set_ylim(ylim)\n",
    "    axPlot.invert_xaxis()\n",
    "\n",
    "    for idx in range(len(stars_ra_toplot)):\n",
    "        if (stars_ra_toplot.deg[idx] > x_min) and (stars_ra_toplot.deg[idx] < x_max) and (stars_dec_toplot.deg[idx] > y_min) and (stars_dec_toplot.deg[idx] < y_max):\n",
    "            axPlot.text(stars_ra_toplot.deg[idx]-0.002,stars_dec_toplot.deg[idx]-0.0007,'{:d}'.format(idx+1),fontsize=10,color='#0000ff',zorder=3)\n",
    "    \n",
    "    props1 = dict(facecolor='white',alpha=1.0,boxstyle='round,pad=0.50',edgecolor='black')\n",
    "    props2 = dict(facecolor='white',alpha=1.0,boxstyle='round,pad=0.50',edgecolor='white')\n",
    "    axPlot.text(0.03,0.035,'{:s} (Field center: {:s} {:s}\\nField size: {:.2f} arcmin x {:.2f} arcmin'.format(a_mid_ra.to_string(unit=u.hour,sep=':'),a_mid_dec.to_string(unit=u.degree,sep=':'),fov_arcmin,fov_arcmin),fontsize=12,transform=axPlot.transAxes,bbox=props1)\n",
    "\n",
    "    axPlot.grid(True,zorder=1)\n",
    "\n",
    "    #axPlot.legend([window_start,window_end],['Observing window start','Observing window end'],ncol=1,scatterpoints=1,handletextpad=0.2,columnspacing=0.,fontsize=13,edgecolor='#000000',framealpha=1,labelspacing=0.4,loc='upper right',bbox_to_anchor=(0.99,0.99))\n",
    "\n",
    "    plt.draw()\n",
    "    plot_filename = finder_chart_name\n",
    "    plot_filepath = output_path + plot_filename + '_gaia.pdf'\n",
    "    plt.savefig(plot_filepath,format='pdf',transparent=True,dpi=400,bbox_inches='tight')\n",
    "    plt.clf()\n",
    "    plt.cla()\n",
    "    plt.close()\n",
    "    output_write_log_entry(log_filepath,'>>> Finder chart written to {:s}.'.format(plot_filepath))\n",
    "\n",
    "    return plot_filepath\n",
    "\n",
    "\n",
    "def generate_finder_chart_gaia_refcat(output_path,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,gaia_results_table,refcat_table,log_filepath):\n",
    "    # Generates finder chart showing star field (from source catalog), object path\n",
    "    #   during the nightly observing window, and close encounters with background stars\n",
    "    output_write_log_entry(log_filepath,'Generating Gaia finder chart for {:s} for {:s}, {:s}...'.format(object_id,obs_code,obs_datetime))\n",
    "    \n",
    "    majorFormatter0 = FormatStrFormatter('%d')\n",
    "    majorFormatter1 = FormatStrFormatter('%.1f')\n",
    "    majorFormatter2 = FormatStrFormatter('%.2f')\n",
    "    majorFormatter3 = FormatStrFormatter('%.3f')\n",
    "    majorFormatter4 = FormatStrFormatter('%.4f')\n",
    "    \n",
    "    finder_chart_name = 'finder_chart_' + remove_extra_characters(object_id) + '_' + obs_code + '_'+ remove_extra_characters(obs_datetime[:10])\n",
    "\n",
    "    ra_toplot  = coord.Angle((ephem_table['ra']/180*math.pi)*u.rad)\n",
    "    dec_toplot = coord.Angle((ephem_table['dec']/180*math.pi)*u.rad)\n",
    "\n",
    "    stars_ra_toplot  = coord.Angle((gaia_results_table['ra']/180*math.pi)*u.rad)\n",
    "    stars_dec_toplot = coord.Angle((gaia_results_table['dec']/180*math.pi)*u.rad)\n",
    "\n",
    "    # initialize and populate star size array\n",
    "    star_size = [0 for idx in range(len(gaia_results_table))]\n",
    "    for idx in range(0,len(gaia_results_table)):\n",
    "        star_size[idx] = 55 * (10**(0.15*(21 - gaia_results_table[idx]['phot_g_mean_mag'])))\n",
    "    \n",
    "    num_lines  = len(ephem_table)\n",
    "    center_idx = int(num_lines/2)\n",
    "    mid_ra     = ephem_table[center_idx]['ra']\n",
    "    mid_dec    = ephem_table[center_idx]['dec']\n",
    "\n",
    "    a_mid_ra   = Angle(mid_ra,u.deg)\n",
    "    a_mid_dec  = Angle(mid_dec,u.deg)\n",
    "\n",
    "    x_min  = mid_ra-(0.8*fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    x_max  = mid_ra+(0.8*fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    y_min  = mid_dec-(0.8*fov_arcmin/2/60)\n",
    "    y_max  = mid_dec+(0.8*fov_arcmin/2/60)\n",
    "    xlim   = [x_min,x_max]\n",
    "    ylim   = [y_min,y_max]\n",
    "\n",
    "    major_xticks = np.arange(int(xlim[0]-1),int(xlim[1]+1),0.05)\n",
    "    minor_xticks = np.arange(int(xlim[0]-1),int(xlim[1]+1),0.01)\n",
    "    major_yticks = np.arange(int(ylim[0]-1),int(ylim[1]+1),0.01)\n",
    "    minor_yticks = np.arange(int(ylim[0]-1),int(ylim[1]+1),0.002)\n",
    "    \n",
    "    fig = plt.figure(figsize=(15,15))\n",
    "    rect_plot = [0.10,0.10,0.60,0.60]\n",
    "    axPlot  = plt.axes(rect_plot)\n",
    "    axPlot.plot(ra_toplot.deg,dec_toplot.deg,ls='--',zorder=4)\n",
    "    axPlot.scatter(stars_ra_toplot.deg,stars_dec_toplot.deg,marker='.',color='#000000',s=star_size,rasterized=True,zorder=2)\n",
    "    window_start = axPlot.scatter(ra_toplot[0].deg,dec_toplot[0].deg,marker='+',s=150,lw=3,color='#00a822',zorder=7)\n",
    "    window_end   = axPlot.scatter(ra_toplot[len(ra_toplot)-1].deg,dec_toplot[len(dec_toplot)-1].deg,marker='x',s=80,lw=3,color='#cc0e00',zorder=7)\n",
    "    axPlot.xaxis.set_major_formatter(majorFormatter2)\n",
    "    axPlot.yaxis.set_major_formatter(majorFormatter2)\n",
    "    axPlot.tick_params(axis='both',which='both',direction='in',labelsize=10)\n",
    "    axPlot.set_xlabel('Right ascension (deg)',fontsize=14,labelpad=10)\n",
    "    axPlot.set_ylabel('Declination (deg)',fontsize=14,labelpad=10)\n",
    "    axPlot.set_xlim(xlim)\n",
    "    axPlot.set_ylim(ylim)\n",
    "    axPlot.set_xticks(major_xticks)\n",
    "    axPlot.set_xticks(minor_xticks,minor=True)\n",
    "    axPlot.set_yticks(major_yticks)\n",
    "    axPlot.set_yticks(minor_yticks,minor=True)\n",
    "    axPlot.set_xlim(xlim)\n",
    "    axPlot.set_ylim(ylim)\n",
    "    axPlot.invert_xaxis()\n",
    "\n",
    "    for idx in range(len(refcat_table)):\n",
    "        if (refcat_table[idx]['right_ascension'] > x_min) and (refcat_table[idx]['right_ascension'] < x_max) and (refcat_table[idx]['declination'] > y_min) and (refcat_table[idx]['declination'] < y_max):\n",
    "            axPlot.text(refcat_table[idx]['right_ascension']-0.002,refcat_table[idx]['declination']-0.0007,'{:d}'.format(idx+1),fontsize=10,color='#0000ff',zorder=3)\n",
    "\n",
    "    print(start_time,end_time)\n",
    "    \n",
    "    props1 = dict(facecolor='white',alpha=1.0,boxstyle='round,pad=0.50',edgecolor='black')\n",
    "    props2 = dict(facecolor='white',alpha=1.0,boxstyle='round,pad=0.50',edgecolor='white')\n",
    "    axPlot.text(0.03,0.035,'{:s} ({:s}, UT {:s} - {:s})\\nField center: {:s} {:s}\\nField size: {:.2f} arcmin x {:.2f} arcmin'.format(object_id,observatory_code_lookup(obs_code),start_time,end_time[-5:],a_mid_ra.to_string(unit=u.hour,sep=':'),a_mid_dec.to_string(unit=u.degree,sep=':'),fov_arcmin,fov_arcmin),fontsize=12,transform=axPlot.transAxes,bbox=props1)\n",
    "\n",
    "    axPlot.grid(True,zorder=1)\n",
    "\n",
    "    axPlot.legend([window_start,window_end],['Observing window start','Observing window end'],ncol=1,scatterpoints=1,handletextpad=0.2,columnspacing=0.,fontsize=13,edgecolor='#000000',framealpha=1,labelspacing=0.4,loc='upper right',bbox_to_anchor=(0.99,0.99))\n",
    "\n",
    "    plt.draw()\n",
    "    plot_filename = finder_chart_name\n",
    "    plot_filepath = output_path + plot_filename + '_gaia.pdf'\n",
    "    plt.savefig(plot_filepath,format='pdf',transparent=True,dpi=400,bbox_inches='tight')\n",
    "    plt.clf()\n",
    "    plt.cla()\n",
    "    plt.close()\n",
    "    output_write_log_entry(log_filepath,'>>> Finder chart written to {:s}.'.format(plot_filepath))\n",
    "\n",
    "    return plot_filepath\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8b20659b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_radec_to_xypix_dss(ra_deg,dec_deg,ra_max_deg,dec_max_deg,ra_total_deg,dec_total_deg,fov_arcmin):\n",
    "    xpix = (ra_max_deg - ra_deg)/ra_total_deg * (fov_arcmin*60)\n",
    "    ypix = (dec_max_deg - dec_deg)/dec_total_deg * (fov_arcmin*60)\n",
    "    return xpix,ypix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d3740ea0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_finder_chart_dss(output_path,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,log_filepath):\n",
    "    # Generates finder chart showing star field (from DSS image) and object path during the nightly observing window\n",
    "    output_write_log_entry(log_filepath,'Generating DSS finder chart for {:s} for {:s}, {:s}...'.format(object_id,obs_code,obs_datetime))\n",
    "    \n",
    "    finder_chart_name = 'finder_chart_' + remove_extra_characters(object_id) + '_' + obs_code + '_'+ obs_datetime[0:4] + obs_datetime[5:7] + obs_datetime[8:10]\n",
    "    sky_field_image_filename_png_final = 'finder_chart_' + remove_extra_characters(object_id) + '_' + obs_code + '_'+ obs_datetime[0:4] + obs_datetime[5:7] + obs_datetime[8:10] + '_dss.png'\n",
    "    \n",
    "    obs_datetime_start = Time(obs_datetime,scale='utc',format='iso')\n",
    "    fov_arcmin = 0.8 * fov_arcmin\n",
    "\n",
    "    ra_toplot  = coord.Angle((ephem_table['ra']/180*math.pi)*u.rad)\n",
    "    dec_toplot = coord.Angle((ephem_table['dec']/180*math.pi)*u.rad)\n",
    "\n",
    "    num_lines  = len(ephem_table)\n",
    "    center_idx = int(num_lines/2)\n",
    "    mid_ra     = ephem_table[center_idx]['ra']\n",
    "    mid_dec    = ephem_table[center_idx]['dec']\n",
    "\n",
    "    a_mid_ra   = Angle(mid_ra,u.deg)\n",
    "    a_mid_dec  = Angle(mid_dec,u.deg)\n",
    "\n",
    "    ra_min_deg    = mid_ra-(fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    ra_max_deg    = mid_ra+(fov_arcmin/2/60/np.cos(a_mid_dec.rad))\n",
    "    dec_min_deg   = mid_dec-(fov_arcmin/2/60)\n",
    "    dec_max_deg   = mid_dec+(fov_arcmin/2/60)\n",
    "    ra_total_deg  = ra_max_deg - ra_min_deg\n",
    "    dec_total_deg = dec_max_deg - dec_min_deg\n",
    "    xlim          = [ra_min_deg,ra_max_deg]\n",
    "    ylim          = [dec_min_deg,dec_max_deg]\n",
    "    \n",
    "    cutoutbaseurl = 'http://archive.stsci.edu/cgi-bin/dss_search?'\n",
    "    query_string = urlencode(dict(v='poss2ukstu',r=mid_ra,d=mid_dec,e='J2000',h=fov_arcmin,w=fov_arcmin,f='gif'))\n",
    "    url = cutoutbaseurl + query_string\n",
    "\n",
    "    # this downloads the image to your disk\n",
    "    sky_field_image_filename_gif = 'sky_field_' + remove_extra_characters(object_id) + '_' + remove_extra_characters(obs_datetime_start.iso[:10]) + '.gif'\n",
    "    sky_field_image_filename_png = 'sky_field_' + remove_extra_characters(object_id) + '_' + remove_extra_characters(obs_datetime_start.iso[:10]) + '.png'\n",
    "    \n",
    "    query_successful = True\n",
    "    max_retries = 50\n",
    "    ephems = ''\n",
    "    output_write_log_entry(log_filepath,'{:s}'.format(url))\n",
    "    for _ in range(max_retries):\n",
    "        try:\n",
    "            urlretrieve(url,sky_field_image_filename_gif)\n",
    "        except Exception:\n",
    "            print('.',end='')\n",
    "        else:\n",
    "            print('')\n",
    "            break\n",
    "    else:\n",
    "        output_write_log_entry(log_filepath,'DSS image retrieval failed. Maximum retries reached.')\n",
    "        query_successful = False\n",
    "    \n",
    "    if query_successful:\n",
    "        im = Image.open(sky_field_image_filename_gif)\n",
    "        im.save(sky_field_image_filename_png)\n",
    "        os.remove(sky_field_image_filename_gif)\n",
    "        \n",
    "        img = mpimg.imread(sky_field_image_filename_png)\n",
    "        \n",
    "        majorFormatter0 = FormatStrFormatter('%d')\n",
    "        majorFormatter1 = FormatStrFormatter('%.1f')\n",
    "        majorFormatter2 = FormatStrFormatter('%.2f')\n",
    "        majorFormatter3 = FormatStrFormatter('%.3f')\n",
    "        majorFormatter4 = FormatStrFormatter('%.4f')\n",
    "\n",
    "        fig,ax = plt.subplots()\n",
    "        plt.figure(1, figsize=(25,25))\n",
    "    \n",
    "        img = mpimg.imread(sky_field_image_filename_png)\n",
    "        imgplot = plt.imshow(img,cmap='Greys')\n",
    "\n",
    "        xpix_start,ypix_start = convert_radec_to_xypix_dss(ra_toplot[0].deg,dec_toplot[0].deg,ra_max_deg,dec_max_deg,ra_total_deg,dec_total_deg,fov_arcmin)\n",
    "        xpix_end,ypix_end = convert_radec_to_xypix_dss(ra_toplot[len(dec_toplot)-1].deg,dec_toplot[len(dec_toplot)-1].deg,ra_max_deg,dec_max_deg,ra_total_deg,dec_total_deg,fov_arcmin)\n",
    "    \n",
    "        window_start = ax.scatter(xpix_start,ypix_start,marker='+',s=75,lw=1.5,color='#00a822',zorder=8)\n",
    "        window_end   = ax.scatter(xpix_end,ypix_end,marker='x',s=40,lw=1.5,color='#cc0e00',zorder=8)\n",
    "    \n",
    "        ra_toplot_pix = ra_toplot.deg\n",
    "        dec_toplot_pix = dec_toplot.deg\n",
    "\n",
    "        for idx in range(0,len(ra_toplot)):\n",
    "            ra_toplot_pix[idx],dec_toplot_pix[idx] = convert_radec_to_xypix_dss(ra_toplot[idx].deg,dec_toplot[idx].deg,ra_max_deg,dec_max_deg,ra_total_deg,dec_total_deg,fov_arcmin)\n",
    "        \n",
    "        ax.plot(ra_toplot_pix,dec_toplot_pix,ls='--',zorder=7)\n",
    "\n",
    "        nullfmt = NullFormatter()\n",
    "        ax.xaxis.set_major_formatter(nullfmt)  # turn off labels\n",
    "        ax.yaxis.set_major_formatter(nullfmt)  # turn off labels\n",
    "        ax.tick_params(axis='both',which='both',direction='inout',top=True,bottom=True,right=True,left=True)\n",
    "\n",
    "        plt.draw()\n",
    "        plt.savefig(sky_field_image_filename_png_final)\n",
    "        plt.show()\n",
    "        os.remove(sky_field_image_filename_png)\n",
    "        \n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "        output_write_log_entry(log_filepath,'>>> DSS finder chart written to {:s}.'.format(sky_field_image_filename_png_final))\n",
    "\n",
    "    return sky_field_image_filename_png_final\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3809d954",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_fscoords_finder_chart_gaia_ephems_refcat(object_id,obs_code,obs_datetime,output_dirpath,time_step_min):\n",
    "\n",
    "    obs_datetime_time = Time('{:s}'.format(obs_datetime),scale='utc',format='iso')\n",
    "\n",
    "    log_filepath = output_dirpath + 'log_fscoords_retrieval_finder_chart_generation_{:s}_{:s}_{:s}.txt'.format(remove_extra_characters(object_id),obs_code,datetime.datetime.today().strftime('%Y%m%d_%H%M%S'))\n",
    "    with open(log_filepath,'w') as logfile:\n",
    "        print('\\n{:s} - Retrieving field star coordinates and generating finder charts for {:s} for {:s}...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),object_id,obs_datetime_time.iso))\n",
    "        logfile.write('{:s} - Retrieving field star coordinates and generating finder charts for {:s} for {:s}...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),object_id,obs_datetime_time.iso))\n",
    "\n",
    "    fscoords_filepath = output_dirpath + 'field_star_radec_coords_{:s}_{:s}_{:s}.txt'.format(remove_extra_characters(object_id),obs_code,remove_extra_characters(obs_datetime[:10]))\n",
    "        \n",
    "    os.chdir(output_dirpath)\n",
    "    horizons_obj_desig = get_horizons_obj_desig(object_id)\n",
    "    ephem_table,distance_traveled_arcsec,start_time,end_time = retrieve_nightly_ephems(horizons_obj_desig,obs_code,obs_datetime,time_step_min,log_filepath)\n",
    "    if ephem_table != '':\n",
    "        fov_arcmin = int(3 * (distance_traveled_arcsec) / 60) + 1\n",
    "        if fov_arcmin < 15: fov_arcmin = 15\n",
    "        output_write_log_entry(log_filepath,'Using FOV of {:.1f} x {:.1f} arcmin'.format(fov_arcmin,fov_arcmin))\n",
    "        gaia_results_table,query_successful,ra_deg,dec_deg  = retrieve_star_field_gaia(output_dirpath,ephem_table,fov_arcmin,log_filepath)\n",
    "        if query_successful:\n",
    "            refcat_source_list_filepath = output_dirpath + 'refcat_source_list_{:s}_{:s}_{:s}.txt'.format(remove_extra_characters(object_id),obs_code,remove_extra_characters(obs_datetime[:10]))\n",
    "            refcat_table = retrieve_refcat_sources(ra_deg,dec_deg,fov_arcmin/2,refcat_source_list_filepath)\n",
    "            plot_gaia_filepath = generate_finder_chart_gaia_refcat(output_dirpath,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,gaia_results_table,refcat_table,log_filepath)\n",
    "            with open(fscoords_filepath,'w') as of:\n",
    "                of.write('Field star coordinates for {:s} on {:s} from {:s}:\\n'.format(object_id,obs_datetime,obs_code))\n",
    "                for idx in range(len(gaia_results_table)):\n",
    "                    of.write('{:>3d}  {:f}  {:f}\\n'.format(idx+1,gaia_results_table['ra'][idx],gaia_results_table['dec'][idx]))\n",
    "        plot_dss_filepath = generate_finder_chart_dss(output_dirpath,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,log_filepath)\n",
    "            \n",
    "    else:\n",
    "        output_write_log_entry(log_filepath,'No valid ephemerides found.')\n",
    "            \n",
    "    output_write_log_entry(log_filepath,'Done.')\n",
    "    return None\n",
    "\n",
    "\n",
    "def generate_fscoords_finder_chart_gaia_coords(ra_coord,dec_coord,fov_arcmin,output_dirpath):\n",
    "\n",
    "    obs_datetime_time = Time('{:s}'.format(obs_datetime),scale='utc',format='iso')\n",
    "\n",
    "    log_filepath = output_dirpath + 'log_fscoords_retrieval_finder_chart_generation_{:.6f}_{:.6f}.txt'.format(ra_coord,dec_coord,datetime.datetime.today().strftime('%Y%m%d_%H%M%S'))\n",
    "    with open(log_filepath,'w') as logfile:\n",
    "        print('\\n{:s} - Retrieving field star coordinates and generating finder charts for {:.6f},{:.6f}...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),ra_coord,dec_coord))\n",
    "        logfile.write('{:s} - Retrieving field star coordinates and generating finder charts for {:.6f},{:.6f}...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),ra_coord,dec_coord))\n",
    "\n",
    "    fscoords_filepath = output_dirpath + 'field_star_radec_coords_{:.6f}_{:.6f}.txt'.format(ra_coord,dec_coord)\n",
    "        \n",
    "    os.chdir(output_dirpath)\n",
    "    output_write_log_entry(log_filepath,'Using FOV of {:.1f} x {:.1f} arcmin'.format(fov_arcmin,fov_arcmin))\n",
    "    \n",
    "    gaia_results_table,query_successful,mid_ra,mid_dec  = retrieve_star_field_gaia_radec(output_dirpath,ra_coord,dec_coord,fov_arcmin,log_filepath)\n",
    "        \n",
    "    if query_successful:\n",
    "        plot_gaia_filepath = generate_finder_chart_gaia_radec(output_dirpath,ra_coord,dec_coord,fov_arcmin,gaia_results_table,log_filepath)\n",
    "        with open(fscoords_filepath,'w') as of:\n",
    "            of.write('Field star coordinates for {:s} on {:s} from {:s}:\\n'.format(object_id,obs_datetime,obs_code))\n",
    "            for idx in range(len(gaia_results_table)):\n",
    "                of.write('{:>3d}  {:f}  {:f}\\n'.format(idx+1,gaia_results_table['ra'][idx],gaia_results_table['dec'][idx]))\n",
    "    plot_dss_filepath = generate_finder_chart_dss(output_dirpath,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,log_filepath)\n",
    "            \n",
    "    output_write_log_entry(log_filepath,'Done.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "73ebba39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_dss_finder_chart(object_id,obs_code,obs_datetime,output_path,log_filepath):\n",
    "    output_write_log_entry(log_filepath,'Starting observation assessment and finder chart generation for {:s} for {:s}, {:s}...'.format(object_id,obs_code,obs_datetime))\n",
    "\n",
    "    os.chdir(output_path)\n",
    "    time_step_min = 5\n",
    "    horizons_obj_desig = get_horizons_obj_desig(object_id)\n",
    "    ephem_table,distance_traveled_arcsec,start_time,end_time = retrieve_nightly_ephems_logfilepath(horizons_obj_desig,obs_code,obs_datetime,time_step_min,log_filepath)\n",
    "    \n",
    "    if ephem_table != '':\n",
    "        fov_arcmin = int(3 * (distance_traveled_arcsec) / 60) + 1\n",
    "        if fov_arcmin < 15: fov_arcmin = 15\n",
    "        output_write_log_entry(log_filepath,'Using FOV of {:.1f} x {:.1f} arcmin'.format(fov_arcmin,fov_arcmin))\n",
    "        plot_dss_filepath = generate_finder_chart_dss(output_path,object_id,obs_code,obs_datetime,fov_arcmin,start_time,end_time,ephem_table,log_filepath)\n",
    "    else:\n",
    "        output_write_log_entry(log_filepath,'No valid ephemerides found.')\n",
    "            \n",
    "    output_write_log_entry(log_filepath,'Done.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "60bc09ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "2023-11-27 09:02:56 - Retrieving field star coordinates and generating finder charts for 2023 RN3 for 2023-11-18 00:00:00.000...\n",
      "2023-11-27 09:02:56 - Retrieving nightly ephemeris table...\n",
      "2023-11-27 09:02:56 - Executing Horizons query...\n",
      "2023-11-27 09:02:56 - Trying Horizons in comet mode...\n",
      "https://ssd.jpl.nasa.gov/api/horizons.api?format=text&EPHEM_TYPE=OBSERVER&QUANTITIES=%271%2C2%2C3%2C4%2C5%2C6%2C7%2C8%2C9%2C10%2C11%2C12%2C13%2C14%2C15%2C16%2C17%2C18%2C19%2C20%2C21%2C22%2C23%2C24%2C25%2C26%2C27%2C28%2C29%2C30%2C31%2C32%2C33%2C34%2C35%2C36%2C37%2C38%2C39%2C40%2C41%2C42%2C43%27&COMMAND=%222023+RN3%22&SOLAR_ELONG=%220%2C180%22&LHA_CUTOFF=0&CSV_FORMAT=YES&CAL_FORMAT=BOTH&ANG_FORMAT=DEG&APPARENT=AIRLESS&REF_SYSTEM=ICRF&EXTRA_PREC=NO&CENTER=%27F65%27&START_TIME=%222023-11-18+00%3A00%3A00.000%22&STOP_TIME=%222023-11-19+00%3A00%3A00.000%22&STEP_SIZE=%225m%22&SKIP_DAYLT=NO\n",
      "\n",
      "2023-11-27 09:02:59 - Start time and position: 2023-Nov-18 04:35 357.622480 3.205050\n",
      "2023-11-27 09:02:59 - End time and position:   2023-Nov-18 11:00 357.623000 3.197180\n",
      "2023-11-27 09:02:59 - Distance traveled:       28.4 arcsec\n",
      "2023-11-27 09:02:59 - Using FOV of 15.0 x 15.0 arcmin\n",
      "2023-11-27 09:02:59 - Retrieving Gaia sources...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/hhsieh/opt/anaconda3/envs/astroconda/lib/python3.7/site-packages/astropy/table/table.py:3409: FutureWarning: elementwise == comparison failed and returning scalar instead; this will raise an error or perform elementwise comparison in the future.\n",
      "  result = self.as_array() == other\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Query finished. [astroquery.utils.tap.core]\n",
      "\n",
      "2023-11-27 09:03:03 - 304 Gaia sources found.\n",
      "Retrieving refcat sources within 7.5 arcmin of RA=357.622670,Dec=3.201060\n",
      "['/users/hhsieh/astro/tools/refcat', '357.622670', '3.201060', '-dir', '/volumes/Fantom12a/refcat_catalog/00_m_16_binary,/volumes/Fantom12a/refcat_catalog/16_m_17_binary,/volumes/Fantom12a/refcat_catalog/17_m_18_binary,/volumes/Fantom12a/refcat_catalog/18_m_19_binary,/volumes/Fantom12a/refcat_catalog/19_m_20_binary', '-rad', '0.125000', '-var', 'RA,Dec,g,dg,r,dr,i,di,z,dz', '-nohdr']\n",
      "Refcat catalog sources written to /volumes/Fantom12a/BackupData/Data2023/ut20231118_ftn/refcat_source_list_2023RN3_F65_20231118.txt.\n",
      "67 refcat catalog sources retrieved.\n",
      "2023-11-27 09:03:03 - Generating Gaia finder chart for 2023 RN3 for F65, 2023-11-18 00:00:00.000...\n",
      "2023-Nov-18 04:35 2023-Nov-18 11:00\n",
      "2023-11-27 09:03:12 - >>> Finder chart written to /volumes/Fantom12a/BackupData/Data2023/ut20231118_ftn/finder_chart_2023RN3_F65_20231118_gaia.pdf.\n",
      "2023-11-27 09:03:12 - Generating DSS finder chart for 2023 RN3 for F65, 2023-11-18 00:00:00.000...\n",
      "2023-11-27 09:03:12 - http://archive.stsci.edu/cgi-bin/dss_search?v=poss2ukstu&r=357.62267&d=3.20106&e=J2000&h=12.0&w=12.0&f=gif\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADrCAYAAACICmHVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACN3ElEQVR4nO29a4xl13Xfuc699a7uZjfZbIukSclWFIl5OSITa2RStJBkLDiyEE8ESP4QIAIzUoDRfAkCeAzDTmJPJoIzQRJ54JmEIkXKMCbxYPzBDmCPFDiC9bAwiKKRIAfD6GVbL1Jkk+xXdVXdetz50P6d+p3/3bf6loeOR5reQKGqzj1n77XX87/W3mffbjqd1q12q91q/99voz9uAm61W+1WW6zdMtZb7Vb7Nmm3jPVWu9W+TdotY73VbrVvk3bLWG+1W+3bpN0y1lvtVvs2aQsZa9d17/mjJuQ/V/tOmkvVd9Z8bs3l+LZoZP2OYWJ9Z82l6jtrPrfmckxbOu7DP/AO76mqBzc2Nqbnz5+v8+fPH9th13U1nU7733n9Zq3ruqqqhZ5t3XszutbX1+vBBx+c/r/dDMIc/3O2Fh/W19frgQcemJmM701a5/HT982T40lpvJmMuN51XXMu+fwfNd8Xme8iejdPLjnOxYsX6+LFizWZTOrg4OA/VNVj0+n0seYziwhiY2Nj+txzz9XS0lLt7+/XeDyuqqqDg4Oqqjo8PKzR6EaQXlpaqq7ran9/v7quq8lkUuPxuEajUXVdV6PRqA4PD/ufrutqaWmpDg4O+n6tJKPRqKbTaR0eHvb0jMfjnlHup+u6Go/Htbe3NxiP5w8ODurNb35zffKTn+xph+7pdNqPPxqNajKZVNd1dXh4WOPxuLquq4ODg5k50OgHQRwcHPRjwzee4zpC5zPPj/kzBvODbj5/5JFH6pOf/ORAcQ4PD2tpaakfC/5AP/zzXMyPvb29Wlpa6nl4eHg4MBaeQQfG43E/lo2cZy1T8wo9gUcPPfRQfeITnxiMS3/+oY/9/f1aWloa9Amt/O85mxbuQRb8dF1Xy8vLtb+/X3t7ezUej3udYizkBR9ME/088sgj9YlPfGIgr/F4PNAL+HF4eFiPPPJIfeYznznWCx0bWWlEUwhj0hgZk4AR/o3iwXgmlUqytLTUM4fG5DAYGyjMNjMtTBhkwVRVPfroo33fzIW/YSB0LS0t9XNDKel/f3+/rly5UrfddlstLS0NjP3w8LAXOJ+trKwM+seBtfjH50mbeUhfjz76aHVdV3t7e72DXF5eHjzDGDYgz9fKhgEjF+bj++xw7OBaxs2zlrcdCTp0eHhY73rXuwZ9YxTWn3Te8Mjy5m/zD5rgATRbV3zN9DE211q6Y52cTqf16KOP9uONx+Pa39/v9Twd66LIZaHI+sADD0w/8YlPDCIBk84fTwwlahkzRDsyOsplsxARoBmNx0Lw9NXyno7E9O2WkdS0oxD7+/u1s7NTa2trvQAZzw6C+c+DlClkK42jQlXV7u5uLS8vD5wQf3/1q1+t973vffWWt7ylfvRHf7SPPObB/v5+7zQsR0dZPsu5o7Qtedm4bQA28kQF9L28vNw7Wu7FsaRS82xV9Y499cay5W87PPpPtGBHzPyNcLiW9MNXOzn00wjPegUdBwcH/dgPP/zwTSPrQgUmKzmdQ4QFiqG0YDH/8+zh4eFAWSwImvOnVlTAUwOlGNtwen9/fxAlYHDCNBhp+hI22bhHo1Ftbm72jsieHppReD+X/cAbaPGY165dq89+9rO9MkObx4Hn9HPq1KlBepHRw3Smk4EPVioMHR7yWeqG/09HyThLS0u1srJSy8vLvV6k4VvG8NARERkZKSB/ImEiFvqAjtQ737O3t1eTyWTmunllftmAMwgwL8sAB8p9y8vLM8FnXlsIBltwqXAMatjjUG/oZwZwr2FdKh5/E1HsBAzJ+ayqBgpto2Zce9B0NgmxDHHM9JxT8oBxeN5IBAcDnfRpyMe1tbW1On/+/EyeZkVBOe+77776x//4H9f6+vqATkfszJk9F/Mhc0Xzwo7L83ZuzxwNx+1wrFP0nXpl2lp/O2+2U7BO5Twmk8lg7qbDc8mUwTwEthudTCaTnu929MvLywPdYjw+N0pcBOEuvCkiPYmjF0RU1cxEMUCan3cRI8cx0+3RDD/8P8/gucwQlMaeLnNpG7/p9XM0w1rPOWnMeaTn9/1WRMZfXl6ue+65pzd4+sAwgMTcv7m5OYDNlo8N3hE/FcaR3det0MvLyzMyafVHxOOa6d3a2qoXXnihz7UdlXg+UxpklGgJ3tmIiKruD5k7mlk+6FrqBPNKvjHOxsbGgHep/1xzf+igab5ZW9hYIWBlZWUwIXs4hOyqqWGTYXQqg++3ASNk/mYMGGvBIgiYhIDMZKKcDduG6DEcxVKRDRd5zvkxNGT1tiUY6MOQjQgwrjQY5t9yJK2ihRXI8uAzG4v7zIgAneaJ+ZsoCjqdRh0eHta/+Bf/on78x3+8rl271o+H7sAT+Lm3t9fPLZ21aTPiMRznejpo35PpkCMp/7uizDPWXSNK6wJ64MDltCPTtHltIRjMhK0A9njptbnmIkVCBE+WyNcq/LT+d1Rzgu/Pq4ZC9DOuMNs4MGYYj5KR53hZykaZsNI8o1mJrZxeBkApnUeSVhg+WvmtMFZG92mnaoRkx2kewBs7GjuuyWQyUHTP389XVR9dGQdD/0t/6S/Vn/pTf6o2NjYGtNmpOXpl9LTx+DPzO/XX0dS66ehpI8vPrF+ZftlG0laMLqy/BwcHg4LfzdrCxmqldu7nKMB9RCz+Pzw87JUmPVc+7yjTYir5DxO/ePFife1rX6vXve51tbq62tOJ4/ByiqEbY6WSJBLw8kvm2J6fldQOw3mr/7fzsGNrKd/y8nL/eRoW/KBfoCYG5jHp307QBRyaI3fV0bIaraVcXkMcj8c1mUx63q6srPT0IJMHHnhgsE5qnbAz8tqmdcpRzuN6jtBlROSxuNYyRut8wmEHEPQk+3CRC/niCHBgVTcq/EYJx7UTbeRPCMBkgamJ2XM5A6NIr1dVffRwTtBiqmFYVdX29nb97u/+bt8nURY6nCNnpCb3ShhlWO2Cgg1xb29vEKlcKLLy0+/S0tJg2YU5+bfzuq7rBkqexS/3Q+Q5ODjoiyjw33UBaGL9mD4wBFc+jYK88cFKlwjBiu3aBX0bJjvaGkZ2XddHbpTcDqa1ccb8gBY70Fa9xQbivtLBW5fScNPRwTv4ZD11n+4nazbHtRPBYEK/meRJwXQml0xO75Vro94oYIOfB1Om02ndfffd9ba3va1fr3PlzoLwDqKqI++8u7vb041CImycUCKAhOYpZFd/LUDooU+ikT2wx2Iu8CXRDONnjl91FG1QqsyVzEtDbGh3OmNYSLSzgtkR2gHBP0dMQ/2cK8twGK71DpqM6hJxcQ0kwnztZFroxE4nYbzps45W1WB92GgkDZl0JHcAHhe8Wm3hpRsYmAvDVj57klznauV1Jja9siMEnyWE5hoMTY+YSy80FM5GgkC51wrs9UrGtyLyG0UDBtJ3wkZXAs03OzCMy/m8HRl5tA3V/WThwoUURygXdFBKmpXaubmhnnXEtQhHpHRoLZkYJTiaWuHTSTmtsdPy3zbSdCams7WG6rnZCSZ6M33QbF3hc/Qu57pogelE1WAz1R4FBrmkbkOjgpwtIRgG461zNp4sKJgeFDAFZgXe29sb5MBANzPLSpxwheuecxp2Cj8dUVYkDSk9H8/DkK6q+hoAf2eUsUOFFvjrSM39zM/OxQqaeT18QHG9Pp7OyRHfaIl+7SSzWGcURRTLZTOcnSNvq0jkYiJzIw1yn3aUjGldgX6vk1uPQWVOIzxP0Mvly5d7vudGk3lt4R1MGCPCN1Rq5QAQZ4VwVdPLJhg0jF9ZWanV1dWqqrlQ1ArR8pwudHEduoFjGQ29YcKw2BXILDYlLY58RgoZnfO5rCrm74S9WQ9g47nhOPehVFZgGy9KmtGwJVvnqtaFTC/sEL2Uhx55PuYbfZoX0I7+tJYDLX/qCU5DDE2ZYzpAp2SmOdMM5pFOw7UCy9Py2t/fr2vXrtUzzzxT29vbvYNJ9NVqC8NgY217C/K+zAc8Oe4lysFsVwO9mSGjlY2jVWDIQg/XrATekJ1M9v0ZQRIuJiS0kK5fvz5QJkcre9vc2QLP4LGjJFFxMpnUysrKoHoKFE6etdCAjSuXbdLAebal2IZyzgWzyo1xmFc845yX+2x8rQidBofuwA/vTrJzzijvqnXWEVIHjJ7Mh0QdyMEokf7tAODx6dOna2Njo65cuVKnTp2qtbW1NLlmWziyGuLlJPC+uSPJod3CtgdyXgQTcpmAcQxXXGQwBDJsdDTMPM3OxnQgAEddnrXnpy/3aRjuOXsvKvNjDHvxFv+yX+dy9MeOnITDzq0MVaG7lcq4YRjOq1w38CYOKzD8gi7fZ8STxu15t/jhdMV645QCaGtjpCWaMc0e1/NPXUZepuW44l1Cf2hbWVmpM2fO1MrKykyOPK8tXA3OQk5W+jKMQ6xfXXPxwIWGqhpU8NKLG+Kw5ueIhgLkvekdocVRmDnQvPiPF86cwhHCsJnc3Ibk9UHmjbc1r2iGyJnT8JnzL57JIp4Nx3O04uQOLWjnN2iE5jFMH3PDwOzo8o0r5rC7uztAWru7u70epExyPl03+xpgogjDfTt4O0n6zbV8fjv9yeBieaVBGybPQznMoeUk57UTbze0MjkvTThjhtjje7MDnj9zLBurDdOMM1MdrVOJiWzuBzq8dML9QJqEy0RD98cc3F/SyrwwAkdGC4m//QYJkSmN1wULw9hUGsYzT7nfDtM0uH/LhDlyLR0i9Dn/svxN+wsvvFDb29u9QU8mk/r1X//1+uhHPzrYT1x1VH214ru/jMAtI2vxO/XEDnw0GtWVK1fqwx/+cF2+fLmPfvThCO3I6aIk/aRcnJ7Rp5HkcW1hYzXTEACTdB7hz02gXyXjmo3SP2Ya/bogQvN4FhQtBZmFFjPV//u1PvfpooINk906GLFhtPNKC9uRGlqpTDqPn06ntbW1NVhjXFlZ6Z/L6Jpv+NiJpidPp4fipDxdqIMPOAvPJx0vemNlxVm7cDUej+vP/tk/W/fff/8AZUAb2/qc8tBswPSVxse8HZ3tmCwL5vWtb32rPvzhD9fFixdn+GWdQT9c08n83Q7PKIk+X3YYjGBgXMLMhFGtKuDBwcGguutCxXQ6Hey1hCkWLsqU+YC9tqOFCyeGTNyTiuGING8Rfd56nPtLiAjEt2IY4lYNX5FzRO66G2cTGcKxw8d5p5eUWjAto5tzTStOpiKOYPDGKCGN3UtljJsnPvAmken7nu/5noHjT6RkHaoa1i8sHxtu0umXHLIwlnr3qle9qn72Z3+2NjY2+mDDurzzd/PAubgDFOkQjpSVDxf+FmkLGysMslFAWK59IQQaSyWZw9k75RofzcqY+Zr7dWQg8njN114OgWDMKK/f+Gdsz8GOwblQQmYckaNhVfWnS2xubg6E5Kq4kQX9WiE8d+hoKY+joB0VcjSvHf0cfdhk7ihkvnRdN4jkjOf91M57GQf6cyOG0xnzvrV2n4gC5JApQxqwdYG5tApvy8vLdfbs2eq6rq80Y4CtpaVciswX3+eliIvmq1UngMFAiRb0ZPKulHkRmclnRcwMZBJmLM1C9A8GbMaZDsNYQ25HSa594xvf6A8eM2T1vtzMScxs32+eOIrv7u7W5cuXm5HM0cvw30LNaAg97DumT/PD0J57TCPRwmvBhuGZ7xkCcy/VceRKrm2Uw2/D30RHWZSxHjhV8Wd+a8lG6r6N0Lxxx7m45WnnaDTJ/87L6YPTL6A1nQL8sNO3fizSTrQ32JA2IZKT+qoh/Mooa2Xmx0bRwvFmniPkyspK790Tgu3t7Q2KNOkU3Pcdd9zRv7zt+dipOArbATg65JIIBoHAzp07N2N0zqG8Bu1xXCW18WTe6X5cFzD/jRDs7e1wnPs52vu5jMbIJ3nroh28bSEaK28WgzxPBwmvO+cz6Gk6bUN9R9aMsMwHlORAYD1LmtA95g2NPo/LgeJlz1kdBTN042URsuGPmWCYlvmelcNGYG/pxfQs9jgv2dnZGQg+DSJzWbz95uZm33dGRUdjL/YjTPrnGRCGHdbGxsaMgSaPoc3KTPOcGdt0OQdGEaqOihitJYoWpLcR+97JZDLjNKz8hoeeB/RRvEp4ah4ZNhLRWmkWNBgGp+Nx30ZnNkr4C8Jw9AYlZBHNy4zM204BO8gqr9GGVxy8vHVcWxgGQ7zL6I4eVcONy3zOvXikzKmm0ykHHA+S/zRWJ+OO7M4LnS9kxIFZdhSp6FaAr3/96/Ubv/Ebg7mh2AjRimfHhFC899SRd39/v770pS/V888/P0ObnRF9ZQpg2kEUNgA7Ditj9o1MfB/3OjoyF5/CAe+tD/TtKO2c2iilZUimPXPfdODIxPR7rp5DSz/NQzvk5HGmWa5226m5b/PfDszRF5r9CubN2rHG2nXde7qu+/TTTz9db3rTm+rxxx8fKA3hPgUCYUzQimqjMuRIA6MPwzZ7RCupYY6rfBmFW9Cx6mj/Kp4Qxb5w4cIAZvEbZ+JFfXtHCzjvY9xXvOIVdfvttw8gqfPFeQjDTtLzI7/k79aeavdn3noM98vvg4ODwekQhveONk5RHKWdhnhM5plwfGlpqdbW1vr94e7PMNv8gJ6WDpqX89I1/k7jNJLx5+ahHaGbeQ/dVIWh9QMf+EA98sgj9YUvfKG6rvt0d8x35Cx8bvBv//ZvD5iBR3MhyRUwfttj51prQpusSqbn5dmE0WYMjDUdrSjr3UcZFQ3tEoamozk4uHE0h3M9w6GEOdm3+0Rxc6nDBp+IJauqViLoMc9skIaQ1gNDT/O8hVqMLGj8z2aCdF7wwfNxju1+oPvw8LBfw6ZKbcSVeTRj2IG03mlOZ+aaSs4RXjm9wgEkVLe+eQ6mg2vj8bgeeuihl+dEfoQKsS7keDJcS++aRQf6TAVKJc6+W0JNJePzquEJjK0Ig8OBjnxp2Xmh6fbYGAvP2MujYERM+rKBJYwyrR6HaGklYq7eeJHLMtlXwmV4YOgGPZkje63Y/WRubeXMvNq0OqL6hX8XoJCl1+yhN43HOpO6kagAfUnj5t7UmbzXY+WSmSN713X9OcS8fG7U4XFv1hYyVkcjGJfEuswPAYa/XMM7++gNC97Gaa+ZEc0MNqOSGWn0Ho/Ps08btCOux8xijSNtC/bb2cEPb7zA0Px8IoOEfW42JheT8Po4J8aDBzgAQ1rz1LlgOkg7s4xWjG0ZZ+5qp2XepzOjb6MLy7BlgHaCnpOvpVx53jRyzXrtWkULIXI/p5fA/0zj0JFWMGu1havB9tqOrvnVBzAbRpixKZiEQmYen/kwMEcHmiNDRuCq2W2RWZCyYPxjeOe+GDON0p9bEPxvY2D83CiPQXmOhvyutntulomdDGPDwzQYGyKysozo08UjrrciDzLm4DqPk7DQTsd9ebmNNGJlZaWvJ5gvRkbwG6NIY7JcuNe0JKIYj8d9REzHb0fkgpONDpSQ+bh5YloXia4LV4MzslipUQRH04RaLjBkP1lsMNQxI9JbWpENmTIfqDraCpZLHAg9TzzwnLmnaviaVlXN5KOGZAjbQprHK+hprYPCx4T4NmJHSy9xMS5j0I/7czSi33nQ1rJ1xBuPx/1XY3gunr//t6PzeBhpfsazXg7xMlXqkg2ntTRl50EQyuW2qqPNQHaG8NKfQSNVc7YUptPN8TOnP66dyFjNTJfteaPFxFUNt6NlZKUgs7u72xPrV8wQBjCNCaNcrurZgHITPLSnAfJM5kWONhZ85iHOv8yLvGeesGieG89mbmrFywJPwmr6SsiJQuX3zGTxKGl1ldkKn6lA8snKn3PPyJLNDt2FrhZ0dkHPjtiobR5UT8hs52AUY0Rmh5y65XlTrzBfnF5YvgnF57UTvSKXipUw158haKAA1xwdr1y5Uv/8n//z+sY3vjGYbJ5EkFDLE/T90OiXCBCQo7xhWMJDmGv0QI6dSpc5tz24YVZCYOZFGT/zQ+eWVsZMITJiORLn2qnRiGWV87H8nJ5wv+VtJ9E6qsW0MWfn4QnD8/6UYx4iYAcCfUQ7v9vs+Vv+mf64PxudEYGdF59nNLecbOh8noHE9x/XFjZW71NNDG7FoTlPNMRhUl3X1ebmZv3QD/1Q/81njAFDDecSPqLQwE2YhlBXV1d7Iftd1oTiNK/1wmTGpEiQVXD+5l48useyAthpuQqdjsg8gH4jAUNj07C3t9ffm3lmHhnjHJ7m6GEHZBk4qtnAvK7r/DWNxZVnQ3DTlujFY6cjgFeWtc9ggk6nFshlXiAxfM6W1Xbzwg6ce+wQMHw7YRvuzdqJtxu6uOKcxxHAibMngMKg3Ovr6/X6179+ECXtECyYVnRlXEdNb1Z3jmcGMz7XbVwZ/TKPszLzmYWa8NCRP19uoLXyT+dBdgL0y5sgKXzntI6y9G00lBHDzoh7rPRErkQrqWw2BuTOvFJ33DKHTtSGg0mUc7OcNQ3SaYEdBH3yv+eDYbVqL4bNmWOT1hgKWx52UjdrJ9puCAN8Aj8EOL+ZlxPwTFZvzRx7Uz9rI2kd38KzGCUvgsNsw5acD9ECRWh5e+deKSjPIyGOS/MJ25hjwkfTZKXknv39/bp+/fqM0RhSJjzz9dZXYbT4iZIaKnocpwGWixUeZ1VVtbW11Y9tqNmKLOaJ6QAB0Yejpg3B/LeDSlhvGqhBoD/wwXrd4q2RCHTxnPu2U+feP5LImgQ5ojmStryE8w0TaqjF3y6IMC5rhc61TEcWP7ys0mKioaphbNXwC5sYyxHctDuKpWMyFOq6rofi5L/0TTEpPTPNKYFz4o2NjcEpiXZE5oNlwJimdzo9+gLhU6dODe6FpkQhyMlw3Ubj+/h7eXm5Njc3B44rkYzTmjRKIxbrU2st2kgkYaqN02jKeuLob9o4/8uOxI5wPD56vzojs2E/8/D66yJt4cjKIBZcCzYgDITs3NKCtrd2dMKLodCZA3CPvWnCHJS+lUOY5szXbPBp2C4eGM6aF/DIiswYS0tLfTm/avgN4xZ2js19GWVzX3RGfDsooxlHI2j81V/91Xr/+9/f5/7mlfvLPs0f5mD+JaKAD8i1xYeqGxGYOeNYkLkdp51e5qmuRCOrg4ODwSkbNiD6MV+5BpzPfD4RD9CfwOLvUkJmRgDQ3cqNW23hHUyspUGkcxALjIknvMmIQR+GTobSCcvsNTO6eXE8i0n2vBkd8zrzsQfPKGxvnfPjmg3ZlVpgtr14zqWFAng21y2JEtD69NNP1/r6er3yla+s5eXlmRez6Qc+YygPPvhgve51r5tZN3REMB02vEQ4WRy0XthJZ+Q1Orrtttt6vmaakfxvITJHLesozye0Np2OzNxnmSFn+ndKgN7nmO7P6/I2/kXaQnclTIDYqiGUbeVrhjP2smZYfj8rz7awvJWnn4Q8uml2ocKRyTmxafLaZAshML6fY15ZqTSPUtlSiVNxLUBHMPPQhlFV9dJLL9U/+2f/rH7jN36jp50oYahpGohMr3zlK+v+++/vx/PccGaJZNJR2ilCJ2N404bn5rmgJ3Yq6cAyuufn1oM0aEfnqhoEhaqj0ytNO30nmjGUdhGypauWewtJpt4e10603ZCWOVq+bZOeOaGyDScVueooClkg+/v7/dlKmU/aY2fBoeuOtttVDQ8mY0wrjyGe6bYi8Dn5pyGsFR3aDYXyfOQ8u8dGMhrdODkxaeS6c66zZ8/WT//0T9fq6mqzAEfxBNl5rdKOkWcy108jNC/NK8smaWwZNHPe29ur1dXVXq+gc2lpqX/p3XThfLwd0vqWjp4x2VU2mUxmxoNG5MpvIzLzlXEceZ2L8xw8Zg4OJC0jn9dOvHRDoSOZ7khpD8VkyCNhshXJSbwN1h4ID+w8tAUVuZfPfY1oY0ZlXgNdZrQjrh1IQv2EjNCWRRDDnqx+O9pQ1DD9Ng4jkJWVlbr33nsHDgVH4F1PCeksK/r1/C2nquEXc9Ev69BVw3OPMQDDZzsd5m9dMbSGvtxxRWNJxHmtDde6g6P00qH5ZyfM83YI9OUggmE7kicKYXzrih0lvF6kneibzxFEVhOd/Jtw5wWOGC3c7uf9d0+oFqNpzi+y0UfmuHh77vFRJVxzRMx8iHE9T8/bBu1nGJ/58lm+p+rczUbjHxSPZRArgZXKfEhDTyfidc+E5fRnGu3sjGaMclLO3O+0xzw3z1wwsjxsjIa1LfnBV3hkFJSQ1DKzXHPuGRhyO+14PO7fbzbac8ByLeCPZJ2VCdvrZ5TkM+cLMNA/VjyXyR3NcgLsILKh5VsreT8Mp28YlQdYmS6ikQXjKI6XNBJIZW05G55NSAkEzMIajeplOgpD2VR8RzZ4XDXc320k1NrxxDP5vB1Izs2O1wrPM/4ibubhF/StP9YXj530IrN5SIRxrQNGM66rtPhs/bBsLV+/eODjW+nH+ujg5gi7SDvxucEMiJHZi6Wg8jd/2yt60zrMsZEnFOYZwzxPOD3faDQarH05qiG4hE6G4Fkk8HOmG0gELYa2tHlrgvawFjJK5gKN50yk8hq2HSl8SwXK6G2+sZDfyv3T0Kzc3G9aPZ6jsPkJrXaSNh7oc35smujbMst6hw3L49Ovc0nrncfI9MuOyxHZ67A4goS75jnyXaSdaAdTVQ3esHGU8dpXbreDKDPPuaR/rBCutHmCVBgNnXJjgw3KUNzRz17SEd3KgIAtZCuGq+Iem7lnrmKHAL0gBubgeRB5nFujPDk3V17TmFIZHWU9N0M38yrnYGO0DGiJRky35+LxrCeegxEM41k3PK/kj6+hO7u7uwMacmeRg4sdqOlL5JAGjHPyl6G1kFsLgc1rJ4LBh4dH3y3KpLww7RwjvUVCvNzRhCB8H326TO7talZUnjejXGnFIFqVZO5xtMQJ8dsCs/JZ2YhyjtYZwfy/IZ1pYRw7ukQEnl/uInK+57GgwXwwpPQGDjtj+vD7uZYZcrPh2zF6bBycDdlr7YlmrB9WahuoN58wP+sNkTqhtftzGmJd9rzsoFwwTUditJVypk+C2knawuusEOKccWtrq65cuTLwsGaYI5cV0opp3A+THO0yMsOEqvZ+3XmL0c4x7e1a63mO0uw6Smfk3ygyCuI9oYxnXsJHF+Y8frZUOPO4JaOkMSHcvEjB/4yZbz15ic7GZ57DT5bZXKHlc8NsR2o7Mt5zdpUZvubSoWnGWafOoCM2WENoaLNxOlWxA2esTF8Mw3Hc6Yh4lmY9v1lbyFgdbexRNzY26syZMwPI1Hp2Xpg3M73g7GZGcJ/3hVrINshkoB1BVfvLj1MYFpiLUnwOfVaa9LBVRzkj91fNvjHE3A4Pj17pM0RzvzZG+GjH5nGyH8vSdQfLC/qQt43NDsu0ew44KDsQO6NUfN9nQ7TDo1+KVG4U6JChkYzzZs6n5nmvApiPmS7Rcqukq9000+olJc/FBwBA0yLtRDAYg/JkUFDvl4XRGE9WD3MN04JMQcyDLhiQDdV/5489IAzF6xMVe6aMhiV8R6OE63yO8qdTcx/+2/mqP2+lCtPp8KsYM6LYSRhqVQ2XHZLHVkZ/7vln/pb92BEm/+1UshprY/DfWUCrqsH+Wa/j2kE5JyQaoidE5eyXz23YpjmdVBqw5Q0KcH98QbTTATtBaLHMj2sLH/L9yCOP1Ac/+MGZTQSGA/aWZlZWQQ01bHxmlmGyjbEnvBHFbTT20q7IWcFsVBiGt+gRIVyw4Bl7VZ5vGeC8PMw5sZ1XRiAUwMjD82xFUuZsiGlDQWl42yaLgp4zdJgm08G1hNjwgrFSSatmvyKx1Xe+Usd9NmTP36c1Igsv+VkX/GPU5YjPuDyX6JJ7WugEWlxzSKT3i7/4i/XII4/U008/Xd3Ldcj3b/3Wbw2El0SmVyMftMDMVJfyDWcPD4++tcvK4N9mUkI5GMH/WfJ3X9DvyFs1fBvHdOzv79eXv/zluueee+rMmTMDpXKByc4gi1bcCw/oPzc0mCZXPh1d3VciAyufI9h0Oq2dnZ2Z0+4tI/PQURkExaK/o2xWZ803oxs7ZXTFEdmydQQn6voLofntL3xCjxy9kU0GDebqWoaXW+z4cF5sQ8z5OTp6funAdnd3a29vrzY3Nwf0vfGNb7zpId8n3sifEQJvY6hjpbdQUI40Xv/dUppksIVuT2nI52KEPSjRZWVlZeZlBOcWFgZzWl5ernvuuafW19dnlNrRP5/nXlckHZlNu1OHNLZWa/HI/6NYLVTigp2dmvnP/8jNqMHXHT14xrTYKfvrEb0LKQs7LR3MglXXdf3Rp1llbvGVZ8yLlHPe4z3VzAO9NyLx6kOuqXpO7A+3QS/SFi4wpbFx3euAMNSRwcaJ92OiCCZ3ETEO/bcKFIapydyqo/wa4RrepkHZ89O/C19eE+aLkDFez8+RyI7HPHNFlaUQ89WKZT4zxzR0e20cpiG/+zZUdO5oJER13xEvc29HNqMfR5LMw7785S/XZz7zmcFpln7307oGzYkeMmXix07WlWZooiZhHTNa8N+uPdiRJNJI5wB9hvpOO+DL8vJyra2tDXJVdqndrJ3o5XN7TRPsQgie09Ehc5Wu6wYvYjtqJ0SxB7fBJeMMcSx4znF1rp15UfaRO2Lcp4s59OUqo2FqenMrE3yyN08vbIVrQc5UPK+9gnI83+zfMjWkNF9MA/f6YABgnVMPPz+d3iiOfeQjH6kPfehDtb29PUBP6EfqmHXNEcv65/vSsI2i4LNptDOAhjzUgB+nKEY8BAH4Bm+MQDPHdUS1Hi4SXf9Qx7rgjZMgMy+Z4xyq6kalzJGWZw1bueaiAkZlZqJorf2XZmw6gXQ+GUUyZ6MvR7aWcnsx3n3zmRWAHAi6EvL676wGVx29eZL3m27TwhjwP9+5ve222/p5e3mLvhPact2IwI4ZtPKOd7yjrl+/XqdOnRrIIdFMGkIuH2X/7iOXbYwM0sBTvywjxrfT7rrhVkrzCHRoNGSdyKBCAQ+6F4XBJ9ob3D8UJWoYREU0j9iwEG1oMD8JztK+m6OWoQR9G84YymZUtbNBEFVHa6IwN2HQaHRjrzFFFu7xkgJGYQVzlHY+l7Qz3u7ubq2trQ0iRcspwF/zlL/JtbzJA/ovX75ct99++6Bv7kUGdjx52Lo33sMrWqKI8Xhc586d6x2BN6JYfnZg+f4q93uJ0DKERv62MfieRCToM5/ndk3mkbk6/eayGQ2em05DZ+vkosZ64jOYmES+O4gQYKYN2UrEb0/e3tCQ2M+4f5Rjf39/cFau+3U+MZlMZr7vxfnidHqjQuq1OITCPA1pq4ZfWcH99pRe+LY3ZieQlxNaEYavS7S3zo0NdlpOD9wwLDugpaWluvPOOweRDH7wufngNMdG5ZzMesKY0GxHbiRkB4ZRTCaT+rVf+7X6/Oc/P1iDd78Z+ZlXvlY3nU5n8kHLw3K23By9mSeOz/qewSfTjEQAdgRGi372uHbil88RCJMzcbnuZwVHqdLbVR0t4xhuGtcbPpqReHfDc7/ZD5MpUgE5DU9RFuBZVpDt9VAS753FEGyoWSSaTo++cIm52CnxvJcEeMbwzHUApxiZD6VzpPnz8XhcW1tbtbu7W+fOnau1tbWaTm9s8+u6rj8kPZXZcnORKiOyn3FE8/z5nPwdWb3wwguDwiF0uzCV87ShpI5evXq1L+rAcxfd7ATtQB000gagI6N2K1Ki+5maWC6LtBMZq5cB7D0NOywYmGU45Uhp+JcLxzbyg4MbX1jstTLgtiOoPRaCIDLYmNi9ZIjlcd0P8/M8XWk2pGNe6QhQxmyOSER8j1FV/et9zDlzcOfwXDesbOWS0LWxsVFra2s9DYeHh/1RJz4Khzk5+s+D5Gks1hn3l1CauW9sbNSjjz46iMx2MPA86wY4jpZzWV9fnynkWY+NXhizangIt6M4iCuN03pP85iWh1OVVrrXaifKWZ0r+BpEoTiHh4cz76m2CEoowLOZx7pqluPZUx9X4HBumQ4kNzI4krmIRvQj6jGOHdHBwUFdv369JpNJra+v99+sBq/sDPgbR+TP6RtYlnmp6WoVyLKfdCCGlRlFrPh2iFZYVz0tW6cjaRy0hH2JMlLO1gPkaQTlcdIZw3uca6YLwOfUF6+FtiKi9aOFdsznVuR3lG/tMW61hXNWEwGx5HOOTs5fIdLvD9JH/lTNHkWCEOydUX6M3AyiDxs+NHiThJWRZ7g3dxxBkzdxw+C8hkJevXq1/u2//bd15cqVfj6JKuAD/HKuZIeVSmujdg7MQWnMg2jlZ+1oDXHTUOGXobsjow0KI3CUcvQxnIef6Xigp2r4vTmJ0vjcDt4OtWp4Cr4hagv5GWIb5Tm1gR9evjGfElFlfcO/uX9/f7+vETjw3awtZNKGIIYCVnYUxzthYJ6FgbDMeLy3ixWGwVZYQxwrDFDI8NfK43VWMykhjCMkObaFbmXkmq9XVZ09e7b+yl/5K/0GCs8VHljhd3Z2ZjYzcK+V0dArkYQPR/McsuDEXFFOf+Gv52Ie5VwtSyCpjR4dyS8jtoMAameksoxSBzAkOw0bkoNGfkY/1EdslP7Mc3QRMVGCnca8sVIX+Xt1dXVQFF30pIiFD/lu/e28gs+SOAiGwTZ2mOsKH80ejTHsMLw0k5DMjXWthOFWNHtqihzb29u1srLSF1pQLLa2GY7m9rGVlZVaX1/v50OEakFSF85Meyot4yf/svDmZZpUPssl+zJvrly5UpPJpM6fPz+gy4bF3y6C8Zmr9tDhzSDcRzOqaTkIy8zOJOEvY9lhoBOGq4wHrxzpDV/npW52bNZP627yK3WQJcJ5DrXVFobBfl/QRQ4bmIk1LM79nPZ+CXctHJoZnYWohIXuGwVAEM5XnM9gpBjE9vZ2vfjii7W9vd1He/YSQ2seoeIcxBCNeSU0S0hsuhyl0rElpLWyOHdHEbwuaX7iJNMpgJAczUyzK/d8+bVlyu/nn39+EPmtC5kWWMeAiEYhVUeOx6iEudGywGRddUS3sfG/9QdeppOzDJAbPHY+nXNzOuL/WzD/uHaiHUwwyFHQRmD4RUvvw/9Ap6qjYkHCBgSbkNUeCuUhj0Qg9sA+GjILPAgHoS8tLdWpU6f6CqIV0HNqeW+aDdARx/mNI7MjpOlmrK47eunaypbOsuWoUnEyGsE/nt/d3a2NjY3a3Nycif4gG563slkB9/b26uzZs/1pEXZWfJ7ph/mG7JkfUYjlHfjGveiB+duC1EZjyQsbr98qMm8cbe1IE52kHNOxQlPLYR3XTlQNtoL5GsbSKo8nJMjfMLUFFejb+YWZlspvxWlFJjPUCT9Q1ss5hpZWCKID9GW1b97L5BZK5kAeM52gaXAF1MUy85klF9OW8M5KY7gK/XauOAobLGcvQ4Mh7mg06lOApME6ZJ5lnofhkaakQbd46uIO1WAbg+nJaI/8jeD8TKYX1j/n1dYpy3oetDYPF2kn3sjvnSj2SBldHXEhzL+tRDRHjKqjyqBzHQoTMBLoZE9sYdrQvE5XNXQYCNlKAApgjnYCjkieL0rG8wnVuM599AldrkBX1eD/5Kn7NX9xWBiXC1etKAZfHf38DPwzTywz89pLX4arHivl7ihkmaytrQ1e+LCzNTKoOkrTaMcdv2OURT82KNPQ0iXkYqO0TeBweM6pgqvhOcbN2sLG6jK2PYQhbcur0LruaHMCQs37EgInnK0avhJGy2WUljLb8xo6Zx6cx9OgnIbsjvaJJgyTDKWYTyoEfGUO8NMCt+AxKl9r7dM18qDh5A4Pj5ZV8p1e5uF529gyEps+8yLHzmjY4u88pbUxpR7wHI7WjteGyL2tJR/TlnUDmlHYPL00Amo90yoGtnR1XjvRGUyZi5oYMxqlyAiSZ6iakY7Uhrcu3BiSMqaXHngm6TEj2FbHded1redaSp2R3nMwNEKBW3CMz4HVjtj0A7/SedCPjd8QMqPBPAjZinTzjIUxXVuAR+7TY6Xz5nNHFztH0+dIDj+NUrKOkc7Pc0qj828jEdcPEilCA6/bjUajwcF2OYZRolOW1OV0Cse1ExlrDmhm5rumEOBXh3JNysziOStbvpJkWGxBudqM0vJjJXJRy7kRz7UEa6djAfse//Yc0zBpvp4GZ6Mxz1pIhGcY3xtCEvYm9EvY5wIXxmI6XcX1+6WpaDZW6PL87FSckxrdVA1fu8RI7ZgTzbiegLEZQpte69Rx6YGdCuNS16B/65rHRcfdr2FwOv2XFQbTMvIRNdPDWaFTqPa2hpe+n6gJ03JTQ0ZYJpyC8FipWLnpIeeQSmQI7khsOiwI5mXFt0Pzcyiijcb8yWJGQtaW8plfrVMNUCgrjtebzc+u6/ojPz3WdDrceuk0wf1n9EzZWSbwOl9tTJ3J6MXfeXpJGiJ/J09w+obBiTocFa1r8Ng1Ausku/gIGnZGdsTHtRMZa3rn3OUDLDDD7GFS4cyAhCK09Hz2/hk1TGeL2Ya55JHeAtczJYTl8n0KwTSaZvgEL7wWZ6XJXDvhl/uDl87XUVDPORFB1fCEC4zRz3scy4nryNobWEgPvPY6zzDTSJ0bprHBuzyuxYgB+hy9qqq++MUv1ic+8Ymen+adI3tVzQQA1xvcv/XVzXxoFeB8aoffYzXfU0ePawsbq+GkPbGZwX1VR/t57eUcYVqQmOftAMx0xvOzeEgXOUwHTIUmC8aQr/X+pfMnCzLzRQvdXtcRwV7cRmOjMspwIQR+5JycSmBEGVEceZCHeZGOsapmjJ/7cw78TXqSCMV90m/u2Ep98nPppBPhmFaeXV9fr7vuuquXJc9i7Lmu3HLSjJ9OzyjHKZWr/4k0GT/fCU6ZLWKsJ/rmcytfCx5YeFw3wRlZsw8YmUUZCw5G+TNyH5TB/Wd0Mj2ZLziyp9CqhrkiNHKdv1kesqHbg3od1Hyhn8xrcjxvLOEZG70dmaMkEK9VJWWN2fzb2dmpjY2NGWOpqsF6uOWVjtdQ3ddcjbeDaZ3GUHWEPvgW+OwT2XVdV/fdd98ApdA/0NPIxUZrPrvibT3lc+RovrQgvj+3Q89Imo54Xls4sjqytSKqvQREOmKhhDDWzxvqtITtvDG/PqHlxcwEFzRgqBnp3JK5pZGbHqKcva4rmgmTDd2oiCPojFz2yjbaltNifK8bJwQ1v4w66Lu1hMUJF37P1by2cbo4Bb/5P9/DTXhtozF/Umk9d+vfvLVny7qlR3mgWQsNZV9ca73KZkdv5wQa8/VEVuhAC2K32rHG2sWJ/I8//vigWubBHdq5Zq9sqNDKW+2xcvcJzPTmh3xLP+GaaaJln3Y6XXd0FIvzo4Rrhqr07dfbPBeE5iNcHAkNlQ2xnD5wLaO7vXkinKoaVFpt5EY9NpLWTiH3lZEMw/ZcjW5aaIXnqG9YDi20kvSkTrV+wy/mCn/SGUOL9dnowDLzSof7t8MyrTQbO+P5/67r6oknnqg3velNL9+J/A8++OD04x//eD/Z8Xg888VDVv58TS4nkUsLuYRiOJdwMXMOj8H/HE3SiqAJJxEGdCVU8nydS+/v7898/WVGVcM8aMi3Z5iH82BHaDsTO5rM6V3kYd7OS91Hfm2mldDKlLmZldmwlei+uro6MOqcD3y0zngulqX5yb1OxZibIW9W0pF37jyzg8iA4DnlvZYJkRYZJJ9ccPKcuGbedF1XDz300E1P5D/xSRGtgQwrUG4znxxgaWmp3zABIz15mqNQ1XB3iCM0UcHCtGEm1LKhZr7K344+9GFF9zyTbudu7ieNyvMxX9PIbUQez0oB3YZTjrg5V8YCmXjceYpuZcXYnT503dFpF/v7N75ewydMMHdDfP/fqvJWDSvluYRmOrOQZDm03oxybSBlZh4bHZh/djrc7z0B7s9yT4e+SKB0W/h9VitYVuDSQ/kzGOXCSL4gbaahIJPJZLA9MZ0BTHSURbFYn7WgHNm+/vWv15133lmbm5sDJc+KpP82LxwZW+9zml4jAh8rkvmqI0zCPdORr5gZxmVOZUOkyGTF4Z6qo/dSDYez6gpfbfRWOviXaYp1w4bveXqtEUNnXCs50NXr1+nczZtWIZNnKKwhS8/L6ZzlblknpLXeJGrwK6YZIBY12oUjq6MlnXsblb0LDK+a/ZrIVBwbGX3Yk/pvK7pp8me8gW9FqBquk21vb/fn8mKsKysrPU1WNvpypLKzyLlZmFx3Hgzdrdw/va7hlQ0iIZkViWtZBPQ1xnAkgB4v+5iHCfnN41a/SZ+dbmu+WVE23a0imwuHtFyaMSxOGpNfvse0ZmrAdRfr4AfySqSSBmnEx/+LtIWPdZnnJatmT9TzWmUWaxyl8WpeikH4XOOl8Cy9c7+NPKFRRigEfP/9988ouB2PkQPPGQp5zqbFwnTU49nV1dUBwuBeTpLIfcggD8uhqmYcmZ2KYR79p8OB//mOLPRzPzyx/EEnhuWGsfDGzpNtf9zbit6taO5oyrJNno/EfaYlU4jU0dTpeQ4sC2yZElnXqoaFRvPd6QQ0GtYvGlkXMmkTXlWDDcwothUHiGHlQlkSaprpaRAwzkzy5Dy2m6ugpj+9s+lJhqG0jMN96ZUzf/J49pw2CgTpqjef2cGlIregZ+bQlonz1Sx+pDd3wQvZ0ae3E1ZVb3yeo5eZDMeB/5n30ad54XnY4U2nN86pcqT0PO3AXXHl+cxL7eDhjR1jLodZzzx+QnOiLV+ElTk7/GdMj7+IwZ4IBmfhxgUUQ7zcgJ/E2kDsqex5DBsdwRjLZXl7rhYca3lCOwRermYu/u1+MkfZ39+fqQJ6bp77ZDKZ+XJd5uJTCeEh/bVQi5chckwrodEE/ydtbjhcF4K4jvLBezvmNAxveOA6zzFX6wl89ZZMO7bpdFqnTp2agez87aiVhSHrFulXGq310Dzx/zs7O7W8vNxXvBMdmQ/wxzTCRzsP76izfs5rJ9puyOCOmhlBTLCZYsITJnLmDoJGqPSZUdLeruro7B4n8a08xHlGCovr0MQbFa1CEWP6t+lwP86JzTMrqqN2GkoqEMLOSrNfyE/jh29GCH6N0dDM/MnqaRqjo8Pq6urgedOPzmBYXOda65woaIJGQ1uaIxqOz3pih56R1wGglYaZx13X1dra2uD7d6z7hvhe8kEGV69erQ9/+MP17LPPzsiDv/9IjBUC8YSGIv4sF9Jt7IYxeBivcRpeZJ6RJXznji42pDCqZjdu25saFgLTLAx/VjXcseUIbiEQjTyWFbrFL/5uwXbvs4Ze+kteGbqZ5/SVR3om8qka7uqCf96W6L5NQ6YrzNcRN/PCTC2QYaIwIzychaMVdObOIfTV/Xp8r1GbdsbKTSM5Vzvt5OfOzk595CMfqWeffXagtxmhb9ZOtDeYSaZyG6IeHh72hZSEWWk8EGwPaPhjxmbhwRAbZhpa5IkPCPXg4GDwbQF2EGx0MBMZDwFlUcLQLKFRS3AJ/+xUkt+ZehihtHJiX3fk9fxc0HPUgteOmua95eJCHI3nWjAYo4Iu08S4djLIykt83oSSuWMr58sCjunJ/zM9yPk4n6+6gQRZVzbtyIZ5Li8v12QyqbNnz9ZP/dRP9Xut0/kvYqhVJ4isrog5OsAYl64hwpGH+y3kFJYjnpncijJVQxji8dKwPV5uBmAubFP0PP185tg+pqXqaP2T+dBawrTDys9TqZKWluIZWvrtmKzkJrTLqG/j99d2YDg+eiZzPXTElc78bQdseTtCOSWw42NOFKYc4Uy35WQ5Z3DhHhfP0sGmEzCfHaU9hmkwP2677bbBKZvpSBZpC0fW3B9ZVYMSPh7TCpwK4LXXlgEmExCEN024aOBobidgWOTI6ciYQnM0cFRJ78nfhp9umfulsXktMJXDhpTKTbNDNF2tvjG6dAaWz3g8HhxCno6GZxjj8PCwj3DJT+7BKbaiXUZUGx3GwEkLGfUzJ/f/5kH2yThZYXf/Wb9YWloanAOWqxvzgpBTFFCc+7K+ek6JrFpt4cjqtxVauRzwyt4FxngN1MKy8beKMxZYruW5uOD8L9dbzVD+T6gzD7ZWzb5N5GiS0dYRwv3a+C1UK0cuvHu+vt/RlCWGzNVd/Mo6AP1zb8oo5Wt5mP6Wk3RVNw01I4mjFs8kfxO9oAuG7/kFykZ+fhb6XYswzbkO6ggO/52X2rmaX/klZNZd05PXF4HCCxsrXspQwFXbVhSoGr607GMvqoYv7To/geFmjhmAIF0UyAhnr2W6XDWeYcZoeCSpFctKAu0I0caOsG183Ac9pjWXLbIKbsUxH+BTohN/KbILMXZyjI88kY1lZOVv8Z//c3kNnbDTA77aIdjZOAd1xMlckTHdHJnSyPnc98zbTunInsffernJ/MDwLc/kh4+KMV/sWBKxzmsLGWvid8ODVjsOTnlva+YGQCcY7GUNGGovT5/27K0D1ZxT+ZorooY29pRW2Iya5g3X7IBoRGc7O6MHDI8+bCg2bnbxtGA78zNdjONIaBrgV86TKGPkkGjC1WsbimXL1k/DfOS6tLQ02NnEc1ZuV7+Ziw3ZBUAcPn25XpLwmN/kkM6Todf8xYl4vEQg1nnSJI8H701vC20e1xbeyG9cbYblTg/uJbqkh/QmgKr5B405GmTOBDNozi3nRdQWNE4FqxruN/Y99s4Wqg3LRmFvnAJy3sn8QBKMb3r53I6Ia9w3mUzq+vXr9YUvfKFGo1G99rWv7fMlHzyXUcgOIh1svjSA/BxdnKehC8zXjsR8N/9au3wYJx2P+ej83Kgh0w6e9z1eC83AkcbuuVTVAOXYgN0v9Ls/02Q+naQtHFm9WO8CihekzRRPqGp44rs9jIs/hmX27DTu414zk99WCHtA32eFtXCrjjZoeKnBggTKmDYvhxBR7KDsZLa3t2fyJgybQgQ0JoqwU0mDgb6rV6/2J+llZZKUxPyEj8zTsjCN9GNU5SqxUY/lDM0u0PDb69+eryOxIaKfM3JIR4iTspzQP29zNGoxFHeAwaG41mBZmCdVw7Vp664r4nYcnvfN2sKRNZXQsKUVxl2MgPls67Mns5cxzud/nvfWLCsLk/fbKVbOVjR2VM9x8HwJf6DVkSLvydzEVUv46LnCP/OTpSU7vixkMZcslJw+fbp+8Ad/cJCHWikNvWwArQqtHZOXY/xyNgUev4BgRaXYUjVc+jPkNp/SQeVyIPNlzdWpArUI61QaOfcasdlw7PTTiVmHEjEmOrCe2llaDr5/0Ui78Fs34HsmnQKAuMT1ZpaNAe9tRcyChidvI8/JQoshjA/XgmZ7VF5GMCwxYshqXRY7EkrisX0ciBfW4SHz9mkaCf9tEKY5vb9TB5bRiHKpsNDj8QzXk/eOHIb4ptlK6KizvLw8OBGjBU2dQ9tI4OHBwUHt7u72349rfviLsnguC4zQ6GKiHQT8wIAz2Bi+po47AJnHdgrmbToqGn+3aG+1havBdOoBTbAZYQ9mrzSZTGbyIMNi959e1kUmL6f4zQUzjOKFaaw6OojZMMeFFubgn6oafFM13rA1Z2i0IOxNv/GNb9TP//zP19e//vUB3K8avn+ZUTPTD34zJlHMUTUV0ZV5lAl+p8IwJ49BczTe3d0dQON0Lj7GxcpvWfK5ofR0Oh3szyZFchS3Q7WjTgRi5+ftkjjtFk2Ovm7WHxu9obufs555HEfddBTz2kLGijexl8rBrMh4eRhsWNqKvEAa/jfkrRouCyREoU+eMT0tz27aq2ZfLTODrcB8C3oqg+m0ornY5PzvS1/6Uv3rf/2v6/d///cH42aBxRA7nYiV1/9Dcyu98Dg0cjjP03zI9VLz3/kjRSwbiwuMdkJAVvTJMjHszPGIqMjA0a3FA+eInrON2qsN1gc7HgcNaEkDNyrkmqO6I7+Dxze/+c26du3ay2usDAJTHTk9Wec7FmR6zMzFMkdLqABjMmeE+SiVFZK/MydMVMCY/oY15kWEctRCCWF+evSEjvQDNH3DG95Qv/Irv1IPPPDAgK/0ZeNOh5BKhjFY6V2M8quKptdevsVPQ/fkkx1KwlfL0g6APvkb/o3HRwfv5Romzt7GUjU8k9q0m+fmGUbjYg7jHR4e9k4Gnhv2pmw8X58qYqdi/XUdI417Z2en/uk//af1yU9+8uXNWXOC9oJMhMn6VajMgbxZwHCyBVWZkM+9zXXYNHCYDJ186W+ObwEzPxTHRu7ohjDyuBf4kC+S06fn2nVdra+v17333tuPn/QTQSxkDC/f/DD0gh4bhYtjqWw2KBsiP/kiRAt5ZDFonjyYQ26TtJHagbjgBC0EBBui5ecoav7AAy/XQF9GQl+3/jqKu3bhdCf1mc+h184Put797nfXHXfcMSjEHdcWPikio6EFy+AYVu4Q8rlKGflspEz22rVr9Uu/9Ev1/ve/v3Z3d/vPW1DJk7cgnCvlEkBLCPRBnpSKBbOdu9qbpgFhdDYSeMTYPJN5Jv3wjCOjlwb4sRGbJjsZHx3r+zIqZjRpOTeaEQD/GzpCt9ce/dt6kHm5dY6lmFZKYHlW1UyUszxS1k4j0nEkD+Cz55cOhrn6XqdcdmBra2v12te+tu68885atB1rrJ0O+X744Yfr8ccfHwjYRglzEnbCTEc4fluoCffuvvvuuvfeewfwms8yB7HS0ofzKNOFp0/PakHCWPqnsGHHYjjsiOOI2jNZgjV/oMdFH5o9uHNfomhuQEin58hnr+/KJDJ04cVKbcRg3vs3vPW8HRXzWSMSeOMcHN46wuVz5l9CdfeL0doZ+DN45/w20Za3RGYa1goc/sy8n4duPvjBD9bDDz/88h3y/cADD0x/67d+azCYFczKb8WAWCtH5nLcQ8s1VoRnz2qImP0m9MqIl7DFNGZl0fdnkYFIyP0J+6EjFdufW9GcS9u5cc1GAD9SAbi3F243u+bYai04nY6WoqF3WbVk6GiSvE7nk7z2M1Z2H7BuZ5L3tebknJL7mIf5b95bn/f29ma+gNuR2YGGvhNdmkfcDz9BEi/bId8ZDRMuVQ2PSjED2JWTxmFGkY958T0N3Izi8yyfHx4ebcI39EqFSA/cWrC2I2KNuZV70Zerji1DgnYrXPI4jXdeTuWiiCOwN4LYAB1lEpK6ekyERJEY03lYQkHmluvc0JL8z5wQB5G1EI/hVQR40DLOFpIw2krHzfxtnOb1dHpj+Ygx+Syhsg0zC3mt/BcH4n4XaSc65Nt43HDI+agjgqNKvv7WMhjG4vOEeu4z+7cRtu7NHMIePL1vXkfhvBhvBfXiummtqkGhypHf92c08fw9r+RdFkyAte7HOaBzWSs7POMev6CfEJt7eOfUcrNcTBN8tLFiULu7u4MjPI0grBdWeOZtJ5QQ2ctLriq3ICt9Wn7T6fC7exKxWQ7Ja/62XWRkt/xazqfVFn753C8T40nsZTyYo6g9Nc2ey8+lwHyfx7Hi+rmMyNBipdra2qqVlZV+PkTJK1eu1Orqav9lSy1onVGvBXcyotrYcjHd/VnRPH//XzU0bPMPfviMZRezQB3m58HB0TcfeL65ZDKdTvtthY6QNiju4X8Xb6pmN7r7t+G3c0x4mjIAxkLr1tZWnTp1qo+CLnyxfGgEVnW0H9tyy3zUDgue+n70mnHtpBJaGzVZbzLAHdcWNlYY5A3uEGPFTIjn7/I0dHPOANEJEVvLAVZoO4s8ZoO/zYiDg4O6fPlynTp1qlZXVwcvN29ubvbzaUG6NHqYzNxyrQwaEr7ZwAz9MCZDeT9XVTPKntHaCp9ev+W5Myrb2cAvwzd0IFOJlKN55LTG/PWxrKbTz2Q6YN7hiCgAuShVNdwyCo2mnXvSAOnDkdE8zEBgROnUwkXBTANoWRW/WTvRGUwpTCtZy6gSVprZ3sOKB8zSP4zK3JX+yLNSgWAiv/lseXm57r777jp79uwMFAXeJxM9T9qXv/zl+r3f+72B1/RzjlKTyWSwDzkNlHsNBemDPpmTl4LMT/jIMzZ2G27Kg3urhoU9fsM7RxjTZeX0pv/Mw9KZ8b9RlPPGFopwKmFHuby8XHfcccfg2w4cGHw4nSvpdirI3/Pnb8spo7z3VacTzj7gpfP0RGE3ayfaGwwzvMG8RVwuZ8BwvGAqWytnMxP4zDkD11trcB4/81tHY3tsr1/6fFgvxyBg/0/D2PI91SzSGN45h3eByv36HvqBLyjeZDKZWTZCUQ8OhmusLdosp0wpmK+NxDDYOZdl6h/Dajthv1rGGOaXI5kdhXUpdyb5/OTd3d367d/+7bp48eLgVAcbHX1mapcI7vBwuFOM9Xbr7XR6tJ+Zvki57MSyFrFoO/GJ/OQmCG91dbW/x2uZ9rRmpqOyBelrAwLjtSv6MQOtYOPxeJCHJW05jp2LoZxhvwsCo9Govud7vmcQFein6ugLuRCIldDz5VlHLvM4va6LKzZKPkPZMic0nPfXMLbWe9NgoQnlz+KK5+bCll/zs0OwUzZ8ZCz6S4eG/IyueGd3Y2Nj8Jyd1XQ6ratXr9bW1ladP39+kH5Bj4ulpjON1ajNqaCdGr/9Kmmisqzgn6QtbKxWcEMGe2GuHR4eznwhEUzIyEZyjhPIiOtKnI3N4/O3CyJmhqug2XJOdiKem+GZFZz/7Tk9fhpeGq9pSpSR80WJHKU9/6rhppSq2ZeeE9lYcRxxjERshNCTio9ztFEhP6rHiYpslMjJ4yV0tlFDj+fiMatuvHzxlre8ZfC8IbmdMI7fQcbQ++rVq/XCCy/UPffcM7MzzQ4SB2DI7rGSxszZj2sLnxThfCSTcCCXIQkJNwTacyNwmiGvr9tInde40mimWZDe3G8HkUJPD+do5R87nOl0OgO9QRN53S9E21CyCGTF6LqjDRfmB17dtBoWVw0RkHPXnFcaSo7v6I3xO+IkpMwCj3XETiAr4lXDQ8bs0KuqP2htPB4P8vrV1dXB+8GGxpYtyMm8p5lXVTVIYdz29/fr/e9/f73tbW+rz3zmMzM6yjj7+/uDF0KMGuCJ6UM+rSDSagvvDcb4ksiqaip15peGSVnWR6HzKysY01DVxmUm2Gm04LShtD25nUBGdJqdhY0wnYfnhAAsEL/O53vyt+fk62kg8AUYmpGf+Zn3lpOrro7E0+nR4ddGBJ53rt9aKf3b91bd+IInjs5JB+BUYnd3d6AHjEOzQfs3zsD/2yly3RG5JT/rT9d19d3f/d31+te/vu68884BTQmL88gZO8tEU2kzN2sL72Ay4/gfbJ6wJb129mHjIwJUDXf+pBd3hOUZDKVV/YQZmRMnhHZE99wSHTBG9m8jzhwuBYFiWKFcrXSf8DHzckcCG4IjiZ1KwvDkq7/31LLjWTtcxiCXR2bwkHl7Wcl83N/fr6985St1++231x133DGgK3M5K7VTg+QBRsJ9loOXnqAhl2RyXZd7LMuqqr/xN/5G/fW//tdrc3OzmYcDf+100sEnT9O53qwtvIPJgoGphg2OpMlcTx6FH4/HfV5rBXKlzIpHX+RNjD8P5tkIzTDTbwXMCl06KHtdmp/hf8/Xxp1e1oqXO18S+iaUatGUeblzSjtGfhxJLTdHdMZO47NyO8+zTuDgnPuNx+P63u/93h5lmM/MDT3DiNA5pzLQTWphQ7cjZQ6ej52vnX7OOWsHVHUNeb3RI522I675SJ/eU+A05rh2opfP6Xhvb68vXTOQBZ1RzkqSipRQNpXGcMfLJV6aIFeAOR7Xhp/5QUZKCzpPnqDKm/e2PGcaRPLD63+mx8qXz2bkN++c86VRw0e/PeJ7bcS+5u2hVDhtsPSRlXSPmbSkg/WcjZ4SRqfDsvxGo6MipefIva0zmOzM/dloNOp3sFkuCe8tl1w6spwzJ00ZoleLwuCFjdXRBxhnIiAcJvh7cJJ451stGMp90+lRrut7UE5ywIQ+jmbQ3lJIC8O5S2vRmuIG/TK2FRGhu6jkCGsDZj3OczRfbbDeyM7Y8BueGHGgMPbmyC03X1y6dKl+93d/t65fv94badd1/SYD+vQJg1ZkQ1WcOH1vb28P+A2veLnDkdSRhpbICf5lncGR0XWQVj6c/3seRnPWjbwHZ5sw1vqVBTfT7XtOkrMubKy7u7v17LPPDkJ7JvhmlJXbb6lwv6Nk5odmfr410YKA3I9ibm9v17Vr1wY7h/xsRiv6xvCJ0vby3J9FqXxh3ApEv17OsGNg7vmiu6MMQrZ3T8NPdOJ7E2kY3r300kv1D//hP6y3vvWt9aEPfWhgXNBr/vA7jQMZ84XDdmZ2VNYTrrfya8s6n7NTp6DTyqurht/g4L24LQRSNSw4pm4gLzsAO1rLKuVo2Vte0N9CDa22cDX4U5/6VL33ve+tr33ta4Ocz/mMPUkqOUrjHTGOoOPxuNbW1mYqZBlVU0kzj8Q4Xnzxxbp69Wo/JgIzDR4r4TH9IwBHKsZIOkwrz1Lh9llUOCb6ZEzvvjFKsHIzF1/nWfOLPo0QLBsM7Ktf/Wq9+OKLdfny5RnF514bmxWTz4gynHTI58DTXtnEb/rLnWCMa4ecY5t/lqWRiaGukZd/O8pyX64m5AYGpwQORFlrcVCxPKwbtEyl5rWFq8H33Xdf/diP/VidP39+xgOmYdoAzASE4Khiw7dx0A+K4HcjrXQwl7Nqu66rjY2NwRsnLmbgkc3ELFykwgKLXSBIaOeCy/PPP1+rq6t15syZwW4v88znMLc8q5XVcJiCBjQ5stkwUwkyDaiquu222+p973tffe5zn6s3v/nNg61ydr6ulPpUCOTso25SJ3JM+m5t6Eg9asmmpeTQYIfrseGdHaqN3JA16wI2ePhsGRg1GUJngdEpxd7eXq2urjaDzXFt4R1Mr3rVq+ree+8dVMCcnzn/dC7GZ8lImOlqYcJDmARjiB6GV1euXKkzZ84MjJLFchtlq6KJwHnW/2fhxMKm2UgODw/r+vXr9dnPfrbe97731R133FE/+ZM/Wa95zWv6udJcJHN0sZfGwZkPOCWjE0caxgGG2hk4Otr4XvOa19R9993Xv5/K9cxN4XvCOhdwzCcjB2SRzskvD9hJWPG9LIODQj8MPenH/Mv0wPzxfR4/UzrrC/psRGh9TJ7NW3LzdszUjePawsbKIPO8ZuYfTLy1zmljyeqYIYudgKuTGM3y8nKdO3du4J34vLVcBGNbMMiez387gvldUXtN5n/16tX6pV/6pfroRz9aXdfVn/tzf67+1t/6W3XmzJkBjE0Fga4UnGm0IeEk/ZyLYvDLhsN8rMBcw7kZsuf9Vsx0pN4owVjemmh47RwynRP9Gx04MqajTJic+XVVDZYH7cAYC5lad61nfo+V+oRrK7aBVhqYhojzQrZ5WMBxbeECkw0sjw6pGkIShJJe6Nq1a/Xcc8/1HsrFmPQy9qgtjG+YlBAk8y1a5n9ueP/MHWxcPvnQVWh4cvbs2XrjG99YFy5cqHvvvbf+4l/8i4OlABeTeMbzzHwpFdI/mQ5gLPbyNkgMy68iJurgc/OZcWyk8J2+GRdltqJnCmLjYv4YueeavxmLcQ173XcqfyIanBe8Ze7mue9jrsBXLw868PhvFyidIlnW9G/neLO2cGRlcryO5WhopTDctLecTqe1trY22ADgfhPKcC2rwWYgeaQ9Os1GZCibMIe/vZjvEw+S1nQUhtorKyv1wz/8w3XhwoW6cOFCvfa1rx28zMDbQJl/t+C1UUjrVbJ0btCbEM2KmGOYF/b4hoqWi19+z2jgtMhOM1cCLCsbRCuHpQ/QhDdIZA5pJJYooJWPul/ri43Qx/g42qNLoBvXM6DJDso8B4k4lVm0nWhThL1mev3d3d1emSeTyaCoAmHLy8u1ubk5UwW2Mjk/cF4Dg70mx72ZY6UzQRhGAWZUOh3usZBhNn0ZliKE7e3t+uVf/uX62Mc+VnfddVdtbGz0dFD8gke5TGAU4msoqQ3bNLlIYnjsfugriyVpHN40AV/4bf7aIL0xJavsli1zTR47FeC5lJuRnJ09fDQ/rJt2tllfafWVcDoNLtOqRBe+16lhK9IbabUQXav9oYw1YdNoNKpLly71JW0zPAWURQt7QRdP0oh9f+67dSRwDuExfY+V0wpo5UUAXk6ZTCZ9nxY0Qrl48WL9/M//fP3Lf/kv6z/+x/84UzwbMH40/CZ0GkaFcTryWClynkSRLOYZ8hpl+Hfm7zZuaLGTyOhlpwtvvbRF5BmNRoPD1/mMuZr/1iPDcxfSqobIyakJaAmDAHW4D9Phl8nTiSSv7fTQbTvfjOLMCR2DrgwaN2sLL91kvpJe+cKFC03Iw/+eqIXnlkxyXy2mmTboSg+fBkFux/18bvjuopahlOGijYrf58+fr7/zd/5OXbx4sb7v+76vpy3zS4/HNWiDb4afHreqBs+3lMh8T5SB7BzlcyuljSRzvp2dnZn104TGiVRwsEZFmSKl0+IZ35/FISMBB4ZMLyzTdHROpzgUzryc5/h9mojp915j640/cz8g0nTmrXbsId/djdPB37O+vv7g/fffX48++mj9zb/5N/msHwSPkWesWknN3MyJWvAHbG/Pc3h4OIi+LqXbw6GkCQPtZNzwwiiUldtw2zDIaYAVite71tbWBpALhUh4TvTOeVgJod08ssIzbguO0Zc33DMXeLy1tVWrq6uDLZTmS9XRMsTv/M7v1H333VenTp0a8NT8YF7Oz9i+6Co2USbnav5kpdZzt6NMA7QRpUPLkxpprWifS4lcS5lglPAvjZLmYlpV1ZNPPllPPvlkPf3007W9vf0fquqx6XT6WNMeF7Ho17/+9dNPfepTA4/BRF3FQxkcmbJAhNCdn9GS4TD7+vXr9bGPfay+//u/v+68884ZgdmoJpNJf9SMFcNjWRksFAskkUPLu2aRyLSk54YO+sTJ+cBuRwNHYhQIXhpiGd55ucYIwgpn+pGHaUKeRgWeuxf0jV5sOAmnU7aONFWzX1VieqDXjsN6mI7QMjNNPHdwMDyilXmyhOPrrZza9EKTA8S8tMb3mybuWeRE/oVPirAnIZJ6st5IbiNJgbP84UVzb/5OyMgzr3/96+u2227rmZZRrupIoVuQNaMZkYXnnGcYDjp3ozkf8/1WRiuzl1TSoH2vlcNCH4/H/aaFVgEmYZfnZr5ajq0oZCfadV2/7OSlJ78mhvJbJlWzJ3/AE2+SSCOgTyMb6xxz413aRG0tp+9qrSO+YatlkXsC3Bw9ablRJtMM08LYrWJoGvi8dqK3brxlsO9AjGTi9s4WqCGYlTrXn3KP7NLSUt19990DQ+Tz3AXjsbPKZk9sBbXggZT8beFwH/R6fcz9eGwLwkpIjmTY5DyJH3+ekcJKbiPe398fnMaQDs15GXTZMNgkYVla7qnIpsdFOxsqzTDQxR30xFHISm3jNW2m246P+7e2turKlSsDZ2AnNc9pWSfMo0SUhteWE81OIvXCiHSRtrCxwjB+W4huRAET6vthqr3hPPgJc3NfpyGVvymOPgw3ablUghGZaTb+9KLmgSNoelJfs/fMyJpOgPtQ2oRJNjI7OvMYXnrdD2doqDrPmVrhErLzYzRiZwB9/LTqFnZWnoPpt4EeB6Gt+Bhawtjd3d36J//kn9Tf/bt/t/+GcfPWyDCrzEmj7+V36geOMuVinjuHtgwWaQtvivAGBHucTKxNkCGoGUDhgZK8jfbg4MZJeEAdJn/p0qU6f/78AGrxuZltJTft5Ie54T0LWF6KMNPtEZkbStmKADYQL2llVTLhP31ZqTw3G6E3kVs5vM3QeWoasfnlPNIOI6uuWbMgMuCgKCSlcV25cqXW19drdXW1GXnok+Y6g2FqRm34Y3jtdOc1r3nNTJGJce2AGD9hvZ1SBhI/74hsJ+aIja7ke7q2p+PaiY4izWURM9bGlQL3NzsbKjgqmcnea9p1XV28eLH+3b/7d/XX/tpf68+JTdpgMM/TMlpnBLfQUHobSkZRcnYbHr/x7pnDQIfnm96afCyjMtdyrlXDb3RjDl5+4P6EuyiYDZz5WCkTvrYW/vndMjR05j/9p/9Uv/zLv1yvfvWr653vfGdtbm4Ooqo3fTjCWU5GcqbBtNiwmePb3/722t/f778exfNjZcHy9+de5fBck292vkYl1m2/Ww1vF/3G817ei9xEtLSnoUrL62iOJDAXoXkxmlPrWNqg0MQ9GIO91YULF+pHfuRHan19vQmpEi4lxOa3oQf305c3OfgsY0cS95e5qJeafN2KkP3lvVXD5RL+zgJEVo7dB4qTkNvRdx5djlZVs1sAExGkPlhxuX97e7ueeOKJeuyxx+rcuXP1Z/7Mn6nv+77v6w0gEYbHJ6K76MVSlzdZZNpjp7y8vNzrWkLadADMwc/bBnwtDc9zp7lGwn2JBhNRHNcWrgY7mjLghz70ofrCF74wUA4+J5fc398fFBAuX77cT8BVPVeEvW7Jz5kzZwZCM4NhVu5CakU1R2zPx1APB2JPTeRzDmLFTGfgnLdV4UQBfJ/pTOU1dMpiioslWXTCWfp6zjEjQ8sJQJMdR6tAZd6MRqNaXV2tH/zBH6xXvepV9dBDD9Vdd901dxOFmyv2aYQey4bF2ImsbJjcm9tW4bvPFsutro70lpfn4hpAi2Z+53r+IgZ7oo38VUff6dJ1Xb3zne+s9fX1ZrGkBVPG43HdfvvtNZ1OB9EWj2mPzDOORIaiVUdvXnj83IrIM0T+XOfyb/fvtVEz2nApldgK4zwuI07rWnrl9OzOGc0/f2ZP7qhnozACyghih0Yz9M/cHL5l2uC2urpab3nLW+pP/+k/XWfOnKkLFy4MeI1BZC3B/SMH8wV+YXD+cX3B6KAlD6c1NORnnXQU9eYSaG/JyfWHVsS2w01n1WoLG6u9GMz7ru/6rgEEMHyDaTYcFqTtAe05rQx4Ni/r2AgsLBSQSdtguRdvllVMRxRawmKac04riSOePSxOworkyrbH93h2Uo52Lv9buC2Yi5IkCnGakBHJ2wjNY9cQcAaZgkA399NH1Q2DfdWrXjVQVusRuuFmJ+ZiWVZ8QWX05cKUHX9G50QF1hHTkLLITT7wdWVlZZBvo+s4jERetqVFDLXqD/FdNy3I60g4T/m4t7UAbIEjmK985StVdeOEiqrhBov0iobHXXe0TuqcxJHLFVnuy80UVqh5jM58z3S2og2K7EjtXDEr0K1ing3dMN6f4xBpS0tL/WHejjjIyIe1W575fUWeY0uWXvc0XZkrW5EzCPh5K7PREc+b94nAbCTmvXntvltzSBk78EAb9LiKTLNjpFn34PnLvs56nCdMA3aEc+Th2ZZnpR+EubGxUaurq/29V65cGXiujAxVs9HXcNU0ZETIPI/+bcxp+DZke36UhXF8PpEhtZWCJSwrr/uqOnoX18UU59nma0Y96MlX7bz04/d5l5aW6sUXX6zf+Z3fGeR/9AFPLDs2eZh++MZzXvpDltDk+Zm/qVPMhbGNolIvnAoY6VCFddrAOD4NJQOR6xi5nOZcteu6Pl20XvOcc/p5AazVTnSsSyo5E/aSQ0ZgJmrmcw9e2wypuiH4u+++u3Z3d6vqxjGoL7zwQm1sbMzAF/dvWORKtJ1G5o5WaENFBO03ifxjz5zQyNEy97/yTObaNkT3w/2tYhZ9OAJZ+VCizMU9tqGroe6ZM2fq1a9+9UCuRh6MD/8w6itXrtTa2lovq+RbRlKnQlZq5+GgHi8J7e7uDoo0GDsIIiNWq29SN2/qZ14YuR0Lz5iO5E/CavMafsFrI9G0q1ZbOLIa+lixrVgtT59RwhO1cvEcny8vL9fq6mr/wvorXvGKgZAZz4pgz+soXFWDyOGcCzq9sJ5RP5XOBpjeFqX1S9ls/eMZ6PA1G4q9dAsG25DhP9AfeGgUY/jrsa24iTiWl5frzJkzgxqCoZufpc+u6+rSpUs1nR4d8p4Qs5U6gMJYNYAe5GiH4qhmg8IxMX/rilGe80/6suN11EMvQAO5L5r3m53PQhe/0QMv39ixZ4pxXDsRDPZAly5dqmeffXZQkcUL2lvYIPjbCp9e1Iplhc4tjL4/80znPR4jI74jhWEVv9PbtZaGLFjTZ+fBUTjmSeZKqfweG0N23z49AfoxJjvEXOPjGryzgjmKoGTQC3TkWp5qsbS0VFtbW3Xq1Kn++B7LO52sf0PL7u7u4PxjV3NtmGxBNbRvGYG3jmYk5DnmyZhZODSURdfNF56zrJxzZzpmfuT1m7UT7Q22If37f//v69/8m3/TTyAhR0ZeJmDGMBEgjK8no7hmqAjjMLrMJ5xD2dk48rpvxrfQ5nnkeYUo516eE3151xHNitsSnqFT6zPyMMP9PLTLcjG9NBfwiCKWaUJhpwPw4Pr163Xu3LmBMSStlhXjgTzW19dnvqYkHX/LoVkehpe+h35SrvDf0NhyM88cePiiqqzLpNO1jrlPO6BFomrVgjkrYdyCf+Mb31h//s//+cFySMKKqqMCgL1IQkIz05PxRBIqENVawkRZYLyVg885ntLMBUYZanqcjIjwJAsENkgigI3azsPChLeZZ1tJWjzjXUyvITvqp3PI+kLmsMzNtCU6ydx3MpnU5uZmbwxEyFbhzI6uqgYRCx0yCjIqsV5YNuwgy7l7zszL8/a3GFoe5LNpxOapawLwzgjEMjdPsy5iFHJcW8hYnRBjiKdPn+63/1mhLBifkp9JPX/TrMhureuepCNxQkuE50hiYTtieowUtnPXpCv/NpQkIqLcwNksPHhDgCEqDWV2lLBT5HkrkecN/10baMFH94FCOwI6r08nNZ1O+wPiEh7aGYE6rKC59xbZ2ZBciGpFMutfrkPbqKEnIyDzNMw3TVlU8iklriYnaklHh/7nCsMibSEYzN7fFtbGE9pjQZwLBi6auPjiCSUEBN7B7NxRYlpyA0bCFuAWnjCXnhLSpqKmAjCnfPOFe/D0foGb4ocLDgmH6BPhoggImb9taEYYNDvWVFy/Vuh8Gx46ClrZqoanOzj6ef/tpUuXZnJ+92Wet5ym9ahFm++Ht4nKTKP75P9Wjpt5Pfe5ma7UeZq3LKYjSb1NvTmuLRRZETCK5Bd4sxmW2GC8BMIknav62Zw8NNAMjx3JDW/sNGAGc3Bks9JaoBgJyMDRkvsy4nr3DP+3cmbGSw9vnhhejUaj2t3dHRgfmwoM940kcG42ZuadUbhVFIFW8zo3fGTRCyW/du1ara6uDtIh6DEtyM9r45Z9pgq5lTR3SiUE9po0dPtNl9Qx5JoIwzrtcaHJNDBPj5XIDt5OJpPa2NioRduJToqwIpjh6S1p3nOZ3s5KYejXit7+yRyBlsoJPV6Qt5fNva72+o5KpsML9oZGNlDGrRp+VUcLcjPvhOHmDePnl0ZZ6AldUw6JSHLOmXPbYdlgE4Imn0ejGwW2u+66qz8HK/NT87l1dK1572cz+tiwjdIcPc0HFwr9HnL26fqC+ZmbNS5dulTPPPPM4FmjB2/caDm+rut6Hvn6cW3ht278NxDX27dgGhNOxUyFNQyqmi3W+HNDvt3d3X4tD+ZgjDDNDLAzSIeRhkNLwWe0dK5hOGNFtpA9bysuPCP38fWEmRa6+YvSzXNgCa+zcOW04LhIalo9b2Tn5506zFNUG7qvEa089zz4znKzo0CO6czgBTyAfufK5rOvoVukbtPpjeLZz/zMz9RP/MRP1OXLl3saMx3MdWPG8PbY5MNxbSFjtWIuLS3V8vJyf04PxpmRNr1UejF7eibgz2h4X36eeeaZ+s3f/M2aTCYDQ8G72mCgJb1WQmQbX1Y5zVQ7EOclOR/o8YHWfJbRH76YJn+DuIVtWeTzNhTgtD+zfBxd6c/r3kkPBmHjTuVKh4ECP/fcc/W5z32url+/PhMBcxkOo/A1G5Bl5E0joCpaa8kJJ5c6Y8SVDtv/O8h0XVd/9a/+1frRH/3R2tzcHERV7jVqgm/zgtciUbXqD3EGkxXIi9uOdHjh9PKZk1jo+dun2DkSXLhwob7/+7+/VlZWBgerQZ/hjaMYn6eSmR5XbLOg4aUe767x887hqqp2dnYGyj9g/Bw4enh4WC+88EJdunRpcA/3mZfQbyWywSSCSSNMuunLnyVMpX+jKp7z2iwo6F/9q39V73nPe+r3f//3BwYDv60fCR/TyEx3Oo7pdPj1k45kBAPrEf16bHS5qgYvCjjfPji48YbNX/7Lf7ne9ra39VtmjSwTIWVQsBNvrUfPawvvDXZhiQlaEaw8/j9zNkei9Kb2ODCWyeCVT58+XadOnRrAxJbRowyt/Mq0J0ypqn7d0vDVSsL9LdiIsImOa2tr/fgt/rnKzP+nTp2a2Tzhoo7HdK5ux2FF9Wc2wtZ91BkyN7cjyghqfhuFra2t1Tve8Y56wxveUN/93d89d+kFvtnZQ5PXJNEHG2mrwuq5OYdkPi5w2vCNJEyjnaPnnYbqyA2N3ntcNVyD39vb6wuFi7RjI2vXde/puu7TTz/9dD3yyCP1xBNPzKzheSIZ/lOozgusIJ4gDS9tI28ZjT2X+2spkiF11dFbJDDbAuYZK0pL8IY6Zvrq6mqdPXt2Jg8y2rDjcHQ8ffp0ra2tzcDQjMQ2Xq9dJsyytwdCOgrZUI0qHJFco+iVZ04UtvHcdddd9dBDDw3OQEJuptfOIFGZeZ0GZQia71NnIOEeUjiuOXCgF4lq3F+iN/TVztq6xbij0fBb9paWlurxxx+vN73pTfX0009X13Wf7m58C0azLXQi/wMPPDD9+Mc/PiOM9Gw2Rgy3xUAzFsVxZLNXtTIlJG4VVQzlzGgEk1D34OCgJpNJraysDHZjWVCGQFZEfnuzgKNROgycAuuvjhxWBL+NlDyHZwjbymsH6Cjo9Vn44Gbl9/92QuY9crUzSGfrPNKRxvfCf8PenKedj3Wq64bfDG+eW4f8BpZllmdDGxlaXpYD0Th5lvrCHBMxonvOxenj4YcffnlO5G8xMnF2Mpn7nQfgQatubE9jT2huaMhm43S/9oCOTLmLxZ7bXg7msyziSqGFN5lMBpDa0dcQxxHBkYtmo64aQmF+QzdRzP97s4QVxYaWjnNeAQWZevE/+7Bjs1Gkwza0dCEqabPRmn6PSzSvOtplBf1EP+tUjmO6XZk1AuPVS8sxDcxo0WODTKAp13Lhsx1A2ozn7PndrC1srEzGXtDro2ZG5q9WsKygGkLRj4tLTN79ICBDCp5//vnn68qVKzORPgWAUPCwVUPGO5rwXC4ReDGf/jPqmXYrZcL6fP2q5QzTWeb9Hg+ecK0lO0dSpwaG1DRHaEcjlJjfNqzJZDLYGw7dnhPpDqjGTgkZufpLQy6WscdOY/JzLupYp7yRIdf+7YwsNzbN2Em0UJEdbaYzL0vOSjOE8CS55i8K9uBeQ8wcA0yPktogEj4yMTPJ3slR3wZnWi0k53+Ocn7Oa4RZuHB/zJPrfnWNufNcrqU6J6aZR+npUWrGcqGEZ115NDTNAko6IubgPoweXBw0zOZ/+rCj4c0U6wL0puPwUlWrb3hhWMxvV5CtQ/Dt8PDojN6Mbobr0+mNN4e2t7cHzgvepzMaj8d9ytKqJ5j+VpRN/bpZO9HSjReGGdBes2VAjq6Z0xhiWJFaCm5onV6Tz7uuqzvuuKP/Ais7izQO9hzbCAylzEhH6VzLzZyTls7LDgBD8JqrI7IVz0LPAp754ajr5xKS0gxXqcj7pQNk5P5yDMNPR07De5y5X0ZwxdyO18Uto48cM2Xizy0T3w9Ndih25ujIlStX+pfKie7p8FIu5lE6NacPqV92rDmHVlv4kG8mVDWEe/ak9pzp+V1Q6LquX8dqweD0TM7lmGBrGcQ0eNHcuRfPe23RdNpDemOB54wCVlXvWS1MOxcLlhcJXFHPZRcLzcsB9JMG5/lbse0k8tiZhMpd19X6+vpAzo7M2bd1AV4ZHrtZpp4v9NN3GlDV8NjbeRAyC1lGRoao1iOe9zzg9blz53r9RIdJd4DJ0OHTNw2bzRcbp1GA0cgiUbXqhEeR+tUeJktLL2IYaEbyN6cnuP9cFqo6OqPHORLNXh3lnE6HXyGRnp9r7p9rzAdmWqldNLDnd46S8Jm+nFOSoxHF7LWTNldczcPW2Gk8REvmYfRj2Nf6bdksLS0N8nrPxbKyQWc0zE0NLd3xXHgGx4LeIa/JZFLXr1+v22+/fTCXlAVONaO1HVbSg6Ha0F0ATF3gd+pH1gQSfpuHi7YTVYMTHhoGm9HG9SbS3/fBYrCTa8NsIAgKmtEBZXW0tcLA3IQYubbnl8MTttvrYQDMxY7IgvPYVlA7Ixsz82n9mPd2NobN5g3jZvRy7mnFNK8yvTAN/j9pQ4aZDhglJKx31O+6oxMgoYk5O9paXszF93n5z2jFAQJdsHOynCwvP2+nZB2eTo+KnIkc6M+ow44y57tIWziyQgRRwa01oCftPqy4uYeTa2aqHUAWA1KB7LUMAT0GDG5tJ7PHswAtEIy2VZ2sOtqGZ+VvoQErT+b9djRehjKdGb2c8yZdKJfXNR21zEufWwSPaGk8yBQeQVfWJrjPMnAEhEdGJ1XVvxY4Gh2dWDkej2ttba3W19f7eSf6cZHPcnXqZuicOmheoLN2gMiZKMx108k4fsZBLeHxIu1EMBjGkKe13t5Po2nlkyY4q7K5KcHGc3h4OIA1hnjOQ6tq8N6ivWBrHdCGCb026JxfFpnsQU2X6U+nQ0t+wWvTjgLCE/dtA8JoKJZl0Yf+PEYaVos/fGZDtBKaxza2vA+DSX6nHLjmzRw7Ozu1urraX08eJsyEXsvPEdVGYgNHz+1c5jlbIz1kvru7W5PJpM6ePTsYxw7exp9LS8e1E70iRyk+MXkLnnpN0tHDAnd+x31EMdbcPD59GkKgXMlsR01HulS6eZCU3Mj5XlYhswjlfnK5imdauTmv/Xk7n2Eevz0Gc7fic80Qkme9yYJr8IfmXC8jahqQoR4wFN4bpuNgW8adfLVzSF06e/ZsjcfjQRXfvLTBWj9afLcueN2f+cEj6xvGmI6U63aI169f7+fg3VauvTh1S97MawtXg/3bBmAPNS9K4fFz4d9bwWC0K5EZGbK6Np1O+83QhoWGPulYstBhg/L2NUdj+slCBLxACH7v0nDStCXPHCmc10C3Dc5Oyzmi7zFNKMh0Op05tcAtc3/zZV6khl/0yTxayyY0yzXrAU4r4IUjqPmDA3Vq5iNMbbQez458f3+//8YHpxum1flnBqfcoN91XR/MOOImEYMdfdXRm2XWp+PaiU6KSOxvw0ulsScejUb9t5k7kuX+YE/aUdKCNNTOPMJ5cDLYhui1Ss8to+9oNKpnnnmmrl27VlXV78jJPAea+HIinmf91sgjlR4e+btnHRkTcfgz89K8Sidn5+nCTsL7XHLzUpWLVImGeN7LF/6MXBMn5OiSiuzr/kFf7PS47rek7FT5H4TEc4eHhwNDZe7WMX/GnJCn+9/a2prJvZ2OseHDOgN/rQ+LRNcTbYpwOZ9mReR3S9nIC+ylPUF7TeeEjqSGIAjeJyw6ytkjWvgYdEJaRxT6XlpaqjvuuKM2NzcH0TahtE8JMAyFH+ZDFs+cszP+wcFB7ezsDN57tIJk/uv++cnTLow68l3UhPGmmYiMkWJ0djrwhus+QsVzrRpWsr0ryFVerjt/5jP6sDFY97IyjlP3PJG/55mnnjhY2Il4LuPxuD9+NfXeKQpjZ60D2rJYOa+daOmGRX0Yg0LnBK3I9rYJKXjeCpV5jKOEFcZ02chgAM/gJAxprEgeK73bdDqtzc3NAXy3wbq63erXEci5Fs0K5KiIE/I3uNm5GUbBf1foM3ragdm4/Xl+xvwMOb3EYn7ZaBk77zP/od+Gy/+Wtb+E2/1bJzBqkBPXgKkYZ0JzO6bMW+E33wJh6Mp8ec7OxPNAJq1lwVb6uEg7EQxO48RbJPytqkGOZCNNGGgDSAhlxlUNPbENNgsFZoSXaAypbUgoZhqs4R0tvb2jNTmLoyh0OX83v6wICfGSJo/LtYRx+QpXIh7ozz4Yz/xL55OFMUcyO2ju5T73Y3lb1jhiK7Z/Q3sWuIwuDF2hzTDUc0/jawWAVs0g0ZrpTD237ts55tq/4f1x7UTfIufCDd46iwqOHlYG8iErm5V4f3+/z1Vz2cNMpS97UsOSZEy+MJ1CbXk6+vA6ZBqJhe45OFoxdzfPrVWwMs1ew2T+NgIvbWTBxXkUc03IlWvWjjrmmf82jG/x2/Pw98TaufE5fXmXG4abTsEFIHji9WLPA354P7KNxsZI/zbKeQZmHUN2mXfbmVs+lod1KQ9wP66duMBkT2smeKI2SiuyoZbhLj+tqlyuQ1lZbaQIJyM3Suzvy/TGgOO8P3Rb2U0zzyWfuG4DscE5otCyiMYYKHMakYVrnjvFMNpIuOZ5J9KhgaTcrxXYSpZOxyjBcnBuazlaNjhi05SbGPjc49rQmJdzTUPoNF76SvSBY7AT8th+ngDWCjbwPxHFooZadUJjNfOM223ATNBQzhEHgWSk4lgM5x2ehMfyu4H2WPz2mxK0NDjPJQsHVbMvBzj6pkHSf8JYPrO3dWSyUtmo7bzSQGxcjAu9SRvXPQfP03w2EiEtcJRyIcx8S1lZARmLz8lBW29vMc+VlZW+qu4+cg52UDkHjqDJ51wYc0Hphz/1aP3QJ9410A8brfNOIz/XSJzuuJZhHcpA06pDHNdOtIMpPUHCKHvvhAcJw6z8FEosOFpG7ZyYPadhWUI1X7eXJiLzf0ZbG6rhD/cBkzPn8Xju1/1n7kZRI3Ol1tKV+WC5mOetKG7nlJHG/VqpLO+MWKbLRppFLSuzUY557eY6RqZOfF5VM+kVRomuOS1IZHfEvyO5ZQELB1Y1/FZCR1s7DPa8p47SbLD0v2g1+EQ5a0KOVBAmnMbHWzAUnVymB9u7bxeiUpkzMqVSOx9NxngOPOtKL9ftXGB2Qhj6RHhZBfUzPbOXlmZQRjqWzCldKMrlGHt4F9CAzc7r7SA8/3Qi0OB5wPuq2a8rgeeZI7tox99Uk7vuaGdP5pQ2lIyq6Ew+R3/ZPNdMyxIVpAP1vZ5/0uNrdhCWF/d7Ocf925iPayeKrJlvpNGS9Pskck82jSthwM2iaka5LCwlhE7DNYNQIO86ajHP0BQD8ectWJmKhCHZ85ouo5BcJ+YzLxO0hJ18s0HAm6yK+qB0yzgVMxU5IawhHZHMtDKPNKh0FolqWqgJ/pkWG7hpog+/bcS1t/6f/3VNp1VdV/XxFz5dVVU//KlH/2Dsqo88/FTPH+jZ2trqj5a1czMkt35/61vfqmeeeab+5J/8k9V1NyA+X5mRfSzSFs5ZUUIrXkIY5w7kJa7GOvLZcFsVw1b/9upAKUfGhMFWIsaEMS3IXHX0DWDQ7+IHwn7mmWfqxRdf7O9xjuqoCe3AspZR59zYJWW+Ze5jeIcRJk+5z+/Nep5pDM4l4U8uj3l7oXM439Ny4n6TyrvW7AT8Y+SQaMZjWTY29lZOPm+MhN8ptzxKpwXr4av1Zm9vrx577LF617veVV/60peqqgZHwDg9yxWLeW3hyIrAMyfNHADvV3X0vaKGRInX7b1T0Bi2/z44OOjXTt2ncxsipiOF77PhGCq3IJ9pHI1GtbW1VU899VQ99NBD9QM/8AMDwdpY/Kw9aEJjWtLKeKQPfq4Fpez5ocFG4uhtlNMycJ6zMTCu1zFRWi9TGdL7d0bbdCBASPedMDavQ2vL6JAj9Yiskv8fP/BkP/Zrfu2ddX1yUB95+EM9D4zMcKwbGxuDOdEM9V0ofOtb31qvfvWr60/8iT9Rp06dGkRd8/tliaydDvl+6KGH6sknnxzAUXt3DMlea2trq65duzZQEnuVLJwkrLH3xWM5r7OXc7SyA4E255YZLe2RvdvFfRAhNjY26m//7b9db3jDG2bSAlpWaO1MzCPP3RGfsaqGLx44CltBGIvnLQ+jHdPliAUtXre04dhgvWW0qgan0ucuLRegkJejmx0ozbDXc7aM7LQswyx48dvzS8g9nU7rcFo1EhIxj10PMBpCV53eZMB44IEH6u1vf3udPn16Ji1i3k899dTLe8j3gw8+OP3Yxz42iBRmAoTzf8JGjNoRFkZmLpGM97MYkz2kBeQNDM6fWnCKMX0ETEZoRyHm5flzD8KhX1oiBebg86cS4tLs0KCRZ3zdSpdKmLJF+Swn4G1CPdNhOSR/MzoD63KLnx2xnUJr7hwrYz5nJLfRcJ+Ld5ad9S0R3sHBQb36V3+sDqfT+r3/6n8bpE1GRB7HKcve3l5vjFtbWzUej2t1dXWgq4zjlMD6u7S09PId8j0PtsEAe1/f60hqIfk+bxZwfwiG/I3xvLvF0dTX/GPFtUIYeiU9ztc8H9NmI0UR8wxiw+qMfFVVV69eHXzHJ3ydTqeDk/8y3/TbPC78mWbGeP755+ull16qruv6va7ul/4YG4WiP6MG57Rcs+OERy247Ll7Hnb03svrSjyfMU/GSTTjA9bcsjZgnRuNRvVfvPTf1X95/acGcBT+ufCY897Y2KjTp0/XaHTjy64fe+yx+vznPz/Qh6wBpH20IPy8duJqsCFLVkZ5Dc4EtOCtCXWzAjmv8v8tuuyBXSn0WPkmSBoeCsuCut/McPR3BEUpvV/aDsI5t721lQUDsnLYqyeCSK+PYaeToq/t7e06derUwLnY+H0/dFoGlllGMJai5r0ob8NO2M/nnmNGMOfCX//61+vs2bN15syZuVEzI73TEGRmXeP6Q997e60szcJuO5CM2CASpy3veMc76tSpU/01p1WuSGdakenCvLaQsSahNkYMAEazx5dr9rhmUDLF0Irrzm1c0EKwCWkYE4VOY0kYmNW9paWlfmnA63nAV/PARpuKaWWvOopIpuPw8LBOnz7df2GTFYJ52xExhoWbOZs9N/ffdddd/e4waEkH6HlxPR2EnQL3EP2zoGe4aIMn+vm69ceOgLny//nz5wdr9xQvLXv6MR04auAndDj6/zePvHLA19RNWqJJO77xeFz33HPPQG4UWjOXNlLhO5YWaQtHViab64Mw09HI3sfw2MUd32cvj4Fn0cFKwLgthYYhPOf+ENZxisLnvu5o56JYRnDTt729XWtra30f5Knu2wUzQ/+qGkDWzDNbfLGTcFTJNzrsCGyoNF7U9gn2RBD3ZcjsvM7oA96x7k6kMe/dF04QHhDVRqMbldjks/lh2aCbXkdnHF4Wse75NU0HheQ3fSfio9nJA8lBHrkhZjQa9bWIl3WdNYXQPzw6+pInQ8J8zpHHuavhCsLDY5oBCW9cvMglAUdJmEWulAbva9BGHxgXp8nnC/E8YyXjGovfRhPch8CzkGZn4tzOBmx4y5hWNPM8awC+x0jB8+CYHZ+2kSdsOP9M+SRkduQejUY9T6DBvHNf+eKCUY1zT+jJ8RNppYOE39zzhp/7WP1PH/1yb2SGzzkn+OCx3Se8o/gELaQ8rmUg20Vh8MKbIqxUTNq508WLF+sXf/EX64tf/GJNJpOBgpiZhlr2hgm5UKiqYaXVBQKUyTDZSpsK4wiaQuMZFCPzMEPDjHRWKlpupzOcYp4JGa1w3lrInD0/P5fzR4myymv+uriT/LIDhhbkjgL7ZQrG5z5o8rfvpZPLVQEbYyIgj+H5049hvPmcPGc85ldVtb1/WFd29mt1fLRk5T69ZAgf4V+mVIw5Ho8HSznwK+fadd3gBfebtRMZqzF/5jiXLl2qX/mVX6kvfvGLPcNtsN6I0Mq1zHBHZEdzVygdYbnP0QThOHr7dTd7ScNAC9JjGFKmF3e+mXAUuh35zdOMHszRdQAXLOC/f6qqzwVBJsvLy4OjVO0QEvq65f3j8bhXqFY6YT7yvHlpp5QR3aggnZflSr+tZb/URadpRga5IgEtl6/fMKLTa0eIMNMu5uc82A445eP+qdxDj/NiB6R58nA7cc5aNTzZHeN65StfWR/4wAfq9ttvHxQ0IN5rbIanNlIbsxUz4SRM8kvaPIeA3NJxcK+VLh2FCxWGXyiFD+myIVhxMle0omVVOaF/LhkgeKKl3901pHLENKIw/SiJld2f0V86oHSo0FU1fO+4hQA851b/doJW4PyB945K6KPRB31ltbrXsb29Gq2u1ovXbyy33b7xB9sgr1+vbm1tpphkufA38/d4dsAUraxDifysUzdrJzowzRDKhkbYp2KXEQbimBxfokyfjirJbCsd99tpABHNJO73BgorrZ2Co7whpZXTc0ABTIfpzRzYc+AHelsVZCMNG6qju3MxaPSRrDzD+C70AeUzzcgo57wxHR20mE7POSOIDc16ZGhoOZk35q+hIisOXoP3HFzMTJkebm/X7739h+tbf+/H66U/MNYLp9fq+X/wE/W1d/5IjfTKpPWWH8ZMQ3W12vSaH+Z5y06OaydeZ8Wzw+wUUhaPXN5HcanSzdt4YMjhtVsrnhXYUZj+uu5oGQBFIBo6B7TCGVKlwdqZWNEtVEcJ0wP/eJ6c3hHPzoY52NkYzvO5K5N+19KOxAUc05E5vumyPGzo3OfqtZ1jyjPhZEZNj+F7M9IY+Ti6VR0hPvqAHvpIXT08PKzR+nqtP/AX6sUP/EKd39qrd7z5PXXuf/kf6sUn/+c69+73Vre21tPqOofn5yIcBmfdsdNJ5GSdzDTsuHaijfyGFumFk4lMCGjGhBNmMREXcAwjLeTMCewgKHqlwAx37WFteBZuGqchohnrrwTMfMnRxAJ25c80tarjVcPjP3MJwlHeOS6fOTVoPWvojlM1D+wg6YfnnRKZZzhDz9n5nvtJiJoR0VEpZXNwcDD4qk0Xvjwu9BL1XNm+8DM/V9Np1UuP/0K98399rLaq6vb3/Ld1599/30AGjoSJEKwXyMg8S8dkuRgxer7HtRMdReq800QnfLRCOuJQ+PDyAIzwu5YIhgnb+xt6p8eyImSkqhrCJHtx92U4bIVBieCBt8ZZ0TIH9Fgez+uYHgs6vS7neZpe+rRiuahDs4PylrxMNbKynLDUdDrn5nM768wTbZhOezwHjMmplPlCAW13d7eXL2M71TBiQBaZWkyn0zr/9//RQMe/62d+buDMGDNRFS2RTkJuBwU7SO84m+cEWm3hdVZ3bOKtyF4ITgHaEKuODMvRNKGtIZaFaqNl/P39/X5N1Iygz8xtM1J7HOizo0iDbvEBlJFOxMK1kD0HG5gjqvnPjwXrv02z7wV5VA3XWFFg1yBsyOmAzGvPx87HRs0YOB/ri/P9efmqK/NGNCsrKwMdyLk7D7QBmd6u6+riz/7kQD7P/vSPz+Tc9Es6lzWZ1KGE4PA0Uz1aOrnj2kLGauLsXazYjjYoh/MFM4C8156ehodlcqnYrVfY+H9paal/nzMLK56L85CETIawKCuf+XnGtXI6WlYNXzlzGmDl9D0YKEra8tqGn5k7eeOGn4UPNrwW+uE584D52DC4x3Ok2aEbtvrHS3AZyRz5JpNJX4yEB17aMvriM9Pmv43wDg4O6tmf/vF66QO/UOfe/d563Tev1e3vfm+99Pgv1HP/4Cd6vnjt23rkQiXj2vm4TpAGbN64ZrJIO9EZTBYeE7JyOCnPAg2MQ7lc/PGSA80L+i1IZDrM0BRO1TBvcREpK3KObo52NhArip+zU7Gi4rzS0TmV8DZHFIR7uJ65XY4FzcmndEYZeSxX7gWCe50XniErXoqHfhupnbfpc06fBaWWUeHMconP9BvNOBrS7Ax6yL+zU7uf/Q917t3vre/6mZ+rruvqjr/3j+pwOq2d/+vTtb+1VZ2cPjqTL2zYBqDFxbfUV+hP3i8CgatO+GXKbumtnctaIb0ma6iV+YWLL2yo9ykJNjiiCI2o5YiXhZz0xC6SZGTxfY4gXiZwUcfM9hyorvq677NC2lHkpgMrto0U2qHHc7XXt7E6H4Z32YxITJeNPaGdHQry4Zp57aji9MHOzGMn7M8IBl/zfCnzKdOGbm2tXvm//3pvkBj8Xf/9/1j7W1t1uLxcY8nC697wLg2MsW245gl0UrNJJ7pIW7jA5H2OhgjeAAEBzgsRivOxfnBVBJ3z2cOTL3DdWxkRtiGQcy6YBiOhy0z1nOyBDbv9upyV18/RPHcbquFtCjQNHt64sMNYNFAH1/wZfZse5uKck88dwVyIgVfuP3M0/sZhpwNiP+zy8vJgE0kLZmeRyPIxivCcjUqgLx1aVtu7rquxXgzweN0fvHzBT6JBO3jz2stsTvdSF0CSODTn6DdrCxeYYIwjmL2Oo5mZloZqz2xYZONzvsVYMAVm5CtZo9Go316Xk4cxjGvDM9RzFY/5oOwZha3kFho8Ya7em5u5PteZM5XCeZ7W+Q33m582BkfBg4OD/kvFsmjF/OGPDYb/oQkep+IxnmVlWl1Q8tgeAxnwDLmqaxI0ryS0Xi/DkHj526kF/LdcLbfxeDw4qsZO3xExg9Du7u6gWDqdTmtnZ6d+8zd/sy5dujRANK5bVB19L9TN2sIFJr/dsr29XTs7OwPF8L2GFhhWK6dESXNTuX87B646isa86uTInDQkcxnPEJQczFEBR5Le2s9zPSOxvTjzdHHI6UBGFz5z/mfhM0/DzXnV59zTbGVNI2gppYs5GRWhfV7ebF6bz/DUOtN1XW1vb8+8OMFnCRet6J53FspaaY3Xlekv6XQwsNNKyM/zNNcX4N3Ozk49++yzgyKZC4rZ783aiU83zCoeDObvrNRBmJcPgBMZUemrBV8yKvA7zzPy81U1gNXZHEWrjrawZR+mx0w3PMXAjTro3+982hm4b5QkiyTO48bj8SAPZg5dd3QSAU4hc7d8CcFystPJOdswOLbGfdKsrJZPVn9tUMiUs3jTaSBT88JpAXNw0dAGnjA+HYUdVqYi9Od3edNR5B5012IODg7qzJkz9c53vnOgNw4kzHPRnPVEBSYmhuDt9c1cMzhzHiZqo7Ch5ZqgoxQMx0BhqqEg4/MWise2QnK/ozb5seeZhRQrmgtrLQNqeeCqoyNaE+KNx+M+P7diJMxMZc00JMc2KrJs+DxlNA+mGqVkHubox//+JnhH5FxiyuiYS15ZyKR/fwMbeuEKNC88uLbRdTfqLN6xlUZqmdiBpzyhJ9MbO1dv5vdnDhQtqN9qJ94b3Ip06aEZOKFE5hcJo13cyA0VVhAXO7ICTb8taM7fNn7G5MVr77G1E3G/zGNetDDzHfGyrxZ/vVbbEqSji+eFobs/rqNUVKgxQBepbKBV1ffl79k1/xydbMTcazRh/cio52d9b8J2O4PUOfqys+H+DCZGYZ7PdDqt69ev19of7AvOpUR0xddwgrxoPhrdeME+ESeOxQ7LxtuqLrfaibYbeukCgefJDva2WVl0UQiBuorGc8AEGxiM9z1eh/WxnGY0SuB5ODIhcO9Moc80Mvc1b15pqOmomFsLmmf6wDUXQxJOA8cc5XnOhgCddgD0yf1s4+u6ri+y+JTFhNVVs8tMKLpl5QKT0VMrotl4bFAZeVzoau0qMrLjmh2pi2Hw8Ctf+Uo9++yz9dWvfrWnl/knwgIxULjzWUoElORF7puGDs/luLawsWaByOVtBIh3c+JOc7QyNDGsSljMBLJcnhPjc78FBM2Gf84nPZYNg/55BqZDL0pmR+U3Pzx/j+/3d604PON5pcKZH9Cb+ZKhL/Lh/xZK8GdVVdvb2wPH6znYUDOiZ544j14rsKGx6bSMd3d3B4jKfErH5zHt+O08jIJaGxleeumlesUrXlEvvvhiXblyZVAMhGYHIvfNFkhQSOpda8nO87FcjmvHGmsXJ/J/8IMfnCHYVUcwup7vCQMeWFAuiEDw7u5uH6nyy5Hpk7XdNC7nZlSJE6457/FnMD7zK0c6X+N+mtdxLSg7JJ7zHmp7WJydjSEhoKMJn6fjo+9MBdJgfc/a2lptbGwM+m3RUHW0JGJDrRpWWz0WBmVEwH2tNWbPz/zPSjstYa/1YjKZ9H97XNr+/n7t7u7WaDSq8+fP13333Veve93rBnSl00jnnrw03dznPNxy/8AHPlAPP/zwy3ci/wMPPDD9+Mc/PgMtIb5qCCkzJ8l80pNBsfb29vqTzCkgObqkl7b3NrTOZB8DtGFaYS0U0+xChaMLzM9XrtKQ/Q1t7te0u/Bh/thh2AHg4e0YrEzOq1r5fAv+GgpWHcE7v+ROfmuY6So7Rup0I1GD6YQ/OC07I+fS8MN8cxoASksH43qKI/nBwUGvY/w9mUzqueeeq3PnzvXP8qKAdTuRA7yctxxlfaG5kIee7+7u1ng8rje96U0vz4n8ZjQCI09yws+ErFwZgYEVOVnK9zamzBNQFN6uSY/P814qMF32yp6Ll1ysfMybnIt+vAxihaLle40oVfLK1cTMkRjXXh1FJc+xcbR4Z0PPpYKMVH7X2LUGHzZggzFyQJ52yswTnr700kt19erVun79ev/1GOZpvg5omWU0SohtPYRu00qfLQNEtnwV4+rq6sCobHBOV9x/Olia54Ve0IcLf3bWx7UTHeuCISQ0PA5imDF+jt8uDHH96tWrAwE4EvGTBaXcgdOqcObYzmUyZ/a9zMEGbiMyb7i3BYMSojGuobl5aCG7cFR1pEiOTkYa6cQwLiuPnSr/e+nMhsO4S0tLfeRxESeVzhF1NBrVqVOnam1trTd+RyNHKiMdF2aMUGzoyS/DVEPxjHRGdOfPn6+Dg4OZSnwWjTI/tz3Y0Xvs/Nu1Cp5fBN1WnXDpht+GsJ5IC7PPi5JZNOLa0tJSnT59utlPQkn+tgLmbxdrDKMMX6qGL0s7giYctWKngbjyZ8/N/c570ghoWS1MhXRl1TJBHpkqGMFYYS275AXXvY7p+5y7IRfm1EIb8GJ5eXmwtGFaWk4VmjN/9d8ew884+tJa6CqLT3lfokWcKzA29TzPHfbfk8mkXnrppbpw4cKgPvKyG6s9ER4nMXp6SeN4F4u8BMHz6e35217Mk8+11qrhd+84f0vBtpQco+d/O5eWRzVkpLG7aDQ6+mqJXKinP6cHRMl0BvAvv4bBYxuqObJNJpN+UwKfGVZbMf3btNjhpOzog3t8rw3ahmaomFB9Xv3BQQCZZtqCkabBJU1Zn0hHS0tk4nESQWaNJPese3OO9aWFfm7WTpSzki+msJm4q4GehGEeuVsauKMXOYqZDiMy1zT8MCRxAcxQMSGRlTcjQ3pSR3JDYfeNIPb39/tKJArnnAUaDZ8zz0zjRQkYg40LuUPJiur7kYWVxQ4OGpBhLg8R5bzxfJ7CGYF5jERalmdGNDtJw3P6QX4t+MvfubZpfWX1gm8KyBTIcnEeb/63UiPTYSdyeHhYZ8+ebe4hWKSdKLLy3S2GifZ+yVAbnSecxun1qTQSVz4zomXEtNBhhPcOm3b6MzPJW/jch2RnpDWMzOJC13U93PMzVUdbGOGN0YUrsVxL/lUN9wOnAVYdfeGRq7TIx/czBrzB2BkPSJdvDtkQ0gDsPJ0W2OlYVnzmyJ0Ohnk6b0a/DMUz+sML88CGbHqznuB0JZ2OHSIBhOZgkBDZx9GY9kyF5rWFI2suJEMEBGZErTraA+v9lY4OCMqMt/L4PguD5qJJemGX802njSeVKw0vhYwD8XwdTdMzZ57u+abjMH3mNXkVArYyb21tzTgnmpWVCG84m1GAZjlmdLORW/42EO7110d4Xvmsnbd3fHkLHrJ3cckyYb7wJ0+OSIMwWkuUmClZrjwgD97TNRyHVoxy3vpyphLW6ePaiY51gblprPae9kJW6DTU9MT2qp44fdkoyAltVM5dzRBDFFouY5CPp8Ox4WRkZ34uv2fRhshuL80cHA1odhRGFQnx6Pu2224byARa/Grizs5OjUY3voUtPX2mBSiUz9IyVDYqYFxDSvfdkl86L8/JES83XNjwPG46A+uWHZNTlaxXWK/sbHgJwBViaLJjpk/zwDS25thKtxaBwwuZ9MWLF3smGup5PdMGknDOa4cQmEsR83IPe3crs480tWCrhl/q5CJC13X1wQ9+cMBQG0ZCNv42bLWDMIwlGmQUt8E7B7T3d8HNQk4I7bFQvieeeGKQq/r0gWvXrtU3v/nNunTp0sDh0VpLLs7lWrTQv52IoyFyd1EMXrUgn/n0xBNPzBgY80/l53NyUiqz3Jvv+SbsRp8dVHLzA/dap+2MjOjQEe5hLtZXz5e/FzXUqgWN9YUXXhgwLOFHTsAMAQa7epyesuqoKHPcRguupeIZcragNLTv7e3VU0891ffD5wjcUM2VQyuZ6XFhqbW+lpHE1VuchCOsI32+9cE9aURPPfXUYO52YJcvX67Pf/7z9fzzz9f169d7Wpm318ENs5mLI2iLD4aj8AO5W5Y2eM+RxnNPPfXUQIa5JJcwf150xLEyR7ZHWi89r67rBidTZLAxH4wG4Zkbuv3kk08OaIXPtORtOtJWW2i7Ydd1e1X1uZve+O3R7q+q//uPm4iXsX0nzef/z3N55XQ6vfO4GxbNWT83nU7/wgkG/v9s67ru098pc6n6zprPrbkc3xatBj/2cg76x9y+k+ZS9Z01n1tzOaYtBINvtVvtVvvjbyfayH+r3Wq32h9fu2Wst9qt9m3SbhnrrXarfZu0W8Z6q91q3ybtlrHearfat0n7fwDyzLT9tkXj9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2023-11-27 09:03:14 - >>> DSS finder chart written to finder_chart_2023RN3_F65_20231118_dss.png.\n",
      "2023-11-27 09:03:14 - Done.\n"
     ]
    }
   ],
   "source": [
    "#generate_fscoords_finder_chart_gaia('358P','G37','2023-10-16 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231016_ldt/',5)\n",
    "#generate_fscoords_finder_chart_gaia('P/2018 P3','G37','2023-10-13 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231013_ldt/',5)\n",
    "#generate_fscoords_finder_chart_gaia('426P','G37','2023-10-14 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231014_ldt/',5)\n",
    "#generate_fscoords_finder_chart_gaia('62412','G37','2023-10-14 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231014_ldt/',5)\n",
    "#generate_fscoords_finder_chart_gaia('358P','G37','2023-10-19 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231019_ldt/',5)\n",
    "#generate_fscoords_finder_chart_gaia('P/2018 P3','G37','2023-10-19 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231019_ldt/',5)\n",
    "\n",
    "#generate_fscoords_finder_chart_gaia('358P','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "#generate_fscoords_finder_chart_gaia('426P','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "#generate_fscoords_finder_chart_gaia('427P','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "#generate_fscoords_finder_chart_gaia('433P','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "#generate_fscoords_finder_chart_gaia('455P','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "#generate_fscoords_finder_chart_gaia('456P','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "#generate_fscoords_finder_chart_gaia('62412','675','2023-10-10 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231010_palomar/',5)\n",
    "\n",
    "#generate_fscoords_finder_chart_gaia('358P','G37','2023-09-22 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20230922_ldt/',5)\n",
    "#generate_fscoords_finder_chart_gaia('313P','568','2015-12-09 00:00:00.000','/volumes/Fantom12a/BackupData/TargetedObjects/MBCs/313P/keck_lris/',5)\n",
    "#generate_fscoords_finder_chart_gaia_ephems('2023 RN3','Z31','2023-11-16 20:00:00.000','/users/hhsieh/dropbox/hhsieh_db/astro_db/00_2023RM3/',5)\n",
    "\n",
    "generate_fscoords_finder_chart_gaia_ephems_refcat('2023 RN3','F65','2023-11-18 00:00:00.000','/volumes/Fantom12a/BackupData/Data2023/ut20231118_ftn/',5)\n",
    "\n",
    "#generate_fscoords_finder_chart_gaia_coords(ra_coord,dec_coord,'/users/hhsieh/dropbox/hhsieh_db/astro_db/00_2023RM3/')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}