{
 "metadata": {
  "name": "",
  "signature": "sha256:32ca49506552149a56b0f2cf4ecd6412feda80beb48569a5cfa3a979ed419e88"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import netCDF4\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "import pandas\n",
      "import datetime\n",
      "import scipy.interpolate \n",
      "import pytz\n",
      "import dateutil.parser\n",
      "import bisect\n",
      "from numpy import deg2rad\n",
      "import matplotlib\n",
      "import openearthtools.physics.tide\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "\n",
      "\n",
      "def interpolate_window(lon_i, lat_i, lon, lat, arr, circular=False):\n",
      "    \"\"\"\n",
      "    interpolate tide with a window (of 4 by default)\n",
      "    lon_i, lat_i are the point to interpolate the tide on\n",
      "    lon, lat are the spatial dimensions of arr\n",
      "    if arr is circular it is split up into 2 components and put back together\n",
      "    \"\"\"\n",
      "    \n",
      "    # lookup position of the point\n",
      "    lon_index = bisect.bisect(lon, lon_i)\n",
      "    lat_index = bisect.bisect(lat, lat_i)\n",
      "    # define a window size (for spline interpolation)\n",
      "    window = 4\n",
      "    \n",
      "    # determine the shape\n",
      "    lon_min = max(lon_index-window, 0)\n",
      "    lon_max = min(lon_index+window, lon.shape[0]-1)\n",
      "    lat_min = max(lat_index-window, 0)\n",
      "    lat_max = min(lat_index+window, lat.shape[0]-1)\n",
      "    \n",
      "    lonslice = slice(lon_min, lon_max)\n",
      "    latslice = slice(lat_min, lat_max)\n",
      "    Lat = lat[latslice]\n",
      "    Lon = lon[lonslice]\n",
      "    A = arr[:,lonslice, latslice]\n",
      "    results = []\n",
      "    for i in range(A.shape[0]):\n",
      "        if circular:\n",
      "            f_cos = scipy.interpolate.RectBivariateSpline(Lon, Lat, np.cos(A[i]), kx=2, ky=2, s=0)\n",
      "            f_sin = scipy.interpolate.RectBivariateSpline(Lon, Lat, np.sin(A[i]), kx=2, ky=2, s=0)\n",
      "            result = np.arctan2(f_sin(lon_i, lat_i), f_cos(lon_i, lat_i))\n",
      "        else:\n",
      "            f = scipy.interpolate.RectBivariateSpline(Lon, Lat, A[i], kx=2, ky=2, s=0)\n",
      "            result = f(lon_i, lat_i)\n",
      "        results.append(result.squeeze())\n",
      "    return np.asarray(results)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def predict(location, dates=None, nodal=True):\n",
      "    \"\"\"\n",
      "    Predict the tide for a location for 1 or more dates\n",
      "    \"\"\"\n",
      "    if dates is None:\n",
      "        dates = [datetime.datetime.now()]\n",
      "    T = openearthtools.physics.tide.nodalconstituents(dates)\n",
      "    T = {key.upper(): T[key] for key in T}\n",
      "    T = pandas.DataFrame.from_dict(T, orient='index')\n",
      "    # julian days\n",
      "    jd = openearthtools.physics.tide.dates2jd(dates)\n",
      "    jh = jd * 24.0\n",
      "\n",
      "    total = 0\n",
      "    for const, row in location.iterrows():\n",
      "        phase, ampl = row['phase'], row['amplitude']\n",
      "        if nodal:\n",
      "            component = ampl*T['FF'][const]*np.cos(T['ospeed'][const]*jh - (deg2rad(phase) - deg2rad(T['VAU'][const])))\n",
      "        else:\n",
      "            component = ampl*np.cos(T['ospeed'][const]*jh - deg2rad(phase))\n",
      "        total = total + component\n",
      "    return total\n",
      "\n",
      "def make_location(lon_i, lat_i, lon, lat, A, Phi, constituents):\n",
      "    \"\"\"create a dataframe for a location with the tidal information\"\"\"\n",
      "    amplitudes = interpolate_window(lon_i, lat_i, lon, lat, A, circular=False)\n",
      "    phases = interpolate_window(lon_i, lat_i, lon, lat, Phi, circular=True)\n",
      "    locationdf = pandas.DataFrame(data=dict(amplitude=amplitudes, phase=phases), index=constituents)\n",
      "    return locationdf\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# amplitude/phases\n",
      "ds = netCDF4.Dataset('http://opendap.deltares.nl/thredds/dodsC/opendap/osu/tpxo/h_tpxo7.2.nc')\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "const = netCDF4.chartostring(ds.variables['con'][:]) # the same as for u\n",
      "constituents = [x.strip().upper() for x in const]\n",
      "# we only need the first row and column\n",
      "lon = ds.variables['lon_z'][:,0]\n",
      "# loop longitude around\n",
      "lon = np.r_[lon, lon + 360]\n",
      "lat = ds.variables['lat_z'][0,:]\n",
      "# for indexing assume sorted\n",
      "\n",
      "assert((sorted(lon) == lon).all())\n",
      "assert((sorted(lat) == lat).all())\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lon_i = 4.766667\n",
      "lat_i = 52.958333\n",
      "A = ds.variables['ha']\n",
      "Phi = ds.variables['hp']\n",
      "location = make_location(lon_i, lat_i, lon, lat, A, Phi, constituents)\n",
      "dates = [datetime.datetime(2013,7,7,0,0, tzinfo=pytz.utc) + datetime.timedelta(minutes=i) for i in range(60*24*2)]\n",
      "fig, ax = plt.subplots(1,1)\n",
      "ax.plot(dates, predict(location, dates=dates, nodal=False))\n",
      "ax.plot(dates, predict(location, dates=dates, nodal=True))\n",
      "ams = pytz.timezone('Europe/Amsterdam')\n",
      "locator = matplotlib.dates.HourLocator(interval=12)\n",
      "ax.xaxis.set_major_locator(locator)\n",
      "ax.xaxis.set_minor_locator(matplotlib.dates.HourLocator())\n",
      "ax.xaxis.set_major_formatter(matplotlib.dates.AutoDateFormatter(locator=locator))\n",
      "ax.set_xlabel(\"hours in utc\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 59,
       "text": [
        "<matplotlib.text.Text at 0x1140ea0d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEPCAYAAAC5sYRSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4VMe5/7+jLgFCqFfUkARIAgkEElWiV2ODjXu5Ntdx\n4jj2da4T5+fkxjhO4sSxYyduwXGNOy4YbHqTAIGEBOod9Y4aKqhr5/fHaCWdPWfVtp1dzed5eKyd\nOWVYo++Z85133pdQSsHhcDgc08LM0APgcDgcjvbh4s7hcDgmCBd3DofDMUG4uHM4HI4JwsWdw+Fw\nTBAu7hwOh2OCaCzuhJDNhJB8QkgRIeRZiX5nQsgxQkg6ISSbEPJfmt6Tw+FwOKNDNIlzJ4SYAygA\nsB5ANYAUAPdQSvNGHLMXgDWl9P8RQpwHj3ejlPZrMnAOh8PhqEfTmftSANcopWWU0j4AXwK4VeWY\nWgD2gz/bA2jiws7hcDi6xULD870AVI74XAUgWuWYfwM4QwipATADwJ0a3pPD4XA4Y6DpzH08ns5z\nANIppZ4AIgC8RQiZoeF9ORwOhzMKms7cqwH4jPjsAzZ7H8lyAH8CAEppMSGkFEAIgNSRBxFCeJIb\nDofDmQSUUqLapunMPRVAECHEjxBiBeAuAIdUjskHW3AFIcQNTNhL1AxQ8s/zzz/P+ybYFxsbq7f7\nyenvbezf5VToM+Xvs7CQwtb2eVy7xj4rFBT33UfxzDO6G4s6NJq5U0r7CSFPADgOwBzA+5TSPELI\nY4P9+wD8GcCHhJAMsIfJrymlzRO5T1xcHO+bYJ+fn5/e7ienv7exf5dToc+Uv88//AG49944BAay\nz4QAf/87MH8+8MQTuhmLOjQKhdQmhBAql7GYAnv37sXevXsNPQyTgH+X2sVUv8/SUmDJEqC4GJg5\nU9j3zDMApcCrr2r/voQQUB3YMhyZMpknPUca/l1qF1P9Pj/6CLjvPrGwA8CTTwIffgh0d+tvPHzm\nzuFwOBpCKRAYCHzzDbBokfQxa9cCP/85cPvt2r03n7lzOByOjrh4EbC1BSIj1R9z//3AZ5/pb0xc\n3DkcDkdDfvgB2LWLLaCq47bbgFOn9GfNcHHncDgcDTlyBNi6dfRjHB2BBQuA+Hi9DImLO4fD4WhC\nRQVQWwssXTr2sdu2AYcP635MABd3DofD0Yhjx4BNmwBz87GP3boVOHpU92MCuLhzOByORnyTfAkl\nCx/Ek0efRG177ajHhocDLS1AdbXux8XFncPhcCbJj4U/4rTLbdgwPxo2FjaIfi8aVW2q6bWGMTMD\nVq0Czp/X/dg0TRzG4XA4U5La9lo8fGAPZh79HntfWgZCAHtrezx88GGcuP8EiJrQmdhYICEBuPtu\n3Y6Pz9w5HA5nErx47kUssX4A64KXDYVA/mblb3D95nV8n/+92vNWrwbOndP9+Li4czgczgSpbK3E\nVzlfwaXgN1i5crjdwswCe2P34qULL6nN2BgRwTz3xkbdjpGLO4fD4UyQ99Pex71h9yL1nLNA3AHg\n1rm3oq2nDYmViZLnmpsDUVHA5cu6HSMXdw6Hw5kA/Yp+vJ/2Pu6f9xOUlbGNSSMxI2bYE7kHH6d/\nrPYaS5YAKSm6HScXdw6Hw5kA8WXxcJ/ujr7qcISFAZaW4mPuW3Afvs37Fl19XZLX4OLO4XA4MuNA\n3gHcPu92pKYCixdLH+M5wxOLPRfjSNERyX6luOsyES4Xdw6HwxknCqrA9wXfY+fcnbhyRb24A8CO\n4B34sehHyT5vb5ZkrLJSRwMFF3cOh8MZNynVKZhpPRMhziFjivv24O04UnQECqoQ9RGie2uGizuH\nw+GMk+PFx7EtaBs6OoDyciA0VP2x/rP84WLngsvV0mExUVFAaqqOBgou7hyZk5EB/N//AW+/DXR2\nGno0nKnO6dLTWBewDmlpULuYOpKtQVtx7Noxyb6FC4HMTB0MchAu7hzZ8tFHwMaNwMAAcPIke42t\nrzf0qDhTlc6+TlypuYKVs1ciPZ1tRhqLNX5rkFCeINm3YAEXd84UJCkJ+PWv2TbtP/8ZOHCA1Z7c\nvRtQiC1MDkfnXKi4gEiPSEy3mo7sbJbhcSxWzl6JlOoUdPeLyy/5+QGtrUBzs/bHCnBx58iQgQHg\n0UeBt94CQkKG2/fuZaFj+/YZbGicKcyZ0jNY67cWAMYt7jOsZyDMNQxJVUmiPjMzdg1dzd41FndC\nyGZCSD4hpIgQ8qyaY+IIIWmEkGxCSLym9+SYNp98Ajg4AHfcIWw3MwNefx148UXuv3P0T0J5Atb4\nrwGlTNzDwsZ3XpxfHOLL4iX7dGnNaCTuhBBzAG8C2AxgPoB7CCHzVI5xAPAWgFsopWEA7hBdiMMZ\nhFLgpZeAP/5Rutjw4sVAdDTwwQf6Hxtn6tLT34PM+kws8VyCykpg2jTAyWl858b6xuJcuXQaSNmK\nO4ClAK5RSssopX0AvgRwq8ox9wL4llJaBQCUUh3nQuMYM2fOANbWLC2qOp56ikXP6HJ3H4czkvS6\ndAQ7BWOa1bQJzdoBINo7Gqk1qRhQDIj65CzuXgBG7rGqGmwbSRAAR0LIWUJIKiHkAQ3vyTFh3nkH\nePxx6Vm7kthY9l995MTmcAAgqSoJMV4xAICsrImJu6OtI9ynuyOvMU/UFxYG5OSwdSZto2klpvHM\nnSwBLAKwDoAdgEuEkCRKaZHqgXv37h36OS4uDnFxcRoOj2NMtLaykMexLBdCgEceAT79dFjoORxd\nklSdhM2BmwEwv33t2omdH+0djeSqZIS5Cp8KM2cCLi5AcTEQHDy+a8XHxyM+Pn7M4zQV92oAPiM+\n+4DN3kdSCaCRUtoFoIsQcg7AQgCjijtn6nHoEBAXB9jbj33s7t3Mf3/77bE3knA4mpJUlYS9sXsB\nMHF/8smJnR/tFY3k6mTsWbRH1Dd/PpCXN35xV534vvDCC5LHaWrLpAIIIoT4EUKsANwF4JDKMQcB\nrCSEmBNC7ABEA8jV8L4cE6P8Rjn+eOEFmK17HkVNoue+CF9fICgIOH1aD4PjTGnqOurQ2t2KIKcg\n9PcDBQVMkCeCUtylmDePibu20UjcKaX9AJ4AcBxMsL+ilOYRQh4jhDw2eEw+gGMAMgEkA/g3pZSL\nO2eI0yWnsfjdKJTW3sDsgC4s/2A5fij4Yczzdu8GvvlGDwPkTGlSa1IR5RkFM2KGkhLA3Z1Fy0yE\nhe4Lca35Gjp6O0R9shR3AKCUHqWUhlBK51BKXxps20cp3TfimFcopaGU0nBK6T81vSfHdChuLsY9\n396D//H4BnHdr+Ef21/G0fuOYs+hPcisHz2MYNs24NgxHjXD0S3pdemIdI8EAOTnA3PnTvwaVuZW\nCHcNx9Xaq6I+2Yo7hzNZKKV49IdH8ZuVv0FVYiw2bWLtUZ5R+PO6P+ORg49Iho8pCQ4GrKxY9AKH\noyvS6tIQ6cHEvaBAuGt6IkS6RyKjLkPUPm8ee2hoe5LCxZ1jMA4XHUZjZyOeXPoUjh/HkLgDwJ7I\nPbCxsMEX2V+oPZ8QYMsW4OhRPQyWM2VJr0tHhDvLEqaJuEe4RyC9Ll3U7ugI2NoC1dWajFIMF3eO\nQaCU4oWEF/D72N+jtMQcvb3C3NiEELy45kW8kPDCqLN3Lu4cXdLa3Yq6jjoEOQYB0IK414vFHdCN\nNcPFnWMQEisT0d7Tjl3zduHECZbaV3Xj0hr/NXC0dcTRa+rVOy6OFTzokq5DzOFoRGZ9JsJdw2Fu\nZg6A2SeTFfdwt3DkNeShb6BP1MfFnWMyvHf1PTy66FGYETPEx6vfFPJ41ON4O+VttdeZPp3N+JOl\no8w4HI1Iq0sbsmSam4GeHsDDY3LXsrO0g6+DL/Ib80V9St9dm3Bx5+id1u5WfJ//PR5c+CAoBS5c\nAFatkj72ztA7kVydjMpW9ZWEY2OBBOl6CByORoyMlFFaMqOlxhgLdb47n7lzTIKvc7/G+oD1cJnm\nguJiwNycbUqSwtbSFjvn7sT+nP1qr8fFnaMrtLWYqiTCjYs7x4T5Lu873Bl6J4DhWftos6G7w+7G\nlzlfqu1fsQK4fBno7dX2SDlTmd6BXuQ15iHcjVXl0MRvV6JuUdXLC7h5E7hxQ7Prj4SLO0evtHa3\nIrEyEVvmbAEAnD8PrFw5+jlxfnGoaK1AcXOxZL+DA4t5T0nR9mg5U5nCpkLMnjkbdpZ2ALQ0cx+0\nZahKUDshwJw5QNHYmTfGDRd3jl45XHQYq31XY4b1DACj++1KLMwssHPuThzIP6D2mJUrgcREbY6U\nM9XJuZ6DUJfh+NyCgsntTh2J23Q3mBEz1HXUifqCg7m4c4yYA/kHsHPuTgDA9etAfb0wvl0d24K2\njRoSGRPDI2Y42iWnYVjc+/uBkhKWrE5T5rvMR05Djqg9KIiLO8dI6Vf041TJKWwL2gaA+eRLl7IF\n1bFY678Wl6svo62nTbI/OpqLO0e75DbkYr4LS/9YWsoShtnaan7dUJdQ5FyXFvfCQs2vr4SLO0dv\npNakwsfeB27T3djnVCAqanznTrOahuU+y3Gq5JRkf0AAi0HW9hZuztQlpyEHoa5s5l5YqLnfriTU\nJZTP3Dmmxcnik9gQsGHoc0oKsGTJ+M/fOmcrjhQdkewjhL0F8Nk7Rxv0DvSitKUUIU5M0a9dYwue\n2iDUNRS5DeKs50rPXVsJxLi4c/TGqdJT2BDIxJ1SJu7jnbkDwMbAjThdqr46B7dmONqisKkQvg6+\nsLawBsDK4AUGaufaypm7asSMkxP7b1OTdu7DxZ2jFzp6O3C19ipWzWahMZWVgJkZ4O09/mvMdZ6L\nrr4ulN0ok+zn4s7RFqqRMtoUd5dpLrA0s0RtR62gnRDt+u5c3Dl6IaEsAVGeUZhmxUrYKC2ZiWzl\nJoQgzi8O8WXxkv1LlwJXruimkjxnapHbkKszcQeYNaNuUVVbvjsXd45eOFd+DnG+cUOfJ2rJKBlN\n3GfNYkmddFHVhjO1yGnIGYqUGRgAysrYor22mO8sHQ6pzVh3Lu4cvZBYmYgVs1cMfU5NndhiqpLR\nxB0AIiOBtLRJDJDDGUFuQ+5QpEx1NSuoYWenveuPNnPntgzHaOju70ZaXRqivaIBAArFxMIgRxLi\nFIKegR6UtpRK9nNx52hK70AvSlpKhiJltG3JAPoJh+TiztE5V2quYK7z3KGUAyUlwIwZgKvrxK9F\nCEGsbywSyqXTQC5axMWdoxm6jJRREuoqHTETFMTCLrURDsnFnaNzEisTscJn2JLJzgYWLJj89Zb7\nLMelykuSfcqZu7aLDXOmDnkNeUN+O6AbcXe2c4a1ubUoYsbBge2CrROnnpkwXNw5OiexMhErZw+n\nfszOBsLDJ3+9Zd7LcKlKWtxdXFh1plJp14bDGZPCpsIhSwbQjbgDQIhzCAoaC0Tt2vLdNRZ3Qshm\nQkg+IaSIEPLsKMctIYT0E0J2aXpPjvFAKUVihXjmHhY2+WsudF+IkpYStXlmuO/O0YSCpgIEOwUP\nfdaVuM91mitZck9bvrtG4k4IMQfwJoDNAOYDuIcQMk/NcX8FcAyABkWqOMZGYVMhpltNh5e911Cb\npuJuZW6FSI9IpFRLJ3DnvjtHE0bO3CnV8cy9STxznzOH+e6aounMfSmAa5TSMkppH4AvAdwqcdwv\nAHwDoEHD+3GMjOTqZMR4xwx97u1lvyya5sUezZrhM3fOZKGUCmbuzc1M4JWpAbRJiJO0uAcGst8R\nTdFU3L0AjKxcXDXYNgQhxAtM8N8ZbOJLXVOI1JpULPEcDmgvLGT1Um1sNLtujHfMqOJ+9apm1+dM\nTZq6WGIXZztnAMOzdk2KYqtjrrO0LaMtcbfQ8PzxCPXrAH5DKaWEEIJRbJm9e/cO/RwXF4e4uDgN\nh8cxNKk1qbh93u1Dn7OyNLNklCzzXoZHf3gUlFIQld+82bNZ+t+6OpaDm8MZL4VNhQh2Ch76N6Ur\nSwYA/Gf5o7a9Fl19XbC1HE4UrxR3SqUfKvHx8YiPjx/z+pqKezUAnxGffcBm7yNZDODLwS/LGcAW\nQkgfpfSQ6sVGijvH+Okb6ENGfQYiPSKH2jT125V4zPDADKsZzB91FibaJgSIiAAyMri4j0ZXF/Cf\n/7AdmDt3sjeeqU5Bo34WUwFWPjJgVgCuNV8bKsINsN2whDBLSMoOUp34vvDCC5LX19SWSQUQRAjx\nI4RYAbgLgEC0KaUBlFJ/Sqk/mO/+Mylh55geuQ25mD1zNuyt7YfatCXuALDUaylSa1Il+8LD2VsC\nR5rmZlZ39tAhVkJuyxbg3/829KgMj77CIJWEOIeIrBlCWB4bTa0ZjcSdUtoP4AkAxwHkAviKUppH\nCHmMEPKYZkPjGDspNSkCvx3QPMZ9JFGeUWrFfcECLu4jyarPwkvnX8JH6R+hs7cL994LrF4N/Pgj\n8Oc/s+Liv/89kCC98XfKUNhcqLeZO6DbRVWN49wppUcppSGU0jmU0pcG2/ZRSvdJHPswpfQ7Te/J\nMQ5Sa1IR5TmcQObmTaC2Vnu/LIs9FuNK7RXJvvBwIDNTO/cxdt5JeQfrP1mPhs4G7M/Zj5BXl6C6\nrQZ/+9uwpxsYCOzbBzz6KFuvmKooPXcluhb3uc5z5SvuHI46VMU9N5fVobTQdKVnkEUei5BWl4YB\nhTiBe2goUFDALIepzNGio/jzhT8jaU8S/r7p7/h+92G0J92Jvjt2YABCFd+xg9kBH31kmLEaGgVV\n4FrzNcxxZPX0urpYVaSJFJSZKCFOYlsG4OLOkTE9/T3IbchFhHvEUJs2/XYAmGU7C27T3FDYJN6r\nPW0a4Omp3YLDxkZbTxv2HNqDT3d+Cv9Z/gCAb78liGj7PwR5eOCVi6+Iznn+eeCll4C+Pn2P1vBU\ntlbCydYJ062mA2AJ7nx9AXNz3d1TmYJANYEYF3eObMmsz0SQUxDsLIeTYGtb3AFgsefo1sxU9t1f\nvfgqNgRuQKxf7FDba68B//tLgje3vIlXL72Kug5hhqplywAvL+DwYX2P1vAUNBUIIq90bckAgKOt\nI2wsbET/H7i4c2RLak0qojyECdu1FeM+kigPvqgqRUtXC95MeRN7Y/cOteXkADU1wNatgK+DL+4L\nvw+vXXpNdO5PfjI1I2cKmwoR7Kg/v12JVMSMtzezhLq6Jn9dLu4cnXCl9goWey4WtBli5j5VF1U/\nzvgYm+dsHrJjAODzz4F77hm2GX614ld4L+09UQK23buBS5fYg2Aqoe/FVCVSETPm5swSKimZ/HW5\nuHN0QnpdusBvb2oCOjrY7lFtsshjEdLr0iUXVaeqLaOgCryd8jYej3p8qI1SJu733jt83OyZsxHn\nF4cvs78UnG9nB9xyC/Dtt/oasTxQFfeSEv2I+1znuZKpfwMDubjrDUqpaOGDI6Zf0Y/chlyEuw4H\ntOfksFm7tnN0ONg4wH26u9rsenV1QHu7du8pd86Xn4eNhQ2W+ywfaktJAayt2c7dkeyJ3IP3rr4n\nusYddwBff63rkcoLQ87c85u0HzHDxX0c1LTX4K5v7oLNn2zg+LIjnjnxDLr6NDDDTJyCxgJ42XsN\nldUDdGPJKFnssVjSdzc3B+bPZw+WqcRXOV/h3vB7BTl3Dh9moY6qD9dNgZtQ014jKta8cSN766kV\nFgoyWbr7u1HTXjNkYw0MAOXlgL//GCdqAXVFOzTdpcrFfQzKb5Rj+fvLEeIUguvPXEfWz7JQ3lqO\njZ9u5AKvhoz6DCx0Wyho07W4X6nhETMAe2v6Nu9b7J6/W9B+5AiwbZv4eHMzc+yevxtf5wqn6dbW\nwObN7LypQHFzMfwc/GBhxjZhVFUBzs6s5J2u8XfwR017Dbr7uwXtfOauQ3oHerHzq514MvpJ/GHN\nHzDTZia87b3x1R1fwdveGz87/DMUFgIPPMAiM+6/n22cmepk1GUI/HZAt+Ie6RGJjPoMyb6ptqia\nUJYAH3sfBDoO+wn19az4w/Ll0ufcPv92fJsnNtg3bwaOHdPVSOWFoSwZALA0t4Sfgx+Km4VKzsVd\nh/zp3J/gbe+Np2OeFrSbETO8v+N9nC5MxJJ7jyAsjO3qW7gQWLWK5+dIr08XzNwp1a24L3RbiPS6\ndMn1kKk2c/8m9xvRrP3YMWD9esDSUvqc5T7L0djZKNoMtnEjcOrU1Njla0hxB4Bgp2DRulFAALOG\nBsSxAuOCi7saqtuq8WbKm3h729uifOEA0FBjh/b9b2HaHb/A08/0YtEi4Fe/Ar74ArjzTu0k2zdW\nVGfuNTUs5YCrq27u5zLNBdOspqG8tVzUpxT3qbAOTinF0WtHcUvILYL248fZLFwdZsQMO+fuxIG8\nA4J2Dw/Azw9ITtbBYGWGocU9xClE9HC1tWUpf6tUk6iPEy7uavhDwh/w35H/DW97cWIJStlGj1/v\n2ohw7zn4T8Z/hvrWrQN+/WvWPxUERZW6jjr0DvQKvjdtZoJUR4R7BNLr0kXtbm6AmRmLmjF1CpoK\noKAKzHMeLmNMKXuTHKvuzbagbThWLPZgpoo1U9BUoNdUv6qoq6eqiTXDxV2C+o567M/dj1+v+LVk\n/8GDLIrgV78Cfrfqd3jpwkvoVwy/uz71FHDjBvDVV/oasXxQztpHvu3o0pJREuEWgYw6se9OyNSx\nZo5dO4bNczYLvvvi4uH84KMR5xeH1JpUdPR2CNrXrwfOnNHFaOWFoWfuwU7BamPdubhrkX1X9uHO\n+XfCyU5cBkWhYMmV/vhH5mGu8l0Fj+keOFQwXH/EwgL461+BvXsn75cZK/qOlFES4R6B9HrxzB1g\n987O1u395YBS3Edy7hzL2z7W/oJpVtOw1GspzpaeFbQvW8YqWnV2anu08qGlqwVd/V1wn87KdlEq\nD1sG0GwjExd3FXoHevGv1H/hF9G/kOz/4Qcm3reMsDV/GvVT7LsiTF+/bh0LpZpqs3fVnamAHsVd\nwpYBpoa4d/V1IbEyEev81wnaExKA2Fg1J6mwKXATjhcfF7TZ2bFIMFP23YuaiwR1U5uamJXn6Ki/\nMbhOc0Wfog9NnU2Cdj5z1yLHrx1HwKwAhLlKq9FbbwFPPy2cCd0x/w5crb2KkpbhRywhwG9+A/zj\nH7oesbzIqM/AQvfhmbtCwfK4h4bq9r6BjoFo7GzEje4bor7wcNMX96SqJIS5hmGmzUxBu3LmPh42\nz9ksEneAnW/KEWCGtmQAgBAimWOGi7sW+TTrUzyw4AHJvqIiID2dbc0eiY2FDe4Pvx8fpn0oaN+y\nhS3kXb2qq9HKi66+LpS0lAgW9EpL2RuMvf0oJ2oBM2KGcNdwSd89NJQ9YBQK3Y7BkJwrP4dYX+EU\nvaKCVb+aO3d81wh3DUdrdysqWisE7atXs4eEqVLQaNjFVCUhzmJrRrlLdTLBGVzcR9DW04Zj145h\nd+huyf533wUefhiwsRH33Rt+L77K+UoQa21uzqJm3nlHVyOWFzkNOQh2Coa1hfVQmz4sGSXqrBl7\nexZSVlqqn3EYgoTyBKz2FU7RL11iG5fGm8+HEILVvqtxvvy8oH3FCpabprdXW6OVF/qum6qOYEfx\noqqTExP25uaJX4+L+wi+y/sOa/zWwNFWbLYpFMCXXwIPPih9bpRnFAboANLq0gTtjzwCfPONZnmZ\njYWMOsMspiqJcI9Qu1PVlH333oFeXK6+jBU+KwTtly8D0dETu9Zq39VIKBd6MDNnAsHBTOBNETnY\nMoB0OCQhk7dmuLiP4Ovcr3F32N2SfRcvAg4O6r1jQgjuDr1blD7VwwNYvHhqVLaRWkzNytJ9jLsS\n5U5VKUw5HDK1JhXBTsEiv/3yZWDp0oldK9Y3FufKxR7MsmVAUpImo5QnlFIUNRUhyDFoqM1g4j5K\nxAwXdw3o6O3A+fLz2DJni2T/l18Cd0vr/hB3h90tsmYAViDhiy+0NVL5YqgwSCXhbuHIb8xH74DY\nPzDlmfu58nMiS6a/H0hLA6Ki1JykhjDXMNTfrEd9R72gPTraNCNmatprMN1quuDBaChxn+M4B8Ut\nxaLaBFzcNeRk8UnEeMeIZj8A+0X5+mvgrrtGv0aYaxgszSyRWS/MVLVrF8vR0dqqzRHLC0qpKFKm\nt5f9oxzvgp6m2FnawdfBV7Ka/FQT95wcVhhlpvif86iYm5lj5eyVotm7qYp7YVOhoG5qZyfQ0sLq\nyOqbaVbT4GLnIkqjMWeOgcSdELKZEJJPCCkihDwr0X8fISSDEJJJCEkkhCzQ9J664FDhIewI2SHZ\nl5TE7JU5c0a/BiEE24O344fCHwTts2axhGKmnD617EYZZljNgLOd81BbYSErFSa1AK0r1C2qzp3L\nfkFMbVFQQRW4VHVJ5LcnJ0/cklGyevZqkbgHBbGiJ6aW372gqUBQN7WkhOXTMTPQtFcqYsYgM3dC\niDmANwFsBjAfwD2EkHkqh5UAWE0pXQDgRQDvanJPXTCgGMDhwsO4JfgWyf7Dh6VzYUtxS/At+LHw\nR1H7rbcChw5JnGAiqM7aAf1aMkoi3KTF3caG/dKaWkrmwqZCONg4wG26m6B9Mn67kpWzV+JC5QVB\nGyGmOXuXWkwdaxKnS6QiZubMYSmbJ4qmz6elAK5RSssopX0AvgRw68gDKKWXKKVKQyIZgDgTl4FJ\nqkqC5wxP+Dr4SvZPRNxX+a5CQVOByLPcvp0lYOrr03S08iSzPtOgfruShe4Lp1TETHJVMqK9xCEx\nmoj7Io9FKGwqxM3em4L2qSLuhvDblUhFzHh5Mavo5k01J6lBU3H3AlA54nPVYJs69gCQnTlxvPi4\nKCeHkooKlrJ2vCFlVuZW2BCwAUeKhH9NDw8WTnb+vJoTjZyM+gwscBM6bgaZubtHIK02TTK3u0mK\ne7VY3Ds6mEgtmKQBam1hjXDXcFypFVa3mgrifu2agcVdImLGzIyV+5tojhkLDccy7n1ThJA1AB4B\nsELdMXv37h36OS4uDnFj5SnVEqdKTuHFNS9K9h05wtKempuP/3rbg7fjUMEhPBz5sKD9lluYNbN2\nrSajlSeZ9Zmi71AfqX5VcZ/uDitzK1S1VcFnpo+gLzycFVUxJZKrk3H/gvsFbenp7EFmZTX568Z4\nxyCpKkke+780AAAgAElEQVSwULt0KZCaypLhTeT3Qa70DvSiorVCULWquJi9ZRsKqaIdwLDvHh4O\nxMfHIz4+fsxraSru1QBG/gb5gM3eBQwuov4bwGZKaYu6i40Ud33R2t2KrOtZWDFb+plz+DBw770T\nu+aGgA14+vjTGFAMwNxs+LfglltY5Mzrr2syYvnR0duB6rZqwQzo5k32xmOIWZByUVVV3E1t5t7V\n14W8hjxEukcK2jMzWVUwTYjxjsH+nP2CNicnlh8/L0//b2S6oLSlFN723rAyH34KGtqWmT1zNho7\nG3Gz9yamWU0bah/pu6tOfF944QXJa2lqy6QCCCKE+BFCrADcBUCwbEgImQ3gOwD3U0onsSygWxLK\nExDjHQMbC3FIR18fS5i0cePErukxwwMe0z1Eu1UXLGChVpNN4SlXcq7nYJ7LvKHiwgDL5RISwjJo\n6ht1O1UDA1mun44OiZOMkLS6NMxzmQdbS2EVZ21sHIvxjsGlqksie2vJEjZ7NwVULZn+fqCyki28\nGwpzM3MEzgpEUXORoH0yETMaiTultB/AEwCOA8gF8BWlNI8Q8hgh5LHBw34PYBaAdwghaYSQy5rc\nU9ucLD6J9f7rJftSU9mX6iRO6z4m6wPW42TxSUEbIaz4wcmTak4yUuTitytRt1PV3JyFRObmGmBQ\nOkDdYmpm5uT9diW+M31BKUVlW6WgffFi00mEpyruFRWAuztgbT3KSXogxDlEKxEzGkdzUkqPUkpD\nKKVzKKUvDbbto5TuG/z5vymlTpTSyME/k1zD1w2nSk9hfYC0uJ85M3l/fH3AepwqPSVq37DB9MQ9\nsz4TC1zlI+6j5XY3pTQEUouplGpn5k4IQYx3DJKrhCuoixaZlriPzAZ57ZphwyCVSC2q6n3mbuxU\ntVWh4WYDIj0iJfvPnAHWrJnctWN9Y5FclYzOPmEJG2XZMlOq0JRZnymLGHclwU7BqO2oRVtPm6jP\nlHz35OpkRHsLxb28nGXB1EahCeWi6kgiI1llJlP491vQVCCrMEglUouqvr5AdfXENuFNaXE/XXIa\na/zXwIyIv4bubhYrPN5CB6rMsJ6BCPcIJFYkCto9PdmfK1fUnGhkUEqRWZ+JcFfhVNGQ4m5uZo4w\n1zBRGgjAdMS94WYDWrpaBOIEaMeSURLtFY1LVZcEbQ4OgKsrq21g7Mgtxl2JVNEOKysW715WNv7r\nTGlxP19xXlTgQElSEjB/vmZFJtYHrMepErE1s3EjcOLE5K8rJypaK2BnaQeXaS5Dbc3NbNFy9mzD\njWuh20LJwh2mIu5pdWmI9IgUTUy0Ke6LPRcjsz5TUPwdYNaMsU9O2nvacaP7Brzsh7flyEbcB1MQ\nqC5mTzTHzJQW98TKRFFODiWa+O1K1vmvw+nS06J2U/LdM+szRYupOTksNfJ4i0ToAnW+u7c3y63f\n2GiAQWmRq7VXRSGQgHZTLNtb28Pb3ht5DXmCdlNYVC1qLkKQU5Dg4SgXcXe0dYSlmSXqbwp3uU90\nUXXKintjZyNq2msQ7ib9m6ANcV/qtRT5jfki73flSjbzMYUCHlJpB7KyDB8HHeEegfR6sbgTYhqz\n97S6NCzyWCRq1+bMHWBFaFJrhLGPprCoqlpaj1L5iDsgHTEz0UXVKSvuFysvItorWhCbraSri+XC\nXr5cs3tYW1hjidcSXKy8KGifMYPNbE1hK7fcwiCVhLuGI+d6jshSAExE3GvTRDP3ri7myYaESJ8z\nGaTEPTKS/X4Yc01aVb+9rg6ws9N9rd/xIhUxw2fu4ySxQr0lk5rKxHfaNMnuCSGVPhUAYmNNo6K8\nlC0jB3GfYT0D3vbekpVtjD0csq2nDdXt1YI85ADbORoUpFnaAVUWeywW5ZhxcWEiaMw1aaXqpsoh\nDFKJVMQMn7mPkwuVF7By9krJvsREVhRYG0jVpARMQ9w7+zpR3louEBlK5SHuAMsQKeW7G/vMPaMu\nA+Gu4aK3Tm1bMgAQ6RGJrOtZ6BsQpjNdvNi4F1XlGimjRCpiJiCAPVDHG4Y6JcW9u78b6XXpohhh\nJdoU9xjvGGTUZYji3VeuZKGWPT3auY8hyG3IRYhTiCA3R20tSzng6mrAgQ2iLrd7aCgTd4nEkUaB\nusVUXYj7dKvp8J3pi9wG4bZeY/bdKaUoaJRnjLsSqaIddnZst3x19fiuMSXF/UrNFcxznofpVtNF\nfQoFK4atqd+uZJrVNCxwWyDa6TdzJvNGjbmifEad2G+Xw2KqEnURM87O7BelSpTizjgYbTFVF1k4\n1S2qGuvMvf5mPawtrOFoO7zTy9CpflUJnBWI8hvlojemifjuU1LcL1RcUOu3FxQwP9HTU3v3W+1r\nmr67lN+uzVA8TVGKu1Rud2P23ZUx7qpkZWl/5g6oF/e0NON8+1G1ZAD5zdytLazhZe+FkhZhlsGJ\n+O5TUtwTKxPVpvjVpiWjxFR998zr8hZ3zxmeUFAF6jrqRH3G6rt393ejqKkIYa7C16P6epbVUJuT\nEiVSi6oeHsx+M8a3H2MQd4AtqmoSMTPlxF1BFbhYeVHtzF0X4r7CZwUuV19G74AwMcSqVWwnrDGW\n3lOmHZCKcZeLuBNC1Fozxiru2dezEeQUJEpRrbRkdLFxLMI9AtnXs0X/fpUhkcZGQaOwKHZrK0s3\n4uY2ykkGQGpRlc/cR6GgsQAzrGcIth2PRBfiPtNmJoKdgkWvto6OrHyWMXqX1e3VsDCzEBRm7u8H\n8vPZgqVcGE3cjdGWuVp7VdJv15UlA7B1o4BZAci+LnwaGqu4FzYXCiK8lLN2Q+6oliLESbPUv1NO\n3EdLOXD9OvujC3GK9Y1FQpnpWDNSfvu1a+x1fbp4ndpgqCvcERrK1lf6xXucZI3U5iVAN5EyI4ny\njMKVGuEsJDLSOCNm5FY3VR3BTsEobJZO/TuetY4pJ+4XKtTHt1+8CMTE6KY+5CrfVThfIa6OvXq1\ncRbNlsrhLidLRom6mfu0aexBNNEc2Ybmap30zF1XkTJKRltUNSb6Ff0obSlF4KxhNb92jW3+khtS\nKQgcHFgxkevXxz5/yon7aDP3ixe1b8koWTl7JS5WXsSAQrgDYdUq4MIF48uPLZV2QI7iHuIUgorW\nCtzsvSnqMzbfvV/Rj+zr2aJ1Dn3YYYs9FiO1Viju/v7Mr25q0t19tU3ZjTK4T3cXlCYsLASCg0c5\nyUB4zfBCR28HWrtbBe3j9d2nlLjXd9Sj4WYDQl2lfwt04bcrcZ3mCo8ZHqIc425urLSXsfm/6XXp\niHCPELTJUdwtzS0xz2Uesq6Lv2BjC4csaCyA1wwvzLCeIWgvKmK5vnVphy10X4i8hjz09A/vujMz\nAyIijGv2nteQh3ku8wRtchV3QgiCnIImHTEzpcQ9sTIRy32WSxbn6OlhFWaW6rAI4KrZpmHN3Oy9\nifIb5ZjvMl/QLkdxB9hOVVPI7a5uMVXXlgwA2FnaYY7jHNFD0th897zGPMxzNg5xB6QjZsab131q\nifsoycLS0tj/YF3OftSJ+6pVwDnxHifZklGfgfku82FpbjnUdvMmUFMjT+/SVMIh0+qkF1N1GSkz\nksWeiyUXVY1q5q4i7k1NLBRZDukypFBXT5XP3FUYbfNSUhJbTNUlyp2qqjsmV69m4m4su/2kIjZy\nclg6BQtxBmWDoy63e3AwqzlqLHn1R5u560PcozyiRJuZjG1RVdWWKSpi/w7kFgapRCo7JLdlVOjs\n60TW9Sws9ZL2XfQh7r4OvrA2t0ZRs7AApa8vYGPDXg+NAant73K1ZABggdsCZNVniRazrazYL0p+\nvoEGNgEUVIH0unTJtAP6sGUANnNXjZiZOxeorGRlFeUOpRS5DbmCmbucLRlAOmImKIiNe6zJ4JQR\n95TqFIS5hsHO0k6y/9IlYNky3Y9jle8qnC+X9t2NxZqRsgfkLO4zbWbCdZorrjWLpzvGYs2UtpTC\n3toeznbOgvbWVlYyMCBA92NY6LYQ+Y35gkVVS0tWazhTXItcdtS018DW0hZOdk5DbXIX92CnYBQ1\nF0FBhyujODuzN42xSkVqLO6EkM2EkHxCSBEh5Fk1x/xzsD+DECKeeuiB0UIga2rYzEMfyfpXz16N\ncxViFV+1yjgWVfsG+pDXkGcUYZAjUbeZyVjEXV2ysOxsFgKpi70Zqtha2iLIKchoF1WNbTEVYHVs\n7a3tUd02nOeXEDbmsd70NRJ3Qog5gDcBbAYwH8A9hJB5KsdsBTCHUhoE4CcA3tHknpNltM1LSktG\nH76bsc/ccxty4efgh2lWw2WqKDUOcU+rFZvDxpKGYLSdqfr83hd7iK0ZY1lUzWswPnEHpCNmdC7u\nAJYCuEYpLaOU9gH4EsCtKsfsAPAxAFBKkwE4EEL0mqJHQRW4VHVJ7cw9KUk/lgwAzHOeh7aeNlS1\nCdPphYQAnZ1sgU/OSM0ga2tZHnxdZCTUFpHukUirEytQeLgRzdwNGCmjRCoNgbEsquY1ChdTFQq2\noCrHCK+RSGWHDAlh6TNGQ1Nx9wJQOeJz1WDbWMd4a3jfCZHbkAsnWydBkquR6GMxVQkhRHL2Tohx\nxLtLzSDT0tjsTa4RBwATpZSaFFGkkp8f0NzMvGs5I7VpDNBfpIwSqZ2q4eFsUbq3V81JMkHVlqmp\nYbUb5FIUWx1SCcTGM3PXNHBtvMF7qr/2kuft3bt36Oe4uDjExcVNalCqjGbJ9PUxv1CXm5dUUca7\n3xN+j6BdKe7336+/sUyUq3VXsSNkh6BNKe5yxmOGB2wtbFF2owz+s/yH2s3M2IJgdrbudidrSn1H\nPbr7uzF75mxBuyHssAVuC1DQWIDu/u6htMN2diwVQU6OvP8dqIZBGoMlA7CImZMlJ4c+x8fH48yZ\neJw/D4yQTBGains1AJ8Rn33AZuajHeM92CZi72gj1YDEykSsnr1asi8zk/3D1OfTe9XsVfgw/UNx\n+yrgHYOsSIwPBVUgoy5DZMukpQF33GGgQU0AZfKrkeIODFszchV3pRVGVF6NKirYpjsnJzUn6gBb\nS1sEOwUjqz4LS7yWDLUrfXe5intzVzM6+zrhNWPYWDAWcVe1ZeLi4rB0aRw++gj4v/8DXnjhBcnz\nNLVlUgEEEUL8CCFWAO4CcEjlmEMAHgQAQkgMgBuU0noN7zshEivUb166dEl/loySSI9IlN8oR3NX\ns6B9wQLmX48n45shKG4uxizbWYLak4C8f6lHIpXZEJB/xIyh0vyqQ2pRVe6+u3LWPvIBaSzi7u/g\nj9qOWnT2dQ612dkBLi6jr9FpJO6U0n4ATwA4DiAXwFeU0jxCyGOEkMcGjzkCoIQQcg3APgCPa3LP\niVLTXoPWnlbMdZ4r2a/PxVQlFmYWiPGOwYWKC4J2c3M2e5Sr7y61qHfjBnsYyX1RChj23VWRvbir\nWUzVd6SMkihP8U5VuUfMSIVBKnenyh1Lc0sEOQYhv1G42y4kZHTfXeM4d0rpUUppCKV0DqX0pcG2\nfZTSfSOOeWKwfyGlVK8RsYkV6pOFAfpdTB3JqtnGFxKZUp2CJZ5LBG3p6Wz2qI84a01RitLIDSHA\ncDikXNM/6Lsg9lgs9hTXVI2IYIn3FAo1JxkYqTDI/HzjEHcACHUNFVXCCg4ePWLG5HeoXqi4gFWz\nV0n2Xb/OdnnNlZ7U6xR1xTvkvJkppSZF4LMCxmPJAICznTMcbR1FO1Xd3VmkT02NgQY2Cm09bahp\nr0GIU4ioz1C2zMhFVSWzZrGdk0VFo5xoQLIbsgVFxbu7WXFvfWxc1AZhLmHIuZ4jaNP5zF3uXKhU\nHymTnAxER7OICX0T7RWNrOtZoiISUVHsf5jcQvMGFAO4WnsVUZ5RgnZjEndA2ncnRL6ecUZdBsJc\nw2BuJnw16u4GSksNMzGxsbBBiHOIqDaBnK2ZrPoshLsNe1iFhSyQwtJylJNkRJhrGLIb+Mx9iLae\nNhQ0FmCxx2LJfkMspiqxtbRFhHsEkqqSBO1WViwsMzHRMONSR35jPtymuxntYqqSKI8opFSLffdF\ni+S5hV6d356Xx2adVlYGGBTYoqqxbGZq6mxCZ18nfOyHg/Zyc1kIrLGgzpaZsjP3pKokLPJYBGsL\na+l+A/ntSkYr3iE33/1y9WVRRs2uLpZ6NCxMzUkyZInXEtEmHEC+4p5ely6rSBklxpSGIOt6FsJc\nwwSRMnl5xiXu/g7+aOxsRFtP21Cbry/Q0KD+HJMW99E2Lw0MAKmpzJYxFMr87qJ2Ge5UTakRL6Zm\nZ7PZg7X0s1OWLPJYhPS6dFH630WLgCtX1JxkQNQtphoqUkaJuoiZq1fltzCdWZ+JcFfhl5WbC8yb\np+YEGWJuZo55zvOQ25A73GYO7Nql/pwpK+7Z2SwXiqOjZLdeWO6zHCk1KegdEO7bjolhUSidnWpO\nNABSM/crV5goGhMONg7wmO4hCisLCADa2+W1x6B3oBcFjQUiYQIMFymjJNwtHIVNhejqG6504uHB\nBKdKdRujgVH12wHjs2UAZs2oLqp+9pn6401W3PsG+pBSk4LlPssl+w1tyQBMaAJnBeJqrdAPsLNj\nv7jJyQYamArd/d3IbcgV5Ta5fNmwbz6TJcozCperLwva5LiomnM9BwGzAmBraSvqM7QtY2Nhg7nO\ncwWLqoTI05rJup4leED29QElJcYTBqkkzCVM5LuPhsmKe1pdGgJmBcDBxkGyX1/FOcZite9q2ce7\nZ9RlIMQ5RFToxFjFfZn3MtFCNiA/a0adJXP9Oivo7qWaok/PLPYQx7vL7QGpoArkNOQIwiCLiwFv\nb8BW/MyUNVIRM6NhsuJ+vvw8VvpIWzKAYSNlRrJq9irJ4h1y8t0vV1/GUk+hJdPeDpSVGddiqpJl\nPstwseqiqH3xYnktqqbVpiHCTX0mSENn4ZQquye3mXvZjTI42Dhglu2soTZj89uVSNkyo2Gy4j5a\nfHtDA1BfLw9hWuW7CokViaJdkytWMFtGDmlUk6uTRX57aiqwcKHxxAmPZKHbQpTdKENrt3Azgdwi\nZtTN3DMy2HdvaEZbVJULWfVZkoupxua3A4CPvQ86ejvQ1Nk0ruNNUtwppaMupiYmMktGDlvm3ae7\nw9nOWfREdnAAAgPl8YuSWClOvHb5sn7TJGsTS3NLLPJYhORq4aJGUBB78Dc3qzlRjwwoBpBRnyEZ\nBikXcQ93DUdRU5FgUTUggG3Aaxqf/ugcVb8dML4wSCWEEDZ7bxjf7N0kxb2gqQB2lnbwmekj2X/h\nArBSvWOjd1bNXqU2JNLQvnt1WzXae9pF29+NWdwBYLn3clysFFozZmYsR4ocbIX8xny4T3cX2AlK\nDL2YqsTawhrzXOYJatOambEHjxy+Q2AwDNIEImWUSKUhUIdJint8WTzi/OLU9l+4IK/c3eryzMhB\n3BMrWeI11VziRi/uPstxqeqSqF0uvntqTapoXwHAbLrCQnlYioD80/+q7vDt72ffnyHSNmiDcLdw\nUdoHdZikuJ8tO4s1fmsk+zo7WYywnIRpte9qJJQniErArVrFLKSBATUn6oHEikRR7dmaGvY9BgQY\naFBaIMY7BslVyaK1jsWLgRRxdgK9k1KTIsrjA7BcIr6+8on0iPGOEUUeycV3b+1uRW17LUKch986\ni4rY/pbp0w04MA2IcI9Aen36uI41OXGnlOJsqXpxT0lhO/vs7CS7DULArADYWtiKvDQ3N8DV1bC5\nxqX89uRk9nA0dLSGJrhMc4HrNFfBjj+AhXbKYX9Bak2qpLhnZMjDklGyzHuZ6A1ILhEzGfUZWOC2\nABZmwwXn5Pb9TZSFbguRVZ8l2mEthcmJe05DDqZbTYevg69kv9z8diUbAjbgZPFJUbshrZmO3g7k\nNeaJRObCBfZWYews9xH77kFBLMyzttZAgwLbmZp1PQuLPMTbfzMz5bGYqiTEOQQtXS24fnN4a++8\neawEYEeHAQcG4GrtVdF3KLfvb6LMtJkJt+luKGoeO7eyyYn72dKzWOu/Vm2/XMV9Y+BGnCg5IWpf\nvRpISDDAgMDi2xe6LRwqhKxErt/hRFnmvUwk7oQYfvaecz0Hfg5+mG4l9g7kNvM0I2aI9o7Gpcrh\n2bulJVuwzByfNawzpMRdbt/fZIhwj0B63djWjOmJ+yh++8AA27wkp8VUJWv91yKxIhE9/T3C9rXA\n2bOG8d0vVFwQ+e03bzKbaIl4rc/oWO6zHImV4tzKMTGGFXepJG1K5BIGORIpa0YOi6qmOHMHgAi3\nKSjuCqpAQnkC1vhLi3t2NvOxXVz0PLBxMMt2Fua7zBfNJL282AJQqjhLrc6RijpKTmbhgnJZ0NOE\nUNdQNHc1o7qtWtAeHc1yDxkKdX57fT2LlvH2NsCgRkGd727IRdXOvk6UtJRgvstwzGNzM9DWBvj5\nGW5c2iDSIxJpdWM/OU1K3DPrM+Fi5wLPGZ6S/efPy9tO2BCwASeKxdbM5s3AsWP6HUtXXxcuV1/G\nKl+huW4qlgzALIVY31gklAt9r6VL2cPUUFFK6mbuckk7oEq0dzSu1l5F30DfUFtUlGGjjjLrMzHf\nZT6szIermWRksGAKuX1/EyXCPQJptWmi6DpVTErcT5ecVmvJAMCZM8zmkCsbAzfiZIl4UXXTJuD4\ncf2O5VLVJYS7hcPe2l7QLvcH5ESJ84tDfFm8oM3Rkb0t5Yw/jYfW6OrrQkFjARa4iY1huVoK9tb2\n8HPwE8RfL1zIEnS1txtmTKZqyQCA1wwvUFDUddSNepxJifvRa0exec5myT6Fgi1MrlGv/QYnxjsG\nRc1FaOxsFLSvXMkspZYW/Y3lTOkZrPUTPgn7+5ktI8c1i8myxm8NzpadFbXHxBjGmsmoZxk4pdL8\nytFvV6K6OG1lxayZy5dHOUmHpNSkiMprmsJiKsDSEES4R4xpzZiMuHf0diC5OlltpExGBvPaPaUd\nG1lgaW6JWN9YUUikjQ0LPTx1Sn9jOVN6BusC1gnarl5lG2gMWeBE24S6hqKlqwVVbcIKE4YS96Sq\nJER7SedRTk+XrzhJ+e7LlrEABkOQVJWEGG9h2lc5PxwnyngWVTUSd0KIIyHkJCGkkBByghAiSp5O\nCPEhhJwlhOQQQrIJIU9qck91nC09iyWeSzDDeoZkv9wtGSXbg7fjh8IfRO36tGbaetqQWZ+JZd7C\nhPenTgHr1+tnDPrCjJgh1i8WCWVC391Qi6oXKy9KFphR1qs1ZGm90VjmIw4rXbbMMN+h8mEd6ho6\n1NbdzRKGmYq4j2dRVdOZ+28AnKSUBgM4PfhZlT4AT1NKQwHEAPg5IUTr2ZSPXjuKLXO2qO03JnE/\ndu2YYHEKGF5U1Ud9yvPl57HUa6nIGjh50vTEHQDifMW++4IFrFycvjNEXqq6JCnuGRksH4pc69WG\nOIWgq78L5TfKh9qUbz/6rql6ufoyojyjBDtTMzNZ5SVTiPICWC3gKzWjV5bRVNx3APh48OePAdym\negCltI5Smj74cweAPABaNUcopaP67X19LMojLk6bd9UNnjM8McdxjiiRWFAQ+8XOyFBzohY5Xnwc\nGwI2CNpu3mTRD7Gxur+/vonzi0N8ebygzcKCidOFC/obR2VrJXr6exA4K1DUd+UKy3sjVwghosgj\nT09g2jSWz0WfJFUlIcZLaMmkpprG3gwlwU7BaOpqEq3PjURTcXejlNYP/lwPwG20gwkhfgAiAWh1\ni0hhUyH6BvoEpbRGcuUK4O8PODtr8666Y0fIDhwqOCRoIwS47Tbg4EHd3ptSisNFh7EteJug/cIF\ntjHFWBMujUaoayhau1sFs05A/6kflJaMagZOQP7iDoCJu4q9ZQhrJqla7LenpLDwTFPBjJhJZuQU\nHDPWRQY99SyJPztGHkdZ0KXaFzBCyHQA3wB4anAGL2Lv3r1Df+Lj48ca2hCHiw5jy5wtkr8UAHD6\ntLyjZFRRirtqHOtttwEHDuj23gVNBegd6BUVODh5EtiwQc1JRo4ZMcPGwI04Xixc1DCUuEthFOLu\nFyt6A4qJ0e+iqoIqkFyVLBL31FTTEff4+Hjs3bsXPad68Nc//lX9gZTSSf8BkA/AffBnDwD5ao6z\nBHAcwP+Mci06WVa8v4IeLjysvn8FpcePT/ryekehUFDf13xpVn2WoL2/n1IXF0pLSnR371cSX6GP\n/fCYqH3BAkoTE3V3X0PzScYndOeXOwVtXV2U2tlR2tamnzEseXcJPVd2TtTe2UmprS0bj5wZUAxQ\n55edacWNiqG2pCRKFy7U3xjyGvKo/+v+graODvb99fTobxz64Oucr+ktn99CB7VTpKma2jKHADw0\n+PNDAL5XPYCw6fT7AHIppa9reD8Rte21yL6ejXX+6yT7m5vZYsrq1dq+s+4ghGBHyA4cyBNO083N\ngR07gO9F37L2OFx0GNuChJZMZSVQXS2vHPjaZlPgJpwpPYPegeGitTY2bLasj5lnZ18nchpyJNMO\nZGYCISFsPHLGjJgN1SZQEhnJNjPduKGfMSRWJGKZjzDKKy2NFTexslJzkpGyxHMJUmrUbwPWVNz/\nAmADIaQQwNrBzyCEeBJCDg8eswLA/QDWEELSBv9Ir3xOgoMFB7E1aCusLaTDCE6eZIuAcv/FUOWu\n0LvwZc6XktaMrsS9tbsVqTWpor0CP/wAbNvGFhlNFZdpLghyChJkNwT0Z82kVKcg1CVUcvOSMVgy\nSlQjj6ysWFipvham48vjEecbJ2gzJUtmJLNnzh41r7tG4k4pbaaUrqeUBlNKN1JKbwy211BKtw3+\nfIFSakYpjaCURg7+0VqmlAP5B7Br3i61/UeOAFu3autu+mOZzzJ09HYg+7qwUsf69Sxipr5ezYka\n8GPhj4j1i8U0q2mC9kOH2BuDqbM5cDOOXRP+09SXuCeUJyDWVzoUyZjEPdYvVhRWGhenn7TVlFIk\nlCUg1k/4PV66xLx/U4MQgrMPiXdXKzHqHaotXS24VHlp1JQDx44BW9SHv8sWM2LGZu/ZXwrabWyY\n0Fy5ZHAAABgHSURBVH71lfbvuT93P+6cf6egra0NuHgR2LhR+/eTG1uCtuDotaOCtmXL2M7cri7d\n3vts2Vm12UyNSdzDXMPQ1tOG0pbSobbYWGAC8RGTpvRGKfoV/QhyDBpqo9S0kt2pMnKjlipGLe77\nc/Zj05xNkkUNAPZL6eRkvCk+7w67W9Kaue8+4NNPtXuv1u5WxJfFY0eIcIp+4gTLJTNDeuOvSbHU\naykqWitQ014z1DZjBtvVqEtbobu/GynVKVg1W1ze6uZNFiduLDsrzYgZNs3ZJIg8WroUyM8HWlt1\ne++EsgTE+cUJoubKy9kkz99ft/eWI0Yt7v/J/A8eXPCg2v7vvwduuUWPA9Iyke6RsDCzwOVqYfal\ndetYGTNtbg45WHAQcX5xmGkzU9D+zTfM558KWJhZYGvQVhzMF24m2LCBrd3oiqSqJIS5hkmmzkhJ\nYbtl5bozVYrNgZsF4m5tzTYQJYrromiV+PJ4kbWlnLUbe5rfyWC04l7cXIyipiK1lgylwNdfA3fc\noeeBaRFCCB5a+BDeT3tf0G5hAdxzD/DZZ9q711c5X4ksmfZ24OhR4/4OJ8quebvwXf53gjZdi/vZ\n0rOioihKLl1i1pAxsTFwI86WnhVEHsXF6daaoZQivixe5LebsiUzFkYr7p9kfoJ7wu6BpbmlZH9O\nDksWZOyr5A9HPIyvc79Ge48wMfZ99wGffMJeOTWltr0WFysv4ta5twraDx5kC4pOTprfw1jYFLgJ\nl6svo7lrOKnM0qVAaaluFrGB0UtDGqO4S0Ue6VrcC5oKQClFiFOIoD0x0bRSVE8EoxT3fkU/Pkz/\nEP8V8V9qj/n2W+D2243/dcxjhgfi/OJEC6uLFwP29mz3raZ8nPExbp93u2jt4vPPgXvv1fz6xsQ0\nq2lYH7AePxQMZ+a0tGTipI3vWpWbvTdxtfYqVswWKxClxhvpoRp5FB0NFBToLhHb0aKjol3qLS3M\nc4+I0M095Y5RivuhgkPwtvdGpEek2mO++cZ07ISfLPoJ3r36rqCNEOBnPwPefluza1NK8UHaB9gT\nuUfQXlfHhGUqhECqsmvuLnyb962gTVfWzNmys1jitUQyKKC4mPnVPj7av6+u2TxnM45cOzL02dqa\nvQXqqibBseJjIos2IYE9GE15f8ZoGKW4v3H5Dfxi6S/U9qenM7/YGGc8UmwM3IimzibRBpt772Ux\n2JWVk7/2+YrzsDCzEOXieP99YPdultVvqrEteBviy+IFVtjGjUzctZ2+9mjRUWydI70RIynJ+CwZ\nJTHeMajrqENJS8lQ2+bNbA1H23T2deJi5UVRcZlTp0w3H9J4MDpxz6jLQGFTIW6fd7vaYz78EHjo\nIcDM6P520pibmeN/l/0vXr74sqB9+nTmve/bN/lrv5b0Gp5Y+oTgdXZgAHj3XfZmMBVxsHFAnF8c\nvssbXlidM4ftMcjMHOXECUIpxZFrR7AlSHojhjH67UrMzcxxW8htgu9QVzUJzpaexWKPxaJ6v6ac\n7G48GJ38vXjuRTwd87TahdTeXuCLL5i4mxIPRz6MxIpE5DfmC9p/8Qsm7m1tE79mUVMRLlRcEK1d\nHD0KuLuzvCBTlQcXPohPMj8Z+kwIs6i0mXK5oKkA/Yp+hLpIb0S5cAFYLp0k0ijYNU9obwUGsgmJ\nNh+QAAvjVc2HVFHB/H25liXUB0Yl7ul16bhYeRGPL3lc7TE//gjMnw8EBOhxYHrAztIOP1/yc/w1\nUZjiMyiIWQaT8d7/funveGzxY7CztBO0v/wy8KROiiEaD9uDtyOtLg2VrcOe144dLBWDtlBaMlKp\nqpuaWISOsexMlWKN/xoUNBaguq16qE05e9cW/Yp+fJ//Pe6YL1xgO32a7Qcxlbf3yWA0f3VKKZ47\n/RyeXfGsSIxG8uabwE9+oseB6ZGnYp7C4cLDyLmeI2h/7jngtdfYbsbxUtpSiv25+/FktFDFz59n\nGSDvuksbIzZebCxssHv+bnyWNbyZYOVKJrhVVaOcOAEOFR4SFUVRkpDAQvgspV9QjQIrcytsD96O\nA/nD2U23bdPuA/J8+Xn4zPSB/yzhFtQTJ0yzJOREMBpxP1hwEGU3yvCzJeqN4PR0oLCQLQSaIg42\nDnhu1XP4zWlhqdrQUBaq99pr47/W8/HP4xdLfwHXaa6C9j/+EXj22akbYTCSBxY8gI8zPh5K/2Bh\nwZLQ/SCuXz5h6jvqkV6Xjo2B0kl7zp41rgIz6rgr9C58nvX50Oe1a1mh6pqaUU6aAN/mfYs75gln\n7X19rJi8MSYM1CayEnd1yZlaulrw5NEn8fa2t2Flrj4p82uvMQ/amGc7Y/GzqJ8htyFXlL3wL39h\nf//xzCpTa1JxovgEfrnsl4L2EydY+J2prVdMFmVVpJH1bLVlzXyX9x22Bm2FjYV0Lur4eNMQ901z\nNqH0RikKGgsAsBTA27Zpp6LYgGIA3+V9h9vnC4MrEhKYXemp1UrNxoesxP2VV8RtlFI8cugR7Jy7\nU+0WbQC4dg04fNh0LRkl1hbW+Ne2f+GxHx9DW8/wKqq/P/D448Azz4x+fr+iH4/+8Cj+tuFvguiC\n/n7gl79k/w+MKY+JLiGE4OdLfo43L7851LZ5M8uSqWnxia9zv8bu+dKvmA0NLLzVFBa0LcwscH/4\n/fg44+OhtttvZ/tQNOVE8Qn4zPRBsFOwoP3gQeDWW9WcNJWQKs9kiD8AqKMjpUVFw2WkFAoFffrY\n03TZe8tod1/3qCWn7rmH0hdfHKMulQmx5+Ae+tCBh6hCoRhqu3mT0pAQSr/8Uv15e8/upRv+s0Fw\nHqWUvvACpZs2UarSPOVp7W6ls/4yi1a1Vg217dxJ6QcfTP6ade11dOZLM2lnb6dk//79lG7bNvnr\ny42s+izq9aoX7R/op5SysoH29pRev67ZdXfv303fSXlH0KZQUOrjQ2l2tmbXNiagpsyewUV9aCAA\nffVVSleuZLVCu/u66U8O/YQufGchbe5sHvUvl5xMqbs7pe3tGn5LRkR7TzsNezuMvpH8hqA9NZXV\nWS0uFp9zrOgY9XzVk9a01Qjak5IodXWltKpKfA6H0icOP0F/d/p3Q5/376d0w4bJX+9viX+jDx14\nSG3/ww9T+s9/Tv76cmTJu0voDwU/DH2++25K33xz8tdrvNlIZ740k7Z0tQjaU1IonTNnak1SjELc\nU6qu0MitV+i2P/yT+r/uT3fv303bukevTtzbS2l4OKWffabhN2SEXGu6Rt1fcaefZQr/8m+9RWlw\nMKUNDcNtlyovUZeXXUQFmMvLKfX0pPT77/UxYuOkoLGAurzsMvRvsbOTUgcHSmtrJ34thUJB5745\nl54vPy/ZPzDAJioj32BNgU8yPqFrP1479PnYMUqjoiZ/vVcSX6H3fXufqP2ppyj9/e8nf11jxCjE\nfcE7C+i8fy6gtvc+SP/46dlx/cV+9StKt26dWk/qkWTXZ1PPVz3pn879aei1l1JKn3uO0tBQJt5f\nZX9FnV92pkcKjwjOLSqiNCCA0tde0/eojY97vrmH/uX8X4Y+P/DA5GbXiRWJNOSNEJEtpuTKFfZg\nNjV6+nuo56ueNL02nVLK3s69vCjNypr4tXr7e6nP331oSnWKoL2vj1I3N0oLC7UxYuPBKMRdyeXL\nlDo7U3r8+Oh/qY8+otTfn9LGxkl+KyZCZWsljf0wli54ZwH94OoHtKipiJbfKKeP/PV7avXIJur6\nYhBNKrs6dHx3N6X/+hf7jvftM+DAjYic6znU9W+utKOng1JK6ZEjlEZHT/w69393P335wstq+//w\nB0r/538mO0p586dzf6IPHnhw6PP/+3+UPvPMxK/zeebndPWHq0Xtk/1/YuwYlbhTSumFC8w7fuMN\n8ax8YIDNNj09Kc3J0cK3YwIoFAp6MP8gvWP/HdTvdT/q+aonXffxOvrct/to3Loe6uTEfOINGyh1\ndGT/TU839KiNi7u+vou+mMBW7fv6Jj7zLGspo45/dRT5xCOJjqb0xAlNRypPWrpaqPPLzjSvIY9S\nymbYLi7M5hovA4oBuuCdBfRQ/iFR3+23M0tyqqFO3AnrMzyEEKo6lsJC4IEHgJ4e4JFHWPKmqirg\ngw9Y8qHPP2f5KjhjU1UFZGWxnxcu5DHAk6G0pRRL/r0EGT/NgJe9F55/nuUveeON8Z3/9LGnYWFm\ngb9t/Jtkf2Ul+39TW2u64ah/ufAXXKm9gq93fw2A7RvYvn38IcyfZX6GNy6/gUt7LgnSNii/u/Ly\nqVHvdySEEFBKRTksZC3uABPxI0dY7GplJeDqyv5B3HYbYG5ugIFypjS/O/M7lN4oxWe7PhsSlMrK\nsVMjV7ZWImJfBDJ/mgkvey/JY157jT2AP/hABwOXCZ19nQh+Ixhf3P4FVvmuQnw8yz6akzN2Hpiu\nvi6Evh2KD2/9UFRO77nnWPqNf/xDd2OXK+rEXRMbxRHASQCFAE4AcBjlWHMAaQB+GOUY3b23cDha\nor2nnQb8I4AeyDtAKaV0+3ZK339/7PMePPAg/e3p3456THQ0iyIxdb7N/ZaGvBFCu/q6qEJBaWTk\n+KK1nj35LN29f7eovbWV2TtTbSFVCdTYMprsUP0NgJOU0mAApwc/q+MpALkA5PGawOFMkulW0/HJ\nzk/w0x9/ipr2Gvz858Drr4+eo/x0yWmcLjmNZ1c8q/aYsjKW+mHtWu2PWW7smrcLoa6h+O3p34IQ\n4IUXgN/+ltURUEdSVRI+TP8Qb2wRe2D//CewaRNLOcAZRhNx3wFAuaf4YwC3SR1ECPEGsBXAewCM\nvKIph8NyzjwZ/SR2fLEDK9fchJmZ+gpDTZ1NePjgw3h/x/uYYa3eDP7oI+DOO007L9JI3t3+Lr4v\n+B6fZX6G7dsBBwdW8F2KmvYa3LH/Drx3y3twm+4m6GtuZlbM73+vh0EbGZP23AkhLZTSWYM/EwDN\nys8qx30N4M8A7AE8Qym9Rc316GTHwuHoG0pZzqPyG+W43+IAPto3E+fOCY+52XsTGz/diNWzV+Ol\n9S+pvVZ/P+Dnx9aWplJxiaz6LKz/ZD3e2PIGZrffiV27gOxswNFx+Jiqtiqs+8867Incg1+v+LXo\nGj/9KVt7e+stPQ5cZqjz3EdN7EoIOQnAXaLrtyM/UEopIUSkzISQ7QCuU0rTCCFxYw1y7969Qz/H\nxcUhLm7MUzgcg0AIwXu3vIcnjz6JV8qWobXnfSQkLEPs4DpfbkMu7vvuPkR5/P/27j/W6rqO4/jz\nlebQaeJdCyVEnMNBSxTSlg2CJJsIEmUDXaHiIIdtouOnq2WOTLH4o625DKR0LSdDlm4iQeqdQJSB\nkBfh8sN5S3MCCRTiDyTf/fH9MI+Hc7iHe88995zvfT227/h8P9/P5/v98rkf3vfwOZ/v93Mp94y+\n57jnevrpbBHsnhTYAS7qcxGrvruKcY+O49rBf+aaiT/i9tubeOQR+DA+ZNnWZcxYOYNZl89i5pdn\nHlN/zZrsDZ1bt3bDzXej5uZmmpub2y3XmU/urcCoiHhT0jnAcxExqKjMT4HJwBGgF9mn98cj4oYS\n5/Mnd2s4EcFjLz/G95+YxXv/PodvjRjM6wdfY8ueLdw96m6mXzq95EpLha64AqZMyab99kT73t3H\nnNVzWLZ1GR/8cygX9D+dfadsou8ZfVn49YWMOG/EMXX27MlWqXrwQb+3vepTISXdD7wVEQskzSOb\nLVP2S1VJI/GwjOXU4SNHGDJuLSMntDFhdB9GDhh53BXDjlq3Lgvq27f3nPH2cva9u48n/raBO+a8\nwz13fI5bJw4s+YvxwIFslaWxY7MvY3u6rgjuTcBSoD/QBkyMiAOS+gKLImJsUfmRwMyIGF/mfA7u\n1tDWroXrrsvmqp91zLdPx4rIgtSkSflfh+BEbNgA11wDM2Zk6xMUrgq2eTNcfz2MGQMLF2YLl/d0\nDfsQk1kjue022L+//MyPQsuXZ7M8Nm3yp/ZibW1wyy3Q2poF+qamLOi/+CIsWODVwgo5uJvVwKFD\n2VjwrFkwdWr5cnv3wrBh2S8Bzxsob/NmePZZOHgQBg3KhmJOP72776q+OLib1cjOnTB8ODz0UPbe\nlGLvv58Fqcsug3vLz5A0q0i54F5Xa6ia5cHAgdkUvalTs6dXC5+83LMnG2bo3Rvmz+++e7T88yd3\nsy6yaxdMm5a9WGz06GxoYeXK7MvT+fM9zm7V4WEZs24QARs3wgsvwKmnwpVXQr9+3X1XlicO7mZm\nOeQxdzOzHsTB3cwshxzczcxyyMHdzCyHHNzNzHLIwd3MLIcc3M3McsjB3cwshxzczcxyyMHdzCyH\nHNzNzHLIwd3MLIcc3M3McsjB3cwshxzczcxyyMHdzCyHOhzcJTVJWi1ph6RVknqXKddb0jJJ2yRt\nlfSljt+umZlVojOf3OcBqyPiQuCZtF/KL4AVETEYGAJs68Q1rULNzc3dfQu54basLrdnbXQmuI8H\nHk7ph4EJxQUknQmMiIglABFxJCL+04lrWoX8D6h63JbV5fasjc4E9z4RsTuldwN9SpQ5H9gr6TeS\nXpS0SNJpJ3qh43UGHyutra2tZterp793o7dlTzjm9qzusXKOG9zTmHpLiW18Ybm0snWp1a1PBoYB\nD0TEMOAQ5YdvyqqnhmyUYw7u1TvmYFTdY27P6h4rR1lcPnGSWoFREfGmpHOA5yJiUFGZs4H1EXF+\n2h8OzIuIcSXO17EbMTPr4SJCxXknd+J8TwI3AgvSn38occE3Jb0m6cKI2AF8DXi50pszM7OO6cwn\n9yZgKdAfaAMmRsQBSX2BRRExNpW7GFgMnAK8Akzxl6pmZl2rw8HdzMzqV4dmy0i6SlKrpJ2S5qa8\nSh9qOqZuLevXG0lLJO2W1FKQ97P00NffJS1PU0pL1XVbluD+WT3unw0sIk5oA04CdgEDgE8Cm4HB\nwP3AnFRmLnBfpXXTsS6vX48bMAIYCrQU5F0JfCKl73Nbun+6fzZ+W9b8Z9eBH/blwMqC/XnAnUAr\n2dx3gLOB1grrzkvpLq9fr1vqwC1ljn0T+J3b0v3T/bPx27KWW0eGZT4LvFaw/3rKK/lQk6S+kp5q\npy5dVT8HbgZWgNuyQu6fteX+Wac6EtxLfQP7sbzIfqVGSr8RaeZMiboqdb5q1m9kkn4AHI6I34Pb\nskLunzXi/lnfOhLc/wWcW7B/bsrbreyhJZQ91LSngrr9Uh41qt8wJN0EXA18p0wRt2Vp7p814P5Z\n/zoS3DcAAyUNkHQKMAl4go8eaoIyDzWVqftkOlaL+g1B0lXAbOAbEfFemWJuy9LcP7uY+2eD6MhA\nPTAG2E72bfadKa8J+BOwA1gF9E75fYGnjle3K+vX+wY8CrwBHCYbY7wZ2An8A9iUtgfclu6f7p+N\n3Za13vwQk5lZDnmZPTOzHHJwNzPLIQd3M7MccnA3M8shB3czsxxycDczyyEHd2so6aGWlvZLdtn1\nF0kaXIXznClpejXuyawUB3czQFJFS05GxLSI2FaFS54F3FqF85iV5OBujegkSb+WtEXSHyX1ApB0\niaS/FCwi0TvlN0v6Qkp/WtKrKX2TpCclPQOslnS2pOclbZLUomxB949J5xqW0m9L+omkzZLWS/pM\nifI/ljSzYL9F0nlk70G/IF1rQTo2V9JL6Xz3Vr3VrEdxcLdGNBD4ZUR8HjgAXJvyHwFmR8TFQAtw\nV8o/3psDhwLXRsRXyV6CtTIihgJDyBaIKFZ4ntOA9RFxCfA8MK2d8oV5c4FXImJoRMyVNAYYD3wx\nne/+MvdrVhEHd2tEr0bESym9ERgg6VPAmRGxJuU/DHylgnOtiogDKf0CMEXSXcCQiHi7nbqHI+Lo\n+8c3ki1qUSkV7Y8GlkR6EVdE7D+Bc5kdw8HdGtH7Ben/kS3JVqwweB7ho77eq6jcO0cT6RfDCLJX\ny/5W0uR27uODgvSHQKlx+8Jrl7p+uXs26xQHd8sDRcR/gf0F4+STgeaUbgMuTelvlz2J1B/YGxGL\ngcVkQzad1QYcHaMfBpyf8g8CZxSUW032v4ZTU9mzqnBt68EqmiFgVmeKx7GP7t8I/ErSacArwJSU\n/3NgqaTvAU8VlC8eix8FzJb0AVnwveEE7qPcuP7jwA2StgB/JXuFLRHxlqR1aVrnijTufgmwQdLh\ndJ8/bOf6ZmX5lb9mZjnkYRkzsxxycDczyyEHdzOzHHJwNzPLIQd3M7MccnA3M8shB3czsxxycDcz\ny6H/A/Izq2NA54AkAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1121a1d10>"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import shapely.geometry\n",
      "import shapely.wkt\n",
      "geom = shapely.wkt.load(open('coordinates.txt'))\n",
      "\n",
      "\n",
      "# convert geometry to multipoint\n",
      "points = shapely.geometry.MultiPoint(geom.coords)\n",
      "\n",
      "# lookup bounds\n",
      "lon_min, lat_min, lon_max, lat_max = points.bounds\n",
      "\n",
      "# expand for interpolation\n",
      "window = 4\n",
      "# lon can be < 0, use modulus\n",
      "lon_min = (lon_min - window) % 360\n",
      "lon_max = (lon_max + window) % 360\n",
      "if lon_max < lon_min:\n",
      "    lon_max += 360\n",
      "\n",
      "# lookup the indices in the wrapped array\n",
      "lon_min_idx = bisect.bisect(lon, lon_min) - window\n",
      "lon_max_idx = bisect.bisect(lon, lon_max) + window\n",
      "\n",
      "# lookup lats\n",
      "lat_min_idx = bisect.bisect(lat, lat_min) - window\n",
      "lat_max_idx = bisect.bisect(lat, lat_max) + window\n",
      "\n",
      "lon_min_idx, lon_max_idx, lat_min_idx, lat_max_idx, lon.shape, lat.shape\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 60,
       "text": [
        "(1418, 1465, 555, 567, (2880,), (721,))"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "arr = A\n",
      "lonslice = slice(lon_min_idx, lon_max_idx)\n",
      "latslice = slice(lat_min_idx, lat_max_idx)\n",
      "Lat = lat[latslice]\n",
      "Lon = lon[lonslice]\n",
      "\n",
      "if lon_max_idx > arr.shape[1]:\n",
      "    # first part\n",
      "    arrs = [\n",
      "        arr[:, lon_min_idx:arr.shape[1]+1, latslice],\n",
      "        arr[:, 0:(lon_max_idx % arr.shape[1]) , latslice]\n",
      "    ]\n",
      "    arr = np.concatenate(arrs, axis=1)\n",
      "else:\n",
      "    arr = arr[:,lonslice, latslice]\n",
      "\n",
      "arr.shape, Lon.shape, Lat.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 63,
       "text": [
        "((13, 47, 12), (47,), (12,))"
       ]
      }
     ],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xy = pandas.DataFrame(np.array(points), columns=['lon', 'lat'])\n",
      "xy['lon'] %= 360\n",
      "xy = xy.sort(['lon', 'lat'])\n",
      "xy.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>lon</th>\n",
        "      <th>lat</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>12</th>\n",
        "      <td> 0</td>\n",
        "      <td> 49.740272</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>31</th>\n",
        "      <td> 0</td>\n",
        "      <td> 49.756716</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>52</th>\n",
        "      <td> 0</td>\n",
        "      <td> 49.773155</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>74</th>\n",
        "      <td> 0</td>\n",
        "      <td> 49.789587</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>97</th>\n",
        "      <td> 0</td>\n",
        "      <td> 49.806014</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 164,
       "text": [
        "    lon        lat\n",
        "12    0  49.740272\n",
        "31    0  49.756716\n",
        "52    0  49.773155\n",
        "74    0  49.789587\n",
        "97    0  49.806014"
       ]
      }
     ],
     "prompt_number": 164
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ampl =  []\n",
      "for i in range(A.shape[0]):\n",
      "    f = scipy.interpolate.RectBivariateSpline(Lon, Lat, arr[0].squeeze(), kx=2, ky=2, s=0)\n",
      "    f.ev(xy['lon'], xy['lat'])\n",
      "\n",
      "f.ev(xy['lon'], xy['lat'])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 174,
       "text": [
        "array([ 1.76363885,  1.77278755,  1.78447549, ...,  2.64448707,\n",
        "        2.64929895,  2.63500672])"
       ]
      }
     ],
     "prompt_number": 174
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "results = []\n",
      "\n",
      "circular = False\n",
      "for i in range(A.shape[0]):\n",
      "    if circular:\n",
      "        f_cos = scipy.interpolate.RectBivariateSpline(Lon, Lat, np.cos(arr[i]), kx=2, ky=2, s=0)\n",
      "        f_sin = scipy.interpolate.RectBivariateSpline(Lon, Lat, np.sin(arr[i]), kx=2, ky=2, s=0)\n",
      "        result = np.arctan2(f_sin(lon_i, lat_i), f_cos(lon_i, lat_i))\n",
      "    else:\n",
      "        f = scipy.interpolate.RectBivariateSpline(Lon, Lat, A[i], kx=2, ky=2, s=0)\n",
      "        result = f(lon_i, lat_i)\n",
      "    results.append(result.squeeze())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "TypeError",
       "evalue": "x dimension of z must have same number of elements as x",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-98-a54eb36e10f4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marctan2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf_sin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon_i\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf_cos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon_i\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m         \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRectBivariateSpline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLon\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mLat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mA\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mky\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon_i\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \u001b[0mresults\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqueeze\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/baart_f/.virtualenvs/main/lib/python2.7/site-packages/scipy/interpolate/fitpack2.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, x, y, z, bbox, kx, ky, s)\u001b[0m\n\u001b[1;32m    958\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'y must be strictly ascending'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    959\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 960\u001b[0;31m             raise TypeError('x dimension of z must have same number of '\n\u001b[0m\u001b[1;32m    961\u001b[0m                             'elements as x')\n\u001b[1;32m    962\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mTypeError\u001b[0m: x dimension of z must have same number of elements as x"
       ]
      }
     ],
     "prompt_number": 98
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lonslice, latslice"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 152,
       "text": [
        "array([ 49.74027236,  49.74027236,  49.74027236, ...,  50.70094289,\n",
        "        50.70094289,  50.70094289])"
       ]
      }
     ],
     "prompt_number": 152
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}