{ "metadata": { "name": "" }, "nbformat": 2, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Basics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "LSL provides readers for the two DP transient buffer outputs: the transient buffer - wideband (TBW) and the transient buffer - narrowband (TBW) through the lsl.reader module. First, download a small portion of a TBW data file:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# This may take a bit...", "import os", "import urllib2", "from lsl.reader import tbw", "", "if not os.path.exists('temp.tbw'):", " fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/firstlight/firstLightTBW.dat')", " fh2 = open('temp.tbw', 'wb')", " fh2.write(fh1.read(tbw.FrameSize*1000))", " fh1.close()", " fh2.close()", " ", "print \"TBW Size: %.1f kB\" % (os.path.getsize('temp.tbw')/1024.)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "TBW Size: 1195.3 kB" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To read in the TBW data:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from lsl.reader import tbw, errors", "", "fh = open('temp.tbw', 'rb')", "for i in xrange(5):", " try:", " frame = tbw.readFrame(fh)", " except errors.syncError:", " continue", " except errors.eofError:", " break", " ", " print \"DP Stand Number: %i\" % frame.parseID()", " print \"-> Frame Count: %i\" % frame.header.frameCount", " print \"-> Time tag: %.6f s\" % frame.getTime()", " print \"-> Mean value: %.3f\" % frame.data.xy.mean()", "fh.close()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "DP Stand Number: 22", "-> Frame Count: 1", "-> Time tag: 1302134255.000000 s", "-> Mean value: 12.954", "DP Stand Number: 21", "-> Frame Count: 2", "-> Time tag: 1302134255.000002 s", "-> Mean value: 14.061", "DP Stand Number: 22", "-> Frame Count: 2", "-> Time tag: 1302134255.000002 s", "-> Mean value: 8.617", "DP Stand Number: 102", "-> Frame Count: 1", "-> Time tag: 1302134255.000000 s", "-> Mean value: 7.439", "DP Stand Number: 101", "-> Frame Count: 2", "-> Time tag: 1302134255.000002 s", "-> Mean value: 3.271" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The try...except blocks help deal with problems in the file as the come up and keep things running smoothly.", "", "One of the most confusing aspects of this interface is that the stand values returned by the \"frame.parseID\" method are DP stand numbers, which are used internally in DP to keep track of signals, and not LWA1 stand numbers. However, the DP stand numbers can be related to LWA1 stand numbers using the lsl.common.stations.lwa1 instance:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from lsl.common.stations import lwa1", "from lsl.reader import tbw, errors", "", "# Get the antenna list", "antennas = lwa1.getAntennas()", "", "# Read in the first 5 frames", "fh = open('temp.tbw', 'rb')", "for i in xrange(5):", " try:", " frame = tbw.readFrame(fh)", " except errors.syncError:", " continue", " except errors.eofError:", " break", " ", " # Calculate the digitizer", " dig = 2*(frame.parseID()-1) + 1", " ", " print \"DP Stand Number: %i\" % frame.parseID()", " print \"-> LWA1 Stand: %i\" % antennas[dig-1].stand.id", " print \"-> Frame Count: %i\" % frame.header.frameCount", " print \"-> Time tag: %.6f s\" % frame.getTime()", " print \"-> Mean value: %.3f\" % frame.data.xy.mean()", "fh.close()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "DP Stand Number: 22", "-> LWA1 Stand: 85", "-> Frame Count: 1", "-> Time tag: 1302134255.000000 s", "-> Mean value: 12.954", "DP Stand Number: 21", "-> LWA1 Stand: 242", "-> Frame Count: 2", "-> Time tag: 1302134255.000002 s", "-> Mean value: 14.061", "DP Stand Number: 22", "-> LWA1 Stand: 85", "-> Frame Count: 2", "-> Time tag: 1302134255.000002 s", "-> Mean value: 8.617", "DP Stand Number: 102", "-> LWA1 Stand: 112", "-> Frame Count: 1", "-> Time tag: 1302134255.000000 s", "-> Mean value: 7.439", "DP Stand Number: 101", "-> LWA1 Stand: 187", "-> Frame Count: 2", "-> Time tag: 1302134255.000002 s", "-> Mean value: 3.271" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interface for TBN data is similar. To download a small section of data:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# This may take a bit...", "import os", "import urllib2", "from lsl.reader import tbn", "", "if not os.path.exists('temp.tbn'):", " fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/TBN/055987_000496599_DEFAULT_TBN_00000_30.dat')", " fh2 = open('temp.tbn', 'wb')", " fh2.write(fh1.read(tbn.FrameSize*1000))", " fh1.close()", " fh2.close()", " ", "print \"TBN Size: %.1f kB\" % (os.path.getsize('temp.tbn')/1024.)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "TBN Size: 1023.4 kB" ] } ], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "from lsl.reader import tbn, errors", "", "# Read in the first 5 frames", "fh = open('temp.tbn', 'rb')", "for i in xrange(5):", " try:", " frame = tbn.readFrame(fh)", " except errors.syncError:", " continue", " except errors.eofError:", " break", " ", " print \"DP Stand Number: %i, pol. %i\" % frame.parseID()", " print \"-> Time tag: %.6f s\" % frame.getTime()", " print \"-> Mean value: %.3f + %.3f j\" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)", "fh.close()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "DP Stand Number: 101, pol. 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.254 + 0.369 j", "DP Stand Number: 191, pol. 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: 0.064 + -0.055 j", "DP Stand Number: 151, pol. 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.371 + 0.217 j", "DP Stand Number: 81, pol. 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.188 + -0.326 j", "DP Stand Number: 171, pol. 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.158 + 0.129 j" ] } ], "prompt_number": 29 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similar to the TBW case, the stand and polarizations returned by \"frame.parseID\" are from the interal DP structure. To relate this to LWA1 stand/polarizations you use a similar method:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from lsl.common.stations import lwa1", "from lsl.reader import tbn, errors", "", "# Download a small bit of TBN data (1000 frames to be exact)", "# This may take a bit...", "import os", "import urllib2", "if not os.path.exists('temp.tbn'):", " fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/TBN/055987_000496599_DEFAULT_TBN_00000_30.dat')", " fh2 = open('temp.tbn', 'wb')", " fh2.write(fh1.read(tbn.FrameSize*1000))", " fh1.close()", " fh2.close()", " ", "# Get the antenna list", "antennas = lwa1.getAntennas()", "", "# Read in the first 5 frames", "fh = open('temp.tbn', 'rb')", "for i in xrange(5):", " try:", " frame = tbn.readFrame(fh)", " except errors.syncError:", " continue", " except errors.eofError:", " break", " ", " # Calculate the digitizer", " stand,pol = frame.parseID()", " dig = 2*(stand-1) + pol + 1", " ", " print \"DP Stand Number: %i, pol. %i\" % (stand,pol)", " print \"-> LWA1 Stand, pol.: %i, %i\" % (antennas[dig-1].stand.id, antennas[dig-1].pol)", " print \"-> Time tag: %.6f s\" % frame.getTime()", " print \"-> Mean value: %.3f + %.3f j\" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)", "fh.close()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "DP Stand Number: 101, pol. 0", "-> LWA1 Stand, pol.: 187, 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.254 + 0.369 j", "DP Stand Number: 191, pol. 0", "-> LWA1 Stand, pol.: 13, 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: 0.064 + -0.055 j", "DP Stand Number: 151, pol. 0", "-> LWA1 Stand, pol.: 90, 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.371 + 0.217 j", "DP Stand Number: 81, pol. 0", "-> LWA1 Stand, pol.: 125, 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.188 + -0.326 j", "DP Stand Number: 171, pol. 0", "-> LWA1 Stand, pol.: 130, 0", "-> Time tag: 1330573596.004160 s", "-> Mean value: -0.158 + 0.129 j" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Generating and Plotting Spectra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although computing the mean I/Q value for each frame isn't particularlly interesting this example shows how data can be extracted from within each TBN frame. A more complicated example that generates spectra would look like:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline", "from lsl.misc.mathutil import to_dB", "from matplotlib import pyplot as plt", "", "fh = open('temp.tbn', 'rb')", "srate = tbn.getSampleRate(fh) # Data sample rate in Hz", "frame = tbn.readFrame(fh)", "cFreq = frame.getCentralFreq() # Data tuning in Hz", "frame.data.iq.shape = (1,frame.data.iq.size) # SpecMaster expects 2-D data", "", "from lsl.correlator import fx as fxc", "freq, spec = fxc.SpecMaster(frame.data.iq, LFFT=512, SampleRate=srate, CentralFreq=cFreq) ", "print freq.shape, spec.shape", "", "fig = plt.figure()", "ax = fig.gca()", "ax.plot((freq-freq.mean())/1e3, to_dB(spec[0,:]))", "ax.set_xlabel('Relative Freq [kHz]')", "ax.set_ylabel('PSD [arb. dB]')", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(511,) (1, 511)" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4FFW6/7+ddHaWgEAiRGTfMQRZlEEnDATcQBRlBERcRhl1cBi9DApzR5i5CO4/N/TOXBUcHcWRUdBBBNS4I/saRZRFDPsWIEDW+v1xfFOnTp+qrk530tv7eZ5+uruquupUddX5nnc553gMwzDAMAzDMBoSwl0AhmEYJnJhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbAmbSJw9exb9+/dHr1690K1bNzzwwAMAgKNHj6KgoACdOnXC0KFDcfz48XAVkWEYJu7xhLOfxOnTp5Geno7KykoMHDgQjz32GBYvXoxmzZrhj3/8Ix5++GEcO3YMc+bMCVcRGYZh4pqwupvS09MBAOXl5aiqqkKTJk2wePFiTJgwAQAwYcIEvPPOO+EsIsMwTFwTVpGorq5Gr169kJWVhUGDBqF79+44cOAAsrKyAABZWVk4cOBAOIvIMAwT13jDefCEhARs2LABJSUlGDZsGD7++GPLeo/HA4/HE6bSMQzDMGEVCaJx48a48sorsXbtWmRlZWH//v3Izs7Gvn370KJFC5/tO3TogB9++CEMJWUYhole2rdvj++//z6g34TN3XT48OGazKUzZ85g+fLlyMvLw4gRIzB//nwAwPz58zFy5Eif3/7www8wDCNmXw8++GDYy8Dnx+cXj+cXy+dmGEatGtdhsyT27duHCRMmoLq6GtXV1Rg/fjwGDx6MvLw8jB49Gi+++CLatGmDN998M1xFZBiGiXvCJhI9e/bEunXrfJY3bdoUK1asCEOJGIZhGBXucR2B5Ofnh7sIdQqfX3QTy+cXy+dWW8Lama62eDweRGGxGYZhwkpt6k62JBiGYRhbWCQYhmEYW1gkGIZhGFtYJBiGYRhbWCQYhmEYW1gkGIZhGFtYJBiGYRhbWCQYhmEYW1gkmJjl178Gbr013KVgmOiGRYLxYdo04KOPwl2K4HnzTeC118JditjgxAlxXzDxB4sE48Ps2cDTT4e7FIHzn/8AN9xgXVZdHZ6yxBrFxcA//xnuUjDhgEWC0ZKcHO4SBM4//gEsWGBdxiIRGsrLgYqKcJeCCQcsEoyWaBQJ3Uy3PA5kaGCRiF9YJBgtSUnhLkHgJGjuZhaJ0FBRwSIRr7BIMFqiUSR0lgQTGtiSiF9YJBgt0SgSqiWhsyyY2lFXIvH440BZWej3y4QOfowYfPstsH69dVksxCTYsggdoRCJFSuARo2sy/7rv3zvPQCorAROngzueExoYJFg8O9/+/YniCSRWL4c2LfP/3bhtiQMA/jxx/o9Zn2wYQMwfLg4v6qq2u+nqEhf8eviRm+/DUyaVPtjMaGDRYLRBiUjyd00dCjwwAP+t1Mth/oWiX/9Czj//NDs69VXga1bQ7OvYFmxwvwcjDXRvLl4V4VCJxKnTokXE35YJJiIFwkASEz0v024LYlvvw1+H2fPivfx490Jo0xZmXDThBqv1/wcjEg0bizeGzWy7kcnEuXl4sWEHxYJBuXlvpVLJLmbAKtI3HwzcPy47zbhtiTcuMR0lJWZlWZaGvDFF+JzoOm7bdsCEybUrgxOyA2GYERCPp/ly/XLCc6mihxYJBitJRHJIjF/vvBvq4Tbkti7t3a/69YNGDXK/E5iE6hI7NsHbNxYuzI4ESpLQo5n7NplftadZ0UFWxKRgtf/JkysI4sEPciRlj5KIkF+ap37qT6zm44cAZo1s1Zwhw6J9/LywER2xw6gpMT8Tte+NkOKeOvgiQ6VJUHn06gRUFpqLmdLIrKJsKqACQfyA0k567Xpqbx/vxgIri4gUaBW9qlTotwLF5rb1KUl8eSTYuhxgiq5zz8H3n1XfD5yxCxboJw5Y352I26vvgq88ILv8mBiSYYBHD5sfk9MFOUKpSVxzTXAH/7gLBL33GN1wTHhhUWCQUWFGZMgkahNK7Z3b6Br19CVCzArkKoqIC/PdOmcPAl8/z0wZYq5bbAxicpK+8DvF1+IoccJuj6XXw6MGCE+Hz4MpKTUTiQoYA2Y5+Ek1HfdBdx5p+/yYCyJf/zDzEACxDkePx46kaiuFsKTkeEsEs88A5w+ze6mSIFFgrG4m6iyqo1I7NsX+g5QVFGcOCHy9ffsEd9PnRJlLSsT/u3du925m06fBgYOFJ8vvhg4cMBc99BDwLPPis+LFlldRt26iXfanq4XXafKSuEyat06cJHweq3Xm8TNSSTsLIZgREJnBVZWWssRTMVdVSXOLSPDeo3k/dN1OHOGLYlIIWwisWfPHgwaNAjdu3dHjx498PTPExgcPXoUBQUF6NSpE4YOHYrjujQWJqTo3E21HWI71Kmz1OI8fVq8//CDeD95UpS1rExU9h06uHM3lZYCmzaJzytXAsuWmeuOHxdiBABff22tpKhyJBGkdVTBHTsmUjwzM0UFeOIEcPSo/XkdPGhe49RU6zo3loTdddYtNwzgrbfs90XI/zkdW7Wu7Crukyf9i6MbS4JiYmfOsCURKYRNJJKSkvDkk09i69atWLlyJZ577jl88803mDNnDgoKCvDdd99h8ODBmDNnTriKGDfo3E0//CBiDIESapGgiocq54MHze9nz4qKJCNDlN+Nu6mqypplI2dJyRaVKpJkYVHsQN3u8GERyG7QQJStoADo1Mn+vLKygN/9TnxOS7OuC7VI/PQTcP319vsi5OtCn9W+F3YiMWiQcDf62z9ZEqWl5rWTj0vHOnuWLYlIIWwikZ2djV69egEAGjRogK5du6K4uBiLFy/GhJ+TvSdMmIB33nknXEWMG+TKkUTi//5PPPiBQi6akhLgjTeCLxu1OEksSCRkdxP1cnZjSVRWWgXgm2/Mz2RRTZ5sDSQD9iJB12v3blMkTp0SHesokG0HzfRGIkGiQCLhZM3ZZU/p3E3+ykHYVdZuROLECWD7dv/7ly0J2VWnloEticghImISu3btwvr169G/f38cOHAAWVlZAICsrCwckJ3GTJ0gu5vkB5NSOmV++sm50klJAdq3F/79MWOs64qLgXvvFRVcZaU1W8gOEgd6pzLJlkTr1mKZm5iEaknIbo/ycuDDD4GnnjKzqGgfJBLk9lIry127RGoniYS/VnByspn2Su4mKheJRSCWBG2risSOHSLgL+/fDlmUaFu3ItG5s/O+af+yJUHH0FkwHJOIHMIuEqdOncKoUaPw1FNPoWHDhpZ1Ho8HHh7Ks87R9ZMA9EM4P/qoeP30k35fycmiYqIAs8xf/iJSSSsqRCX/5pv6ivDbb80hLqjlLotEZqYZkwCEMAHAd9+J95MnhUi5cTfJnysqgDVrzM8EWS1yedQK7MgRUfm5FYnsbPP4JBIk0G6G1lBFgsq1eDEge2jl/9BfmYKxJM45R7yTiNrtnyyJU6f8WxIsEpFBWDvTVVRUYNSoURg/fjxGjhwJQFgP+/fvR3Z2Nvbt24cWLVpofztjxoyaz/n5+cjPz6+HEkce9ODJrFwJXHSR+33IMQn5gZXTMokzZ4CnnwYeflhfwVNLllwoJ08CpP1ygJb2XVbmG7jt1UssNwxRtqQka0yCRIL2QdbABx+I959+Al56yV4kqqutqbWEbEXJFdTGjWZlaycSR48C6emmSNhV9KdPi2PSdTp40Dx/2if91o0l4fEIQW3a1Fz3wAPA/fdb90X7V6+1jE4k3MYk6PqcPi2ugw41cO3PkqD/429/E/dEv372ZXdDdbW4XvHU7iwsLERhYWFQ+wibSBiGgdtuuw3dunXD5MmTa5aPGDEC8+fPx9SpUzF//vwa8VCRRSKe8XpFf4H27cX36mqR7UNBQkA8bF6vfb8B2d0kVwjy52eeAUaOtLZMKyt93RtUYVPsYN8+8VAahl4kzpzxrbjkclZUCDcOWRLHjgHnnSdcNbILqGdPYPNm8Z2sBbuYhPxuJwzy8m3bxLGaNLF3Nx05IirHhg3NsuriA1dfLUQnPV3sr7g4OJEAxP/fpo35XY5XyOX01zKX3U2BWhIknmosR8YucK2zJOTA9cSJoj/KkiXO5fdHaiowdSrw178Gt59oQm1Az5w5M+B9hM3d9MUXX+DVV1/Fxx9/jLy8POTl5WHp0qW4//77sXz5cnTq1AkfffQR7qcmEeMDtYhk14/q1wZEx6t//9t+P3buJpl77hHBbFkkdNnJ5Gcnn/6+fUB+PtC9u1UMnCoVuQIsLxciQdtVVYkA8bFjVpEglxNtYycSctaO/E7H0n0md1OTJs7uprQ0M7tJPQ/iwAHhMquuFudx5IivSOj+QxVZCBIThZusbVvxXc6WKi83y3HffXoLZ8kS37ki7GISZ87o96G643SQJdGkifj/aD92lkRZmZhXAjBHkAWEy/LYMd/9V1Y69/ivqNBPcMQ4EzZLYuDAgai2Sd9YIQ9gz2iRKxA5+EoPHj2QgEhlpfx/HeRuuvNO4Be/sN8uMdG/j5vWkwWxebMQsQMHrJWXU6WitoSVUBWaNRPZRLJIyL+prBSVjSwchJNIyOcjL3crEhdcIESC+nLoMpCys8U1qa4WlsZll5nrqF+F/B/aIQuQ1wusXi06Ce7c6SsS/fuLONG8ecCsWUDLltZ9XXmlWO/kbqJzGzMGeO454LPPrPsIxJJo0EC4x3bssB6LtpH3c+214l2e0e7ee0UM5KabrPt/5BFg+vTaDSnD2BP2wDVTO9yIBFFS4uxqIHfTCy8AH39sv11CgjVOQfs8dUo8uDL79wO33SY6cdEDLldsgYiE6uMmS0L2g+tEwi67ST6+fD52rid/IpGaarqbKCahni8gZq37OXEPVVW+68kilN1N8v9cVgZ88olwf8mZZ16vEOWOHcV3WSQqKsS1oWPZuRyPHTOvzY4dZs9ysiSuuw5o104s+/xz39+7sSTk+FmPHqZ70EkkCNmSsDsPXbKESjzFI0IFi0SUIotAsCIhu5ucxjvyevUt7717xYBzavmGDhX9EEgk5DI4tTypwu/fX6SWZmRY159zjnB10W9VkSgrs7qbtm8H3nlHuIFICOg8dKKnfqZMKjUmQZVx48ZWkSCXm3wtt28X/TkoD0MnEuQqoYry88/F3BnEq68K193/+3/Ali3m8qQk8b83aCC+k/tqwgRhYSQn248sSyJ09Ki5rn17Mz2ZRMLrde4oeeaMSCiQ/0+1x7ls3XbrZo0hEXYioc6NrbtPdYkWKpE2unE0wJcsSpFbmPJwCDqROHHCOa1Szm5yeojs3E2lpfrUxxYtRAVLrUD5t04tT3ITbdwoWt+qJZGWJios6q9RWmp1LZWXW0Xi1lvF6KONGgFz55plSUiwj0no3E1Nm1otCaqUVZGgVrh8bnR9SViOHbMGtjMyTJGQ/yud51XnOqysNMWUjvHKKyITLSlJHySm8wBEooFcWdM9JYuELMTkKiLU6/P112ZaLCH/J+3aObub1AqfrjWhswjciARbEoHDIhGlqJYCEYy7CXAWCbVSlUVCtmaIpk3F9lR5yZUwPdDHjwNffSU6sBFUGZWVicpLtSQSEkSrnoLjdpYEVQgUIwDMa1VWJioeO0uCypqUZO9uooqLWrlpaSJ+cuCAqCBPnzbFnFrhVAlWVlpb5uee6+tukrcHTFFRB1Gk8ZVUkQDENUpONu8H9T6g7+PGWa8TWSMUk/B6rULcvr0QAsLOHWcYorwUrCdLonVrUyR0loS/viK6+1TXr4cJHhaJKEUWATnLSBUJwxAtz9q6m6ZONT/buZvsOlBRpUmVtSwSVJlcfz0wYIAYCoOQK/xDh3wtCcqQobGl1Owm1d0kTytKqblnz4qy+ctuatLEKhKnT4v9TZxotSQAqyXRuLE4vjoQoGxdyCLRsqUpEkuXmsvlypIqWFUkZs+2BqzVgH1ysnXAPhn5nOX/kURCtiTU/coW7JkzVkuCBO3ECWD4cCAnx/qfnH++KLNaJrvsOtVNVlt30wcfWMWN8Q+LRJSiE4mePUVPY3k95aPrWmbHj5tWhN0D9sor5mcnd5OO1FSrL7msTIhOdrZ5PN34PHLlefCgXiTOOcd0z/iLScj8+KO5TWqq2I6ujSykVL6mTa3uptJSM3BrJxKy9aRaB3IFL7ubOnY04wzvv28ul/83O0vigw9EzKNJEyG6uvnK/VkSdE0IeTgSO5EARJmpYyNZEh4P8J//iPX79ome8OXlVkvi/PNNkdFlValUVVnTXs+eFUkR6jJ/nD3rm5nFOMMiEaXIIkGVxpYtZgVD68l/XVEhKlXZpdCkCXD33aYPH/CttPPyzFao7G5q2NC/SJAlIVeQaWniWE4PtFzh24lEixamJWEYzjEJ6o/ZoIE5aVFZmRCj1FTx+YknRKc0QCyjcyJLoqxMiNvx42b5KDVXFgkSxYwMMRkR9fmkayC3vmUxvOgiU8BkdCJhl86clCRmfVP/w6QkvSXx3HNiVr1zzxUNDPl3VPk6iYRhiN/NmiXEQXY3rV4t3vfts/b/oP8kM9N6jgcOiAaOkyXRtKlZ/rfeEoIo9zR3IxJUDsY9LBJRBqUqyiLx+edmRUiVAa0nH3xFhciM6dBBZJ1QQNRu3gSiTRtrSqYsErTcnyVBZaAOZ9XV7kTC4xFlVfsbkEgAooKTf0PnIMckmjUT79TZDBDHl0XivvvMdRkZZvmaNBHiVl0tZm2Ty6PGJNLTTcFITQV++1vfYLRsBcgi0b27/lq4sSTk9cnJvv+hnSWxfj2wdatY7/X6Fwm1ZzztMzFR3BsNG5oiQWXct8/cr24IGTrH2bNFHwe7ClzteKe7f9yKhJuxsRgTFokoY8IE0ZdBFoniYjGTGmCNRQBmq7Oy0nyIZswQ8x0A4iGW3TJyRXHFFdYKvaLCFIlGjfzHJFJSRMVBIlFSYloSTvn0VAnTUBOqSCQkmCLRoYPvNmVl5oijgCkScquzrMxsHauVi2y5NG0qrAdycR09au5XHb31/PPNSrCyUpSRYiC0jWxJyO4mu6G/5QqNRMWfSKjzQ9vFJCoqxH+XlOQbb6LtTp92djcBopNedbX4b+l/pfvuyBHr3Bs6kaiqMi1gWSTkqWnpvlSHaJdhS6JuYJGIMk6fFg+iXW9cJ0uCHnJZCA4ftlZW8joa0pseSPmzLBJ2loTXa7UkTpywupvs0hGpnDRftpqfT4FrwLQO1MC1DImEHB8hS+LUKWDBAuv2cjZVgwbimicmCsE4etTcP13rfv1EZzPqKEdlkEXCX0xCV3nSMahSo+tld70TE8V1KC8X9wjtU06BVa1GEonERPO/l8vlJBIkBHv3AjfeKLahoYHoPEtKrJaE3XhaNIKvXIHLSRPqUPZuLYl163x7YLMlERgsElEGjUtkN7CjG5GgSo4qWrlSKCszW7UpKVZhkC0Jp5iEXPnbWRJnz/p2kCLo+F26WL8TiYnmsalCVy0J+RroRIIsCV1vcdmSSE0VFVNionAlnThhWk50ra+6CvjXv6z7KCsTxz14UAjZ2rViuZ27KTHRty8AYO3B7a8FLLubKEZAy+0sidJSX3eTHC9wEondu83PV15pnS+cBGTtWrPcdpbE4cPiPTvbeo7y/0BlI0tFJ5Q66+LCC4UbS74n2ZIIDBaJKKOqSgRsr7vOeTs1cF1Z6WtJUOWpWhJUeSUnWzOfzp419+skEnJrUZ70/vhxMyZRWuo7JpP6e5rIRmdJ9Opl7eClBq4Bs2KkTl2qSMj7lcssV04kal6vOC71g5DRuYpIhADRa1wnEqoloY6K/+abQiQyM8X/GIhInD5tCmhVlXktdu0SwWrA3pKguEqjRs4iIQ8smZoq5rGg3548Ka4pTSyZlGRvSWzZIvpNVFdbz1EdagXwnV9Exm4mO3VQPxaJwGCRiDIqK52nddRZEuQaoodO7igGWCsreTuyJMjvT63OiRPFOruYhFwRJCRYy0KV7qFDZtBZhbZ3Eolf/cpqHeksCboWTpZEYaE56Bwhu5uo7wG1gJs2NftdUDl1w1XI/1GnTmafALlyk1vVCQm+08Vef71ZrhMnnAf8A6wxCXleh4oK87fLlolhPWi5LiZBx2zYEPjyS1HR6wLXcv8cij/R4IGVlab77R//MFONdZbE/v1Aq1bmmFt0PeRt6V4jkZA7kNqhG/5e951xhkUiyqiqsvqV1TkL5MnlZ88Wrb1zztG7m+R9UMWuczcRpaWi0nzhBVGxvPmmSIF0siTkzxTcrKoSrfFWray/+/574RaorhZxApq9TRe4JuicdJYEXQudJUExiV/+UpRJvqaqJQGYFVbjxmacoaIC+POf9SJBld2JE8Dvfy9a8IBVUOXzSEwU11WeVU5myRJg1CjrMvW/p9Y+uZuo7LJIFBebllB5ufjvSCTUxoN8TXWWhFxRk4DIZSKRaNVKXNMTJ/SWxJEj4j+SRUL9z6lsdP10Q4WrMS66L9VYBVsSgRHWmemYwFFFIiPD+rBS63nJEmDaNNGbmR5A1ZKgh8rrFS3kw4d93U2ySJw6ZVYUSUmmK4Eq+44dRacuO5EATHfTgQPCZSSzfbt5jgkJpjtKZ0kQdE46S4LKnp4uyq2zJADh5qGWPmC1JNTKLyXFdBmdPWsGamU2bDCFpmFD0dtY9t8Tqkio4yPJ5/C3v/n+Pi3N130lu5uoDJWVVpEga4a2a9LEek3ps1zhk7tNRrYkdCIhpyk3aiQqdp0lcfSo2ZCxEwnVktAhX8/bbtP3cgdYJAKFRSLKqKqyujLsRILy8w8fFoOpyW4k1R3j9Yr9+BOJ0lKzMpAr7iNHRFpu06ZAbq61IlBbd6mpYtn+/b7zGqhzYZDbQxe4JnTuJmo5yoPqpaZaK38nkZAr+JQUq+sjJcWM8wwdCi25udbvLVvq5zjQVcx2rhF5GOyGDYU49Owp3EEEjdSqEwn5vqDKWHU3EVTZyss8HmdLgtbpLInsbOHyO3DAPotLtSTUhoEauNYh32svvWR+Vn/D7qbAYHdTlKHGJNTeyNRiJCGgweZkl4PaY9frFa1dQOxbdjfJVktJiVnRypXBkSNiSlESECdLgnL2T5/2HSWUUh6pjwOJhOpW8ScSsjvoiSeEeKWk+E56RBURZQG99pp4p3Ns3twM6qoi8ec/W8ebcsJuzmfVkgDsRYLOCTCv2yuviEmACHIber2i1U7nUVFhikR5ufhMDQI5cK2WRb7uR4+a/y9VxseP+/YZkX9D5WzcWAjG3r32A0j6czcFaknI11yNmbElERgsElGGzt0kQ5UBiURJiagk5fRVSjkkvF4RA+jWzd6SuOwyUTnSwye39A4fFg85VS46kaB3cvk0berbWvzmG/McExKsc3TL6NxNJBZpaebwFhUVYr5vQFRicuBVtiTonUSJxOScc8xWtuxuOnHCvmOZDruOcm5EQjcwI3UKTE+3dhCkMlZWAmPHmuInu5sAUckfPOibAquWSy6fHJOg9+PHzWPoRKJtW5G55PGYImFnSdB+6HzVe8Nfx021vJSsALAlESwsElGGKhJqbr1qSQDiAZT7O5CbgFpUXq9oNV9wga8lUVkpKoDLLzc7wwHWh/jMGdFapArgF7/wbXVSxUJ+arViUs9BrkzUwKMucE1lPu880zVTUWFuq1oScgosHYv2QdcwLc0sp2pJBCISdts6iQSJvU4k5NRl+TrR9RwyRLyTgGRnW91dnTuLMbxkS0LXsY+O/cILYj4OOg/6b0tKzD4VOndTSoo53EhWln0KLCDuH6/X/K+DtSScRELut+EkOoyARSLK0MUkZGS3AtG4sXjIaJncwQ4wH+yEBKslQSJRXi5arSUlekuiSRMhBlS53HGH+WDSg0vf5Wk07URCrUzUbCwnS6J1a+tc0SRSqam+IkHHp/0lJophT6iiTkoyXTG0bXJy6CwJNQWWyixz/vm+vyM3jmzhyPt78EHx3qSJsBgeeshqSeTnA2vWWFvturLQ+pYtrTEJeq+qcrYk5GtEjQM7S6JRI2sarl12k1zh2wlOdbU7kXj6ad/nh/GFRSLKqKwMzN0EmCJBs7HJ+wLsRYJy7gFRCdi5m6ii0Pmy7R5kueJVUQeCowqC9uUUkzjvPHOdbEl06WLtlyHHJOT9zptnVioZGcJSC9aSCMbdVFgoLDwZshDk1GX6DljHqmreXBx/+HBzOxIJOeVV91+onetSUsS28jHpv9dZErogtt39kJpqFQm1o6PO3UTZb4YhhvUgIdy71+paVK0xusaUTcc4wyIRZajuJsoWInTupsxM092kDtEMWFvSsrspOdnqjjhzRi8S5PIKRCScLAm5cl+xAvjNb8RnOk8nkZBbkLIl8eabIkWX0FkSdEzKr3/lFdHBTQ1cy+9uqI27icjMNPuLELIlIbuR6HyaNxfvcrxi4ULRXwMQAycePWqKgOr6kyeJKioCLrlEfM/IEIIh3290P8mWBJ2XnDFGsSidJfHXv4oh6ZOS9O6mpCS9u4nuO8MQLjTqAX7vvcDixb7HIciS0GWcMb6wSISJ+fP1Pad37RJ9HOxQRSIpyfpA6URCtiSo5QfoLQm1Mx1N6EPb6ERCbtmq6+wG8ZMrXpWyMnPd4MHmMXWWhNpPQnYpyWKjO4adSFA/AoqzqIFr+d0NbtxNdiIB+AZxe/Qwf0MVncdjlr9JE3OKVxm5lS5bpKq7iSgvF4Ms0n/YoIFILlAtCdkqJKvrkUeA0aPN7dROiYC5nz/9SQifnSVx/vl6kSArWq3s7UbJJdS4D+MMi0SYmDbN7MsgM3myGCzNDjUFlkSCWlW0rqzMrMhkkZAtCeo7IVeW1dXWmASJBC3TpXO6cTc99JB1DB0nS0J1i6n70lkScnYTUVlpLxKygKgiMW6cyA4igrUk7MRQZ0n89rfA//yPdTu5wpw0SQxJAlgtCfWan3OO1ZIAzFa612u9j1R3E4mCXdBctSTUntlerxjmm4ZVAXyFXv1Mv1VjEocPA3/8o7NIvPee/jz/93/FxE8qn30GLF++XYXAAAAgAElEQVTOIuEWFokwQZPXq9hVaoCY3+HQIb0lQZU/7fPsWbPCzMiwdzelpVktCdonYLqbdJaEPMidG5HIybH2sHaKSfgTCafsJlXE7CwZWSRUsRg50uwzQeejioQ6jpETdmXQiUTXrmLUUhnVPy//XzqRAMT4WjTmFkEVcGKi1ZJITra3JPyRmWm9FnYWopMlQehE4pxzrBMpyTGJvDzxThNuESQSOTm+1wAQwjN0qHntNmzwTQtnTFgkwoQagCbowbnjDt91NLSDzpKQK3/A6rIhn64sEmlp9iJhF5MAzEpYtoLUzBWdu0mtKP1ZErqKxsndpLMk5N+oOFkSKjp3EwVzg8GpVS0jt3gTE/XuK/Va/vWvvmVULQlZJOxiEjrk/zI7W29JqARqSdgFrmVLYuZM/fhVdJ5JSc4WH13XvDzg9tvtt4t3WCTChD+R+PvffdfRQys/vL1724uEbBmQSMj9HKgPhOqbVwd4k91N9PvzzxcP/iefmD2lnSwJtUIItbtJF5MA7FvxsivKn0jo3E3q9a4NTq1qGbmHcEKCOEeKm9hZEjqoApYrUkBkCaki8dFH5hS3KnJZe/e2WhK6MZ4A95aELnAtD2Mui4TuOMnJ5ja6sbBkZPF1M6psvMIiESZUkXjkEeFDtfNfA75TOLZvD9x8s71IyIHkM2esDwUFwP1ZEvSuWhJz5oget5de6i67KRBL4uxZ3+tw/fXArbdajwOIStvj8S2fenwVnSVhd+117qZQWxJ2YgZY+zhQOdSgbaAiIbtt5I6QxKBBIlVWh1zWLl3EaMNEIJaEesxjx8yBDOXK3c6SULO7AGuGlG7kWhn5urqd+jQe4QH+woQak5g6VQQaL7vM/jeqSBD+3E0pKWJ8JXnK0fJy8fCpgWvA15LQuZvkYTMInUhQhaLbNhBL4s03xXAbTz3l626S0y7dWhK6mISbPh3BWBIej7VSc7IeZHQiQdRWJOTKVhUJJ8GS1x86JK7Hr39trgvGkpArajf9JHTnrFoSumHcCfm/0M1qxwjCaknceuutyMrKQs+ePWuWHT16FAUFBejUqROGDh2K4/J4xFHMxo2+UyjSjU83a6NGzhWHPI2o/LtbbjHTIuVt5TkQqqvNLCfaB7X67ALXOneTehwZnSDUxpKwi0mo5aQyyu4gNaBsdz1JJAF3MQm1PKoYuUF3Ddwgi4STkPmDWuc6S8LfhEYyVAa5TwphZ0m4EQm58aO6m2idPDilnUjIMYkjR+y3lUWCLQl7wioSt9xyC5YuXWpZNmfOHBQUFOC7777D4MGDMcduFpYogwadA8zJ7enGJx2kKR/l7QCRsrd4sSkSakDxvvt8ZzUrKzPneUhKEu6JRo2sv6WWlmpJ6NxNF14oUlj79AnsvJ1EQq3Y/vAHIUJ2/Rt0bqHUVGuev7pPJ0tC/Y2bmARVJv5a2zrsLC9/hMqSeOMNYNs2321VkXBrSeiwsySokeEkeHJr3s6S2LHDeix/7ibqYKf7b1kk3BFWkbjkkkvQROnxs3jxYkyYMAEAMGHCBLxDM9tEObqJ2MndVFEhesmWllpvVnporrkGuPpq04zWBbz79RMjucq//a//EvM2AKIikN1NQGCWREKC7yRBduh87W4C1wkJvkNlyOiEIDNTDF1hV9G7iUnoAuIy8jVyGmDOH3VhSQQiEs2bi6lU5Y5vQOhFQlcW+o0sBHbnAtgHrk+eFPOjyOWXkdNlvV4RWO/cmUUiGCIucH3gwAFk/TzQS1ZWFg6os85HKfKDJYsDYI5V1LixtbcofabMC3XwPrUVpY53lJpqjpmTmSlEoksXcxuqqO0C13JMIhDkc7WzJHQxCTkArauw7ayF3r3tRcJNTCIQd1MoRMLf8VScLAnCjUio21Lwu1GjwOZYaNvWed9OFpKTSMjL5TiHbEkAImEDsPY4J+R7NSlJWKfffutfJOxiEvffr880jCciOnDt8XjgsXnKZ8yYUfM5Pz8f+XapGBGCfJPaiURGhq9ING/u2+nOLn9dfjirqqwPTGamEKHHHxfZKOnpotLSuZvUwdoCqYAA38wU+V1erhMJJ0vCjYD4syQ+/lic/1df1S4FduBA6xhQgSCLhNOw2SqhsiQIOUOqpCRwS+Ktt+xb3naWBCFXxnZi0rKlGG7ebrs2bcS77vrJ954uy07GLrvp5puB224D+vcHHn5YpHtHaz+KwsJCFBYWBrWPiBOJrKws7N+/H9nZ2di3bx9aUE8tBVkkogG6Se+7z3eugMpKcUNnZpr57wDwf/8nhrMARODPyd0E+D5MskiQuykx0QwiUstdrSRJJPylhurYssU6Emsgnelkd1MgQuC0Tv2eny8Gu/vss8BSYKmsQ4aIEUdrg2pJ1HdMglAztdLTA7MkGjTwncdE3newloScUqvbTrZknCwJ+ZrIz5Xut7JIzJ8vXg8/LL5H80x2agN6pm5Sdj9EnLtpxIgRmD9/PgBg/vz5GKn2uY9SqIJ44gkxjj1gWgh27ib5YdENfezkbgJ8LQkaiZOgSlltodOxAnWLAOYkM/IxdPvQuZtk4XAbkyDsWr+65ZQtE0gKbCBCaUdt3U1qZzodNGy2G+j4zZuLMYw8HusxahOUJ/xZEpRMAejPpXVr4MYbrcvcuhBldyWVhXjySRGjk6FzJstOhTKjqqrEHB12jbNYJ6wiMWbMGAwYMADbtm3Deeedh5dffhn3338/li9fjk6dOuGjjz7C/fffH84ihgz5xqZWmL+YBGAKgXwT292s6sOkczepZdLFJKiFaRdPCIRAUmBld1MgMQkZVTjtxCaQHtf+WsdusTuev8rHjbvpgw/cl4MmiEpONmexC8Td5ITTtSovB4YNM7///vdiQEMZXec32l+vXmKcJbt0XfWeku//yZNFAogMpQHbdbiTG3FPPGFOZhVvhNXd9Prrr2uXr7AbDyCKkSuvhg2FL9hJJO6/X/TC1uWHuwlcA9aHZMwY39amKhKqJUEV0kUXuTtHHYFkN/mLSdRGJHQVnl1Wk5vOdMHQurUYals9nj+Lol074PPPxWf1/Oh7oP021I5m9WFJqAkQ997ru41uGA26Ph07Arm5wH/+47tNx45imH353lDLod435IJKTRWCMWmS7yjJgDkbZLx2uHO89RcuXAiPxwPDYUzdtLQ0XHHFFSEvWKwhP4RUWcsxCRIJijtkZYkH/9gx3335czeR+Sw/lDRxjIwauNbFJJYsAQoK3J2jjkAsCX8xCZozIRSWhNO7SigsiUOHgM2bgV/9yvd4/irlv/1N9DRv0sReJAJFHdeoPiwJNzhZEvSuO+dzzxV9kexiEoDvvUANL2oUPfusGB2WkJ/Pigpg0yYRbwvm+kQjjiJxxx13YIRuQPafMQwDn332GYuEC+QMJXI3yeas12t1ByUni++HDoksFJotDfAfuKYB4PylrtqlwMrupssv939u/o4hv8tlpfLec4+I08g+ZadWvVPL240lYdeSr0t3U7NmZsqpfLyEBP+VTkqK+Z+o5/eXv/i6UdzgJBLB4C8m4Q+dSKj/j04k6L6Qp/N1siRSU81OrPJ16NQJWLZMfJaHsDEMMQXsBx9YhSQecPw7L7vsMrz88suOOxg3blxICxSryCJBAeRJk4RPmNxNcmBZFom0NKtIEHaWRGqqEAmnETAB8UDq3E1qTCIY3FgSkyebIuHkbgJ8h7VWMQzgiy9EZsrixe4yoerL3aSmFJNIBENWlvN4X05lqSt3U20F9bXXrBMVEaoloRM0upaySKiNJLlczZqZloS8nfyckXtJ7hvjZo6NWMPxFn1NnnkliG0Yq0jIFc7nn5siIVfqNGjfoUP+K3tCtiQA/5YEjcnkL3AdDG5iEvLxndxNALBypbUi0DFggDmukBtLoj7cTbQf9fiBXuNQzaamWhKBzLTnb7+1FdSxY8XwLyrq/6S7BnQt5fNQ/zNVJKiTqvycyKmyuuHDQ9FYiDb83qK7d+/G4Z+nbfrqq6/w6KOP4u23367zgsUacktNbo0kJ5sxCVUkyJKwq+ztLAm3ImFnSVA5QuF7dbIk1MrZXwos4DzAIGBeEzcB8PpOgdXFQALdb6hEIjHRen888ojIHAKC+99Dda3UfcrvOpciuZtkkXDqM9OsmfkcyuWVRaKoyLcs8SgSjqf8l7/8pabPwpgxY7BixQrk5+djyZIlKCwsxFNPPVUvhYwFZEtCFQmKSegsiYMHxWd1mGkdtbEkGjTw9ZUHOgyHE04xCWr9yQFctSNfoNA1ov0EYkk4xSRC6W6KREuicWOROQTUbT+J2lBbS0JFtSQI2X1FIpGWZh2UU97H9u0idhEvc2Q7/p2vv/46ioqKcPr0abRu3Rr79+9HRkYGKisrkUt3FKPl0CHg9ddF3MHjsYqE2uuU3E1y5ZySIirv48fFw9y8uRAMJ9QUVjci8d//7ZtlE8qWoFOPa8DqyvF4zMyv2rq66MGlCjBUMYlQupvUwHUghFIkQtkYkPcbaksiEJFwcs3S/1hVZRUJ2conkWjXTkxr+uqrvvvRiUcs43iLpqamIiUlBU2aNEGHDh2Q8XOT0+v1ItmtozxOee450VmIgl527qYbbxS9XnXupuRkEUhLShIioRKsuyk5WQzJQK0valGFMsVP526aOlWIJ+CbXUXB+9pWNG7iKnaxiLqOSYTC3RQqVEtCJpj//9JLgbvvrv3vdagi7tbdpNsPrT/nHHO5zpJIT9dbRPIQ//GCoyVRUlKCf//73zAMo+YzgJrvjD0kBOqw4PI6QNxwBw/ai8SxY+K9SRNg61bnY9ZGJGTqYowanbtJniLk978XLq8LLhBzYqxe7bu9W774Arj4YvHZKa5SmxTYSMluqit3U6ho2VK8Qokq6k7ZTU7nlJAgrOyWLcXzRMj3PWU3ZWTon5+KCt8BN2Mdx1v/0ksvxbvvvuvzGQB++ctf1m3Johy6kXUiofbcLC31dTeRSFB/h1deAa691qxEdcjuJjcVkJ1IhNLXauduImgAw40bxfu2beK9NiIxYID52cmSCJe7SSdOkehuirTOYm4C125jEl6viCm89JK5vLpaXIuKCjN1vEMHfcNAnZs+HnAUiXnz5tVTMWIPEgl6t7MkADEkgNyyO/dc0bOT3E3JyUBOjhjBVBYJu97FaWnu/M31IRJuh54gKCYRbKUciCXhz900cmTtpitV0QWuw5XdVFfuprqgNoFr3TbNmgETJ4rPsphUVQn3UmmpEIkhQ8QcEr//ve8+KiutLuRwuQvrE0eRePzxx23ncwCAe3WDrzAAfC0JpwH6SkuFeUsP7dy5QhTIkiDT2F9FK1sSbkRCbXXVpbvJbcVDvdGDrahqY0nYHVOe2yAYVJHw13tcR6QHruuCQFJgndxNqamihzrg29s8PV18LikxEz90lkRFhemSKiszfxfLOIrEyZMn4fF4sG3bNqxevRojRoyAYRh477330K9fv/oqY1QSqLupUSPfjmRJSWIdzS6njuGvPix0zKSk2lkSVMa6cDe5rQz9dZRzS6AxCY+n7lvQun4h4RIJNQYmE6mWBL3r7hG6lp06udunTiQSEoRIOM3IWFlpBrfLy1kkaib2ueSSS7Bu3To0/NkXMHPmTB6vyQ+BuJsoJkE3rtypjbKbADEefuvWwM9TgPtAgiT3XHaiPgPXbiueUIlEoJZEKHqX+yMU7qZQliVaLAnV4psyBRg1ylwvxyQGDXIXM5DvfXI30fVwsiQqK62WRDzg6tE4ePAgkqQ7KikpCQf9Je3HOeo8EP5EQvYRqyJBy1NTxZj6dtAx5DGQnIjEmER9WBJqmeorFZWO0bWrefxIzG6KNEtCtcBSU4Fu3Xy3kfve+ENnSdC959TPiILbgFUkNm60DjMeS7i6RW+66Sb069cPM2bMwIMPPoj+/ftjgl1zNo754gvzxlEtCbmVXlUlxqmZPl18V7Ob6CanmITdvL1qhdGunZjExa0lYReTyM0V/TxCgb/sJpWmTUNzXDfupmDcPrWB/ru5c8XgcuGOSURLVyfVkrDbJpBrqQ5umJ5uunPVOd5lHnnEHK5DbuwdOSJcVfJggLGCq8s6ffp0vPzyy8jMzETTpk0xb948TJs2ra7LFnUMHCiyIgDnmARgfUid3E1nz/r2xLYjLQ14/nn/ATwAGD9eCJUMlTU5GbjrLuffuyVQS6JpU3PGsGAINLupPkQiMRFYtUr8nw0bht+SiJYUWOodvWePfr3sbnKLaklkZJgi4eRu+u474L33xGfZkqDfrFnjvgzRgusuQhdeeCEu1A3RyFhYtw6YNs0UiWefFX0BdCJBD6mTSMjvgHWua7sKw40l8corvssiISYBhCbd1On8dSJRX7GBvn3Nz4HGJIYPF69QEG3upuRkc/4Hu21qKxJkSdB5ux3WRrYk6HkvLnZfhmghDsc0rFto+o077hDvTz4JfPSRmAB+zBjR4WvSJGsv3upqe3cT4N6SINy6m1TqQiQCzW4KFU5+6XBZEiqBtn4XLw7dsRs08J3ONpI5dMj5Pw3WkqDsJkDvbsrJEe7cTz81l8mWBIlELHa0C8OjER/IQwds3CimAe3XT3SSA3xveLvAtfwO1K1ItGsX+G/clAWo/9ZpIJZEfcUkdOUIx3EB4Ra1m3Uw0iwJQKSIO6WbBiMSVVWi4TZ+vPiusyRSUswhXwDRd0n+TlZ9LE5KxCJRR+jGl5EHivN6re4inbuJblI7SyIYd5OOu+8OTTxALQsQ2SIRTksiXCmwGRmx01tY7kznFtWS6N4daNtWfNfFJNR5XWgbshzYktBw++23h7IcMYfOdaOm6clCIrubnCwJN5V/bUXC4wlNPEAti/xeXzidf7hSYFUaNbJOWRspRKIl4Y9gLYmEBFMU1M5055wDtGhhFQ26X2icU3qW2ZKQmEiDoMQx8ixWKrpUOK/XKhKqJaEOD6ETCTcPcG1Foi4INAU2VLixJIIZHiMUdOoEFBbW/3GdGDoUuO66cJcicIKNSdDgf4CvJbFmjRgzTRYJ6lBHwfRYtiQCClyfOHECHo8HDRs2RJ8+feqqTFHBli1Az572Lh+dgKjuJtmS8HrNipSW041s54sNtbupLohESyJSYhJyWSKFDz4IdwkCJ1h3k2GIfagiQfdQkyZi5j5ZJE6cEO9kScR9TGL16tXo2bMnevbsiR49eiA3NxdrYjEhuBYcO6ZfTi0NGdWSkJErC1UkdBMOORGJIhGJlkS43U1M6AjGkqCxu+ieUbOb1Hks1qwx3cmqJaGKxPPPAwcOuC9XJOLqst56662YO3cudu/ejd27d+O5557DrbfeWtdliwrsegjrREKNSajuJsKtSESDJRGuFNgWLezXRUrgmgkNtelMJz8fasNNtSRoOfV1kruLUe9rO3fTXXcB8+e7L1ck4uqyer1eXHLJJTXfBw4cCG+oZzuPMpxmp0pNtRcJO3eTLBJyz2fAucLTEUkiES5LolEjexGNJHcTExoC/Q/l581OJFRLQk5GadNGvN9zD3D4sHPgOlqGP7HDsaZfu3YtADEL3cSJEzFmzBgAwIIFC+J+ZjqnjmeNG+tjEgcPik51gLgBc3LMdXQjTpgg0vEAs6KPBXdTJFXCbEnEFrWJScjIDTfA15Kg9bKV8O23Yp6RtWuFy5kaJLrAtSwSX31l7V8RDTiKxH333WeZdGjmzJkAxBzXTpMRBcvSpUsxefJkVFVV4Te/+Q2mTp1aZ8eqLVVVwuxcu9YMfBGNG4u8apXycmur5cYbxZSkDRqYv5cnA/QnEtHkboqktMpISYFlQkcw1qBqUarDzNO77D1ISTETSo4dc29JDBggGouBNvzCiaNIFBYWoqqqCm+99RZ+/etf10uBqqqq8Lvf/Q4rVqxAq1at0LdvX4wYMQJdaXzlCKGyUlTECQlmCh3RuDGwa5d1++++EybqunXiO2Uz0fDEdimzgLte1jJuhwqvD6LBkkhPD33/EKZ+CYVIAOKZI0tCbdio3oN33xVTnT7zDPDqq2KZbEksXCjeb79d7IvCuNHmqfd7WRMTE/HII4/UR1kAAKtWrUKHDh3Qpk0bJCUl4YYbbsCiRYvq7fhuofltvV7f+ETjxr4tinPPNUUF8L1RdDGM7Gy9ReIPN6PA1hfhikk4ofaTaN8+8vorMO4J1t1kJxLq/nTPedu2wBtvmMs++gjYtk18lmd3/uADMaJzNOLqshYUFOCxxx7Dnj17cPTo0ZpXXVBcXIzzaIAjADk5OSiOwKEVq6rMQfrUm0fXi9ZuXmVCJxKAOUyyDrspPRo3Dt28DMESruwmJ9T/AojdCWPihWAsCbkBk5RkWu7+LAlA9KGQn/8ffjA7I8pD3BiG6S3QDdkTybgyfN544w14PB48p8xEs3PnzpAXyG2sg6ZWBYD8/Hzk5+eHvCxOVFbaWxK60TXVTAk3loQTl18OPP64fh2NQBsJRLIlEUnCxdSe2qTAqr8nHnhADMOhLgf0GY1NmvguI0tETl6pqjItifoUicLCQhQGaSa7EoldqoO9DmnVqhX2SLOL7NmzBzlyGtDPyCIRDpwsCXnOB0KtmIIVCafu/5FUIUdDTIKJfmrjbho8GPjwQ+vzIufIqM/Rgw/6zuchi0RKihg+PDVVCIHsXnr7bTPuVZ8ioTagKfkoEFyHULZs2YKioiKclc78pptuCviA/ujTpw+2b9+OXbt2oWXLlliwYAFef/31kB8nGP70J5G+qrMk5IHCCGrpAKGxJAYODN3kM3VNJGY3sUjEFrW1JFasEL+1uzfV5dnZwBVXWJfJbkpZJHSJKP/8p3iPSXfTjBkz8Mknn2Dr1q248sor8f7772PgwIF1IhJerxfPPvsshg0bhqqqKtx2220Rldn0+efArFnAoEHmcMuVlWY6akqK780qC4LOkhg40H5sfx2ffVa7socDdjcx9UGo3E1ulsuolgRg35mWiEmReOutt7Bx40b07t0bL7/8Mg4cOIBx48bVWaEuv/xyXB5IrVnHnDolbsD0dEDqeG5xN1FQKynJNyjt1LsTiK5KP1Ai0d0UiWVigiOY7CY7MdB4uX2QRYLczGlpsSUSri5rWloaEhMT4fV6UVJSghYtWljiBrHOlCnWTm6E7G5yEgndOPTRlitdWyLZ3cQd6GKDUAauZc4/377DKmFnSThN3hVtIuGqqurbty+OHTuG22+/HX369EFGRgYGDBhQ12WLGNLT9X86WRJVVWamUXJy4JZELBOJrXZ2N8Ue4Rp/S2dJxKW7ae7cuQCA3/72txg2bBhOnDiB3NzcOi1YJKETCerAQ5bE9OliuT9LIl5FIhItCRaJ2CDYznTB3JuySFAH1pSU2BIJx8u6b98+n2Vt27a1CIRum1jDn0hs3GguZ3eTFbYkmPqgLtxNbmjc2PxMLufq6thyNzle1iuvvNLvDtxsE+1kZDi7m8aONZfVJnAdy0RyTIJFIjaoq5iEG+Rnm/ouVVbGliXhWFVt3LgRDXXdhyUaReJM7iHGLiahGz1UF5Ngd1NkVcgsErFHuERCpmdPMQlRXIlEldOkCXEEiYR6OShoLePPkohXd1MkWRKRKFxM7fF4gK5dnV08dckbbwA33CDiE/PmiV7chw4BkyeLqUvVvsAxJRKMgERCHgqDRoF1IxLxbElEywB/THTTv7941YZgU6G7dRPvNI9LZSWwc6eYPEw3+Vi0iQQ/Ji7QiURFhV4kfvc7FgmZSLQkuJ8EIxPssyhbppTtuGMH0K6dviHCIhGDkEjIYzSVlendTePGsbtJJhJdO2xJxBbBNkCCbSzoRGLnTr1ItG0bfSLhqqratGkTvv32W3g8HnTt2hU9evSo63JFFDpLorxcb0kA7iyJeGnFcnYTE+mESiRoRsiKCmDfPqBlS+s99txzwIsvxphIlJSU4Oqrr8aPP/6I3NxcGIaBzZs3o3Xr1li0aFFcZDYBZgpsbUWCLYnIqpBp5M9IKhNTeyLRkjh7VozhJN9jV1wBvPxy9ImE42Pypz/9CX369MH333+Pt99+G++88w6+++479O3bF9Opi3EckJ4uUtpUkdC5mwDzpqObl2MSkWVJAGJimXix5hhngn0WZcvU6xWuaJpvRhYJmr6YRGLLltpNT1zfOF6eFStWYNOmTUiQzjQxMRGzZs1Cz54967xwkUIggWvAKhKGwZ3p5PdIYedO8dAy0U8kWRJJSaJBmZJijspAJCdbRWLaNOC224Crrw7u+HWNoyWRnJyMJM2TlJSUhBQa8jAOaNhQpLLJedhuLAmdIMSjuynSBAIAGjQIdwmYUOFvpFZ/hNrddOqUOdifzpIoKRFWxOnTogEa6ThWVWVlZVi3bh0Mw7DMPW0YBsrKyuq8cJFCcjLQujXwzTfmMjcxCXY3Ra5IMAwRysC112taEvI6wBSJyZOBbduAiy4S8c5Ix7Gqys7Oxn333addd+6559ZJgSKVrl2BTZvM74GIhHwT6oQjluEAMVPXBNsICWU/iaQkYUlQ5a8TCRoTNSYsicLCwnoqRuTTuTOwdav5ndxNct8JQg3WqjchjR4bD7AlwUQ6oXY3lZaKxAh5HWCKxIkT4ntpaXRYEo5tvNWrV1uGAp8/fz5GjBiBe+65B0ePHq3zwkUSDRqYfy7gPnAtfyeaN4+foGliYvycKxMeIilw7fWK9Fedu0m1qk+fjgGRuOOOO2oC1J9++inuv/9+TJgwAY0aNcIdd9xRLwWMFBITxZ8vU5vANQDs2mVOUBLrpKYCa9aEuxQMY0+wIiE3CqlBpBMJ9XtpaQy4m6qrq9G0aVMAwIIFCzBx4kSMGjUKo0aNiquZ6QDx55aVAb/4hTAlFy92b0morYU4SgwDIOI5DFMXfPhh8JlqoR67CXAnEjFhSVRVVaHi584BK1aswKBBg2rWVeqc8TFMQoKwJDIzRXwCcC8Szz9fP2VkmHjjVycR5HUAABVaSURBVL8C+vULbh+hdjcB+hRY3fdocMU6auiYMWPwy1/+Es2aNUN6ejouueQSAMD27duRmZlZLwWMFBIThSVBwSfAv7upoAD44QfRz4JhmMgk1J3pAF9LQmdZRIMVAfgRienTp+NXv/oV9u/fj6FDh9b0vDYMA88880y9FDBSIHeTPF9EYqIY1XHnTuu2tP7554GsrPotJ8MwgRHqfhKAO5GIhngE4MfddObMGaxcuRIrVqzAq6++WuNi6tSpE3r37l0vBYwUKHCtWhLr1onhwdVtAe4fwDDRQF24m1RRGDXK+h0QAwBGA47V2IQJE7B27VpccMEFWLJkiW3HuniALAmv12pJZGYCzZpZt2WRYJjoIVTZTbK7SY1JPPus9XsojltfOLqbvvnmG2zevBkAcNttt6Fv3771UqhIRHY3qXNCqMEnFgmGiR5CaUmoFoRaV0SjSDhWY14pN8wbL12EbUhMFB3o5JgEXRKvV6TGfvmluS3AIsEw0UCoUmDlfhIUp1TrArlOoF7ZkY5jNbZp0yY0bNiw5rV58+aaz8FMOPSvf/0L3bt3R2JiItatW2dZN3v2bHTs2BFdunTBsmXLan2MUCO3COgzdYhLSgLOPRe4+GJzG/k3DMNELsG26HXD4W/bJt5VcZDrhObNgztufeGooVW6/M4Q0LNnT7z99tuYOHGiZXlRUREWLFiAoqIiFBcXY8iQIfjuu+8s81mECzkOobqZkpKsLicWCYaJHkItEpMnm/WBndsJiBGRqCu6dOmiXb5o0SKMGTMGSUlJaNOmDTp06IBVq1bhoosuqucS+uLPktDNGcEiwTCRT6hiAzSvxZNPmstogEsSELlOaNEiNMetayKqGtu7dy9ycnJqvufk5KC4uDiMJTLRWRKySLAlwTDRSV2GWxMSrPUAfZ42DfjNb+ruuKGkzi5PQUEB9u/f77P8oYcewvDhw13vx2MzxOOMGTNqPufn5yM/Pz/QIgaEbEno3E1sSTBMdBJqS0JGFQli1qzQHNMfhYWFQU/5UGcisXz58oB/06pVK+zZs6fm+08//YRWrVppt5VFoj6QfYuqu6lfPxG4JlgkGCZ6qMtUVFUkaH7r+kJtQM+cOTPgfYS9GjMk+R0xYgTeeOMNlJeXY+fOndi+fTv6BTt6V4hwClxfdBFw/fW+27JIMEzkE8siEQrCUo29/fbbOO+887By5UpceeWVuPzyywEA3bp1w+jRo9GtWzdcfvnlmDt3rq27qb5xClzbbRshRWcYxoG6jknIIqRzSUU6Ycluuuaaa3DNNddo102bNg3Tpk2r5xL5xylwrcLiwDDRwf/8D3DDDaHZl5uYRDRaEvHdjToAdJaE3VjwLBIMEx1Mnx66fbkRiWi0JNhr7hJddhNbEgzDOBELlgSLhEtkdxNbEgzDuCEWYhIsEi5hS4JhGCd0AiA3KgG2JGIaOa2VLQmGYdzA7qY4QrYkqMVg1w+CxpRnGCa+4cB1HCGLhL/BcbOzge+/r/syMQwTOdhlN8kxCbYkYhg5cO1mBPX27eu2PAzDRBacAhvnBGJJMAzDAByTiCtkS6KyMrxlYRgmOmBLIo5gS4JhmEDhmEQcwSLBMEygsLspjpD7SUTjH80wTN3Cges4hy0JhmGciNVRYFkkXMKBa4ZhAoXHboojZEti0CAgLy+85WEYJvJhSyKOkEWiTx9g3brwlodhmMiCYxJxjuxuYhiGcQOPAhtHyJYEwzCMis5KyMsDHn3UeZtIh0XCJWxJMAzjhG5+mdRU4NJLze/RaEnwHNcuIUvCbnhwhmHil1WrgJ49/W/HIhHDsCXBMIwdffu6247dTTEMxyQYhgmWaLQkWCRcwiLBMEywZGSEuwSBw+4ml7C7iWGYYFm4EDh2LNylCAwWCZewJcEwTLA0by5e0QS7m1zClgTDMPEIi4RLOAWWYZh4JCxV3pQpU9C1a1fk5ubi2muvRUlJSc262bNno2PHjujSpQuWLVsWjuJpYXcTwzDxSFhEYujQodi6dSs2btyITp06Yfbs2QCAoqIiLFiwAEVFRVi6dCnuuusuVEdIzhi7mxiGiUfCIhIFBQVI+Llp3r9/f/z0008AgEWLFmHMmDFISkpCmzZt0KFDB6xatSocRfSBLQmGYeKRsHvYX3rpJVxxxRUAgL179yInJ6dmXU5ODoqLi8NVNAtsSTAME4/UWQpsQUEB9u/f77P8oYcewvDhwwEAs2bNQnJyMsaOHWu7H4/HU1dFDAi2JBiGiUfqTCSWL1/uuH7evHlYsmQJPvzww5plrVq1wp49e2q+//TTT2jVqpX29zNmzKj5nJ+fj/z8/KDK6w8WCYZhoo3CwkIUFhYGtQ+PYdT/kFNLly7Ffffdh08++QTNmjWrWV5UVISxY8di1apVKC4uxpAhQ/D999/7WBMejwdhKDY8HmDvXuDcc+v90AzDMEFTm7ozLD2uJ02ahPLychQUFAAALr74YsydOxfdunXD6NGj0a1bN3i9XsydOzdi3E0MwzDxSFgsiWBhS4JhGCZwalN3hj27Kdpgw4ZhmHiCRSJAOHDNMEw8wSIRABs2RN8IjgzDMMHAMQmGYZg4gWMSDMMwTEhhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbGGRYBiGYWxhkWAYhmFsYZFgGIZhbAmLSPz3f/83cnNz0atXLwwePBh79uypWTd79mx07NgRXbp0wbJly8JRPIZhGOZnwiISf/zjH7Fx40Zs2LABI0eOxMyZMwEARUVFWLBgAYqKirB06VLcddddqK6uDkcRw0phYWG4i1Cn8PlFN7F8frF8brUlLCLRsGHDms+nTp1Cs2bNAACLFi3CmDFjkJSUhDZt2qBDhw5YtWpVOIoYVmL9RuXzi25i+fxi+dxqizdcB54+fTr+8Y9/IC0trUYI9u7di4suuqhmm5ycHBQXF4eriAzDMHFPnVkSBQUF6Nmzp8/r3XffBQDMmjULP/74I2655RZMnjzZdj8ej6euisgwDMP4wwgzu3fvNrp3724YhmHMnj3bmD17ds26YcOGGStXrvT5Tfv27Q0A/OIXv/jFrwBe7du3D7iODou7afv27ejYsSMAEYfIy8sDAIwYMQJjx47Fvffei+LiYmzfvh39+vXz+f33339fr+VlGIaJV8IiEg888AC2bduGxMREtG/fHs8//zwAoFu3bhg9ejS6desGr9eLuXPnsruJYRgmjHgMwzDCXQiGYRgmMom6HtfPPPMMunbtih49emDq1Kk1y2OpE97jjz+OhIQEHD16tGZZtJ/flClT0LVrV+Tm5uLaa69FSUlJzbpoPzdi6dKl6NKlCzp27IiHH3443MUJmj179mDQoEHo3r07evTogaeffhoAcPToURQUFKBTp04YOnQojh8/HuaSBkdVVRXy8vIwfPhwALF1fsePH8d1112Hrl27olu3bvj6668DP7/aBpzDwUcffWQMGTLEKC8vNwzDMA4ePGgYhmFs3brVyM3NNcrLy42dO3ca7du3N6qqqsJZ1Frz448/GsOGDTPatGljHDlyxDCM2Di/ZcuW1ZR56tSpxtSpUw3DiI1zMwzDqKysNNq3b2/s3LnTKC8vN3Jzc42ioqJwFyso9u3bZ6xfv94wDMM4efKk0alTJ6OoqMiYMmWK8fDDDxuGYRhz5syp+S+jlccff9wYO3asMXz4cMMwjJg6v5tuusl48cUXDcMwjIqKCuP48eMBn19UicT1119vfPjhhz7LH3roIWPOnDk134cNG2Z89dVX9Vm0kHHdddcZGzdutIhELJ2fYRjGv//9b2PcuHGGYcTOuX355ZfGsGHDar6rmXqxwNVXX20sX77c6Ny5s7F//37DMISQdO7cOcwlqz179uwxBg8ebHz00UfGVVddZRiGETPnd/z4caNt27Y+ywM9v6hyN23fvh2ffvopLrroIuTn52PNmjUARCe8nJycmu2itRPeokWLkJOTgwsuuMCyPFbOj3jppZdwxRVXAIidcysuLsZ5551X8z1az8OOXbt2Yf369ejfvz8OHDiArKwsAEBWVhYOHDgQ5tLVnj/84Q949NFHkZBgVoWxcn47d+5E8+bNccstt6B37964/fbbUVpaGvD5ha3HtR0FBQXYv3+/z/JZs2ahsrISx44dw8qVK7F69WqMHj0aO3bs0O4nUrOinM5v9uzZFp+84ZBTEInnZ3duDz30UI2/d9asWUhOTsbYsWNt9xOJ5+aPaCyzW06dOoVRo0bhqaeesgypA4jzjtZzf++999CiRQvk5eXZDscRzedXWVmJdevW4dlnn0Xfvn0xefJkzJkzx7KNm/OLOJFYvny57brnn38e1157LQCgb9++SEhIwOHDh9GqVSvLSLI//fQTWrVqVedlrQ1257dlyxbs3LkTubm5AMQ5XHjhhfj666+j5vyc/jsAmDdvHpYsWYIPP/ywZlm0nJs/1PPYs2ePxUKKVioqKjBq1CiMHz8eI0eOBCBan/v370d2djb27duHFi1ahLmUtePLL7/E4sWLsWTJEpw9exYnTpzA+PHjY+b8cnJykJOTg759+wIArrvuOsyePRvZ2dmBnV9d+MLqihdeeMH485//bBiGYWzbts0477zzDMMwg59lZWXGjh07jHbt2hnV1dXhLGrQ6ALX0Xx+77//vtGtWzfj0KFDluWxcG6GIYKC7dq1M3bu3GmUlZXFROC6urraGD9+vDF58mTL8ilTptTEkWbPnh3VgV2isLCwJiYRS+d3ySWXGNu2bTMMwzAefPBBY8qUKQGfX1SJRHl5uXHjjTcaPXr0MHr37m18/PHHNetmzZpltG/f3ujcubOxdOnS8BUyRLRt27ZGJAwj+s+vQ4cORuvWrY1evXoZvXr1Mu68886addF+bsSSJUuMTp06Ge3btzceeuihcBcnaD777DPD4/EYubm5Nf/b+++/bxw5csQYPHiw0bFjR6OgoMA4duxYuIsaNIWFhTXZTbF0fhs2bDD69OljXHDBBcY111xjHD9+PODz4850DMMwjC1Rld3EMAzD1C8sEgzDMIwtLBIMwzCMLSwSDMMwjC0sEgzDMIwtLBIMwzCMLSwSTMSSmJiIvLw8XHDBBbj22mtx6tQpx+1nzJiBxx9/3HGbRYsW4Ztvvqn5/uCDD1p6gNeWm2++Ge3atUNeXh7y8vLw7LPPBr1PmcLCQjRu3BhXXXVVzXca6kRm3rx5mDRpkmVZfn4+1q5da7vvcePG4ZxzzsHChQtDWmYmNmCRYCKW9PR0rF+/Hps2bUKjRo3wv//7v47buxlj5+2330ZRUVHN95kzZ2Lw4MFBl9Xj8eCxxx7D+vXrsX79evzud7+zrK+qqgr6GJdeeinee+89v+XQLXO6Nq+99hpGjBgRtWMUMXULiwQTFVx88cX44YcfAAA//PADLr/8cvTp0weXXnoptm3b5rP93//+d/Tr1w+9evXCddddhzNnzuDLL7/Eu+++iylTpqB3797YsWMHbr75ZixcuBAffPABRo8eXfN7uaW+bNkyDBgwABdeeCFGjx6N0tJSbRnVfqn5+fn4wx/+gL59++Lpp5/G2rVrkZ+fjz59+uCyyy6rGQxx7dq1yM3NRa9evTBlyhT07NkzoGuzevXqmvOxwzAMvPvuuzWWTufOndGuXTvH8jMMwCLBRAFVVVVYtmwZevToAQC444478Mwzz2DNmjV49NFHcdddd/n8ZtSoUVi1ahU2bNiArl274sUXX8SAAQMwYsQIPPbYY1i3bh3atWtX08oeMmQIvv76a5w5cwYAsGDBAowZMwaHDx/GrFmz8OGHH2Lt2rW48MIL8cQTT/gczzAMTJkyBXl5eejduze2bNkCj8eDiooKrF69GpMmTcKkSZOwcOFCrFmzBrfccgumT58OALjlllvw3HPPYcOGDQGPOvrll1/izjvvxOLFi9GuXTsYhoEFCxbUiEFeXh7WrFkDj8eD4cOH11g6JEgM44+IGwWWYYgzZ84gLy8PxcXFaNOmDX7729/i1KlT+Oqrr3D99dfXbFdeXu7z282bN+NPf/oTSkpKcOrUKVx22WU163Qt5sTERFx22WVYvHgxRo0ahSVLluCxxx7Dxx9/jKKiIgwYMKDmWPRZhtxNNEox8etf/xoA8O2332Lr1q0YMmQIACF8LVu2RElJCUpKSjBw4EAAwPjx4/H++++7uj7ffPMNJk6ciOXLlyM7O7umHDfccEPNVKMAMGjQIMvvHnnkEaSnp+POO+90dRwmvmGRYCKWtLQ0rF+/HmfOnMGwYcOwaNEiDBkyBJmZmVi/fr32N9QKv/nmm7F48WL07NkT8+fPt8wXYNdSv+GGG/Dss8+iadOm6Nu3LzIyMgCIeTL++c9/+i2vTnxoH4ZhoHv37vjyyy8t69X5hd26fDweD84991yUlZVh3bp1NZM4+dvHihUrsHDhQnz66aeujsMw7G5iIp60tDQ8/fTTmD59Oho0aIC2bdvirbfeAiAqxE2bNtVsSxXkqVOnkJ2djYqKCrz66qs1wtCwYUOcOHHCsn/6zaWXXop169bh73//O2644QYAQP/+/fHFF1/UxENKS0uxfft212WnfXfu3BmHDh3CypUrAYh5GoqKipCZmYnMzEx88cUXAEQQ2e1+MzMz8d577+GBBx7AJ598Yjmejt27d+Puu+/Gm2++iZSUFNfnwMQ3LBJMxCK3+Hv16oUOHTrgzTffxGuvvYYXX3wRvXr1Qo8ePbB48WKf3/z1r39F//79MXDgQHTt2rVm/Q033IBHH30UF154YU2gl36TmJiIq666CkuXLq1JNW3evDnmzZuHMWPGIDc3FwMGDNAGytXyqsuSk5Px1ltvYerUqejVqxfy8vLw1VdfAQBefvll3H333cjLywvo2ng8HrRo0QLvvfce7r77bqxatco2pmEYBubPn4+jR49i5MiRyMvLqzlHhnGChwpnmAhi9+7duOqqq7B582bL8sLCQjz++ON499136+S4N998M4YPH45Ro0bVyf6Z6IUtCYaJIAzD0FoCKSkp2LJlS520/seNG4fPPvsMaWlpId83E/2wJcEwDMPYwpYEwzAMYwuLBMMwDGMLiwTDMAxjC4sEwzAMYwuLBMMwDGMLiwTDMAxjy/8HJpGbYLWw4qcAAAAASUVORK5CYII=", "text": [ "" ] } ], "prompt_number": 37 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The \"SpecMaster\" function used above also works for TBW data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Post-Acquisition Beam Formings with TBN" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For post-acquisition beam forming, you need need an azimuth (in degrees) and elevation (in degrees) to point the beam towards. For planets, this can be accomplished using the pyephem package that is required by lsl. For example, compute the location of Jupiter at LWA-1 on 12/17/2010 at 21:18 UTC (JD 2,455,548.38787):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import math", "import ephem", "from lsl.common import stations", "", "lwa1 = stations.lwa1", "lwaObserver = lwa1.getObserver(2455548.38787, JD=True)", "jove = ephem.Jupiter()", "jove.compute(lwaObserver)", "print \"Jupiter: az -> %.1f, el -> %.1f\" % (jove.az*180/math.pi, jove.alt*180/math.pi)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Jupiter: az -> 112.4, el -> 24.4" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "For fixed positions, use:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "lwaObserver.date = '2010/12/17 21:18:34'", "", "cyga = ephem.FixedBody()", "cyga._ra = '19:59:28.30'", "cyga._dec = '+40:44:02'", "cyga.compute(lwaObserver)", "print \"Cygnus A: az -> %.1f, el -> %.1f\" % (cyga.az*180/math.pi, cyga.alt*180/math.pi)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Cygnus A: az -> 10.0, el -> 83.2" ] } ], "prompt_number": 35 }, { "cell_type": "markdown", "metadata": {}, "source": [ "After TBN data have been read in and a pointing position has been found, a beam can be formed through phase-and-sum beamforming:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from lsl.misc import beamformer", "", "antennas = [lwa1.getAntennas()[0],]", "", "beamdata = beamformer.phaseAndSum(antennas, frame.data.iq, sampleRate=1e5, azimuth=10.9, elevation=83.2)", "print beamdata.shape" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(1, 512)" ] } ], "prompt_number": 39 } ], "metadata": {} } ] }