Jake's_Cleaned_Up_Work.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating Synthetic Populations \n",
    "\n",
    "- Inherited from Holmans' PS1_debias.ipynb notebook \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import math\n",
    "import random\n",
    "from collections import defaultdict\n",
    "import MPC_library as MPCL\n",
    "import healpy as hp\n",
    "import os\n",
    "import sys\n",
    "\n",
    "\n",
    "from collections import Counter\n",
    "\n",
    "\n",
    "import glob \n",
    "\n",
    "\n",
    "import warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import kepcart as kc\n",
    "\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Centaur Samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# M == Mean Anomaly\n",
    "# E = Eccentric Anomaly\n",
    "# e = Eccentricity (shape of orbit)\n",
    "# f = True Anomaly\n",
    "def M_to_E(M, e, E1=None, iter=200):\n",
    "    if np.any(E1==None):\n",
    "        E1=M\n",
    "    for i in range(iter):\n",
    "        E0 = E1\n",
    "        E1 = M + e*np.sin(E0)\n",
    "    return E1 % (2*np.pi)\n",
    "\n",
    "def E_to_M(E, e):\n",
    "    M = E - e*np.sin(E)\n",
    "    return M % (2*np.pi)\n",
    "\n",
    "def E_to_f(E, e):\n",
    "    f = 2.0*np.arctan(np.sqrt((1.+e)/(1.-e)))*np.tan(E/2)\n",
    "    return f % (2*np.pi)\n",
    "\n",
    "def f_to_E(f, e):\n",
    "    E = 2.0*np.arctan(np.sqrt((1.-e)/(1.+e)))*np.tan(f/2)\n",
    "    return E % (2*np.pi)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_vec(nsamp=10):\n",
    "    u = np.random.uniform(-1,1, size=nsamp)\n",
    "    incl = np.arccos(u)\n",
    "    longnode = np.random.uniform(0,2*np.pi, size=nsamp)\n",
    "    r = np.zeros((nsamp, 3), dtype=np.float)\n",
    "    r[:,0] = np.sin(incl)*np.cos(longnode)\n",
    "    r[:,1] = np.sin(incl)*np.sin(longnode)\n",
    "    r[:,2] = np.cos(incl)\n",
    "    return r\n",
    "\n",
    "def random_ortho(nsamp):\n",
    "    rvecs = random_vec(nsamp)\n",
    "    pvecs = np.zeros((nsamp,3))\n",
    "    pvecs[:,0]=rvecs[:,1]\n",
    "    pvecs[:,1]=-rvecs[:,0]\n",
    "    p = np.linalg.norm(pvecs,axis=1)\n",
    "    pvecs=pvecs/p[:,None]\n",
    "    cvecs=np.cross(rvecs,pvecs)\n",
    "    theta = np.random.uniform(0,2*np.pi, size=nsamp)\n",
    "    fvecs = pvecs*np.cos(theta[:,None]) + cvecs*np.sin(theta[:,None])\n",
    "    return rvecs, fvecs, np.cross(rvecs, fvecs)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creates a random sample of objects of size nsamp.\n",
    "#Most attributes are chosen uniformly, but eccentricity is via a normal distribution defined by muEc, sigmaEc.\n",
    "\n",
    "def centaur_sample(nsamp=1000000, e_max=0.99):\n",
    "    r0, pvecs, _ = random_ortho(nsamp)\n",
    "    x, y, z = r0[:,0], r0[:,1], r0[:,2]\n",
    "    xn, yn, zn = pvecs[:,0], pvecs[:,1], pvecs[:,2]\n",
    "    \n",
    "    \n",
    "    #normal distribution centered at muEc with std dev sigmaEc\n",
    "    muEc, sigmaEc = .4, .2\n",
    "    lowerEc, upperEc = 0, e_max\n",
    "    TRN = stats.truncnorm((lowerEc-muEc)/sigmaEc, (upperEc-muEc)/sigmaEc, loc=muEc, scale=sigmaEc)\n",
    "    e = TRN.rvs(size=nsamp)\n",
    "    \n",
    "    meananom = np.random.uniform(0,2*np.pi, size=nsamp)\n",
    "    E = M_to_E(meananom, e, iter=10)\n",
    "    #muQ, sigmaQ = 30, 10\n",
    "    #lowerQ, upperQ = 5, 25\n",
    "    #radiusRandomizer = stats.truncnorm((lowerR - muR)/sigmaR, (upperR-muR)/sigmaR, loc=muR, scale=sigmaR)\n",
    "    q = np.random.uniform(2, 35, size=nsamp)\n",
    "    a = q/(1.-e) \n",
    "    incl = np.arccos(zn)\n",
    "    longnode = np.arctan2(xn, -yn)\n",
    "    \n",
    "    cosnode=np.cos(longnode)\n",
    "    sinnode=np.sin(longnode)\n",
    "\n",
    "    rcosu = x*cosnode + y*sinnode\n",
    "    rsinu =(y*cosnode - x*sinnode)/zn\n",
    "    u = np.arctan2(rsinu, rcosu)\n",
    "    trueanom = E_to_f(E, e)\n",
    "    argperi = (u - trueanom) % (2*np.pi)\n",
    "    \n",
    "    \n",
    "    return a, e, incl, longnode, argperi, meananom\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Given a set of state variables, this method runs the nbody-r and probe executables to integrate \n",
    "#Change nbodyDir to wherever the probe and nbody-r files are.\n",
    "#Returns a dictionary with times as keys and a position/velocity array as output. \n",
    "\n",
    "def createFilesAndIntegrate(states, timestep, steps_per_output, num_outputs) :\n",
    "    nbodyDir = \"/Users/jake/Dropbox/centaurs/nbody_safe/\"\n",
    "    outputDir = \"/Users/jake/Dropbox/centaurs/nbody/jakeTests/\"\n",
    "    filestem = \"DoIExist_Year_\" + str(int(time.time()))#_%fZ\")\n",
    "    # to set up ICs\n",
    "    # - this is *NOT* the direct input to the n-body integrator\n",
    "    # - it should more be regarded as a header file / info file \n",
    "    #with open(nbodyDir + filestem + '.txt', 'w') as file:\n",
    "    with open(outputDir + filestem + '.txt', 'w') as file:\n",
    "        outstring='# These are initial cartesian ecliptic states for a sample of orbits\\n'\n",
    "        outstring+='# These are a replication of the test-particles in %s\\n' % nbodyDir + filestem\n",
    "        outstring+='# Intended to demonstrate tracklet creation at a number of distances\\n'\n",
    "        outstring+='# Generated using pyTrax_Data_Ingest_Example_Synthetic_3.ipynb notebook\\n'\n",
    "        outstring+='#\\n'\n",
    "        outstring+='# x y z xd yd zd\\n'\n",
    "        outstring+='# units are AU and AU/day\\n'\n",
    "        outstring+='#\\n'\n",
    "        file.write(outstring)\n",
    "        for state in states:\n",
    "            outstring='particle\\n%.16f %.16f %.16f %.16e %.16e %.16e\\n' % (state.x, state.y, state.z, state.xd, state.yd, state.zd)\n",
    "            file.write(outstring)\n",
    "\n",
    "    # THIS is the main input file \n",
    "    # need to add in the planetary positions, etc, at the top of the main input file \n",
    "    with open(outputDir + filestem, 'w') as file:\n",
    "        outstring='time 2454466.5\\n' # We should set this to correspond to the Julien Date at the start of 2010\n",
    "        outstring+='hill_sphere_factor 0.25\\n'\n",
    "        outstring+='dt '+str(timestep)+'\\n' # timestep in days \n",
    "        outstring+='nupdate '+str(int(steps_per_output))+'\\n' # how many timesteps to take before producing output \n",
    "        outstring+='noutput 1\\n' \n",
    "        outstring+='ntotal '+str(int(num_outputs))+'\\n' # how many sets of output to produce: we want to run for ~10 years total \n",
    "        outstring+='planet_states\\n'\n",
    "        outstring+='particle_states\\n'\n",
    "        outstring+='G 2.959122082322128e-04\\n'\n",
    "        outstring+='Msun 1.0\\n'\n",
    "        outstring+='em_index 0\\n'\n",
    "        outstring+='gr_constant 0.0\\n'\n",
    "        outstring+='quad_constant 0.0\\n'\n",
    "        outstring+='barycentric\\n'\n",
    "        outstring+='sun\\n'\n",
    "        outstring+=' 1.561402529627539e-04  4.953005432015677e-03 -7.105631884623408e-05\\n'\n",
    "        outstring+='-6.389749013486683e-06  6.495964171502198e-07  1.100598619353550e-07\\n'\n",
    "        outstring+='m_planet 1.660114153000000e-07\\n'\n",
    "        outstring+='2.602148208137757e-01 -3.292244628782728e-01 -5.123727708368509e-02\\n'\n",
    "        outstring+='1.659532456067476e-02  1.864190471122602e-02 -3.929871874523593e-07\\n'\n",
    "        outstring+='m_planet 2.447838288000000e-06\\n'\n",
    "        outstring+='-7.178008846563310e-01 -3.564142618684388e-02  4.080959356442847e-02\\n'\n",
    "        outstring+='1.013934063822650e-03 -2.028288042635279e-02 -3.364175614938229e-04\\n'\n",
    "        outstring+='m_planet 3.003489615000000e-06\\n'\n",
    "        outstring+='-1.675444337270498e-01  9.738318399949668e-01 -8.848710649717524e-05\\n'\n",
    "        outstring+='-1.723965496695509e-02 -2.998089479566788e-03  2.014628834321812e-07\\n'\n",
    "        outstring+='m_planet 3.227156038000000e-07\\n'\n",
    "        outstring+='-1.650127782726454e-01  1.579950905388059e+00  3.698226487000504e-02\\n'\n",
    "        outstring+='-1.339404168689282e-02 -2.715126682046571e-04  3.232197353626530e-04\\n'\n",
    "        outstring+='m_planet 9.547919384000000e-04\\n'\n",
    "        outstring+='1.483366733392910e-01 -5.240607574939024e+00  1.838922497888381e-02\\n'\n",
    "        outstring+='7.449937516689820e-03  5.707878757173034e-04 -1.691124368854220e-04\\n'\n",
    "        outstring+='m_planet 2.858859806661000e-04\\n'\n",
    "        outstring+='-8.261542977968610e+00  4.193055697215763e+00  2.558059889927936e-01\\n'\n",
    "        outstring+='-2.822018174472759e-03 -4.987848378405353e-03  1.990803802348707e-04\\n'\n",
    "        outstring+='m_planet 4.366244000000000e-05\\n'\n",
    "        outstring+='1.964046375645850e+01 -4.238756428479231e+00 -2.702098023085751e-01\\n'\n",
    "        outstring+='8.010336315728508e-04  3.661284926651336e-03  3.230110758793685e-06\\n'\n",
    "        outstring+='m_planet 5.151389000000000e-05\\n'\n",
    "        outstring+='2.346585183697751e+01 -1.714470518227547e+01 -7.601640065764197e+00\\n'\n",
    "        outstring+='1.939090244950700e-03 2.304786228786757e-03 8.950850954553373e-04\\n'\n",
    "        file.write(outstring)\n",
    "        for state in states:\n",
    "            outstring='particle\\n%.16f %.16f %.16f %.16e %.16e %.16e\\n' % (state.x, state.y, state.z, state.xd, state.yd, state.zd)\n",
    "            file.write(outstring)\n",
    "            \n",
    "            \n",
    "        # to run nbody code ...\n",
    "    #(i) make the command as a string ...\n",
    "    os.chdir(nbodyDir)\n",
    "    command = \"%snbody-r %s \" % (nbodyDir, outputDir +filestem) \n",
    "    #(ii) execture the command ...\n",
    "    os.system(command)\n",
    "\n",
    "\n",
    "    # N.B. the  output files produced by nbody will look like ...\n",
    "    #filestem.pb # <<-- Binary file with posns of planets\n",
    "    #filestem.tb # <<-- Binary file with posns of test particles (what I am interested in)\n",
    "\n",
    "\n",
    "    # Generate detailed output files using \"probe\" \n",
    "    # probe : can look at the file while running or when finished\n",
    "    # \"case 4\" outputs barycentric posn & veloc of all test \n",
    "    outfilestem = filestem + \"_out\"\n",
    "    command = \"%sprobe 4 %s %s \" % (nbodyDir , outputDir + filestem , nbodyDir + outfilestem)\n",
    "    os.system(command)\n",
    "    \n",
    "    nbody_dict=defaultdict(list)\n",
    "    # Find all the output files ...\n",
    "    # ... then read them all into a dictionary  ...\n",
    "    # ... key-ed on times \n",
    "    for name in sorted(glob.glob( '%s/%s*' % (nbodyDir, outfilestem) )):\n",
    "        i = 0\n",
    "        for v in open(name):\n",
    "            if v.startswith('#'):\n",
    "                t = float(v.split()[1])-2454466.5 # not permanant solution. FIX THIS! This is a problem!\n",
    "            else:\n",
    "                nbody_dict[t].append(tuple(map(float, v.split())))\n",
    "                i+=1\n",
    "    # Sorted keys (times)\n",
    "    sim_times = np.array(sorted(nbody_dict.keys()))\n",
    "    # Make the data into arrays \n",
    "    for key in sim_times:\n",
    "        nbody_dict[key] = np.array(nbody_dict[key])\n",
    "    for fn in range(2001):\n",
    "        command = \"rm \"+ filestem +\"_out_\" + str(fn*2+2454466.5)\n",
    "        os.system(command)\n",
    "    return nbody_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chebyshev Approximations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fits a piecewise chebyshev polynomial to the first 3 columns of nbody_dict (which should have x,y,z positions \n",
    "#and velocities). The maximum error at any point will be less than maxError. Based on some testing, \n",
    "#I recomend order 15 and pieces = 1, but feel free to mess around and change this up. \n",
    "#The output will be a list of objects. Each object is a list of coordinate functions, and each coordinate function \n",
    "#is a dictionary which maps startpoints to an array of polynomial coefficients. For example, if an object is: \n",
    "#[{0.0: [1,0,0,0...]}, {0.0: [1,0,0,0...], 2000.0:[2,0,0,0...]}, {0.0: [1.0]}]\n",
    "#then its X coordinate is always 1, its y coordinate is 1 from time 0 to 2000, and 2 from 2,000 on, \n",
    "#and its z coordinate is always 1.\n",
    "\n",
    "def getPiecewiseChebs (nbody_dict,maxError, order = 15, pieces = 1):\n",
    "    sim_times = np.array(sorted(nbody_dict.keys()))\n",
    "    Xchebs = []\n",
    "    Ychebs = []\n",
    "    Zchebs = []\n",
    "    allYs = []\n",
    "    errors = []\n",
    "    # loop over the objects in the dictionary ...\n",
    "    nSamp = len(nbody_dict[sim_times[0]])\n",
    "    for samp in range(nSamp):\n",
    "        # Get the positions at each of the times\n",
    "        ys= np.vstack([nbody_dict[t][samp] for t in sim_times])\n",
    "        allYs.append(ys)\n",
    "        x = ys.T[0]\n",
    "        y = ys.T[1]\n",
    "        z = ys.T[2]\n",
    "        Xchebs.append(getGoodApproximation(sim_times, x, order, maxError))\n",
    "        Ychebs.append(getGoodApproximation(sim_times, y, order, maxError)) \n",
    "        Zchebs.append(getGoodApproximation(sim_times, z, order, maxError))\n",
    "    return ([Xchebs, Ychebs, Zchebs], allYs)\n",
    "\n",
    "\n",
    "#This gets a piecewise chebyshev polynomial (each piece of some order/degree) which fits the points (x_1, y_1), (x_n, y_n)... with some maximum error.\n",
    "#If a single piece does not fit the data with the prescribed accuracy, it will bisect itself and try fitting a polynomial to each side.\n",
    "#Often, when an object has a close encounter with a planet, this ends up repeatedly bisecting to cover one particular small region (near the encounter) very accurately.\n",
    "#Otherwise, one piece usually fits. The average number of pieces for my centaur sample was about 1.3\n",
    "def getGoodApproximation(x, y, order, maxError):\n",
    "    chebCandidate = np.polynomial.chebyshev.chebfit(x, y, order)\n",
    "    if len(x) <= order:\n",
    "        return {x[0]: chebCandidate}\n",
    "    quickEval = np.polynomial.chebyshev.chebval(x, chebCandidate)\n",
    "    if error(quickEval, y) < maxError: #error threshold\n",
    "        return {x[0]: chebCandidate}\n",
    "    else:\n",
    "        half = int(len(x)/2)\n",
    "        dict1 = getGoodApproximation(x[:half], y[:half], order, maxError)\n",
    "        dict2 = getGoodApproximation(x[half:], y[half:], order, maxError)\n",
    "\n",
    "        dict1.update(dict2)\n",
    "        return dict1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Finds the maximum difference between any pair of nth elements of two lists. \n",
    "def error(a, b):\n",
    "    err = 0\n",
    "    for i in range(len(a)):\n",
    "        err = max(err, abs(a[i]-b[i]))\n",
    "    return err\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Evaluates a piecewise chebyshev polynomial at a given time. Returns a list of 3 values. \n",
    "def evaluate3DPiecewiseCheb(evalTime, cheb):\n",
    "    times = [list(i.keys()) for i in cheb]\n",
    "    for i in range(3):\n",
    "        times[i].sort()\n",
    "    iterators = [0,0,0]\n",
    "    for i in range(3):\n",
    "        while iterators[i] < len(times[i])-1 and times[i][iterators[i]+1] < evalTime:\n",
    "            iterators[i] +=1 \n",
    "    vals = [np.polynomial.chebyshev.chebval([evalTime], cheb[i][times[i][iterators[i]]])[0] for i in range(3)]\n",
    "    return vals\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creates a sample of objects, then simulates their motion, and returns chebyshev approximations for the objects' orbits.\n",
    "def getChebsFromSimulation(nSamp, sim_timestep,num_outputs,total_time,maxError = .0000005):\n",
    "    errors = ([],[],[],[],[],[],[]) \n",
    "\n",
    "\n",
    "    steps_per_output = int((total_time/num_outputs)/sim_timestep) \n",
    "    time_per_output = int(total_time/num_outputs) \n",
    "    sim_times = np.arange(0, total_time+1, steps_per_output*sim_timestep)  \n",
    "    \n",
    "    \n",
    "    a, e, incl, longnode, argperi, meananom= centaur_sample(nsamp=nSamp, e_max=0.99)\n",
    "    states = kc.cartesians(len(a), MPCL.Constants.GMsun, a, e, incl, longnode, argperi, meananom)\n",
    "    nbody_dict = createFilesAndIntegrate(states, sim_timestep, steps_per_output, num_outputs)\n",
    "    chebs, ys = getPiecewiseChebs(nbody_dict, maxError, 15, 1)\n",
    "    chebs = np.array(chebs).T\n",
    "    \n",
    "    return (chebs)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Same as getChebsFromSimulation, but instead of returning the objects, writes them to files. \n",
    "#Default is 10,000 objects per file, computed 1000 at a time.\n",
    "def makeChebFiles(start, end, sampleSize=1000, numSamples = 10):\n",
    "    ender = ''\n",
    "    for fileNum in [start, end]:\n",
    "        print()\n",
    "        print()\n",
    "        print()\n",
    "        print('File: ' + str(fileNum))\n",
    "        print()\n",
    "        \n",
    "        #This section is the only important stuff. Runs a simulation of sampleSize objects numSamples times, and \n",
    "        #writes the objects to a file named chebs<number>.txt. \n",
    "        with open('/Users/jake/Dropbox/centaurs/nbody_safe/chebs'+str(fileNum)+'.txt', 'w') as file:\n",
    "            for sampleNumber in range(numSamples):\n",
    "                totalTimer0 = time.time()\n",
    "                chebs = getChebsFromSimulation(sampleSize, 2.0, 2000.0, 4000.0,.0000005)\n",
    "                for cheb in chebs:\n",
    "                    file.write(str(cheb) + '\\n')\n",
    "                \n",
    "                #Just printing out progress. Not important. \n",
    "                totalTimer1 = time.time()\n",
    "                if sampleNumber % 10 == 1:\n",
    "                    ender = 'st'\n",
    "                elif sampleNumber % 10 == 2:\n",
    "                    ender = 'nd'\n",
    "                elif sampleNumber % 10 == 3:\n",
    "                    ender = 'rd'\n",
    "                else:\n",
    "                    ender = 'th'\n",
    "                print('through the ' + str(sampleNumber) + ender + ' part in ' + str(int(totalTimer1-totalTimer0)) + ' seconds')\n",
    "        \n",
    "        file.close()\n",
    "        \n",
    "def getChebsFromFiles(numFiles):\n",
    "        obs = []\n",
    "        for fileNum in range(numFiles):\n",
    "                with open('/Users/jake/Desktop/savedChebs/chebs'+str(fileNum)+'.txt', 'r') as file:\n",
    "                        fullStr=file.read()\n",
    "                        fullStr=fullStr.replace('\\n', '')\n",
    "                        fullStr=fullStr.replace(' ' , '')\n",
    "                        fullStr=fullStr.replace('\\t', '')\n",
    "                        abc = fullStr.split('][')\n",
    "                        abc[0] = abc[0][1:]\n",
    "                        abc[-1] = abc[-1][:-1]\n",
    "                        for particle in abc:\n",
    "                                piecewiseChebStrings = particle.split('}{')\n",
    "                                piecewiseChebStrings[0] = piecewiseChebStrings[0] [1:]\n",
    "                                piecewiseChebStrings[-1]= piecewiseChebStrings[-1][:-1]\n",
    "                                obj = []\n",
    "                                for coordinatePiecewiseCheb in piecewiseChebStrings:\n",
    "                                        pieces = coordinatePiecewiseCheb.split('),')\n",
    "                                        pieces[-1]= pieces[-1][:-1]\n",
    "                                        newDict = {}\n",
    "                                        for piece in pieces:\n",
    "                                                key, value = piece.split(':')\n",
    "                                                key = float(key)\n",
    "                                                value = value[7:-1]\n",
    "                                                coefficients = value.split(',')\n",
    "                                                for i in range(len(coefficients)):\n",
    "                                                        coefficients[i] = float(coefficients[i])\n",
    "                                                coefficients = np.array(coefficients)\n",
    "                                                newDict[key]=coefficients\n",
    "                                        obj.append(newDict)\n",
    "                                obs.append(obj)\n",
    "                file.close()\n",
    "        return obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XdYFOfax/Hvs3QBUcSCImJH7IhdE40x9pJoLNGoMdH05JyUk/qmnfTezTGJNWpiS2I3xtg79oIgCgqKiiJd+rx/zJpoggKy7Gy5P9e1l8uyzPxY8Gb2mWfuR2mahhBCCMdhMjqAEEIIy5LCLoQQDkYKuxBCOBgp7EII4WCksAshhIORwi6EEA5GCrsQQjgYKexCCOFgpLALIYSDcTVipwEBAVpISIgRuxZCCLu1e/fuC5qmVS/peYYU9pCQECIjI43YtRBC2C2l1MnSPE+GYoQQwsFIYRdCCAcjhV0IIRyMFHYhhHAwUtiFEMLBSGEXQggHU+7CrpSqq5Rap5SKUkodVko9aYlgQgghbo4l5rEXAE9rmrZHKeUL7FZKrdE07YgFti3s0KWsPPYnphJ3IYuMnAI83UzU8PWkVZAf9QO8UUoZHVEIh1buwq5pWhKQZL6foZSKAuoAUtidSGGRxoqDSfy46xRbj1/kekvp1qnixdC2tRnTsR61q3hZN6QQTkJZcjFrpVQIsBFooWla+t8+NxmYDBAcHNzu5MlSXUAl7MCmY8m8sfQIx85nUq9aJYa0rk2nhtVoWtOXyl5u5OQXcjr1MntPpbL68Fk2HbuAi1KMbF+Xp3o3oaq3u9HfghB2QSm1W9O0iBKfZ6nCrpTyATYAb2matvhGz42IiNCkpYD9y8kv5I1lR5i74xT1A7x5+o4m9G8RiMl046GWxEvZfL3+OD/tSqCKlxuvDW7OoNa1rZRaCPtl1cKulHIDlgGrNU37uKTnS2G3f+fSc5g0K5KDp9OY1L0BT/VugqebS5m2EZWUzvOLDrA/MY0xHYN5ZVAYHq5l24YQzqS0hd0Ss2IU8D0QVZqiXh5FWhEJ6QkVuQtRCieSMxny5RaOn89k6r0RvNi/WZmLOkCzwMosergLD97agDk7TjHyf9u5mJlbAYmFcC6WmMfeFbgXuE0ptc9862+B7f7Dq1tfZcKqCeQU5FTE5kUpnEjOZNTU7eQXFrHgoS70DqtZru25uph4oV8zpowJJyopnbv/t43ES9kWSiuEc7LErJjNgFXmrw1qMIhfYn9hYcxCxoaNtcYuxVXOpF5m9LfbKSzSmDupE01r+Vps2/1aBhLg68HEGbsYPmUbPz3YiXrVvC22/QpRVAgXYyE5Gi5EQ0ocZCXrt+wUKMyHogLQCsHVCzx8wMMXvGtAlbrgVxf860Otlvp9mQYqLMSis2JKqzxj7Pevvp/jqcdZOWwlXq4yXc5aMnMLuPubbSSkZLPgoc40C6xcIfuJSkrnnm+34+3hyoKHOhPoZ0M/46IiOB0Jx9fBqW2QGAl5GX993qcW+NYE7+pQqRq4uIPJFZQJCnIgNwNy0yHjHKQlQF7mX1/r6Qe1WkFwZ6jfHYI6gJun9b9HYdOsPiumLMpT2Pec28P4VeN5ut3TTGgxwbLBRLGKijQmzYpkfUwy0ya059YmJS7gUi4HE9MY/e12alT2YP6DnQnw8ajQ/d1QURHErYdDiyFmlX40joKazaFuRwhqDzVCoVoj/Wi8tDQNLl/Sj/jPHtRvZ/bC2QOgFYGLB9TrDKED9VvlwIr6DoUdcdjCDvDgmgeJuhjFymEr8Xaz8bfrDuCrdbF8sDqaN4Y0Z1znEKvsc1d8Cvd+v4Pmtf2Y80DHmzo5Wy5pp2HPLNg3Rz+69qgMjW6H0AHQqBd4Va2Y/eakwcltEL9J/0NyMVZ/PKg9tBoJLYdX3L6FzXPown4g+QBjVozhibZPMKnVJAsmE3+3My6FUVO3MaBVbT4f1caq7QBWHEzikTl7GNqmNp+MtNK+k6Nhy2dwYL4+Pt6gB4TfC00HWH9oRNP0PEeXwuFf4dxB/Ui+2UBoe6+eTcblnYpDF3aAR9c+yr7z+1g1bBW+7pY7iSf+cikrj36fbcLTzcTSx7vh6+lm9Qxf/nGMD3+L4eneTXi8V+OK29HF47D2dTjyq36iM3wcdHpYP7lpK5L2w945cOAnyEmF6s30jK1Gyni8k7DaPHajPNrmUdLz0vnhyA9GR3FYry89zIXMXL68J9yQog7waM9G3Nm2Dh//HsOGmGTL7yDrAqx4Fr7qAMd+h1uehX8fhv7v21ZRBwhsred6OhqGfgMurrD0CfikOax/Tx/GEQI7Luxh1cK4re5tzI6aTebVswuERfx+5By/7DvDoz0b0aKOn2E5lFK8fWdLmtTw5d8/7SMp7bJlNqxpsPcH+KId7PpeP0J/Yi/c9jJ4V7PMPiqKmye0GQ0PboLxS6FOO1j/NnzaEjZ8ADnpJW9DODS7LewAk1tPJiMvgx+jfzQ6ikNJu5zPS78cJLSWL4/2bGR0HLzcXfh6bDi5+YU8Nncv+YVF5dvgxeMwcxD8+ijUaAYPb4WBn+hTFe2JUlD/FhgzHx7cCPW6wro34bNWsPkTyJcL+ZyVXRf25tWa07V2V2Yfmc3lAgsdyQneWRHFhcw8PhjeGndX2/gVaVjdh3eGtWL3yUt8+Fv0zW1E02D3TPimGyQdgEGfwYQV+nRFexfYGkbPg0nr9Bk0v78GX7XXp2kacB5NGMs2/teWw+RWk0nJSWFRzCKjoziEfQmp/LgrgYldQ2gZZNwQTHEGt67N6A51mbrxBDtOXCzbF1++BAvG62PSdTvAozug3QQw2f1/gWvVCYcxC2Dcr/oUzYX3wbS+cHq30cmEFdn9b3V4zXDa1WzH9EPTySvMMzqOXSsq0nj110PU8PXgiYqcgVIOLw8II9i/Ek/N309GTn7pvihpP3zTHY4uh95vwNifHf+CnwY99OGZQZ9BynH4thcsf0ZOsDoJuy/soB+1n798nl+P/2p0FLu2YHcC+xPTeKF/qGGzYEri7eHKxyPakJR2mTeWlmKRroML4fs++tWc9/8GXZ90vKP06zG56O9KHt8DHR+EXd/BVx31KZ0yPOPQHOI3vHNgZ1pUa8H3B7+noKjA6Dh2Ke1yPu+viiaiXlWGtqljdJwbalevKo/0aMSC3YmsOnS2+CcVFenjzIvuh9ptYPJ6ffaIM/KsDP3eg0lrwTsA5o+DeaMh/YzRyUQFcYjCrpRicqvJnM48zcq4lUbHsUtfr4slJTuP1wY3t4vFpp/o1ZjmtSvz8i+HSM3+2xBcQS4sfkCfGdJuAoxbAj41DMlpU+q0g0nr4Y43IW4DfN0JDiyQo3cH5BCFHeDWurfSuGpjvj34LUVaOafDOZkzqZeZvjWeu9oGGTpnvSzcXU28P7wVl7LzeHtF1F+fyM2AuSPg0CK4/XUY+Cm4ypqqf3JxhS6Pw0ObIaCp/gdwwQS9zbBwGA5T2E3KxKSWk4hLi2Ndwjqj49iVT9bEAPDUHU0MTlI2zWv7MfmWBsyPTGRL7AX9KtIZAyBuEwydAt3+Jb1UrqdaQ5i4Cnq9qp9U/roTxPxmdCphIQ5T2AF61+tNHZ86TD803egodiP6bAaL9iQyrlM96lSxod7npfRkr8bUD/Dm3UWbKZoxEJJj9Pncbe4xOprtM7lA96dg8jq9h/zcu2H1S1Ags8vsnUMVdleTK+PCxrE/eT97z+81Oo5d+GD1UbzdXW3iCtOb4enmwgf9a/NB1ssUXjwB9/wITfoYHcu+1GoJD6yF9pNg25cwvS9cijc6lSgHhyrsAEMbDcXPw49ph6YZHcXmRcan8HvUeR7q0ZCq3nY6Dp11gYj142ngcp77cp7msGdboxPZJzdPGPAhjJgFF2Lhm1vgyBKjU4mb5HCFvZJbJUaHjmZ9wnpOpJ0wOo5N+3hNDAE+HkzsamNdDEsrJx1m3wkpceSNmMcRr3Be/fUwRUUyy+OmhQ2BhzbqY/Dz74WVz+trtwq74nCFHWB06Gg8XDyYeXim0VFs1q74FLYev8hDtzbAy93KqxNZQn4O/HgPnD8CI2fj06wXz/cLJfLkJRbvPW10OvtWNQQmroaOD8OOKTBrKGRWQMtkUWEcsrD7e/oztNFQlh5fSnK2/EIW5/O1xwjwcWdMx3pGRym7okJYPElfPm7I19C4NwDDw4MID67COyuiSLssR5nl4uoO/d6FO6fqC3hP7QGn9xidSpSSQxZ2gHFh4yjUCpkTNcfoKDZn98lLbDp2gcm32OHRuqbpC2NELYE73oLWI//8lMmkeGNICy5l5/HxzXaAFNdqPVI/elcmvZnYvnlGJxKlYJHCrpSappQ6r5Q6ZIntWUJw5WB6BfdifvR8svKzjI5jUz5be4xq3u6M7WSHR+vbvoLI76HLE9DlsX98ukUdP+7tVI/Z209y+Iw0vLKIKy0ZgjvCLw/p4+5FhUanEjdgqSP2GUBfC23LYia2mEhGfgYLYxYaHcVm7D11iY0xyUy6pQGV3F2NjlM2Mathzf9Bs8H6VaXX8dQdTalayZ3/++WQnEi1FO9qelfMTo/o4+7zRutX+QqbZJHCrmnaRsDmrkluEdCC9rXaM/vIbPKLZMwV9KN1f2937rW3o/XzUbDwfqjZAu785oYdGv283Hi+Xyh7TqWyZL80urIYF1fo+w4M+Bhif4dp/SAt0ehUohgOO8Z+xYTmEziXfY5VcauMjmK4Q6fTWB+dzP3d6uPtYUdH61kXYe5IcK8Eo38Ed+8Sv2RYeBCtgvx4b9VRLufJsIFFtb9fX8wj9aTe5/2MXAxoa6xW2JVSk5VSkUqpyORk681U6VanG/X96jP7yGw0J+9iN3XjCXw8XO1rbL2wQF/5KPMcjJoHfqVrKWwyKf5vYBhJaTlM3SjXM1hco176SVUXd5jeH6KWGZ1IXMVqhV3TtKmapkVomhZRvXp1a+0WkzIxttlYolKi2H3OeZcHS0jJZvnBJO7pGIyfl20uolGsdW/p0xoHfgJBZeun3j7EnwGtAvlmw3GS0mRNXIurGab3eK8RBj+NhR1TjU4kzBx+KAZgUMNB+Hn4MfvIbKOjGOa7TScwKezrKtPolbD5Ywgff9NNvZ7vG0qhpvHBKpn+WCF8asCEZdC0P6x8Ftb+V/q72wBLTXecB2wDmiqlEpVS91tiu5bi5erFiCYjWJewjoT0BKPjWF1KVh4/RSYwtE0davl5Gh2ndFLi4OcHoVYr6Pf+TW+mrn8lHuhWn8V7T7MvIdWCAcWf3Lz0HjPh42HTh7DkMX0ITRjGUrNiRmuaFqhpmpumaUGapn1vie1a0qjQUbiYXJhz1PkuWJq5NZ6c/CIevLWB0VFKJz9HH1cHvWC4le+P0SM9GxHg48F/lx1x+vMsFcbFVV84+9bnYO8P+tBMXrbRqZyWUwzFANSoVIO+IX35+djPZOQ5z/zb7LwCZm2L5/ZmNWlUw9foOKWz5hVI2g9DvwH/8g8d+Xi48p8+Tdl98hLLDiRZIKAollLQ80UY8BHErIJZQ2RlJoM4TWEHuDfsXrILsll8bLHRUaxm/q4ELmXn83APOzlaj/kNdv5Pb0AV2t9imx3WLoiwwMq8u/IoOfky/bFCtX8A7p4BSfv0NgSyaLbVOVVhD6sWRrua7ZgTNYeCIscfAywoLOLbTXFE1KtKu3r+RscpWeZ5+PURqNEcbn/Nopt2MSleHtiM06mXmbUt3qLbFsVoPhTGLob00zC9nyzcYWVOVdhBP2pPykpi7am1RkepcMsPJnE69TIP3drQ6Cgl0zT45RG9x/qw78o9rl6cLg0D6Nm0Ol/+EUtqtiz/VuHqd4dxS+Byqn6VanKM0YmchtMV9h5BPQjyCXKKqY/TtsTTIMCb20JrGB2lZDv+B7Fr4I439fnRFeS5fqFk5Bbw9frjFbYPcZWgdjBhORTl60fuZw8ancgpOF1hdzG5MDZsLPuT93Mg+YDRcSrMnlOX2J+QyoSuIZhMyug4N3buiH7CtHEf6DCpQncVWqsyw8ODmLElnoQUmbVhFbVawH2rwNUTZgyAhF1GJ3J4TlfYQV8X1cfNx6GP2qdtjsPX05Vh4UFGR7mxwnx9vrqHLwz5Sp9ZUcGeuqMJSulLAworCWgEE1eCl78+WyZuk9GJHJpTFnZvN2+GNR7GmpNrSMp0vOlvSWmXWXnoLKPa17X9Zl+bP4WzB/SWAT7WaTUR6OfFxG71+WXfaQ6dlp7tVlMlGCau0v+dMxyOrTE6kcNyysIOcE+ze9DQmHfU8VaEmb3tJJqmMa5ziNFRbuzsIdjwHjS/C8IGW3XXD/doSBUvN95bddSq+3V6vrX0MffqTfU1a2N+MzqRQ3Lawl7bpza9gnux8NhCsvMdZ6z1cl4h83aeondYTer6VzI6zvUV5utTGz39oP+HVt99ZU83HrutMZuOXWBjjKyLa1Xe1WDcr+bmYWMgWlpqW5rTFnaAMc3GkJGXwfK45UZHsZhf9p3mUna+7Tf72vypfnXpwI/1/+gGGNspmLr+Xryz8qistGRtXlX14l6zhd5+4OgKoxM5FKcu7OE1wgn1D2Vu1FyH6CGiaRrTt8QRFliZDvVt+IKkc4evGoIZYlgMD1cXnu0TSlRSOr/sO21YDqflVQXu/RkCW8H8cdLT3YKcurArpbgn9B5iU2PZddb+p2BtPX6RmHOZ3Nc1BGWF2SU3pTAffnnYPATzgdFpGNgykJZ1/PhwdbS0GjDCn8W9td74LWqp0YkcglMXdoB+9ftRxaMKc6Lsv+vjtM1xBPi4M6h1baOjXN/WL/QhmAEfgXeA0WkwmRQv9A/lTFoOM7fGGx3HOXn66cW9djgsmABHfjU6kd1z+sLu6erJ8CbDWZ+4ntOZ9vt2PP5CFn9En+eejvXwdHMxOk7xUk7oQzChA/VeIjaiS8MAejStztfrj5OWLYueG8KzMoxdBHXawYL74PDPRieya05f2AFGNh2JQvHT0Z+MjnLTZmyNx9WkGNsp2OgoxdM0WPYUmNxsYgjm757rG0p6Tj5fr481OorzulLcg9rDwvvh8C9GJ7JbUtiBWt61uC34NhYdW8TlAvtbGzMjJ5+FuxMZ2Ko2NXxtdIWkgwvgxDro9QpUtr2homaBlbmzbR2mb43ndKr9/Q44DA9fGLtQL+6L7oejjjNjzZrsq7Af+x1WPl8hayqOaTaG9Lx0lp+wv1+k+ZGJZOYW2O4Ux+wUWPU81ImA9ja1auI1nr6jKQCfSKsBY3n4wpgFENgG5o+HmNVGJ7I79lXYT22DHVNgxbMWL+5Xpj7OiZpjV1MfC4s0Zm6NJ6JeVVoG+Rkdp3i//Z/eunXQZ2Cy0fF/oE4VLyZ0CWHRnkSOnk03Oo5zuzIsU7O5Ps899nejE9kV+yrst70MnR+DXd/qR4AWLMD2OvXxj6PnOZWSzX22erQetwn2/QBdHte7/Nm4R3o0xNfDlfdWSqsBw12ZClm9Kfw4Bk6sNzqR3bCvwq6U3q+706Ow4xtY/aJFi/uVqY9zj8612DYr2rTNcdT286RP85pGR/mn/BxY9i+oUk9f5NgOVKnkziM9G7EuOpltxy8aHUdU8od7fwX/BjB3FMRvNjqRXbCvwg56ce/zlr4m5vav4beXLVbcr0x9XJewzi6mPkYlpbPtxEXGdQnB1cUGf5SbP4aLsXrnRncb7lvzNxO6hBDo58m7q47a1bCcw/Kupq/EVCUY5oyAU9uNTmTzbLAalIJS0Pcd6DAZtn0Jv79qseJuT1MfZ2yJx9PNxKj2dY2O8k/J0bDpY2h5NzTqZXSaMvF0c+HfvZuwPyGVFQfPGh1HgN7SefwSqBwIPwyHxEijE9k0ixR2pVRfpVS0UipWKfW8JbZZip1Cv/ch4n7Y8hmsfcMixb2Wdy161O3B0hO2fWnzxcxcft53mrvCg6hSyd3oONcqKoKlT4K7N/R5x+g0N2VYeBBNa/ryweqj5BcWGR1HgN7yd/xS/Yrl2XfBmb1GJ7JZ5V6FQSnlAnwF9AYSgV1KqSWaph0p77ZLsXO95atWpL/tVyb9BGs5+6TU8anD1jNbSc9Lt9mWvjO2xpHPJQaFN+Fslo0dVR5cAKd36udDVCHYWr5SevA2f55ecIBvt+zhrnYlr0SlsNH+PI7ExQVGTNebhv1wJ9w9A2o0MzpVmVTxrIKHi0eF7kOVdwxRKdUZeE3TtD7mj18A0DTtuodqERERWmSkBd9KFRXBsidhzyz9JF3PF8u1uY8iP2LG4Rm4KlcKtAILhRRCCJhy+xS61el2U1+rlNqtaVpESc+zxLppdYCEqz5OBDpaYLulZzLBwM/0I/cN74FygR43PwvjypFXgVbAyKYjaeZvW0cE+xIu8eOuBCZ0CSG0VmWj41xrz0z9LXKPF8DHBmfqlNGplGy+Xh9Lr9Ca9A67/vejISdZrS7rgj4MqxXp02l9axmdqFQaVWlU4fuwRGEv7v3nP37LlVKTgckAwcEV0M/EZIJBX+hH7+vf1j++5dmb29ZV31HX2l3pGdzTMhktZNaaLdR1bcQL3W7FZLKht//H/4CjG/R3TeGPGJ3GYo4d3836Pcn8t1dPqvtW7FtoUUZ1e8GMAbDuc5iwQl80W1jk5GkicPW0jCDgzN+fpGnaVE3TIjRNi6hevYIWLTaZYMiX0GoU/PEmbPro5jZz1ctia33N95y6xP6EVCZ0DbGtop5/WW/y5d8Quj1ldBqLeuaOpuQWFPH52mNGRxF/F9BYnwpZVAgzB+kdRIVFCvsuoLFSqr5Syh0YBSyxwHZvjskFhn4NLUfoM2U2f1rmTVxdzG3thNj0LfH4eroyLLzkk3lWtekjuBSnz1l3s9FGZDepQXUfRneoy7ydp4i7kGV0HPF3NUL12TIFOTBzMKQmlPw1Dq7chV3TtALgMWA1EAXM1zTtcHm3Wy4mFxg6BVoM0+e4b/2iTF9+dTG3pSP2pLTLrDiYxMiIunh7WGIUzUKSo/U/oK1GQYNbjU5TIZ7s1QR3VxMfro42OoooTs0wGPcL5KTrR+7pSUYnMpRF5rFrmrZC07QmmqY11DTtLUtss9xcXOHOqdD8Tv3q1G1flfpLbfWIffa2k2iaxvguIUZH+UtRESz9F3j46FcEO6jqvh5M6t6A5QeT2JeQanQcUZzA1nrjsKxkmDVY7yrqpOzzytPScnGFu77VF0xe/SJs/6ZUX2aLR+w5+YXM23mK3mE1qetvQ5fn75sDp7ZC7//axFJ3FWnSLQ0I8HHnnRVR0mrAVtVtD/fM11tZ/PKIvsauE3Lswg7g4gbDvodmg2DVc7BjaolfYlJXnTy1kSP2n/ee5lJ2vm11ccy6AGv+D4K7QNuxRqepcD4erjzRqzE74lJYH51sdBxxPSFd9avSY1bqi3UUOt+1KI5f2MFc3KdB0wGw8lnY9d0Nn25rR+yapjF9SxxhgZXpWN/f6Dh/+e1lyM3UT5jawOtkDaM7BBNSrRLvrjxKYZEctdusDpOgz9v6wtg/P6jPmnEizlHYAVzd9cuPm/SD5U9D5PTrP1ddfdf4grUl9iIx5zKZ2K2+TfyhAeDEBtg/D7r9S5+V4CTcXEw806cp0ecyWLwn0eg44kY6Pwq3vw6HFsKvjzpVcXeewg56cR8xExr30fuE755Z7NNsbR779C1xBPi4M6h1oNFRdHnZepOvqvWh+9NGp7G6AS0DaR3kx8drYsjJd55iYZe6/Qt6vqwfhCx9Uj/Z7wScq7ADuHrAyNnQqDcsfQL2zP7HU2xpVkzchSzWHj3PmI718HC1kWXl1r2lz1kf/AW4eRmdxuqUUjzfrxlJaTnM3BpvdBxRkluf1a+G3jsbVjxdIWsm2xrnK+xgLu4/QMNesORx2HftiklXnzy9+r4RZmyJw81FMaZTBbRhuBmJu/UFTtrdB/W7G53GMJ0bVqNH0+p8tS6W1Ow8o+OIkvR4Abr9GyKnwcrnHL64O2dhB/3qyFFzoEEPfVrU/h+NTvQPaZfzWbA7kUGta1PD1wau5izI1ccqfQOh9xtGpzHcc31Dycgt4Ov1x42OIkqiFPR6VV8zeef/LLrymi1y3sIO+jDCqLn6kecvD8OB+cDfZsUYOBSzIDKB7LxCJtrKFMdNH0FyFAz8VF9F3sk1C6zMsPAgZmyJ5+RFaTVg866smdzhQfPKa685bHF37sIO+lqco3+Cel31aVEHF147j92gk6eFRRoztsbTIcSfFnX8DMlwjbOH9MLeaiQ0ucPoNDbj2T5NcXVRvL0iyugoojSUgn7vQcRE2PKp3izQAYu7FHbQi/s9P+kX2iyehDp76M9PGTXGvubIORIvXWZitxBD9n+Nwnz49RHwqgp93zU6jU2pWdmTR3o0ZPXhc2yNvWB0HFEaSkH/jyB8HGz60KJrJtsKKexXuHvrxb1uJ9TBBX8+bNRQzLQtcdSp4kXvMBtYPGDDe5C0HwZ8DJVs6AIpG/FA9wYEVfXijWVHKJD1Ue3DlcV5rqyZvOoFhyruUtiv5uEDY+aj/OqW/NwKdOh0GjvjUpjQJQQXo3uun9qhD8G0GQNhg43NYqM83Vx4sX8zjp7N4Mdd0jLWbphMMOAj6PQI7JgCy59ymHnuUtj/zsMX1W7Cnx+ajq+zeoTpW+Kp5O7CiPbG/oEhNxN+ngx+QTIEU4J+LWrRob4/H/0WTVq2czaesktK6a0Huv5Lnwq59HGHuEJVCnsxTO5/dU9U694udVdIS0jOyGXp/jMMbxeEn5eb1fZbrNUvwKWTevtjmQVzQ0opXh0URurlfD6TlZbsi1Jw+2tw6/Ow9wd9hpydNw6Twl4CFdJN7wr5x1tWGYObs+MkeYVFTDC65/qRJbBnln5Jdr3OxmaxE81r+zGqfV1mbYsn9nym0XFEWSgFPV+A2/4PDvykd4UssN8Lz6SwF+Oaeey3vwZt74WN7+vNwyrwbVpuQSE/bD9Jz6bVaVDdp8L2U6KUE/qFSHXaQY8Xjcthh56+oylebi68tfyI0VHEzbjlGbjjLTjyC8wbqQ9H2iEp7MW4Zh67i5veE6XrkxD5PSx6oML+ki/dn8SFzDxje67n58D88aBMejdMV3fjstihAB8X2910AAAeEElEQVQPnujVmHXRyayLPm90HHEzujwGQ77SO5jOGgxZF41OVGZS2IvxjytPldIvoe/9BhxeDD/cBZcvWXSfmqbx3aYTNK3pS/fGBq5EtPpFOHsA7vwGqthIfxo7M75LCPUDvPnvsiPky/RH+9R2rN5P6txhmNbH7hbIlsJejGu6O1595WnXJ/UTiae2w/d3QEqcxfa56dgFjp7N4IHuBvZc3ztHf1fS5Qlo2s+YDA7A3dXES/2bcSI5i1nbThodR9ys0P5w78+QeV7//37+qNGJSk0KezGCfIL+vP+PC5Raj4Rxv+oL5n7XS5/nbQHfbjpBDV8PBrepbZHtldnJbXq/6vq3Qq9XjMngQHo1q0H3xgF89nsMKVn2exLO6dXrAvetAK1QP3KP22R0olKRwl6MDoEd/rxf7NFzSFe4/3fw9IOZg+DQonLtLyopnU3HLjC+S4gxPdcvnYSfxuhDLyNm6ksJinJRSvHKwDCy8gr5eE200XFEedRqAff/Bj41Yfad+pRIGyeFvRiuJleqelQFICMvo/gnBTTSi3udcFg4Eda8ctNzX7/ddIJK7i6M6WjAmPblSzBvFBQV6C0VvKpaP4ODalzTl7Edg5m74xRRSelGxxHlUTVEL+4h3fQZY2tetemrVMtV2JVSdyulDiulipRSEZYKZQtaVm8J3KCwA3hX04dlIibq/SZ+uAuyytYI6mxaDkv2nWFERF2qVLLyDJS8LJg7Ei4cgxGzIKCxdffvBP7duwmVvdx4dclhNAfqReKUvKrAmAV/dYacf6/+f8gGlfeI/RBwF7DRAllsytMRTxNeI5wOtTrc+ImuHjDwE3161Knt8L9b4fTuUu9nxtZ4ijSN+7tZeYpjQR7MHweJu2DYd/qCI8LiqlRy59k+TdkZl8KS/WeMjiPKy8VNb4bX5x04uhy+6w0XbW+hlXIVdk3TojRNc8gBxAZ+DZjZbyaV3CqV/GTQp0fdv1qf/z2tL2z9osS3apm5BczZcZJ+LQKp61/K/VhCfg4sGA+xv+uLZjQfar19O6FR7YNpWcePt5ZHkZlr35eqC/Tpz50fgbELIeMMTO0BUcuMTnUNGWO3pNpt4cEN0PgOfemt2UMh/fpHafN3JZCRU8AD3a14tJ6XrY+pR6+Afh9Au/HW27eTcjEp3hjSnPMZuXwufWQcR6Pb4cGNUK2RPvmgHOfZLK3Ewq6U+l0pdaiY25Cy7EgpNVkpFamUikxOTr75xLaukr9+YcOgz/Vhjild4ODCf/SZKSgs4vvNcbQPqUrbYCudsMxM1v/YxG3Qh446TrbOfgVtg6syMqIu0zbHcezcDc7bCPtSJRgmrvqrr/v0vjYxNFNiYdc07XZN01oUc/u1LDvSNG2qpmkRmqZFVK9e/eYT2wOl9CPhBzeBfwO9odCc4XAp/s+nrDp8ltOpl3mgewPrZDp3GL69TV8wY/h0fehIWNV/+jalkruLnEh1NK4eMPBjGD4NLsTAN90gcrqhC3fIUExFCmgE96+Bvu/pJ1a/6gSbPkbLy+bbjSeoH+DN7c1qVmwGTdOvKP2uNxTlw30rZUzdINV8PHi2T1O2Hr/I8oNJRscRltZiGDy8DYLaw7J/wdwR+jUiBijvdMc7lVKJQGdguVJqtWViORCTC3R6CB7dAY16wdrXyfs0nMZJS3iga3DFrpCUmQwLJujrldZuC5P+0OfdC8Pc07EezWtX5q3lUWTJiVTH41cH7v1FP5iL3wJfdYTNn1i9BbAy4i1hRESEFhkZafX92oS4TRyf9zQN86IpqtYYU9cnoOUIcPO03D4K82HX97DubcjPhtte0vu/mAy4qlX8w+6Tlxg2ZSsP92jIc31DjY4jKkpqAqx6Ho4ug4CmcMd/9YkV5egFpZTarWlaidcMyVCMlR10a0Wv9FdYHfYOJldPWPI4fNoS/ngTLsSWb+N5WbBjKnweri8OEtQOHt4K3f4tRd2GtKtXleHtgvhu0wmOJ9tnv29RClXqwqg5MPpHfRh07giYMRDOVXyvftcK34O4xpQNsfh6uNF58CTweFifobL1S33B6I0fQO1w/a96w576/ZL6oWenwMmtELVUv+VnQd2O+iK9jXuX6+hAVJzn+4Wy+vBZXltymFkTOxjX0VNUvKb99KmRu2fAhvf19h0VTIZirOhEcia9Pt7Aw7c25D9/fwuengSHFsKhxXBmL6CByVWfI1utkT6N0qMyaEX68EraaX2loxTz1CpPPwgbCm3ugeBOVv/eRNnN2BLHa0uPMGVMOP1aBhodR1hDQa4+i+YmlXYoRo7Yreh/G07g7mIqfoWkyoHQ5XH9lp0C8Zv0qYnno+BiLCSmQk6afkmzqwf4Bupd51qP1hsT1Wknqx3ZmbGd6jE/MpHXlx6he5Pq+HjIf0eHV46iXqbdWGUvgqS0yyzem8io9sFU9y3hh1vJH8KG6DfhsFxdTLx9V0vu/HoLH/0WzauDmhsdSTgIOXlqJd9viqNIg8m3WOmCJGEX2tStwr2d6jFzazwHE9OMjiMchBR2K7iUlcfcnacY3Lq2dZt9CbvwTJ+mBPh48OLPBykskitSRflJYbeCmdviyc4r5KFbGxodRdigyp5uvDIojIOn05i1Ld7oOMIBSGGvYFm5BczYGs/tzWrQtJav0XGEjRrQMpBbm1Tnw9XRJKVdNjqOsHNS2CvYD9tPkpqdz8M9GhkdRdgwpRRvDm1Boabx+pKKv4BFODYp7BUoO6+AqRtP0L1xAO3qyVqi4sbq+lfiiV6NWXX4LL8fOWd0HGHHpLBXoLk7TnExK48ne8laoqJ0JnVvQJOaPry65DDZedIkTNwcKewV5HJeId9sOEGXhtWICPE3Oo6wE24uJt6+syWnUy/zyZoYo+MIOyWFvYLM3XmKC5m5crQuyiwixJ/RHYL5fnMc+xNSjY4j7JAU9gqQk1/INxuO06mBPx0bVDM6jrBDL/QPpbqvB88tOkBewY0XRRfi76SwV4Afd54iOSOXJ+RoXdykyp5uvDW0JUfPZjBlvfFraAr7IoXdwnLyC5my4TgdQvzpLEfrohxuD6vJ4Na1+XLdMaLPygLYovSksFvY/MgEzqXn8uTtjaXHtii3VweF4evpxn8WHZB2A6LUpLBbUE5+IVPWH6ddvap0aShH66L8qvl48OqgMPYnpDJ9S5zRcYSdkMJuQXN3nCIpLYenejeRo3VhMYNb16ZXaA0+/C2akxezjI4j7IAUdgvJyi3gq3WxdG5Qja6NAoyOIxyIUoo372yBm8nE84sOYsSqZ8K+SGG3kBlb47mYlcczfZoaHUU4oEA/L17o34xtJy4yb2eC0XGEjZPCbgFp2fl8s+E4vUJrSE8YUWFGd6hLl4bVeGv5ERJSso2OI2yYFHYLmLrpOBk5BTx9hxyti4qjlOL94a1QSvHMgv0UySwZcR3lKuxKqQ+UUkeVUgeUUj8rpapYKpi9SM7IZdrmeAa2CiSsdmWj4wgHF1S1Eq8MDGNHXArTZJaMuI7yHrGvAVpomtYKiAFeKH8k+/L1+ljyCot4qncTo6MIJ3F3RBC9Qmvw/upoYs/LhUvin8pV2DVN+03TtCu9RbcDQeWPZD9Op15mzvZTDAuvQ4PqPkbHEU5CKcU7w1ri7e7CU/P3k18ovWTEtSw5xj4RWHm9TyqlJiulIpVSkcnJyRbcrXE+NbdVlZ4wwtpq+Hry5tCWHEhMk14y4h9KLOxKqd+VUoeKuQ256jkvAQXAnOttR9O0qZqmRWiaFlG9enXLpDdQVFI6C/ckMq5zPYKqVjI6jnBCA1oFMrh1bT5fe4xDp9OMjiNsiGtJT9A07fYbfV4pNR4YCPTSnOjKiXdWHsXXw5XHbpO1TIVx3hjSnO0nLvLU/H0seawbnm4uRkcSNqC8s2L6As8BgzVNc5qJtRtjktkYk8wTvRpTpZK70XGEE6tSyZ33hrci5lwm7606anQcYSPKO8b+JeALrFFK7VNKfWOBTDatsEjj7RVR1PX34t7O9YyOIwQ9m9ZgQpcQpm+JZ93R80bHETagvLNiGmmaVlfTtDbm20OWCmarFu9J5OjZDJ7tE4qHq7ztFbbh+X6hNK3py7ML95OckWt0HGEwufK0DC7nFfLRbzG0DvJjUKtAo+MI8SdPNxc+H92WjJwCuSpVSGEvi2lb4jibnsOL/ZtJW15hc5rW8uXlAc3YEJPMjK3xRscRBpLCXkrn03OYsv44tzerKQtUC5s1tlM9bm9Wg3dXHuXImXSj4wiDSGEvpfdWRZNbUMhLA5oZHUWI69IbhbWmSiU3nvhxL5fzCo2OJAwghb0U9p66xKI9idzfrQH1A7yNjiPEDfl7u/PxiDbEns/k9aWHjY4jDCCFvQRFRRqvLTlMDV8PuRhJ2I1ujQN4pEdDftyVwOI9iUbHEVYmhb0Ei/Yksj8xjef7heLjUeKFukLYjKd6N6FjfX9e+vkQx85JF0hnIoX9BjJy8nlvVTRtg6swtE0do+MIUSauLia+GN0Wbw8XHp6zh6zcgpK/SDgEKew38MUfsVzMyuW1Qc0xmWR6o7A/NSp78vmothxPzuSln2UhbGchhf06Ys9nMm1zHHe3C6J1XadbGEo4kC6NAvj37U34Zd8ZftwlC2E7AynsxdA0jZd+Poi3hyv/6RtqdBwhyu2xno3o3jiAV5cclha/TkAKezEW7TnNjrgUnu8XSoCPh9FxhCg3k0nx6cg2+Fdy5+E5u7mUlWd0JFGBpLD/zaWsPN5eEUVEvaqMjKhrdBwhLKaajwdTxoZzLi2Xx+ftpUCW1HNYUtj/5p2VUaRfzuetO1vKCVPhcNoGV+XNoS3YHHuB91dHGx1HVBCZmH2VnXEpzI9M5KFbG9K0lq/RcYSoECPa1+XQmTSmbjxB89qVGSJTeR2OHLGb5RUU8eLPBwmq6sWTsji1cHD/NzCMDiH+PLfoAIfPyMlURyOF3WzK+uPEns/kjSHN8XKXBTSEY3NzMfHVmHCqVnJn8qzdpMjJVIcihR04ejadL9cdY3Dr2twWWtPoOEJYRXVfD/53bzuSM3N56Ifd5BZIJ0hH4fSFvaCwiGcXHMDPy43XBjc3Oo4QVtUqqAofDG/FzrgUXlgsV6Y6Cqc/efq/jSc4eDqNKWPC8fd2NzqOEFY3pE0d4i9k88nvMTQI8Oax2+Qck71z6sJ+7FwGn/1+jAEtA+nXUtYwFc7riV6NiL+YxYe/xVCvmjeDWtc2OpIoB6cdiikoLOKZhQfw8XTl9SEyBCOcm1KKd4e1pH1IVZ5esJ/dJy8ZHUmUg9MW9qmbTrA/IZXXBjeXtgFCAB6uLvzv3ggC/TyZPCuSUxezjY4kblK5CrtS6r9KqQNKqX1Kqd+UUnbx/u1gYhof/xbDgJaBDGolQzBCXOHv7c60Ce0pKNIYN20HFzJzjY4kbkJ5j9g/0DStlaZpbYBlwCsWyFShLucV8uRPewnw8eCtO1uglLQNEOJqDav7MG1Ce86m5zBh+k4yZYEOu1Ouwq5pWvpVH3oDNj9X6q0VR4i7kMXHI1pTpZLMghGiOO3qVeXrMeFEJWXw4OxImeNuZ8o9xq6UeksplQCMwcaP2NdGneOH7aeY1L0BXRoFGB1HCJt2W2hN3hvWii2xF3lq/n6Kimz+uE2YlVjYlVK/K6UOFXMbAqBp2kuaptUF5gCP3WA7k5VSkUqpyOTkZMt9B6WUnJHLfxYeoFlgZZ6+o4nV9y+EPRreLogX+oWy/EASry89LBcw2YkS57FrmnZ7Kbc1F1gOvHqd7UwFpgJERERY9bejsEjjXz/tJTO3gHmj2uDhKr1ghCitB29tyIXMXL7dFIeXuyvP9W0q56ZsXLkuUFJKNdY07Zj5w8HA0fJHsrwv/jjGltiLvDesJU1qSjteIcrqxf7NyM4r5JsNx/FwNfHv3vKu15aV98rTd5VSTYEi4CTwUPkjWdaW2At8tvYYd7WtwwhZEUmIm6KU4r9DWpBfWMRna4/h7mri0Z6NjI4lrqNchV3TtGGWClIRzqfn8OSPe2lU3Yc3ZWqjEOViMineuasV+YUaH6yOxsPVxAPdGxgdSxTDYXvFFBQW8fi8vWTlFjJvUjiV3B32WxXCalxMig+GtyKvoIg3l0dhUoqJ3eobHUv8jcNWu/dWHWVHXAof3d2axjKuLoTFuLqY+HRUGwqLNN5YdoScgkIe6SHDMrbEIXvFLN6TyLeb4hjfuR7D2gUZHUcIh+PmYuLLe9oytE1t3l8VzUe/RctUSBvicEfsBxJTeX7xQTo18OflgWFGxxHCYbm6mPhoRBs83Vz44o9YsvMKeXlAMzmXZQMcqrCfz8hh8qzdVPfx4Osx7XBzccg3JELYDBeT4p27WuLp5sL3m+PIzivkzaEtcDFJcTeSwxT2nPxCHv5hD2mX81n4cGdZDUkIK1FK8eqgMLw9XPhq3XEuZuby+ei2eLrJhYBGcYhD2qIijWfMiwN8eHdrmtf2MzqSEE5FKcWzfUJ5bVAYa6LOcc+320nJyjM6ltNyiML+/upolh1I4vl+oQyQ/upCGGZC1/p8fU84h86kM3zKVlmswyB2X9hnbz/JNxuOM7ZTMA/eIhdLCGG0fi0DmfNARy5m5XHXlC3sPSXL7FmbXRf2tVHnePXXQ9wWWoPXBjWXs/FC2Ij2If4sergzXu4ujJy6nUW7E42O5FTstrDvjEvh0bl7aF7bjy9Gt8VVZsAIYVMa1fDl10e70S5YXyD7zWVHKCgsMjqWU7DLanggMZWJM3ZRu4oX0+9rj7eHw0zuEcKh+Hu7M+v+DozvXI/vNsdx34xdpGXnGx3L4dldYY85l8H4aTvx83JjzgMdCfDxMDqSEOIG3FxMvD6kBe/c1ZLtJy4y4ItN7EtINTqWQ7Orwn7yYhZjv9uBm4uJuZM6EujnZXQkIUQpje4QzPwHO6NpMHzKVr7bdELaEFQQuyrsn609Rn5hET880JF61byNjiOEKKO2wVVZ8UR3eobW4M3lUUyatZvUbJnvbmnKiL+YERERWmRkZJm/Lie/kISUbOnWKISd0zSN6VvieWdlFP7e7rw7rBU9m9YwOpbNU0rt1jQtoqTn2dURu6ebixR1IRyAMvdx//mRrvh5uXHf9F28sPgAmbkFRkdzCHZV2IUQjqVFHT+WPNaNB29twE+7EujzyUY2H7tgdCy7J4VdCGEoTzcXXujXjAUPdcbd1cTY73fwxLy9nE/PMTqa3ZLCLoSwCe3q+bPyye482asxqw6fpddHG5i+JU4uaroJUtiFEDbD082Ff/duwup/3ULbelV5fekR+n++ibVR52RqZBlIYRdC2Jz6Ad7MvK8934wNJ6+giPtnRjJy6nZpKFZKdjXdUQjhfPILi/hx5yk+/f0YF7Py6BVag0d6NqJdvapGR7M6q053VEo9o5TSlFIBltieEEJc4eZi4t7OIWz4T0+e6t2E3acuMWzKVkZN3cbGmGQZoilGuY/YlVJ1ge+AUKCdpmklzlWSI3YhxM3Kyi1g3s5TfLvpBOfSc2lQ3ZsxHesxPDwIv0puRse7IU3TytVevLRH7JYo7AuB/wK/AhFS2IUQ1pBbUMiy/Un8sOMke0+l4ulmon/LQIa0qUPXhtVsqpV3Qko2i/YksmhPIt+OiyC0VuWb2k5pC3u5+t0qpQYDpzVN2y+LXAghrMnD1YVh7YIY1i6IQ6fTmLPjJMv2J7F4z2n8vd3p37IWvcNq0bG+vyELayekZLP68Fl+O3yOnfEpKAVdGwaQV1Dx0zdLPGJXSv0O1CrmUy8BLwJ3aJqWppSK5wZH7EqpycBkgODg4HYnT54sT24hhPiHnPxCNsQks3T/GdZGnedyfiEeriY6NqhG90YBhNerQvPafhVS6C9l5bEzPoUdJ1LYevwCR89mABBay5cBLQO5q10QdaqUryNthQ/FKKVaAmuBK6vVBgFngA6app290dfKUIwQoqJdzitkR9xFNsQkszEmmePJWQC4mhTNAisTWsuX+tW9aRDgTbC/NwG+7vhXcr/hEE5BYRHnM3I5k3qZM2k5HD+fydGz6Rw9m8FJ88LdHq4mwoOr0qtZDe4Iq0VwtUoW+56sNsZ+1Q7jkTF2IYSNOp+ew96EVPYnpLIvIZXY85mcz8j9x/P8vNzwcnPB1UXhalJoQHZeIVm5BWTnFV7zXKWgfjVvQgN9CQusTMcG1WgV5IeHa8UM/VhljF0IIexFjcqe9Gleiz7N/xpZzswtIP5CFqdSsrmYlUdKZh4Xs3LJzS8iv6iIgkINpaCSuyve7i54e7hSs7IngVU8CfTzJNi/EpXcba+MWiyRpmkhltqWEEJYg4+HKy3q+NGijp/RUSzKduYDCSGEsAgp7EII4WCksAshhIORwi6EEA5GCrsQQjgYKexCCOFgpLALIYSDkcIuhBAOxpAVlJRSycDNdgELAEpsW2AAyVU2kqtsJFfZ2Wq28uSqp2la9ZKeZEhhLw+lVGRpeiVYm+QqG8lVNpKr7Gw1mzVyyVCMEEI4GCnsQgjhYOyxsE81OsB1SK6ykVxlI7nKzlazVXguuxtjF0IIcWP2eMQuhBDiBuyqsCul+iqlopVSsUqp5w3Yf7xS6qBSap9SKtL8mL9Sao1S6pj536rmx5VS6nNz1gNKqXAL5pimlDqvlDp01WNlzqGUGm9+/jGl1PgKyvWaUuq0+TXbp5Tqf9XnXjDnilZK9bnqcYv+nJVSdZVS65RSUUqpw0qpJ82PG/qa3SCXoa+ZUspTKbVTKbXfnOt18+P1lVI7zN/7T0opd/PjHuaPY82fDykpr4VzzVBKxV31erUxP261333zNl2UUnuVUsvMHxv3emmaZhc3wAU4DjQA3IH9QJiVM8QDAX977H3gefP954H3zPf7AysBBXQCdlgwxy1AOHDoZnMA/sAJ879VzferVkCu14BninlumPln6AHUN/9sXSri5wwEAuHm+75AjHn/hr5mN8hl6Gtm/r59zPfdgB3m12E+MMr8+DfAw+b7jwDfmO+PAn66Ud4KyDUDGF7M8632u2/e7lPAXGCZ+WPDXi97OmLvAMRqmnZC07Q84EdgiMGZQM8w03x/JjD0qsdnabrtQBWlVKAldqhp2kYgpZw5+gBrNE1L0TTtErAG6FsBua5nCPCjpmm5mqbFAbHoP2OL/5w1TUvSNG2P+X4GEAXUweDX7Aa5rscqr5n5+840f+hmvmnAbcBC8+N/f72uvI4LgV5KKXWDvJbOdT1W+91XSgUBA4DvzB8rDHy97Kmw1wESrvo4kRv/J6gIGvCbUmq3Umqy+bGamqYlgf4fFahhftzaecuaw5r5HjO/FZ52ZbjDqFzmt71t0Y/2bOY1+1suMPg1Mw8r7APOoxe+40CqpmkFxezjz/2bP58GVLNGLk3Trrxeb5lfr0+UUh5/z/W3/VfEz/FT4D9Akfnjahj4etlTYVfFPGbtKT1dNU0LB/oBjyqlbrnBc20hL1w/h7XyTQEaAm2AJOAjo3IppXyARcC/NE1Lv9FTrZmtmFyGv2aaphVqmtYGCEI/amx2g30Ylksp1QJ4AQgF2qMPrzxnzVxKqYHAeU3Tdl/98A32UeG57KmwJwJ1r/o4CDhjzQCapp0x/3se+Bn9F/7clSEW87/nzU+3dt6y5rBKPk3Tzpn/MxYB3/LXW0ur5lJKuaEXzzmapi02P2z4a1ZcLlt5zcxZUoH16GPUVZRSrsXs48/9mz/vhz4kZ41cfc1DWpqmabnAdKz/enUFBiul4tGHwW5DP4I37vUqz8kCa94AV/STHPX56wRRcyvu3xvwver+VvRxuQ+49gTc++b7A7j2xM1OC+cJ4dqTlGXKgX5kE4d+8qiq+b5/BeQKvOr+v9HHEAGac+2JohPoJwEt/nM2f++zgE//9rihr9kNchn6mgHVgSrm+17AJmAgsIBrTwY+Yr7/KNeeDJx/o7wVkCvwqtfzU+BdI373zdvuwV8nTw17vSxWaKxxQz/LHYM+3veSlffdwPyi7wcOX9k/+tjYWuCY+V//q37JvjJnPQhEWDDLPPS36Pnof+Xvv5kcwET0EzSxwH0VlGu2eb8HgCVcW7ReMueKBvpV1M8Z6Ib+lvYAsM9862/0a3aDXIa+ZkArYK95/4eAV676P7DT/L0vADzMj3uaP441f75BSXktnOsP8+t1CPiBv2bOWO13/6rt9uCvwm7Y6yVXngohhIOxpzF2IYQQpSCFXQghHIwUdiGEcDBS2IUQwsFIYRdCCAcjhV0IIRyMFHYhhHAwUtiFEMLB/D9toRBHXbP2qAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#chebs = getChebsFromFiles(100)\n",
    "#counter = 0\n",
    "#minNum=4001\n",
    "#baddies = []\n",
    "#for chebnum in range(len(chebs)):\n",
    "#    keys = list(chebs[chebnum][0].keys())\n",
    "#    for keyPair in range(len(keys)-1):\n",
    "#        if keys[keyPair+1]-keys[keyPair] < 17:\n",
    "#            baddies.append(chebnum)\n",
    "#            print(keys[keyPair+1]-keys[keyPair])\n",
    "#            break\n",
    "\n",
    "i=2\n",
    "for j in range(i*3, i*3+3):\n",
    "    x = baddies[j]\n",
    "#print(chebs[x])\n",
    "    evalPoints = [evaluate3DPiecewiseCheb(t, chebs[x]) for t in range(4000)]\n",
    "    ran = range(4000)\n",
    "    evalX, evalY, evalZ = [evalPoints[i][0] for i in ran], [evalPoints[i][1] for i in ran], [evalPoints[i][2] for i in ran]\n",
    "    plt.plot(evalZ)\n",
    "#plt.plot(evalY)\n",
    "#plt.plot(evalX)\n",
    "#print(chebs[x][0])\n",
    "#len(baddies)\n",
    "#baddies\n",
    "#print(chebs[x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[12393, 22357, 22517, 24518, 30984, 45537, 55445, 70481, 79337, 80593, 84011, 91647, 103809, 104353, 117199, 118243, 120378, 124732, 125107, 128582, 132494, 147951, 148905, 164998, 176452, 187853, 189389]\n",
      "21\n"
     ]
    }
   ],
   "source": [
    "#actualBadIndexes = []\n",
    "#actualBads = [3, 4, 9,10,12,21,37,42,45,54,55,64,65,67,70,72,77,84,91,93,94]\n",
    "#for i in actualBads:\n",
    "#    actualBadIndexes.append(baddies[i])\n",
    "print(baddies[:27])\n",
    "print(len(actualBadIndexes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Absolute positions to Earth-based observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Given a set of objects and a time, this returns a set of vectors from Pan-STARRs to the objects as they appear in\n",
    "#the sky at that time (speed-of-light delay factored in). \n",
    "\n",
    "def getObservationVectorAtTime(atTime, chebyshevs):\n",
    "    c = MPCL.Constants.speed_of_light\n",
    "    baseTime = 2454466.5\n",
    "    offset = 100 #100 days before for speed-of-light buffer\n",
    "    t0Base = baseTime + offset \n",
    "    t0=MPCL.EOP.jdTDB(t0Base)\n",
    "    nSamp = len(chebyshevs)\n",
    "    d = np.zeros(nSamp)\n",
    "    xs = np.zeros(nSamp)\n",
    "    ys = np.zeros(nSamp)\n",
    "    zs = np.zeros(nSamp)\n",
    "    t1Base = t0Base + atTime\n",
    "    t = MPCL.EOP.jdTDB(t1Base)\n",
    "    obsPosition = np.array(MPCL.Observatories.getObservatoryPosition('F51', t))\n",
    "    #print(obsPosition)\n",
    "    totalTime = 0\n",
    "    for objNum in range(len(chebyshevs)):\n",
    "        lightDelay = 0.0\n",
    "        for i in range(2):\n",
    "            delayedTime = t - lightDelay\n",
    "            centaurPosition = evaluate3DPiecewiseCheb(delayedTime-t0+offset, chebyshevs[objNum])\n",
    "            obsToCentaurVector = centaurPosition - obsPosition\n",
    "            d = np.linalg.norm(obsToCentaurVector)\n",
    "            lightDelay = d/c\n",
    "        x, y, z = obsToCentaurVector\n",
    "        xs[objNum] = x\n",
    "        ys[objNum] = y\n",
    "        zs[objNum] = z \n",
    "    return (xs,ys,zs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Given a set of piecewise chebs representing objects, this creates a dictionary which maps healPix to objects at a given time.\n",
    "#Returns a list of lists of objects. The nth list has the integer numbers corresponding to the objects in the nth pixel.\n",
    "#For example, if chebs is 5 objects and there were only 3 pixels, the dictionary might be: [[3], [2,4], [1,5]]. \n",
    "def getHealPixDict(time, chebs, nside):\n",
    "    nObjects = len(chebs)\n",
    "    dicts = [[] for i in range(hp.nside2npix(nside))]\n",
    "    xs, ys,zs= getObservationVectorAtTime(time, chebs)\n",
    "    for i in range(nObjects):\n",
    "        pix = hp.vec2pix(nside, xs[i], ys[i], zs[i])\n",
    "        dicts[pix].append(i)\n",
    "    return dicts\n",
    "\n",
    "\n",
    "#Given the healPixDict (with healPixes numbers as keys and an array of objects in that pixel as values), \n",
    "#this returns a list of objects in healPixes \n",
    "def getObjectsInDisc(healPixDict, nSide, vector, radius):\n",
    "    pixes = hp.query_disc(nSide, vector, radius, True)\n",
    "    obsInTheDisc = []\n",
    "    for pix in pixes:\n",
    "        obsInTheDisc.extend(healPixDict[pix])\n",
    "    return obsInTheDisc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get list of exposures from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Used to load exposures from amalgamated_unique_exposure_reads.txt into python, \n",
    "#and get a list of exposures with tuples containing (JD, ra, dec, and two parts of the filename). \n",
    "def getExposureTimesAndAnglesFromFile(filePath):\n",
    "    exposures = []\n",
    "    with open(filePath,'r') as file:\n",
    "        fullText = file.read()\n",
    "        lines = fullText.split('\\n')\n",
    "        for k in range(len(lines)-2):\n",
    "            values = lines[k+1].split(',')\n",
    "            values[1].strip()+values[2].strip()\n",
    "            for value in range(len(values)):\n",
    "                values[value] = values[value].strip()\n",
    "            if len(values) > 6 and not values[6] == 'None':\n",
    "                exposures.append((float(values[4])-2454466.5, float(values[5]), float(values[6]),values[1], values[2]))\n",
    "\n",
    "    return exposures\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate list of plausible objects "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Same as getHealPixDict, but reads from a premade file instead of from cheb variables. \n",
    "#The n'th element of what is returned is a list of objects in the n'th healpix. \n",
    "#The objects are represented by a number, so you'll need to read both the HPdict and the chebs. \n",
    "def getHealPixDictFromFile(filePath, fileNum):\n",
    "    allObjects = []\n",
    "    with open(filePath + 'HPD'+str(fileNum)+'.txt', 'r') as file:\n",
    "        fullText = file.read()[2:-2]\n",
    "        dicts = fullText.split('], [')\n",
    "        for dictStr in dicts:\n",
    "            obs = dictStr.split(', ')\n",
    "            if obs == ['']:\n",
    "                obs = []\n",
    "            for objNum in range(len(obs)):\n",
    "                obs[objNum] = int(obs[objNum])\n",
    "            allObjects.append(obs)\n",
    "        file.close()\n",
    "    return allObjects\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Searches exposures and returns the set within 12 hours of a given time. \n",
    "def getExposuresOnDay(exposures, day):\n",
    "    exposuresOnDay = []\n",
    "    for exp in exposures:\n",
    "        if round(exp[0]) == day:\n",
    "            exposuresOnDay.append(exp)\n",
    "    return exposuresOnDay\n",
    "\n",
    "\n",
    "# Reads a healpix dictionary and set of exposures in from files, then writes new files with a list of the objects which \n",
    "# might be in each exposure. \n",
    "def calculateObjectsNearExposures(exposuresFile, healPixDictPath, outputPath):\n",
    "    exposures = getExposureTimesAndAnglesFromFile(exposuresFile)\n",
    "    for day in range (1, 4000):\n",
    "        HPD = getHealPixDictFromFile(healPixDictPath, day)\n",
    "        with open('outputPath' + str(day) + '.txt', 'w') as file:\n",
    "            for exp in getExposuresOnDay(exposures, day):\n",
    "                theta = exp[1]*2*np.pi/360.\n",
    "                phi = (exp[2]+90)*2*np.pi/360.\n",
    "                x,y,z = np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)\n",
    "                generousObs = getObjectsInDisc(HPD,256, [x,y,z], .0332)\n",
    "                strictObs = []\n",
    "                for ob in generousObs:\n",
    "                    pos = evaluate3DPiecewiseCheb(exp[0], chebs[ob])\n",
    "                    pos /= np.linalg.norm(pos)\n",
    "                    if pos[0]*x + pos[1]*y + pos[2]*z > .999450:\n",
    "                        strictObs.append(ob)\n",
    "                print(len(strictObs), len(generousObs))\n",
    "                file.write(str(exp[3]) + ' ' + str(strictObs) + '\\n')\n",
    "                file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calculateObjectsNearExposures('/Users/jake/amalgamated_unique_exposure_reads.txt', '/Users/jake/Desktop/healPix/', '/Users/jake/Desktop/exposureObjects2/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get an entirely new set of objects, just call makeChebFiles. As is, when makeChebFiles is called with parameters a and b, it will create b-a files, each with 10,000 objects. To read chebs in from files, run getChebsFromFiles (will get chebs0.txt up to chebs<x-1>.txt where x is the number you pass in). Both getChebs methods return a list of objects. The HealPix dictionary created from this list will represent each number by its index in the list, so don't modify the list of objects after creating the dictionary.\n",
    "    \n",
    "    \n",
    "To make a HealPix dictionary from a set of files, use getHealPixDictFromFile. HealPix dictionary files should be named HPD<day>.txt, where day ranges from 0 to ~4,000. I used 2454466.5 (Jan 1 2008) as day 0, but I believe the first exposure is early april 2009 and the last exposure is about 3,600 days after, so about 20% of the HealPix dictionaries don't need to exist. \n",
    "    \n",
    "To get a list of likely objects for each exposure, pass an exposures file (amalgamated_unique_exposure_reads.txt), a directory with HealPix dictionaries in it, and an output directory to calculateObjectsNearExposures. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}