{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Tutorial on the Analytical Advection kernel in Parcels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While [Lagrangian Ocean Analysis](https://www.sciencedirect.com/science/article/pii/S1463500317301853) has been around since at least the 1980s, the [Blanke and Raynaud (1997)](https://journals.ametsoc.org/doi/full/10.1175/1520-0485%281997%29027%3C1038%3AKOTPEU%3E2.0.CO%3B2) paper has really spurred the use of Lagrangian particles for large-scale simulations. In their 1997 paper, Blanke and Raynaud introduce the so-called *Analytical Advection* scheme for pathway integration. This scheme has been the base for the [Ariane](http://stockage.univ-brest.fr/~grima/Ariane/) and [TRACMASS](http://www.tracmass.org/) tools. We have also implemented it in Parcels, particularly to facilitate comparison with for example the Runge-Kutta integration scheme.\n", "\n", "In this tutorial, we will briefly explain what the scheme is and how it can be used in Parcels. For more information, see for example [Döös et al (2017)](https://www.geosci-model-dev.net/10/1733/2017/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Most advection schemes, including for example Runge-Kutta schemes, calculate particle trajectories by integrating the velocity field through time-stepping. The Analytical Advection scheme, however, does not use time-stepping. Instead, the trajectory within a grid cell is analytically computed assuming that the velocities change linearly between grid cells. This yields Ordinary Differential Equations for the time is takes to cross a grid cell in each direction. By solving these equations, we can compute the trajectory of a particle within a grid cell, from one face to another. See [Figure 2 of Van Sebille et al (2018)](https://www.sciencedirect.com/science/article/pii/S1463500317301853#fig0002) for a schematic comparing the Analytical Advection scheme to the fourth order Runge-Kutta scheme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the Analytical scheme works with a few limitations:\n", "1. The velocity field should be defined on a C-grid (see also the [Parcels NEMO tutorial](https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_nemo_3D.ipynb)).\n", "\n", "And specifically for the implementation in Parcels\n", "2. The `AdvectionAnalytical` kernel only works for `Scipy Particles`.\n", "3. Since Analytical Advection does not use timestepping, the `dt` parameter in `pset.execute()` should be set to `np.inf`. For backward-in-time simulations, it should be set to `-np.inf`.\n", "4. For time-varying fields, only the 'intermediate timesteps' scheme ([section 2.3 of Döös et al 2017](https://www.geosci-model-dev.net/10/1733/2017/gmd-10-1733-2017.pdf)) is implemented. While there is also a way to also analytically solve the time-evolving fields ([section 2.4 of Döös et al 2017](https://www.geosci-model-dev.net/10/1733/2017/gmd-10-1733-2017.pdf)), this is not yet implemented in Parcels. \n", "\n", "We welcome contributions to the further development of this algorithm and in particular the analytical time-varying case. See [here](https://github.com/OceanParcels/parcels/blob/master/parcels/kernels/advection.py) for the code of the `AdvectionAnalytical` kernel." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, we will show how this `AdvectionAnalytical` kernel performs on one idealised time-constant flow and two idealised time-varying flows: a radial rotation, the time-varying double-gyre as implemented in e.g. [Froyland and Padberg (2009)](https://www.sciencedirect.com/science/article/abs/pii/S0167278909000803) and the Bickley Jet as implemented in e.g. [Hadjighasem et al (2017)](https://aip.scitation.org/doi/10.1063/1.4982720).\n", "\n", "First import the relevant modules." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "from parcels import FieldSet, ParticleSet, ScipyParticle, JITParticle, Variable\n", "from parcels import AdvectionAnalytical, AdvectionRK4, plotTrajectoriesFile\n", "import numpy as np\n", "from datetime import timedelta as delta\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Radial rotation example\n", "\n", "As in [Figure 4a of Lange and Van Sebille (2017)](https://doi.org/10.5194/gmd-10-4175-2017), we define a circular flow with period 24 hours, on a C-grid" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def radialrotation_fieldset(xdim=201, ydim=201):\n", " # Coordinates of the test fieldset (on C-grid in m)\n", " a = b = 20000 # domain size\n", " lon = np.linspace(-a/2, a/2, xdim, dtype=np.float32)\n", " lat = np.linspace(-b/2, b/2, ydim, dtype=np.float32)\n", " dx, dy = lon[2]-lon[1], lat[2]-lat[1]\n", "\n", " # Define arrays R (radius), U (zonal velocity) and V (meridional velocity)\n", " U = np.zeros((lat.size, lon.size), dtype=np.float32)\n", " V = np.zeros((lat.size, lon.size), dtype=np.float32)\n", " R = np.zeros((lat.size, lon.size), dtype=np.float32)\n", "\n", " def calc_r_phi(ln, lt):\n", " return np.sqrt(ln**2 + lt**2), np.arctan2(ln, lt)\n", "\n", " omega = 2 * np.pi / delta(days=1).total_seconds()\n", " for i in range(lon.size):\n", " for j in range(lat.size):\n", " r, phi = calc_r_phi(lon[i], lat[j])\n", " R[j, i] = r\n", " r, phi = calc_r_phi(lon[i]-dx/2, lat[j])\n", " V[j, i] = -omega * r * np.sin(phi)\n", " r, phi = calc_r_phi(lon[i], lat[j]-dy/2)\n", " U[j, i] = omega * r * np.cos(phi)\n", "\n", " data = {'U': U, 'V': V, 'R': R}\n", " dimensions = {'lon': lon, 'lat': lat}\n", " fieldset = FieldSet.from_data(data, dimensions, mesh='flat')\n", " fieldset.U.interp_method = 'cgrid_velocity'\n", " fieldset.V.interp_method = 'cgrid_velocity'\n", " return fieldset\n", "\n", "fieldsetRR = radialrotation_fieldset()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now simulate a set of particles on this fieldset, using the `AdvectionAnalytical` kernel. Keep track of how the radius of the Particle trajectory changes during the run." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: Particle initialisation from field can be very slow as it is computed in scipy mode.\n" ] } ], "source": [ "def UpdateR(particle, fieldset, time):\n", " particle.radius = fieldset.R[time, particle.depth, particle.lat, particle.lon]\n", "\n", "class MyParticle(ScipyParticle):\n", " radius = Variable('radius', dtype=np.float32, initial=0.)\n", " radius_start = Variable('radius_start', dtype=np.float32, initial=fieldsetRR.R)\n", "\n", "pset = ParticleSet(fieldsetRR, pclass=MyParticle, lon=0, lat=4e3, time=0)\n", "\n", "output = pset.ParticleFile(name='radialAnalytical.nc', outputdt=delta(hours=1))\n", "pset.execute(pset.Kernel(UpdateR) + AdvectionAnalytical,\n", " runtime=delta(hours=24),\n", " dt=np.inf, # needs to be set to np.inf for Analytical Advection\n", " output_file=output)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the trajectory and calculate how much the radius has changed during the run." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yV5fn48c91MskgE8IIBEKCylTCkqFB0apFcdZVi3611FXt+tb1tfqrWu23/bbVOlocVasVrYpSZwGJLIMaZK+EQCAJAhmMAJnn+v3xPEkjhpFxck6S6/16nVfOuZ91nZPkuc59389z36KqGGOMMQAefwdgjDEmcFhSMMYY08CSgjHGmAaWFIwxxjSwpGCMMaaBJQVjjDENLCmYTkdE7hWR505gvRdF5OH2iOkox/9QRGb46/iN4jihz8t0DWL3KZj2JiLbgCSgDjgIfAD8WFUrWrCvTOAVVU1uwbYvAoWq+j8t2FaBdFXNa+62bUlEHgTSVPX7/ozDdB5WUzD+cqGqRgGjgDFAS07MwW0eVRsJ5Nga6yhxmvZjScH4laoWAR8CwwBE5AYR2SAiB0QkX0R+VL+uiGSKSKGI3CUiXwOvudv2EZEK99FHRB4UkVcabTdJRJaJyF4R2SEi1zcVi4hME5GV7nrLRGTEUdZb5D5d5R7zyiZi+5uIxInIeyKyR0TK3efJjfaTJSI3NXr9X+57LxeRj0UkpdGyoSIyT0TKRGSX2+RzHnAvcKUbxyp33T4iMtddN09EfthoPw+KyJsi8oqI7Aeub+LzGt/o81rl1sbql13v/l4OiMhWEbn26L9d0xFZUjB+JSL9gAuAr9yi3cA0oDtwA/BHERnVaJNeQDyQAvwAOB8oVtUo91F8xP774ySOPwM9gFOBlU3EMQp4AfgRkAD8FZgrImFHrquqZ7hPR7rHfL2J2Gbi/H/9zX3dHzgMPHmUz+FinBP8pW6ci3GSHiISDcwHPgL6AGnAAlX9CPgN8Lobx0h3d68Bhe66lwO/EZGzGx1uOvAmEAu8ekQcfYH3gYfd9/IL4C0R6SEikcATwPmqGg1MaOqzNB2bJQXjL++IyF5gCfApzskNVX1fVbeo41Pg38DkRtt5gQdUtUpVD5/Aca4F5qvqa6pao6qlqtrUieyHwF9Vdbmq1qnqS0AVML4Z7+kbsbnHektVD6nqAeAR4MyjbPsj4FFV3aCqtTifx6lubWEa8LWq/p+qVqrqAVVd3tRO3CQ7CbjLXXcl8BxwXaPVPlPVd1TV28Rn+H3gA1X9wF0+D/gSJ3HXv8dhItJNVXeq6rpmfD6mA7CkYPzlYlWNVdUUVb21/uQkIueLSLbb9LEX52SU2Gi7Papa2Yzj9AO2nMB6KcDP3SaTve6x++F82z5R34hNRCJE5K8iUuA21SwCYkUk6CjHf7zRscsAAfo24z3gxlvmJqF6Be5+6u04xvYpwBVHfA6TgN6qehC4ErgZ2Cki74vIyScYl+kgLCmYgOE21bwF/B5IUtVYnCuTpNFqR14ud7zL53YAg07g8DuAR9xEVf+IUNXXTjD8pmL5OXASME5VuwP1zU7Ct+0AfnTE8bup6rLjvIcjj1kMxLtNTvX6A0XH2ObIOP5+RByRqvoYgKp+rKrnAL2BjcCzx9iX6YAsKZhAEgqEAXuAWhE5Hzj3ONvsAhJEJOYoy18FporI90QkWEQSROTUJtZ7FrhZRMaJI1JEvnvEyfXI46YeJ7ZonH6EvSISDzxwjHX/AtwjIkMBRCRGRK5wl70H9BKRn4hImIhEi8i4RnEMEBEPgKruAJYBj4pIuNtZfiNH9B0cwyvAhSLyHREJcveRKSLJIpIkIhe5fQtVQAXOZcWmE7GkYAKG2+RxB/AGUA5cA8w9zjYbcTpW893mjj5HLN+O0wT1c5wmmZXAyCb28yVOv8KT7rHzgOuPcegHgZfcY37vKOv8CegGlADZOB3FR3sfc4DfArPdpqa1OJ3o9Z/LOcCFwNdALjDF3fSf7s9SEVnhPr8aGIBTa5iD088x7xjvpXEcO3A6ou/FSc47gP/GOVd4cD7HYpzP8kzg1hPZr+k47OY1Y/zEvbT1OVV92d+xGFPPagrG+IGIROA0P231dyzGNGZJwZh2JiI9cZqBPsW5JNeYgGHNR8YYYxpYTcEYY0yDDj8YVmJiog4YMKBF2x48eJDIyMi2DagNWFzNY3E1j8XVPJ01rpycnBJV7fGtBaraoR8ZGRnaUgsXLmzxtr5kcTWPxdU8FlfzdNa4gC+1iXOqNR8ZY4xpYEnBGGNMA0sKxhhjGlhSMMYY08CSgjHGmAZtlhTcERW/EpH33NcDRWS5iOSKyOsiEuqWh7mv89zlAxrt4x63fJOIfKetYjPGGHNi2rKmcCewodHr3wJ/VNV0nFEnb3TLbwTKVTUN+KO7HiIyBLgKGAqcBzx9lMlIjAloOQXl3DdnDS+trSSnoNzf4RjTLG1y85o4k5F/F2e6wZ+JiABn4Qx9DPASzlDDz+AMy/ugW/4m8KS7/nRgtqpWAVtFJA8YC3zWFjEa4wt1XqXsYDWlB6soOVDNF9vKeHJhHnVeZ/iYRX9dxv3ThjBhUCIJkaHERoQS5Glqjh1jAkObjH0kIm8Cj+JMKvILnHHos93aQP28sR+q6jARWQucp6qF7rItwDicRJGtqq+45c+727zZxPFm4kyMTlJSUsbs2bNbFHdFRQVRUVEt2taXLK7maeu4ar3KgWplX5Wy3/15oFrZV63sb1S2v1o5UH38qd8aEyAqFLqHCt1DhWj3Z/cwITrE+dm4vFswON+Z2k5X+T22lc4a15QpU3JUdfSR5a2uKYjINGC3quaISGZ9cROr6nGWHWubbxaqzgJmAYwePVozMzObWu24srKyaOm2vmRxNc+JxFVZU0dJRRUlFdWUHKhyn7uvj3i+91BNk/sID/GQGBVGYlQYJ/cMo0d0aMPrxKgwEqJC2b2/kl/8czXVdV4AQoKE//nuEBKiQimtqKb0YDWlFVXuc+eYG/dWsb+y6WOGBnlIiAolPjKUhKgwEiNDSYhynidEOsePd8sSo8IIDzl+i2tH/j36Q1eLqy2ajyYCF4nIBUA40B1nxqlYEQlW1VogGWe2JoBCnInIC0UkGIjBmcWpvrxe422M+YacbWUszi0hPSmK4vI6Ktd+3ejk7jTllFRUUXrQSQIHqmqb3E90WDCJ0WEkRoWS1iOK8anx3zjRNz7xR4QGndC39r5xEby9opDioiJuv3AcGSlxx92mutbb0AxVnzBKK6opqXCTyEEnoeTvqaCkoorKGm+T+4kMDXISRlQoCZFO4qhPIolu2dKiGr78eBNTTu55QrGZrqXVSUFV7wHuAXBrCr9Q1WtF5J/A5cBsYAbwrrvJXPf1Z+7yT1RVRWQu8A8R+QPQB0gHPm9tfKbzKKmoYmleCe98VcTCTXu+uXB5TsPT2IgQ90QeytA+3d2T+39OivVJ4ES/WTdXRkocGSlxZGWVnvBJNzTYQ6+YcHrFhJ/Q+oeqa92k0SiJHKx2nrtJpGjvYVYX7qXsYDW13iMr3Xk8tTCP747ozWUZyYwbGE9EaIcfH9O0AV/+FdyFM9/sw8BXwPNu+fPA392O5DKcK45Q1XUi8gawHqgFblNVmxS8C6uqrSNnWzmLcktYnLuHdcX7AacZp54AE/sEcfdlp9Mj2mlKCQnq/LffRIQGExEfTL/4iOOu6/Uq+ytrKD1YzXOLtzL78+0oTtvsh2u/5r3VOwkJEkb1j2NyeiIT0xIZkRxrHeJdVJsmBVXNArLc5/k4Vw8duU4lcMVRtn8E5wom0wWpKrm7K1i0eQ+Lc0tYvrWUyhovwR4hIyWOX5w7mMnpPaiu83Ld88upqfUSEuwhs18Iw/rG+Dv8gOXxCLERzpVPl2ck89aX26lTCAn28Lfrx1DrVZbklbAkt4Tf/3szv//3ZrqHB3P6oAQmpfdgUloiAxIi2rzD2wQmqy8avyqtqGJJXgmL3drArv1VAKT2iOSqMf2ZnJ7IuNQEosK++af66k3jyc4vZXxqAge2rvJH6B1SRkocvxwTTlVsCuNTExqatyan94Dznd/Hsi2lLMktYUleCR+v2wVA39huTEpLZJJbk4iPDPXn2zA+ZEnBtKujNQnFRoQwMS2RM9ITmZTeg76x3Y65n/p2e4CsrT4Pu1NJiwsiMzOtyWUJUWFcOLIPF47sg6qyrfSQW4vYwwdrd/L6lzsAGNqnO5PSE5mUlsiYAfE+6Zsx/mFJwfiUqpK3u6IhCSzPL+NwTR3BHmFUoyahYX1jrA07wIgIAxMjGZgYyXXjU6it87KmaF9DLeKFJVv566f5hAZ7GDsgnolpiUxOT2RI7+547HfZYVlSMG2ucZPQktwSvt5fCThNQleO6XfUJiET2IKDPJzWP47T+sfx47PTOVhVy+fbypwkkVvCbz/ayG8/griIECakObWISWmJJ9QZbgKH/VeaVquqrSOnoLyhX2BtUcuahEzHEhkWzJSTejLlpJ4A7D5QydK8EpbklrIkbw/vr94JQEpCBJPcWsTpqYnERIT4M2xzHJYUTLNZk5BpSs/ocC45LZlLTktGVdmyp4LFuSUN95a8unw7HoHhybFMSktgUloPRqXEEhYcRE5BOe9tqSZ6YLndUOdnlhTMceUUlPPJxl0U7ajm/T2rWGxNQuY4RIS0ntGk9YzmhokDqanzsnLH3ob+iL98ms9TC7fQLSSIk3pFs7ZoH3Ve5b1t2bx603hLDH5k/8HmmOasKOLn/1xJ/Q2xkaE7yTypJ5PTncsTk+OsvdgcX0iQhzED4hkzIJ6fnjOYA5U1ZOeXsSR3D/9aXdxwx3VljZcP1uy0pOBHlhRMk3J3HeBPC3Ib2oXBuXv4lsxB3H5Wuv8CM51CdHgI5wxJ4pwhSVx0al+ueTabqlpnPKcXlmylpKKKH5+VTlrPwBudtLOzpGC+IW/3AR5fkMd7q4uJCAnislF9eW/1TmrrvAQJnD4o0d8hmk4mIyWOf/xwPK/N/4ILJo5k+dYy/v5ZAXNXFXPhiD7ccXYaaT2j/R1ml2FJwQCQt7uCJxbk8i83Gdxy5iB+ODmVuMhQrhmXQnZ+KWF7C6xab3wiIyWOA4NCyTw5ibNOTmLm5FSeXbyVlz/bxr9WF/Pd4b254+x0BidZcvA1SwpdXP4eJxnMXVVMeEgQPzpjEDPPSP3GMAb/GfWz0I+Rmq4kISqMu88/mZlnpPLc4nxeWraN99fs5ILhvbnjrHRO6mXJwVcsKXRRW0sO8ucFubyzsoiw4CB+ODmVmWekkhAV5u/QjGkQHxnKL887mR9OTuW5Jfm8tKyA91fv5ILhvbjj7HRO7tXd3yF2OpYUupiC0oM8sSCPd1YWERIk3DhpID86cxCJlgxMAIuLDOW/v+Mkh+eXbOVvS7fxwZqvOX+YkxxO6W3Joa1YUugitpce4s+f5PL2V0UEe4QbJgzgR2cOoke0JQPTccRGhPLzc0/ixkkDecFNDh+u/ZrvDE3ijrPTGdrHhlBvLUsKndyOskM8+Ukeb60oJMgjzDh9ADdnptIz+sRm+DImEMVGhPKzc0/ixkmpvLB0Ky8s3crH63Zx7hAnOdj8Gi3X6qQgIuHAIiDM3d+bqvqAiAzEmYozHlgBXKeq1SISBrwMZAClwJWqus3d1z3AjUAdcIeqftza+LqqHWWHeGphHm/mFOLxCN8fn8KtmYPo2d2Sgek8YiJC+Ok5g/mvSQP529KtPL9kK/9ev4uppyTxk6mWHFqiLWoKVcBZqlohIiHAEhH5EPgZ8EdVnS0if8E52T/j/ixX1TQRuQr4LXCliAzBmZpzKM4czfNFZLBNydk8RXsP8+QnebyZswNBuHZcf27JTDvhuX+N6YhiuoXwk6mDuWHiQF5cuo3nl+Qz7c+7mHpKT+48ezDDky05nKhWJwVVVaDCfRniPhQ4C7jGLX8JeBAnKUx3nwO8CTwpzjx/04HZqloFbHXncB4LfNbaGLuC4r2HeWphHm986SSDq8b059Ypg+gdYyOTmq4jplsId05N54ZJA3hp6TaeW7KVC59cwlkn9+TOs9MZ2S/W3yEGPHHO6a3ciUgQkAOkAU8BvwOyVTXNXd4P+FBVh4nIWuA8VS10l20BxuEkimxVfcUtf97d5s0mjjcTmAmQlJSUMXv27BbFXVFRQVRU4N1G35y4yiq9vJdfw6IdtShwRnIw01JDSOjW9pPXd4bPqz1ZXM3ji7gO1yrzCmr4eFsNB2tgRI8gLh4UQmrsic8U11k/rylTpuSo6ugjy9uko9lt4jlVRGKBOcApTa3m/mxqLGU9RnlTx5sFzAIYPXq0ZmZmNjdkALKysmjptr50InHt2l/J0wvzeO3zHXhVuWJMf26bMsinA9R15M/LHyyu5vFVXOcDFVW1vLRsG88tzufX2ZWcObgHd05NZ1T/49+h39U+rza9+khV94pIFjAeiBWRYFWtBZKBYne1QqAfUCgiwUAMUNaovF7jbYxr9/5Kns7awj8+347Xq1yekcxtU9JsditjjiEqLJjbpqQxY8IAXv5sG88uyufSp5dxxuAe3Hm2M8Bjdn4p41MTuvxQLm1x9VEPoMZNCN2AqTidxwuBy3GuQJoBvOtuMtd9/Zm7/BNVVRGZC/xDRP6A09GcDnze2vg6iwUbdvF01hZWF+7Fq3DZqL78+Kx0SwbGNENUWDC3ZqYx4/QB/D27gFmL8rnsmWXUzwUVGuzp8vM5tEVNoTfwktuv4AHeUNX3RGQ9MFtEHga+Ap53138e+LvbkVyGc8URqrpORN4A1gO1wG125ZHjH8sLuHfOWgCCRHji6lP57og+fo7KmI4rMiyYm88cxHXjU7j11RV8unkPANW1XrLzSy0ptIaqrgZOa6I8H+fqoSPLK4ErjrKvR4BHWhtTZ/LuyiLuf3ddoxJlW+khv8VjTGcSGRbMHWenk51fSlWtF6/C4Zqu/V3U7mgOUKrKk5/k8X/zNnNKr2jySw5SW+clJNjD+NQEf4dnTKdRP5/Dp5t2s3DTbp78JA+vV/nFuSfh6YJzjFtSCEC1XuUX/1zNWysKueS0vjx22XDWFu23jjBjfKR+ePgfn53Or95dx9NZW8jfc5A/XDnS36G1O0sKAWbvoWp+90Ulm8oL+enUwdxxdhoi0vBHa4zxnZAgD7+5ZBiDekTyyAcbuPKvh/mvdK+/w2pXlhQCyLaSg9zw4hfs2Ovl8atOZfqpff0dkjFdjohw0+RUBiREcsfsr/h1iZf04fu6zDhKbX/bq2mRL7aVccnTS9l7qJpfjg23hGCMn00dksSbN09ABK74y2d8vO5rf4fULiwpBIB3viri2meXExcRypxbJzI47sRvwTfG+M6QPt351enhDO4Vzc2v5PCXT7fQFkMDBTJLCn6kqvxp/mZ+8vpKRqXE8vatExiQGOnvsIwxjcSGeXh95nguGN6bxz7cyC/fXE11beftZ7A+BT+pqq3j7rfWMOerIi4blcyjlw4nNNhytDGBKDwkiD9fdRqDekTxxIJctpcd4i/fzyAuMtTfobU5Owv5QfnBaq577nPmfFXEL84dzO+vGGEJwZgA5/EIPztnMH+68lS+2r6Xi59eypY9FcffsIOxM1E7y99TwSVPL2Vl4V6euPo0bj8rHWc6CWNMR3DxaX15beY4KiprueSppSzNK/F3SG3KkkI7Wp5fyqXPLGN/ZS2v/XAcF4208YuM6YgyUuJ557aJ9IoJ5wcvfM4/lm/3d0htxpJCO3l7RSHff345CZGhvHPrRDJS4v0dkjGmFfrFR/DWLROYlJbIvXPW8NB766nzdvwrkywp+Jiq8od/b+Jnb6xizIB43r5lIv0TbLhrYzqD6PAQnp8xmusnDOD5JVuZ+fKXVFTV+jusVrGk4EOVNXXcOXslT3ySx/dGJ/PiDWOJiQjxd1jGmDYUHOThwYuG8tD0oWRt3sPlzyyjaO9hf4fVYpYUfKS0oorvP7ecuauK+eV5J/Hby+wKI2M6s+tOH8AL14+hqPww059cylfby/0dUovYWcoHtuyp4JKnl7GmaB9PXTOKWzPT7AojY7qAMwf34O1bJ9At1MNVs7L516qON6Nwq5OCiPQTkYUiskFE1onInW55vIjME5Fc92ecWy4i8oSI5InIahEZ1WhfM9z1c0VkRmtj84fPtpRy6dPLOFRdy2szx/PdEb39HZIxph2lJ0Xzzq0TGZEcw49f+4rH5+d2qKEx2uKO5lrg56q6QkSigRwRmQdcDyxQ1cdE5G7gbuAu4Hyc+ZfTgXHAM8A4EYkHHgBGA+ruZ66qdog6WE5BOc8tzuff678mNTGKF64fY/MnG9NFJUSF8cpN47jn7TX8cf5mviwoY3RKHJPSewT8EPhtMR3nTmCn+/yAiGwA+gLTgUx3tZeALJykMB14WZ3UmS0isSLS2113nqqWAbiJ5TzgtdbG6Gs5BeVc+dfPqPUqHoH7pw2xhGBMFxcWHMT/XTGS8JAg/rF8O4tzS3jm0y28etP4gE4M0pbVGhEZACwChgHbVTW20bJyVY0TkfeAx1R1iVu+ACdZZALhqvqwW34/cFhVf9/EcWYCMwGSkpIyZs+e3aJ4KyoqiIqKatG2jb24toqsQucyNA9waXoI0wa1fEyUtoqrrVlczWNxNU9njeu9LdW8mVsDgACXtfL80FZxTZkyJUdVRx9Z3mYD4olIFPAW8BNV3X+MjtWmFugxyr9dqDoLmAUwevRozczMbHa8AFlZWbR023o1dV5+tfxToJYggZBgD1dPHdOqbwJtEZcvWFzNY3E1T2eNK3pgOe9ty6ayxhlZ9ZIzRzExLdHvcR1Nm1x9JCIhOAnhVVV92y3e5TYL4f7c7ZYXAv0abZ4MFB+jPKDNWpTP9vJD3HP+yfzs3JMCvmpojGlfGSlxvHrTeK4Z2x8FPtm4+7jb+FOrawriVAmeBzao6h8aLZoLzAAec3++26j8dhGZjdPRvE9Vd4rIx8Bv6q9SAs4F7mltfL5UUHqQJxbkct7QXvzozEH+DscYE6Dq51gXgb8t3cpFI/swsl/s8Tf0g7aoKUwErgPOEpGV7uMCnGRwjojkAue4rwE+APKBPOBZ4FYAt4P5IeAL9/Hr+k7nQKSq3DdnLSHu3YzGGHM8d51/Mj2iw7jrrdXU1AXmRD1tcfXREpruDwA4u4n1FbjtKPt6AXihtTG1h3dWFrEkr4RfTx9Kr5hwf4djjOkAuoeH8ND0Ycz8ew6zFuVz25Q0f4f0LXZHcwuUH6zmofc2cGq/WK4dl+LvcIwxHci5Q3txwfBePL4gl60lB/0dzrdYUmiBRz/cwP7DNTx66XCCPDZ8hTGmeR68cChhwR7ufms13gAbbtuSQjNl55fyxpeF3DQ5lVN6d/d3OMaYDqhn93Duu+AUlm8t440vd/g7nG+wpNAMVbV13DtnDf3iu3Hn2en+DscY04FdOaYf41Pj+c0HG9i9v9Lf4TSwpNAMTy/cQv6egzx88XC6hQb5OxxjTAcmIjx66Qgqa708+K91/g6ngSWFE5S3u4JnsrZw0cg+nDm4h7/DMcZ0AgMTI7nz7HQ+WPM1H6/72t/hAJYUTojXq9w7Zw3hIR7unzbE3+EYYzqRmWekcnKvaH717lr2V9b4OxxLCifizZxCPt9axr0XnEKP6DB/h2OM6URCgjz89rIR7DlQxf9+tNHf4VhSOJ6Siioe+WADYwfE873R/Y6/gTHGNNPIfrHcMHEgr2Rv54tt/h3IwZLCcTz83noOVdfym0uH4bF7EowxPvLzcweTHNeNu99aTVVtnd/isKRwDItz9/DOymJuyUwjrWe0v8MxxnRiEaHBPHLJcLbsOchTC7f4LQ5LCkdxuLqO++asJTUxklszbQRUY4zvnTm4B5ec1pdnsvLY9PUBv8RgSeEonvgkl+1lh3j4kmGEh9g9CcaY9nH/tCFEh4dw99urqfPDEBiWFJqw8ev9PLson8szkpkwqPUzJBljzImKjwzlV9OG8NX2vbySXdDux7ekcASvV7n37TV07xbCfRec4u9wjDFd0PRTnZtk//ejjRTtPdyux26r6ThfEJHdIrK2UVm8iMwTkVz3Z5xbLiLyhIjkichqERnVaJsZ7vq5IjKjLWJrrlc/386K7Xv5n++eQlxk6yfXNsaY5hIRHr54GF6FH/9jBU8tzCWnoLxdjt1WNYUXgfOOKLsbWKCq6cAC9zXA+UC6+5gJPANOEgEewJmicyzwQKOpOdvF/PW7eOhf6xnetzuXnNa3PQ9tjDHf0C8+gqvG9mPF9r38/uPNXPtcdrskhjZJCqq6CDjyjovpwEvu85eAixuVv6yObCBWRHoD3wHmqWqZqpYD8/h2ovGZnIJyZv79S6rrvGzeVcGK7Xvb69DGGNOkBLe1QoGaWi/Z+aU+P2arp+M8hiRV3QmgqjtFpKdb3hdoPIB4oVt2tPJvEZGZOLUMkpKSyMrKalGAFRUVDdv+a0s19R39NbVeXpv/BQcG+af5qHFcgcTiah6Lq3ksrm8L31eHR8Cr4BEI21tAVlahT+PyZVI4mqZuC9ZjlH+7UHUWMAtg9OjRmpmZ2aJAsrKyqN92f2wRb+WuRIDQEA9XTx1DRkq7tl41GVcgsbiax+JqHovr2zKBuvgtPPrhRh64aBjfH/+f6X99FZcvrz7a5TYL4f7c7ZYXAo0HEUoGio9R3i7qs8+MCQN49abxfksIxhjT2CWjnAaTypr2GfrCl0lhLlB/BdEM4N1G5T9wr0IaD+xzm5k+Bs4VkTi3g/lct6xdrC3aR2iwh/u+e4olBGNMwOgZHU7vmHDWFO1rl+O1SfORiLyGU9NJFJFCnKuIHgPeEJEbge3AFe7qHwAXAHnAIeAGAFUtE5GHgC/c9X6tqu02XODaov2c0iuakCC7dcMYE1iG941hdWEHSgqqevVRFp3dxLoK3HaU/bwAvNAWMTWHqrKueB/TRvZp70MbY8xxjewXy7/X72Lf4RpiuoX49Fj2tRgoLD/M/spahvWJ8XcoxhjzLSOSnXPTmnaoLVhSwOlPABjWt7ufIzHGmG8b0TcWgFWFvr9/ypICsLZ4H8EeYXCSzYW1IVkAACAASURBVJlgjAk8MREhpCREWE2hvawt2k96UrQNkW2MCVgjkmNZbTUF31NV1hbtY1gfazoyxgSukckxFO+rZM+BKp8ep8snhV37qyg9WM2wvtbJbIwJXCOSnX4FX9cWunxSsE5mY0xHMLRPdzyCz+9XsKRQvA8ROKW3JQVjTOCKDAsmrWeU1RR8bW3RflITI4kI9cfYgMYYc+KczuZ9OPcA+0aXTwrri/dZf4IxpkMYmRxD6cFqn07R2aWTwv5qpXhfpd3JbIzpEIa7nc2+vF+hSyeF7fudoWiHWiezMaYDOKV3NCFBwipLCr6xbb8XgKFWUzDGdABhwUGc3Ku7Tzubu3RSKNjvpX98hM9HHTTGmLYyIjmGNYX78Pqos7nLJwW7P8EY05GMSI7hQFUtuw9ZUmhT+w7XsPuQWtORMaZDqb+zOX+f1yf7D7ikICLnicgmEckTkbt9dZz1xfsB7HJUY0yHkt4zivAQD1v3+WbO5oBKCiISBDwFnA8MAa4WkSG+ONa6Yqf3fqgNhGeM6UCCgzwMSIgkZ1cdOQXlbb7/gEoKwFggT1XzVbUamA1M98WBFuXuITwICkoP+WL3xhjjEzkF5eTurqCsUrn2uew2Twziy9ulm0tELgfOU9Wb3NfXAeNU9fYj1psJzARISkrKmD17drOOk1dexyPLK1GUUI/wyzHhpMUFzlwKFRUVREVF+TuMb7G4msfiah6L68S8t6WaN3NrAOdb/aXpIUwbFNrs/UyZMiVHVUcfWR5oA/5IE2XfylqqOguYBTB69GjNzMxs1kHWLcxD2QQIdQpVsSlkZqa1IFzfyMrKornvqT1YXM1jcTWPxXViogeW81buMhQIDfFw9dQxZKTEtdn+A635qBDo1+h1MlDc1gcZn5pAsMfJP8FBHsanJrT1IYwxxieG9e2OxwOD4zy8etP4Nk0IEHhJ4QsgXUQGikgocBUwt60PkpESx4MXDQXgZ+cMbvMP1RhjfCV3VwV1XpjaP8Qn566ASgqqWgvcDnwMbADeUNV1vjjWpaP64hGoqKr1xe6NMcYn6icGS+num9N3oPUpoKofAB/4+jgRocH0jfKwcofvJ8I2xpi2srZ4H9FhwfSIaKoLtvUCqqbQ3lJjPKzasdenE1YYY0xbWlu0nyF9uuMRSwptLjXGw/7KWrbZvQrGmA6gts7Lhp37fToSQ9dOCrHOvQmrrAnJGNMBbNlzkKpa3w7k2aWTQp9IoVtIEKt8PBG2Mca0hfpOZl/OFtmlk0KQRxjeN8ZqCsaYDmFt8T7CQzyk9vDdHdZdOikAjOwXw9ri/dTU+WYYWmOMaSvrivczpHd3gjy+6WQGSwqM7BdLda2XTV8f8HcoxhhzVF6vsr7Yt53MYEmBke6EFXa/gjEmkBWUHaKiqtan/QlgSYHkuG4kRIZav4IxJqDVdzIP9fEUwl0+KYgII/vF2hVIxpiAtrZ4H6FBHtJ7Rvv0OF0+KYDThJS7u8LGQTLGBKx1RfsZ3CuK0GDfnrYtKeBcgaQKawr3+TsUY4z5FlVlbfE+n/cngCUFAEa4nc3WhGSMCURFew+z91ANQ3185RFYUgAgPjKU/vER1tlsjAlIa4v2AzCsj287mcGSQoOR/WItKRhjAtL64n0EeYRTeltSaDcjk2Mo3lfJ7gOV/g7FGGO+YW3xftJ6RBEeEuTzY7UqKYjIFSKyTkS8IjL6iGX3iEieiGwSke80Kj/PLcsTkbsblQ8UkeUikisir7vTcbabU/s5/Qqrd1hnszEmsKwt2ufz+xPqtbamsBa4FFjUuFBEhuDMrzwUOA94WkSCRCQIeAo4HxgCXO2uC/Bb4I+qmg6UAze2MrZmGdonhiCPWGezMSagzF+/i90HqojpFtIux2tVUlDVDaq6qYlF04HZqlqlqluBPGCs+8hT1XxVrQZmA9NFRICzgDfd7V8CLm5NbM3VLTSIfnHd+NeqYnIKytvz0MYY06ScgnJufTUHgFezt7fLuclXczT3BbIbvS50ywB2HFE+DkgA9qpqbRPrf4uIzARmAiQlJZGVldWiICsqKhq2zSuvY3tZJV6Fq/66jLvGhJMW5/v2u+PFFUgsruaxuJrH4vq2f22pprrOmS64ts7La/O/4MCgUJ/GddykICLzgV5NLLpPVd892mZNlClN10z0GOs3SVVnAbMARo8erZmZmUdb9ZiysrKo33bdwjzqKz21XqiKTSEzM61F+22txnEFEoureSyu5rG4vm295kLuZgQIDfFw9dQxZKTE+TSu4yYFVZ3agv0WAv0avU4Git3nTZWXALEiEuzWFhqv3y7GpyYQFuKhssaLAgMSItrz8MYY8w37K2t4cVkBAxIiuCwjmQmDEhsSgi/56pLUucBVIhImIgOBdOBz4Asg3b3SKBSnM3quqiqwELjc3X4GcLRaiE9kpMTx6k3jueXMQUSEBPH37AKcsIwxpv397qNNlFRU8fhVp/Hjs9LbJSFA6y9JvURECoHTgfdF5GMAVV0HvAGsBz4CblPVOrcWcDvwMbABeMNdF+Au4GcikofTx/B8a2JriYyUOO46/2Tuv3AI2fllzP5ix/E3MsaYNpZTUM4rywu4fsJARrqXy7eXVnU0q+ocYM5Rlj0CPNJE+QfAB02U5+NcneR3V43px7sri/jN+xuYclJPesWE+zskY0wXUV3r5d6319C7ezg/P3dwux/f7mhugojw2KUjqK7zcv+7a60ZyRjTbp5dnM+mXQf49fRhRIb56gLRo7OkcBQDEiP52TmDmbd+Fx+s+drf4RhjuoCtJQd5fEEuFwzvxdQhSX6JwZLCMdw4aSDD+8bwwNy17D1U7e9wjDGdmKpy35w1hAV7ePDCoX6Lw5LCMQQHeXjssuGUH6rh4fc3+DscY0wn9taKIpZtKeXu80+mZ3f/9WNaUjiOoX1iuPnMVN7MKWTR5j3+DscY0wmVVlTx8PvrGZ0Sx9Vj+vs1FksKJ+DHZ6WT2iOSe+es4aDN42yMaWOPvL+Bg1W1/ObS4Xg8TQ3w0H4sKZyA8JAgfnvZCArLD/N//97s73CMMZ3I4tw9vP1VETefOYjBSdH+DseSwokaMyCe68an8LdlW1mx3UZRNca03uHqOu6bs5aBiZHcNsU/Y60dyZJCM/zyvJPo1T2cu99aTXWt19/hGGM6uCc+yWV72SEeuWRYu8yqdiIsKTRDdHgIj1wyjM27Kng6K8/f4RhjOrANO/cza1E+V7iD3QUKSwrNdNbJSUw/tQ9PLcxj864D/g7HGNMB1XmVe95eQ2y3EO694BR/h/MNlhRa4FfThhAVFswv31xNndeGwDDGNM8r2QWs3LGX+6cNIS6yXaejPy5LCi2QEBXGgxcNZeWOvby0bJu/wzHGdCA79x3mdx9vYnJ6ItNP7ePvcL7FkkILXTSyD1NO6sHvPt7EjrJD/g7HGNNBPPDuOmq9Xh65eDjO9PSBxZJCC4kID18yHI/AvXPW2Eiqxpjj+mjt1/x7/S5+MnUw/QN0dsfWTrLzOxHZKCKrRWSOiMQ2WnaPiOSJyCYR+U6j8vPcsjwRubtR+UARWS4iuSLyujszW0DrG9uNu88/mcW5Jby1osjf4RhjAtiByhoemLuWU3p358ZJA/0dzlG1tqYwDximqiOAzcA9ACIyBGeqzaHAecDTIhIkIkHAU8D5wBDganddgN8Cf1TVdKAcuLGVsbWLa8elMGZAHA+9t549B6r8HY4xJgDlFJTzgxc+Z9f+Kh69dDghQYHbSNOqyFT13+4UmwDZQLL7fDowW1WrVHUrkIczq9pYIE9V81W1GpgNTBenYe0s4E13+5eAi1sTW3vxeIRHLx3B4eo67nztK55amEdOgd3xbIxx5BSUc/WsbL7avpcgjwT8FYvSVm3hIvIv4HVVfUVEngSyVfUVd9nzwIfuquep6k1u+XXAOOBBd/00t7wf8KGqDjvKsWYCMwGSkpIyZs+e3aKYKyoqiIqKatG2R3pxbSVZhXUAhHrgl2PCSYtr2R2KbRlXW7K4msfiap7OGtdfVlWSvdM5N3iAS9NDmDao9a3jrY1rypQpOao6+sjy4871JiLzgV5NLLpPVd9117kPqAVerd+sifWVpmsmeoz1m6Sqs4BZAKNHj9bMzMyjrXpMWVlZtHTbI63x5pJV6AyWV+OFwzH9ycxM93tcbcniah6Lq3k6W1yqyl8X5ZO9cyMizskvJNjD1VPHkJES57e4jue4SUFVpx5ruYjMAKYBZ+t/qh2FQL9GqyUDxe7zpspLgFgRCXaboxqv3yFMGJTIUyF5VNV4UeCTjbuZMWEA3cND/B2aMaadVdd6uf+dtbz+5Q6mjejNteP6s2L7XsanJrRJQvClVs0KLSLnAXcBZ6pq44v15wL/EJE/AH2AdOBznBpBuogMBIpwOqOvUVUVkYXA5Tj9DDOAd1sTW3vLSInj1ZvGk51fyoHDNTy3ZCuXPb2MF64fQ7/4wLz0zBjT9vYdquGWV3NYtqWUH5+Vxk+nDsbjEU4PoPGNjqVVSQF4EggD5rk3YWSr6s2quk5E3gDW4zQr3aaqdQAicjvwMRAEvKCq69x93QXMFpGHga+A51sZW7vLSIlr+BZwxuAe3PLqCqY/tZS/XpfBmAHxfo7OGONrBaUHueHFL9hRdoj/u2Ikl2UkH3+jANOqpFDfMXyUZY8AjzRR/gHwQRPl+ThXJ3UKE9ISeee2idz44hdc82w2j146gss74B+IMebEfLmtjJl/z8Gryis3jmNcaoK/Q2qRwL1YthMYmBjJnFsnMnZgPL/45yoe+3Aj3gC/HM0Y03zvrizimmeXE9MthDm3TuywCQEsKfhcTEQIL94wlmvH9ecvn27h5ldybJ5nYzoJVeVP8zdz5+yVnNY/ljm3TmBgYqS/w2oVSwrtICTIw8MXD+PBC4cwf8MuLv/LZxTvPezvsIwxrVBZU8dPX1/Jn+bnctmoZP5+4zhiIwJ+dJ7jsqTQTkSE6ycO5IXrx1BYdojpTy1l5Y69/g7LGNMCpRVVfP+55byzspj//s5J/P6KEYQGd47Taed4Fx1I5kk9efvWCXQLCeLKv37G3FUd6nYMY7q8vN0VXPL0MlYX7ePJa07jtilpATkEdktZUvCD9KRo3rltIiOTY7njta/4w7zNNvS2MR3AsrwSLn16KYeqa5k9czzTRgTeJDmtZUnBT+IjQ/n7TWO5PCOZJxbkcvtrX1FZU+fvsIwxR/HGFzv4wQufk9Q9nDm3TmRU/8C+M7mlWnvzmmmFsOAgfnf5CNJ7RvHYRxspLDvEsz/41vhUxhg/8qry2Icb+cunW5icnshT147q1MPXWFLwMxHhR2cOYmBiJD95fSUXPbmUW5ocG9YY094OV9fx9Moqvty1hWvG9ef/XTQ0oOdCaAud+911IOcO7cWbN0/AI/DI8ko+Wvu1v0MypkvbfaCSq2Z9Rs6uOv7nu6fwyMXDOn1CAEsKAWVIn+68c/tEkqM83PxKDk9n5VkHtDF+sPHr/Vzy1DI276rgx6eFcdPk1E51hdGxWFIIMD2jw7l7bDgXjezD/360iZ//cxVVtdYBbUx7WbhpN5c/8xm1Xi//vPl0RiV1rVZ2SwoBKDRIePyqU/nZOYN5e0UR1z67nNIKm//ZGF97+bNt3PjiF/SPj+Cd2yYyrG+Mv0Nqd10rBXYgIsIdZ6eT2iOSn7+xiulPLeXn55xE8b7DHWKiDmM6ki+2lfHoBxtYsX0vU0/pyeNXnUZkWNc8PXbNd92BTBvRh35xEVz/t8/56RsrESAsxMOrN423xGBMK6kqsxbl89hHG1GFII9zNWBXTQhgzUcdwsh+sVw5xpnFVIHKGi8vL9tmw3Ab00KqyoINu7joyaU8+qGTENwFfL61zK+x+VurkoKIPCQiq0VkpYj8W0T6uOUiIk+ISJ67fFSjbWaISK77mNGoPENE1rjbPCFdpav/BJ0zpBfhIR5EnDlN311VzPmPL+a91cXUWXIw5oSoKp9s3MX0p5Zy40tfsvdwNbdlDiI8xEOQQEiwh/EdeC6EttDaOtLvVPV+ABG5A/gVcDNwPs68zOnAOOAZYJyIxAMPAKNxvvTmiMhcVS1315kJZOPMzHYe8GEr4+s0Gs8BPXZgPMV7D/PnT/K4/R9fkdYzlx+flca0EX0I8lguNeZIqkrWpj38af5mVhXuIzmuG7+9bDiXjkomJMjDWackkZ1fav11tH46zv2NXkbinOgBpgMvq3ORfbaIxIpIbyATmKeqZQAiMg84T0SygO6q+plb/jJwMZYUvqHxHNDg9Dd8uHYnf16Qx52zV/L4/FxuPyuNi0b2IbgL3GRjzPGoKp9u3sOf5ueycsde+sZ247FLnWTQeKjrI/+3ujJp7c1RIvII8ANgHzBFVfeIyHvAY6q6xF1nAXAXTlIIV9WH3fL7gcNAlrv+VLd8MnCXqk47yjFn4tQqSEpKypg9e3aLYq+oqCAqKqpF2/pSc+PyqpKzq465W2rYccBLUoQwLTWE0/sEE9yGNYfO8nm1F4uredoyLlVlbUkd7+TVsGWfl4Rw4cJBIUzq2/z/ic76eU2ZMiVHVb812NpxawoiMh/o1cSi+1T1XVW9D7hPRO4BbsdpHmrqU9cWlDdJVWcBswBGjx6tmZmZx3wPR5OVlUVLt/WllsR1FvBzrzJvwy6eWJDL82v3M684mNumpDZUkf0RV3uwuJqnM8elqizOLeFP8zezYrtTM/jNJWlcnpHc4klwOvPn1ZTjJoX6b+8n4B/A+zhJoRDo12hZMlDslmceUZ7llic3sb5pBo9H+M7QXpw7JIlPNu7m8QW53PXWGp5YkMdtU1r3j2FMIFNVluSV8Kf5ueQUlNMnJpxHLhnGFRn97G++mVrVpyAi6aqa6768CNjoPp8L3C4is3E6mvep6k4R+Rj4jYjUN96dC9yjqmUickBExgPLcZqj/tya2LoyEeHsU5I46+SeZG3ew+Pzc7l3zhqe/CSXWzIH8b0x/QgLDvJ3mMa0mqqyNK+UP83fzJcF5fSOCefhi4dxxehk+xtvodZeffSYiJwEeIECnCuPwLl66AIgDzgE3ADgnvwfAr5w1/t1faczcAvwItANp4PZOplbSUSYclJPMgf3YHFuCY8vyOX+d9fx1MIt3HxmKleN7U94iP3jmI5HVVm2xUkGX2xzksFDFw/je5YMWq21Vx9ddpRyBW47yrIXgBeaKP8SsJkEfEBEOGNwDyanJ/LZllL+tCCXB/+1nqeztvCjMwdxzdj+dAu1fyQT+FTV+Ruen8vn28ro1T2ch6YPtdpvG+q693J3QSLChLREJqQlkp1fyuPzc3novfU8k7WFH52RyrXj+xMRan8SJjAt2+L0GXy+tYyk7mH8evpQvje6n9V225idAbqo8akJjJ+ZwOdby3hiQS6PfLCBZz7dwg8np/KD01O69NgvJrB85jYTLXeTwf+7aChXjrFk4Cv2n9/FjR0Yzys3jSOnoIwnFuTx2482MmvRFm5yk0N0J56L1gS27HwnGWTnl9EzOowHLxxi/WDtwJKCASAjJZ6X/mssK3fs5YkFufzu403MWpTPjZMGMmPCAGK6WXIwvpVXXse6hXl0Dw/mgzVf81l+KT2iw3jgwiFcbcmg3VhSMN9war9YXrh+DGsK9/H4glz+MG8zzy7O578mDiR4Xy3rFubZ+DCmzX2+tZTHPq+kVjcBEBsRwq+mDeGacZYM2pslBdOk4ckxPDdjNGuL9vHkJ3k8vqD+dpRNhAQJf776NM4b1tuvMZqOS1UpKD3E4tw9LM4tIWvTHmrdMQwEuGHCAP5r0kC/xthVWVIwxzSsbwx/uS6DB+eu48Vl2wCoqVNufmUFAxIimJiWyKS0RE4flEBsRKh/gzUBbd+hGpZuKWFxbgmLc/dQWH4YgL6x3ThzcCILN+5GcYavnpTew7/BdmGWFMwJuXBkH/6RvY06heAgD9eO68/2skO8u7KYV5dvRwSG941pSBIZKXFW7e/iqmu9fLW9nCV5JSzKLWFN4V68ClFhwZw+KIGZZ6QyOb0HAxIiEBGem7OAqtgUa570M0sK5oRkpMTxyzHh3/qnranzsrpwL0tyS1maV8Kzi/J5JmsLYcEeRg+Ia0gSQ/vE2FwPnZyqsmXPQZa4TULZ+aUcrK7DI05f1e1npXNGeiIj+8U2OUBjWlwQmZlpfojcNGZJwZywpv5pQ4I8ZKTEk5ESz51T0zlYVcvnW8tYklfC0rwS/vejTfwvm4jpFsKEQQkNSSLF/XZoOrayg9UszXOag5bkllC8rxKAlIQILj6tL5PTe3D6oAS7eq0DsaRg2lRkWDBTTu7JlJN7ArD7QCWfbSllSa6TJD5c+zXgtCNPSktkYnoiEwYlkBgV5s+wzQmqqq0jp6CcJblO38Da4n2oQvfwYCYMSuS2sxKZnNaD/gkR/g7VtJAlBeNTPaPDmX5qX6af2hdVZWvJQZbmlbAkr4QP1+7k9S93AHBK7+5MSnNqEmMHxttwGwFCVcndXdHQObw8v4zDNXUEe4TT+sfy06mDmZSeyIi+MTbbXydh/3mm3YgIqT2iSO0RxXWnD6DOq6wt2tfQ1PTSsgKeXbyVkCBhVP+4hpqEnXDaV0lFFUvzSli0uYQleXvYtb8KgNTESL43OplJ6T0Ynxpvd7t3UpYUjN8EeYSR/WIZ2S+W26akcbi6ji8L/tMf8Yf5m/m/eZuJDgtm/KAEJg5KYFJ6IoN6RFl/RCvlFJQ3TFQ/tE93vtxW3nDPwPqdztTrsREhTExLZHJaIpPSE0mOsyahrsCSggkY3UKDmJzeg8nuNeplB6ud/gg3ScxbvwuApO5hTExLJDm2G1V1XuIP135jOj/TNK9XKT9UzbKiGl6cl01NnRcRJznX1CkhQUJGShz//Z2TmJxuV4x1VZYUTMCKjwzluyN6890Rzp3T20sPsXSL0x8xb93XHKiqa1j3z6s+omf3cBIiQ4mPDCUhKqzheXxkKAmRYSREhZIQGUpcZGibzFkdCKpq6yg7WE3JgWpKKqrcRzWljZ7X/yw7WIX3iJnPVeG05Bhun5LO2IHxNjquaZukICK/AH4H9FDVEnHq9o/jzL52CLheVVe4684A/sfd9GFVfcktz+A/M699ANzpTtZjDAD9EyLon9Cfq8f258lPnHGZ6k9yaT2j6BPbjdKKavL3HOTLbeWUH6r+1kmwXvfwYBKjwv6TNKK+mTwaJ5P4yNB2m+dXVTlYXUfJgW+e1Evrfx6s+kYC2F9Z2+R+uoUEkRgdSmJUGMlxEZzWP5aEyDASo0JZtSGX97fWUef1EhLs4d4LhtjNYqZBq5OCiPQDzgG2Nyo+H0h3H+OAZ4BxIhIPPACMBhTIEZG5qlrurjMTyMZJCudhU3Kaozh9UCKhC/OoqfUSJHD/tKHfOrHVeZW9h6opO1hN6cH//CytqPpPWUU1BaWHWLG9nLKDR08i0eHBjWoezsm1cUKpTx71CSUsOKhh1M9xA+MZmBhJ6cFqSg5UsafxSb7im9/wSyqqqKr1NhlDbEQICZHOif6UPt1JdJ8nRDnxJEaHkRgZRmJ06DGv3sqqKeD754xs6FOwhGAaa4uawh+BXwLvNiqbDrzsftPPFpFYEekNZALz6udlFpF5wHkikgV0V9XP3PKXgYuxpGCOIiMljldvGk92filhewuaPLEFecRpRooKI/0E9un1KvsO1/wngVRUNTxvnFAKyw+xqnAv5QerqT1KFukWEsThmjpg01GPF+SRhpN8QlQog3pEkRgd9o2yxKgwekSHERfRtrWVjJQ4SwamSa1KCiJyEVCkqquOuBqkL7Cj0etCt+xY5YVNlB/tuDNxahUkJSWRlZXVovgrKipavK0vWVwnbqhARchhn8TVDUgGkoOA7u6jQTCqQRyqhQPVyv5q5YD72F+trN5TR95exRnzE0b2CGJCn2C6hwrdw4SYUCEiBDwiQB1w2H3g1KEPOI8SnEdbCsTfI1hczeWruI6bFERkPtCriUX3AfcC5za1WRNl2oLyJqnqLGAWwOjRozUzM/Noqx5TVlYWLd3Wlyyu5gnEuHIKyrn6r8uoU2fUz19dPi5gvpkH4ucFFldz+Squ4yYFVZ3aVLmIDAcGAvW1hGRghYiMxfmm36/R6slAsVueeUR5llue3MT6xnRIRxtA0JhA1+JGSlVdo6o9VXWAqg7AObGPUtWvgbnAD8QxHtinqjuBj4FzRSROROJwahkfu8sOiMh498qlH/DNPgpjOpy0uCBum5JmCcF0KL66KPkDnMtR83AuSb0BQFXLROQh4At3vV/XdzoDt/CfS1I/xDqZjTGm3bVZUnBrC/XPFbjtKOu9ALzQRPmXwLC2iscYY0zzdY7bOo0xxrQJSwrGGGMaWFIwxhjTwJKCMcaYBtLRx5wTkT1AQQs3T6TtbxhtCxZX81hczWNxNU9njStFVXscWdjhk0JriMiXqjra33EcyeJqHoureSyu5ulqcVnzkTHGmAaWFIwxxjTo6klhlr8DOAqLq3ksruaxuJqnS8XVpfsUjDHGfFNXrykYY4xpxJKCMcaYBp0+KYjIL0RERSTRfS0i8oSI5InIahEZ1WjdGSKS6z5mNCrPEJE17jZPyBHTzDUznofc464UkX+LSJ8Aiet3IrLRPfYcEYlttOwe9xibROQ7jcrPc8vyROTuRuUDRWS5G+/rIhLairiuEJH/3975xEZVRXH4O6EVjFhpFbWhJNjEGCEaNQ3RaEwjRrEQqjsSV+JK4r+4MOgkGBcuwIUsNOnCmEgCIhSNEWOwGhtdCDUCxWoDtsWEakNjENENihwX93Tmzfhm2nkzb2Zozi+5mfvOfXfuN++e9t53z5u5P4jIJRHpKiirG9ccuGMZUmzvHRGZFpGRiK1NRAbs8w7Yz9Un8rWETMtF5EsRGbU+fK5BuBaJyJCIDBvXo+0iAwAABIVJREFUq2aP9Q8RWWjHY1a+IvJesT5YId8CETkqIgfqwqWq8zYRNvo5SPhy23Vm6yH8LLcAdwOHzd4GTNhrq+VbrWwIuMfqfAo8UgFTSyT/LNDXIFwPAU2W3wZss/xKYBhYSNhUaRxYYGkc6ASusHNWWp29wEbL9wFPVcB1K3ALYTOmroi9rlyzMBdlSNHX7wfuAkYitu3AFstvifRp2b6WkKmdsMcKwNXASeu3enMJsNjyzcBhay/WP4DN5P5ONwLvl/LBKvTlC8Bu4EApv02La77fKbwBvEj+1p69wE4NOgQsEZF24GFgQFXPqurvwACw1spaVPUbDVd8J/BoUiBVPR85vCrCVm+uz1T1oh0eIrcTXi+wR1UvqOopwh4Zqy2NqeqEqv4N7AF67W7lAaDf6r9bIdeoqp6IKaor1yyKZUipLQBU9SvgbIG5l/A5If/zluVrFTBNqeoRy/8JjBL2Xq83l6rqX3bYbEkp7h9R3n5gjflTMR9MLBHpANYBb9txKb9NhWveDgoisgH4RVWHC4qWAacjx5NmK2WfjLFXwvaaiJwGHge2NgpXRJvIbXJULte1wLnIAFNNrqgalasUW611g4ZdDbHX681e7rWrWLa0cSdhVl53LluiOQZMEwaZcYr7R7Z9K/+D4E9pXK8dhInsJTsu5bepcKW181pNJCKfAzfGFGWAlwlLIv+rFmPTBPZEXKr6kapmgIyIvAQ8DbzSCFx2Tga4COyaqVaknbgJRWpccdXS5qpAtWwriarmU3NqTGQxsB94XlXPS/HQV824VPVf4A4JsbMPCcuUxdqoCZeIrAemVfU7Eemepe3UuC7rQUFVH4yzi8hthLW0YXPADuCIiKwmjJrLI6d3AL+avbvAPmj2jpjzy+aK0W7gE8KgUHcuC+CtB9bYkhQluChi/41w299ks5dqXq+oUueqQKXYaqkzItKuqlO2DDNt9nJ9LbFEpJkwIOxS1Q8ahWtGqnpORAYJMYVi/jHDNSkiTcA1hKW6avfzvcAGEekBFgEthDuH2nJVGhS5HBLwM7lA8zryg1lDmgtmnSIEslot32Zl39q5MwHdngpYbo7knwH6G4RrLfAjsLTAvor8oNUEIZDaZPmbyAVTV1mdfeQHxjZXoQ8HyQ80NwRXEdaiDCn7+QryA82vkx/Q3Z7U1xLyCCHWtaPAXm+upcASy18JfE2YDMX6B2Fr4WhAd28pH6xSX3aTCzTXlCtVJ22URP6gIMBbhDXE78n/R7OJEJQZA56I2LuAEavzJvZN8IQs++29jgMfA8sahGuMsA55zFJfpCxjbZwg8oQT4WmRk1aWidg7CU9GjZlDL6yA6zHCzOcCcAY42Ahcc+COZUixvfeAKeAfu15PEtaXvwB+steZyUTZvpaQ6T7CssXxiF/1NADX7cBR4xoBtpbyD8KsfZ/Zh4DO2XywCv3ZTW5QqCmX/8yFy+VyubKat08fuVwul6t8+aDgcrlcrqx8UHC5XC5XVj4ouFwulysrHxRcLpfLlZUPCi6Xy+XKygcFl8vlcmX1HyOBz9P/J+BfAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Particle radius at start of run 4000.000000\n", "Particle radius at end of run 4002.483887\n", "Change in Particle radius 2.483887\n" ] } ], "source": [ "output.close()\n", "plotTrajectoriesFile('radialAnalytical.nc')\n", "\n", "print('Particle radius at start of run %f' % pset.radius_start[0])\n", "print('Particle radius at end of run %f' % pset.radius[0])\n", "print('Change in Particle radius %f' % (pset.radius[0] - pset.radius_start[0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Double-gyre example\n", "\n", "Define a double gyre fieldset that varies in time" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def doublegyre_fieldset(times, xdim=51, ydim=51):\n", " \"\"\"Implemented following Froyland and Padberg (2009), 10.1016/j.physd.2009.03.002\"\"\"\n", " A = 0.25\n", " delta = 0.25\n", " omega = 2 * np.pi\n", "\n", " a, b = 2, 1 # domain size\n", " lon = np.linspace(0, a, xdim, dtype=np.float32)\n", " lat = np.linspace(0, b, ydim, dtype=np.float32)\n", " dx, dy = lon[2]-lon[1], lat[2]-lat[1]\n", "\n", " U = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)\n", " V = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)\n", "\n", " for i in range(lon.size):\n", " for j in range(lat.size):\n", " x1 = lon[i]-dx/2\n", " x2 = lat[j]-dy/2\n", " for t in range(len(times)):\n", " time = times[t]\n", " f = delta * np.sin(omega * time) * x1**2 + (1-2 * delta * np.sin(omega * time)) * x1\n", " U[t, j, i] = -np.pi * A * np.sin(np.pi * f) * np.cos(np.pi * x2)\n", " V[t, j, i] = np.pi * A * np.cos(np.pi * f) * np.sin(np.pi * x2) * (2 * delta * np.sin(omega * time) * x1 + 1 - 2 * delta * np.sin(omega * time))\n", "\n", " data = {'U': U, 'V': V}\n", " dimensions = {'lon': lon, 'lat': lat, 'time': times}\n", " allow_time_extrapolation = True if len(times) == 1 else False\n", " fieldset = FieldSet.from_data(data, dimensions, mesh='flat', allow_time_extrapolation=allow_time_extrapolation)\n", " fieldset.U.interp_method = 'cgrid_velocity'\n", " fieldset.V.interp_method = 'cgrid_velocity'\n", " return fieldset\n", "\n", "fieldsetDG = doublegyre_fieldset(times=np.arange(0, 3.1, 0.1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now simulate a set of particles on this fieldset, using the `AdvectionAnalytical` kernel" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "X, Y = np.meshgrid(np.arange(0.15, 1.85, 0.1), np.arange(0.15, 0.85, 0.1))\n", "psetAA = ParticleSet(fieldsetDG, pclass=ScipyParticle, lon=X, lat=Y)\n", "\n", "output = psetAA.ParticleFile(name='doublegyreAA.nc', outputdt=0.1)\n", "psetAA.execute(AdvectionAnalytical,\n", " dt=np.inf, # needs to be set to np.inf for Analytical Advection\n", " runtime=3,\n", " output_file=output)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then show the particle trajectories in an animation" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "output.close()\n", "plotTrajectoriesFile('doublegyreAA.nc', mode='movie2d_notebook')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can also compute these trajectories with the `AdvectionRK4` kernel" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gr/T/parcels-504/cbcb799aee13754c9eae2649e4002775_0.so\n" ] } ], "source": [ "psetRK4 = ParticleSet(fieldsetDG, pclass=JITParticle, lon=X, lat=Y)\n", "psetRK4.execute(AdvectionRK4, dt=0.01, runtime=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we can then compare the final locations of the particles from the `AdvectionRK4` and `AdvectionAnalytical` simulations" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2dfZgcdZ3gP9/umZ4E4RQm7KmEvOhFgWUMhABbApcyo7zEF1TgDmV3ADkjLNkT9273IcfC5SHIcL48yyo5ZdQg4yIsBOWJCCvYoVjMFGyCCMPrQpDIyO6JI+yahUxNun/3R9X09Ezmpbunquulv5/nmWe6q6q7v1NT/a3v7/sqxhgURVGU9JOLWwBFURQlHFShK4qiZARV6IqiKBlBFbqiKEpGUIWuKIqSEdri+uAFCxaYJUuWxPXxiqIoqeTRRx/9rTHm0Kn2xabQlyxZws6dO+P6eEVRlFQiIrun26cuF0VRlIygCl1RFCUjqEJXFEXJCKrQFUVRMoIqdEVRlIygCl1RFCUjxJa2qCiKUo3rguOA3TmINXw32DZYVtxipQpV6IqixI7rQnc3eCOGQvndFHM/xurYCMWiKvU6UJeLoiix4zjgeVAqCx7tOOVT/A2OE7doqUIVuqIosWPbUChAPmcoMIqde8jfYNtxi5Yq1OWiKErsWJbvXXEcwe7chTX8YbC/rO6WOlGF3iiVCI6tF52ihIBljX2VuoIfpV5UodeJ64LTvxt783qs0s/8ZaEGbhRFSQCq0OugEonfezgFcw9FurG8Hb6lrgpdSRgV44MHsXqW6TXaAqhCr4NKJN7k/Ei8rMYqPK6BGyVxuC50f6CEN3IYBc6muHkNltOrSj3jaJZLHVQi8XkodOSwP/feOblbXBd6L9mNe0m//0RRQsI3PoQSbb7xMXqSpgC2AGqh18F4JB5sO49l9TT8XmpBKRMIOcjuGx8Gb2TUTwNs3w5275zfV0k2qtDrZDwSPzfGLag8HgZn9CQs9cW3JG7fIM6l92KXt4VWHWlZUHwgj9M/FPjQ1VhoBVShx4RaUAoEK7V1R+Dtu4oCl1McOTW0G7tvfCwGGl9JKulCFXoMuH2DOHcOc/3nD2f439rUgmphHAe8UhslxF+p5VZjaZBdaRBV6E3G7Ruk+3PvxuNICvd5FG/chbVWLahWxbah0CF+U6o82DecA5YW1SiNoVkuTca5cxiPwnj2wZ3DcYukxMhYoH3jNULxwXastarMlcZRC73J2Gd1UrjPwyNoQnRWZ9wiKTETVqC9pXFd3P7ncViF3bO4Zc+nKvQmY63toojvQ7fP6lSLTFHmiuvi2uvp9u7Bo0DhphLFB/ItqdRVoceAtbYLa23cUihK8mioXYHj4IyeNO7K9Eot241DFXqNtFpzxVb7e5X4abjYzrax29dT8AJXZiFXfzeOjFzwqtBnI/DNdd90Ht6+fEs0V6w0IfOg0FaieOEt2twpRjKia2al4WI7y8Jyein2b2nMh151wbv5k3E+c3Nq/fCq0Gci+Ec7e7+AZwwlxqdipfGfXSuVJmQl8EplnBufw7r54uzfyRJIxWr1hELBZNo3PKdiO8vCsiwaOjXBBe+Wjqe7dA/ejfMo3JzOy13TFmci+EfbZhsFPPJSaompWJUmZFLyv1hmm853jAmnfzfeiPE7fI6Ucfp3xy1SZIy1K9h48SsUL97SvN5GwQXvyGrfD29yqb3cVaHPRPCPtvI7KBbWsPFzQw3ftd2+QXpPc3D7BkMXM2wqudGfG6JYWIOV36HzHWPC5kHfmCCwWnkwbpEixbJg/TcWY32jp3nmcXDB2597L4UO8buppvRyF2NMLB+8cuVKs3Pnzlg+uy5CcGCOV4cWKDBWHZqSdMVWceAmlSAlzxk9Cbt9u3bkjJg0XO4i8qgxZuWU+1ShR0/vaQ5X3ncyJdrIM8rGU7ez/id23GIpaSENWiZCdPLSRGZS6BoUbQJaHarMiRYuJdW5AfWhCr0JaHWoojRGoucGJHDlpAq9SWh1qKLUT2LnBkyIbaxPzKpBFbqiKPXTJOs0qZOX3P7nx3vHeB7F/i1YCZBLFbqiKDVR0eGdg1iXjZUSR186ncTJSw6rqtpgGxxWNVbUFDI15aGLyOki8pyIvCAil0+xf5GIPCAij4nIEyKyJnxRFUWJi7Hq+Cuv9EfmuSMrglLilFbgzBG7Z7Gfsy4lCh057J7FcYsE1KDQRSQPbALOAI4CPiUiR0067K+A240xxwLnAv83bEEVpRVxXejt9X/HycR2EG04udWkugJnjlSqWr+YT1Q7hlpcLicALxhjXgQQkduAM4Gnq44xwH8IHr8VeCVMIRWlFUlSH5exdhC+l0Wwrz8HhucnKsOj2SQxm7QWhX4Y8HLV8yHgxEnHbADuE5E/A94CfHCqNxKRtcBagEWLFtUra2xM8B0O353uiziBqVbK1Ph9XA7zU/ZGRnH6hwJfcvMZawcxful0AZp+mzRqUegyxbbJ5aWfAr5rjPmqiFjA90TkaGNMecKLjOkD+sCvFG1E4GZT6aw5YiiU300x92Osjo3pbMWW0FQrZWr8Pi5njxek8SC1BgajuG8n0SJVJlKLQh8CDq96vpD9XSoXAacDGGNcEZkHLAB+E4aQcVLxHZbFH+pcPgXLeziVPXSTmmqlTI3Vs4zi5jXjfVx6asu/dvsG6V53BF6pjUKHpNL2UBqjliyXHcAyEVkqIgX8oOfWScf8CugGEJEjgXnAq2EKGheVVrK5wErKPZTaQNDEVKt2HFbFLZIyE8HghvVfPHB8NeW6uJf003vJ7qkDpa6Lc+kdeKPiGyEjphWTUFqWWS10Y8w+EVkH/ATIA5uNMU+JyNXATmPMVuB/AN8SkS/gu2MuMHF1/QqZcd+hYHfuwhr+MNhfDt/kaYJv2+5ZTOGmEp5X8sd0JSTVSpmBaj9HLcOQHQe7vI0Cl/uumjzYdnssoivNp6bCImPMPcA9k7ZdVfX4aeCkcEVLDuPfqYgCQYGj3h1ZgZN7E3vTgQ31e5ntnlCpupvhGCXB1DIM2baxOjZSHDkVJ7ca+4ZzggCm0gpopWgScBzckRV0l+/DKxcorDMUu+pTuBOCt/l9FG94dsqbgga2Ukwtw5CDJaXlOFi2DXNR5poRlTpUoScB28bJvYlXDiyvkqk75uo4vjIvlQWvDM6ld2B17dEvYpaodRhyGHftkFaNrYzb1/wOq6rQk4BlYW86kMI6g1cyFDqk7pirbUMhvw+vjB+8LW8DZ74q9Kwxl2HI9RDCqrGVGZ9SdiSF+zyKDDZFqetM0YRgre2i+GA7G69pLM3MsqB4w7NsbNtIMXcqVsfPU5mJoyQE28bJrR7315faNFumDpw7hydmlN053JTPVQs9Qcx1pWyt7fLdLM78aDJxlNah3lWj+tsnENeUMp0pqswd/TJnlpr+tZWIfHPa6aaFqHzoOlNUiQ79Mvtk9KY206qx8if/6nmsSitGL5VV1FEQx5QyVehKzUyps4LeCG7peJy9q7H7n2+9dgIteFOb8Ce3nUcxvxmLn6W2ijorqEJXamLaVq62jZs/me7SPXimQOEmodiTeX02kSAjxCmfgj3yUHKGGEfIhP7o5HE+ezPWou9nboWSNlShh0BGV9sTmLaVq2XhfOZmvBvnUTI5vH2tt+J2Oz9Cd/nzfjl+2aPYuSsR48iiZGJ/dL+tBNb6uMVqeVShV9OAZk7SEIIomamVq92zmMLNVV9uO1ZRm44z3IWXC4q6cnmc4a7MK/T9+6OH+OatYCFFhCr0MRr0gyZpCEGUzNTKNdIvdwqwbSh0yPg0H3vSARlVUFG0kXD7BnEuvRe7vC29cwdiRBU6wfdtwwj2yAqs8va6IvVzGUKQKoKyc8txwN5/MEYr94iZ8YamQ0VqxnX9AdTevqsocDnFkVNbIh4RJi2v0MebWq2iUL7Pr7Is1F5l2egQglTSylp7FqY7NTpUpHYcxx9AXULwMDi51X6DMaVmWr70f8JEotx8nA9eU98yb6ohBBHg9g3Se5qD2zcYyfsr0aBDRWpnzHWVzxkK7WBvOkcNiDppeQt9v2nmG2zqjmhFbLnG1ehHmTs6VKR2JgyTsdu1j3sDtLxCT0NAz2/0c2Rg5RmcO4ebXoEWFhmND06LDhWpD/XqzY2WV+iQ/IsorkY/YdMqKZ6TSfr1pdRJgq2SlvehpwFrbRfFG3ex8dTtFG/clVp3i5/iafwCpJEyTv/uuEVSlPoIspZ6r9iDa69n6knd8aEWekqIo9FP2LRMiqeSWZKetaQWutI0rJ5lFAtr2CgbKBbWYPUsi1uk1sB1obc3cdZkGkl61pJa6ErzmKU4SYmAFuwEGSVJz1pShT5HXNf3Dds86Fuc+mWZGY0QNpeqtojuyAq/InqD/gsaJelZSzqxaA5UsjZGDAU8342gpd1NIcGJBskisNArA59z8yl0NDa3VkkGM00sUh/6HPCNHxn3p42eRGom6abYrzp2I73yijLdHyil8U9oHkGhhfPBa/By8/2KaC89l2nWiPprpwp9DvhVpoY8o37WRvv2UHrHRq5rE556NRua/lgnloW9wfbL6vN1tDie64WYYqMhCsbCGVde6f+O4rSoD30OVPxp/UOBD33u7ha3b9DvOFdqi2xpnPTUq9nQ9Mf6qbsieq7BVA3G7seEKU8RjV7NjEKPy6fqx/gWE4pCcV2cS+/F23eV33FuxOA4EvrfMzH1yuCwKlUDGVqqw2WI1BWPnqv2aYb2Shn7TXmyw/+MTCj0zBgDjoNd3kaBy33rMw+23R76xyQ99WpWNP0xeqq0j5s/GedXn8Z26zjVYWuvDETBm9I3yhgTy89xxx1nwuLaa43J540B//e114b21s1lYMCY+fPNQO4kc23blWbgxici/ahrr/V/K8qUDAyYgYtvNvM79pl83pj58+u8XsK6yAYGzEBhlblW/pcZKKxq+YsW2Gmm0auZsNBtGwptJbwyFNrAtvNxi9QYwS3cchy/sX+E7UM1HVyZFcvCcSy8fQ16TiwLF8u3SGn8ekt7zGc/IlxtZEKhW7gUzXocTsI227Hopf6m5gkhyZo2A8tepT7m4jkJyxWa9pjPBCL2D2dCoeM4WKWfYZkHoZTXAEwUZCZQodTDXPy+YcVFUx/zqSbiYHE2FHozwsetjmYttCyNLhrD+lrOtdw+UQvLiHVVTQpdRE4H/gbIA982xlw3xTH/BdgAGOBxY8ynQ5RzZqrNiM7O8TK42P97GUJvmqkhKQoszKyORm8qiVtYRpzqMqtCF5E8sAn4EDAE7BCRrcaYp6uOWQasB04yxrwmIn8QqpS1YFm4gwfiXHoHdnkbVsfGBPz3MkQaZvUpiZsKFXdIKJELywhPSi2l/ycALxhjXjTGeMBtwJmTjvkssMkY8xqAMeY34Yo5O64L3euO4Mp9V9Fdvg93ZIU2rAgby/JbBjiWVnMnFG2LMJGxhWU+ZyjkRrE7B+MWKVJqUeiHAS9XPR8KtlXzHuA9IrJdRB4OXDT7ISJrRWSniOx89dVXG5N4GhwHvFLbeKOs3Gp1C4RMM3pRtDJu3yC9pzm4fY0rHb8tgjfeX4gHQ5QwfVgWFK8fZGPuf1MsfQDrshMzfeHW4kOXKbZN7rnbBizDTzddCDwkIkcbY16f8CJj+oA+8Nvn1i3tDNg2FDr8cvlCHuwbzok0j3sm3L5BnDuHsc/qTO38z6lI5PI1I7h9g3R/7t14HEnhPo8igzNeO9P5ybUtwv5Yw3djmWuhXAIv21lwtSj0IeDwqucLgVemOOZhY8wo8EsReQ5fwe8IRcoaGHfxCrbdjhWjMq/ni5kmNC4aHc6dw3gcOZ5rfefwtDNkZwz0aVuE/WmhC7cWhb4DWCYiS4FfA+cCkzNY7gI+BXxXRBbgu2BeDFPQWog7AAP1fTHThsZFZ6fR1Zl9VieF+7zxDpJndU57rOPg+8nL0zRwS8IXIUm00IU7q0I3xuwTkXXAT/DTFjcbY54SkavxewpsDfadKiJPAyXgL4wxw1EKnlTq+WKmEdUV0zOX1Zm1tositd0M7M5BCuV349FOoTyK3bkLyMYqMDJa5MKtKQ/dGHMPcM+kbVdVPTbAnwc/LU09X8xmk1XfflKY6+rMWttV0/HW8N0Ucz/GKZ+CnXsIa/jDqEJXICuVogmj1i9mM8mybz8pNG11ZttYHRuxvIcDn/CXo/kcpWHiMp5UobcIWfbtJ4Wmrc5axSeclJLXOonTeFKF3iJk3befFJq2OkuLT9h1cfufx2EVds/i2kVOXM1+7cRpPKlCbxGS7NtXMkowjLzSy/ymUu2tCIKiB7d0PM7e1dj9z6emB3qcxpMq9JRTz6o0ib59JcM4Ds7oSeO9zL1S7TU9to2bP5nu0j14pkDhJqHYkw4jPU7jSRV6inH7BuledwReqY1Ch6RpVaq0AraN3b6eghdYq4Vc7TU9loXzmZvxbpzn96XZl64Cz7iMp1p6uShJxHVxLr0Db1SqCkziFkrJJK4Lvb3190AJqlaLF29h48Wv1N350e5ZTGFejnw+8wWeoaEWehMJNWjvONjlbRS43Ld+8mDb7SFIqShVzDU4aVlYltXQyLhWSeYJE1XoTSJ090iQi1wcORUntxr7hnNi61+jZJiYO7KlJZknKahCbwaui3PpvXj7rqLENP036iUwXyzHwbLt2DpLKhmnhRpbZQFV6M0gKveImi9K1KjfI1VkT6GH5KgOtXRX3SNKmplkOKS0gLMlyJZCD6m6LPTSXXWP1IQqiuST4gLOliBbaYtTBXAaeZs7h8eLIWjHuTOETsCWBevX69U/DWPDja+8okz3B0rZmBLWYLpfGKPooiKkr5gSEem00Kcz5UIK4Gjfk+bjDzc+jBJ5vJFRnP4hLGtx3GLVTcVVd8zrWF//NO7ICpzcm9ibDqxplZf0rpgaI0026VPoM635LAv3+kfGfd8Nuja070njNDyxhwcpcPb4TZQHgZ7oBI2Aycr4ev6Yy7ger1ygsM5Q7Jp9gZb0rpgaI0026VPo/f2wdy8Ys19erOtC92Vdvq5/iJq+QNOhfU/qZ1brcobOe1kYbjxZGd8pZ+GZwHVXMjWlcKdhdajJVcklXQrddWHzZl+ZA7S1TVjz6VT6aKg1WDmjdTlb570MDDeerIzP+vQ8Hrrd4JUMhQ6pyT2hq0NlLqRLoTuOr60BRODCCyd88dW/Fz71VLjOaF3W0nkv5abf/sp4FV2X1u+e0NWh0ijpUuiTNXbPRB9rmvx7qUjRq7PCdUbrci6d91LEZGWc8nuUkjLEjLkvmszKlSvNzp07635dFgYdpyaXt7cX969+THf5Pn/CfDsUH2xvXNZGp9coTScVBkeNZEFnVCMijxpjVk61L1UWephBzzhJja8/7ArXOXTeU5rHWE2A5wmFgqm77W2SSHoaaNikqrAoE0UNrov9q34KbaXk93keq3C95sOs/4czMv1FUMbxawKMP1hipIzTvztukRomkiLBBJMqC33MhT4y4sdEO5OX0TUzga/F8jyK+c04n725btdD05fC6gRuObJQEzBGGtJAwyRVCt2y4PrrYd0630q/7DLoSpPbpWqJYfEzrEXfB2t9zS+v+N5HDIX8Poo3PKtWsxI6UdQExOWTb7U00FQpdIDhYSiX/Z9E+5+nYo55lY7jK/NSWfDK4Fx6B1bXnhSdACVWag1Kh1wT0FSf/BR3jlZKA02dQk91rvkc8yptGwr5fXhl/OVjeRs481WhK7MzW2HXZEJ0tTXUp6cRkz416WPRkaqgKIzrxI0bU/r/mkPXRcuC4g3PsrFtI8XcqVgdP0/ZHU2Jjf0Ku5qXVOD75D3yjFb55GcguPn0XrEH115fe7dKx8EdWUFv6S9wR1akNGtibqTOQofWjtNZa7t8N4szH+wv138ispRgrNROjIVd9frk3f7nx1cSnkexfwtWDdeq2/kRusuf919X9ih27mq5FNlUKvSWp9E72qQlqXv9IzjDXarbW4HAL17s39L8wq46ffIOq6pSDQ0Oq2pSzM5wF14uiDHl8jjDXarQlQxTlWXjjqzwe7SUW9bd2HrEWdhVhxFi9yymcFMJzyv5K4me2vri2zYUOiSwV2prhpY1VKG3ElURZUdW45XaKKUxW0jJNJYFxQfy9Tc1S1Evp6hQhd5KVF3xdudHKFwm6cwWUjJPo17FVo6vgSr01iO44i38XjhZsGY0zqsoPjUpdBE5HfgbIA982xhz3TTHnQ3cARxvjKm/laLSVLJgzWSpkZSizJVZ89BFJA9sAs4AjgI+JSJHTXHcQcB/Bx4JW0hFmY4sNZJSlLlSS2HRCcALxpgXjTEecBtw5hTHbQS+BOwNUT5FmZG6i1YUJcPUotAPA16uej4UbKsgIscChxtj7p7pjURkrYjsFJGdr776at3CKspkrJ5lFAtr2CgbKBbWYPUsi1skJSW4fYP0nubg9g3GLUpo1OJDlym2VcYciUgO+GvggtneyBjTB/SBP7GoNhGThwbhEkQGhksrzSergy9qUehDwOFVzxcCr1Q9Pwg4GnBEBODtwFYR+VgWA6MahEsgWYjuRkVgfbidH9Gq4Cr8wRdHjlej3jmciY6MtSj0HcAyEVkK/Bo4F/j02E5jzL8CC8aei4gD/M8sKnNosHOcUj86f3TuBK0e3JEVfo+TnKHQIVoVTHYHX8yq0I0x+0RkHfAT/LTFzcaYp0TkamCnMWZr1EImibinubSEu6feVq/K1AStHpzyKX5vlLJoVXBAVgdf1JSHboy5B7hn0rarpjnWnrtYySWKaS61Eru7p1l3k/1avZZUCTVC0OrBHnmIQtnDy+XD73GSYgsji4MvtFK0XmIMwsXq7mnm8IAYW71mirEh345DsXNX+D70qmvCzZ+M85n6Z+Qq4aIKvRFiCsLF6u6p6tQY+bo9zlavWaOq1UPopzC4JtzS8XSX7sG7cR6Fm7VzZ5yoQp8Fty85frZY3T2dH8GRN7Fz27AKTZiUFGerV6U2ApeOs3c1nin41brqo48VMSaedPCVK1eanTuTnQgznqtaoIBH8cZdoSv1ul2QDfgsKy/pHMQavrtuf2dlZT1iKOT3Ubzh2dhvbkpCCLKRum86D29fXnvrNwERedQYs3KqfWqhz0DUuaoNuaXrdPdMUMbld1PM/Riro76BrBVvS1nwpL0lJ8Eo0xCspIo9qY2NZorUDYluJvZZnRP7hIScq+r078bbW57glg6bCcqYdpzyKXV/2NhcjHxee6crUzOH2edKiKiFPgOR5qq6Lvbm9RTMPXi0U2jLYdv58N4/oDKkaMRQKI9i5x6qWyvrJBglEaQ4RbJZqA89Lnp74corcUvH48hq7M+9F+sb0WSszNWHHjtaNdryuH2DOJfegV3ehtXx85Z21M/kQ1eFHhfNzOuOizAsqslVox2iVaMthutC96pRvFHxkxNyp2Jd82Hfx9OCaFC0Sbiu7xe3edBv4zqT1sm6HyNQxH6K5Xosp8EirOmqRlGrPZU0cJN3HPyB5oifnJBbjaWBnKkxxsTyc9xxx5ksMTBgzPyOfSbPqJnPv5uBwip/Y4sycPHNZj7/Pn4+Lr65wTcaMAOFVcF7eWZ+xz4zcOMTVdtG/W2te6rTw8CAMfPnm4HcSebativNwI1P1PMyk8+Vzfx2r+bXZRX8HlpT6lW10EPCzyYRvywfgzN6kt8eoEVNR4dV41Y1BodVjaU6TlU16txNr/Z6SR+OE3R+vA+vXKCwzlDsmv3/Nr6YFWy7HcvSGojpUIUeEn42icEbCVIc27f7vV5aFLtnMYWbSnheye/F0jOHnjP7VY1qr5dUYts4uTfxysGNuGRqvhFry/vaSGdQNKHpS3X50FuASP9NmvkSKs26dt2+QbrXHYFXatPe7A2SrSyXVsgOUZQmUmnLPGL8LJLCmsaD2DV+XgLtsdQwk0JPX6XoVF3/FEVpmPH4T5tfTTx6UqTfK60qjY70KXStQ1cawXVxL+mn95LduG7cwiSLsfhPpcVF+3b9XqWU9AVFLQv3+kfGy/E14q3Mho60mxHLguIDeZz+ocCH3tzBLUp4pE6huy50X9blu9Afoqa0J6XF0ZF2s+JnkSymmfNxlfBJnctFXehK3dg2dvv28c6Z6qlTMkrqLPRK90BPXehNI+1pCTrSTmkR0pe2SPr1S6qY0JNle6TpbErzSdKIRaU2MtecS6vGmofb//x4MNHzKPZvwdKTnwnGRyweSeE+jyKDqtRTTup86EpzmdiTpR2HVXGLpISEP2Kx6n9753DcIilzRBW6MiN2z2IKHUJeShQ65tiTRUkUUY9YBHz/aG8vmvzfHFLpclFCJAhIuJ0fwRnu2i8uUclRdjRmkTUiHbEI2qYjBlShtzLBF85vafp5vJyZsmGSxiyyi7W2C2ttRG8etMt1yqdgjzzU0u2km4Uq9ITR1AyeIKnfKZ/i+1LLUsnt1++dMlfczo/4hgIFCmWPYueuxnriKzWjCj1BNH2FGiT12yMPUSh7eLk8hYJobr8SCs5wF17O+IZCLo8z3KUKPWJUoSeICVWwe0s4/UNBOXZEBKNgLMeh2LlrSh+6ojSKbUOhQwIDRQ2FZpDKwqKsMt6XukyB0cj7UitK1GgRYPhkrrCo6TTpqrQsKF54C86Nz2GbbVilHerQVlKNBtSbiyr02ZhQ+r4+covZ6lmGdfPFkTSrUWtJUbKNKvRZaHrp+/iI81A1r6YEK0r2qalSVEROF5HnROQFEbl8iv1/LiJPi8gTIlIUkcyUE8ZS+h7BjK79A667Q3tvRVGSwawKXUTywCbgDOAo4FMictSkwx4DVhpj3gdsAb4UtqBxkZXSd9uGQlvJL/M2I9ibz9dybEXJGLVY6CcALxhjXjTGeMBtwJnVBxhjHjDGvBE8fRhYGK6Y8TFW+r7xi/lUjy0bC7hulA0U6cYq/UyngyhKxqjFh34Y8HLV8yHgxBmOvwi4d6odIrIWWAuwaNGiGkWMn6xE6qMMuCqKEj+1KHSZYtuUyesi8sfASpja0WyM6QP6wM9Dr1FGJSwiCrgqipIMalHoQ8DhVc8XAq9MPkhEPghcAawyxoyEI54SOllZbqh1qHIAABAnSURBVCiKsh+1+NB3AMtEZKmIFIBzga3VB4jIscCNwMeMMb8JX0wljbgu9F6yG/eSfg3AKi1B3O3fZ7XQjTH7RGQd8BMgD2w2xjwlIlcDO40xW4EvAwcCd4gIwK+MMR+LUG4l4Yy3MTiMAmdT3KxtDJRsk4Raj5oKi4wx9wD3TNp2VdXjD4Ysl5Jy/Lx3oUQeD4MzepL2w1YyzYRaj5jaUOsIOiUS/M68Zny8Wfv2urNq3L5Bek9zcPsGI5FRUcIk6EZNPh9fEpmW/iuRUBld1z+EzYNYPfW5W3QivZI2LFyK5z+PwyrsnsWxLEZVoacE1wWnf3egHJelwnXhJ9QsBnrqfq0/kf7IoOWCwblzOLpRaYoyVwIHuuV5WIUC9BQhhnEeqtBTQFYCjKOjowwNDbF3795Zjz1j48Gs+MJzGATB8B8POZhnnnmmCVKmh3nz5rFw4ULa29vjFkVJggMdVejJparXreNYoQcY47D4h4aGOOigg1iyZAlBNtSM7Hn1DX7/WomDDs5z4KEHRC5fmjDGMDw8zNDQEEuXLo1bHGXMgR5zFbYq9ATi9g3iXHovdnkbVsdG7OsfoVA4Cm+kOsDY2/j7x2Tx7927t2ZlDnDgoQdw4KGRipRaRITOzk5effXVuEVRIDFV2KrQE4brQve6I/D2XUWByymOnIo1fDfFB7oaDjBOJs6UwlqVuTI7ei4TRgKqsFWhJwzHAa/URgnxlW1uNZZtzynAOJmxlMKwLH5FUZKB5qEnjLFJ6fmcodAO9qZzQr/rV1oCX/wKxYu3pDLA2ij5fJ5jjjmGo48+mo9+9KO8/vrrALz00kscffTRleO+9a1vsWLFCl577bXKtq985SuICL/97W+bLrei1IJa6Alj3BUn2HY7lhVN7nWYFn+khDwIdf78+fziF78A4Pzzz2fTpk1cccUVE4753ve+x9e//nW2bdvGwQcfDMDLL7/M/fffn6q2z0rroQo9gczmiktjTnpDRNwcw7IsnnjiiQnbbr/9dq677jqKxSILFiyobP/CF77Al770Jc4888zJb6Mo+xPTRHZV6CkjKznpNRFhbm+pVKJYLHLRRRdVtu3evZt169bx2GOP8fa3v72yfevWrRx22GEsX748lM9WMk6MXbrUh54yxjNUgqHVoydld5RcBM0x3nzzTY455hg6Ozv53e9+x4c+9KHKvkMPPZRFixZx++23V7a98cYbfPGLX+Tqq6+e82crLcJUhkiTUIWeMsJoepUaxgIKGzeGZuWM+dB3796N53ls2rSpsu+AAw7g3nvv5Zvf/Ca33HILALt27eKXv/wly5cvZ8mSJQwNDbFixQr+5V/+Zc6yKBklxi5d6nJJGXNtepU6Isrtfetb38rXvvY1zjzzTC655JLK9kMPPZS///u/x7ZtFixYwGmnncZvfjM+s2XJkiXs3Llzgn9dUSYQY5GRKvQUkpoMlYRz7LHHsnz5cm677TZOOeWUyvalS5eydetW1qxZww9+8ANOPHGmmeiKMgUxFRmpQldaij179kx4/qMf/ajy+Mknn6w8Xr58Ob/+9a/3e/1LL70UmWyKMlfUh64oihISiZ8pqiiKosxOEmaKtoSFrtPnFUWJmhizFStk3kJvqUIcRVFiIwkt0TNvobdUIY6iKLERQdlE3WTeQtdWsYqiNIu4W6Jn3kJv5VaxytT88Ic/RER49tlnG36PCy64gC1btsx4zLXXXjvh+fvf//6GPmvDhg185Stfaei1SmuReYUOvv5e/43FWN/oUWWeMqJIA7v11ls5+eSTue2228J70ymYrNAHBgYi/TxFaQmFrqSTsTSwK6/0f4eh1Pfs2cP27dv5zne+U1HojuNg2zZnn302RxxxBOeddx7GGACuvvpqjj/+eI4++mjWrl1b2T5GsVjkE5/4ROX5/fffzyc/+Ukuv/zySiOw8847D4ADDzywctyXvvQlurq6WL58OZdffjngD9U4/vjjWb58OWeddRZvvPHG3P9gpaVQha4klijSwO666y5OP/103vOe93DIIYfw85//HIDHHnuM66+/nqeffpoXX3yR7du3A7Bu3Tp27NjBk08+yZtvvsndd9894f1Wr17NM888UxnWfNNNN3HhhRdy3XXXVRqBjTX6GuPee+/lrrvu4pFHHuHxxx/nL//yLwH45Cc/yY4dO3j88cc58sgj+c53vjP3P1hpKbKr0OMu2VLmTBRN62699VbOPfdcAM4991xuvfVWAE444QQWLlxILpfjmGOOqZT4P/DAA5x44ol0dXWxbds2nnrqqQnvJyL8yZ/8CX/7t3/L66+/juu6nHHGGTPK8NOf/pQLL7yQAw44AIBDDjkE8FsPnHLKKXR1dXHLLbfs91mKMhvZzHIJ1uruyAqc3JvYmw7EWrv/KLeYhoooNRJ207rh4WG2bdvGk08+iYhQKpUQEdasWUNHR0fluHw+z759+9i7dy9/+qd/ys6dOzn88MPZsGEDe/fu3e99L7zwQj760Y8yb948zjnnHNraZv5aGWMQkf22X3DBBdx1110sX76c7373uziaXqvUSTYtdMfBHVlBd/k+rtx3Fd3rjtjPUJ/gn/1ASatIE4plwfr14dxwt2zZQk9PD7t37+all17i5ZdfZunSpfzsZz+b8vgx5b1gwQL27NkzbVbLO9/5Tt75zndyzTXXcMEFF1S2t7e3Mzo6ut/xp556Kps3b674yH/3u98B8Pvf/553vOMdjI6O7uemUZRayKZCt22c3Go8Cn5BUaltP//rBP/sSBnnxufCi7wpieTWW2+dEMAEOOuss/j+978/5fFve9vb+OxnP0tXVxcf//jHOf7446d97/POO4/DDz+co446qrJt7dq1vO9976sERcc4/fTT+djHPsbKlSs55phjKimJGzdu5MQTT+RDH/oQRxxxRKN/ptLKGGNi+TnuuONMlAzc+ISZ3+6ZfK5s5s83ZmBg0v4BY+bPNyYv+8x8/t0M8EfG5PPGXHttpHK1Mk8//XTcIkTGpZdear797W83/XOzfE6VqQF2mmn0ajZ96IC1toti1/T+14p/tn8Ie/P5WKUd8TVgUFLNcccdx1ve8ha++tWvxi2K0uJkVqHD7GW4lck/Pb0aHVUa5tFHH41bBEUBalToInI68DdAHvi2Mea6Sfs7gH7gOGAY+K/GmJfCFTVC4m7A0EKYaTI8lPoxk4qcFGXWoKiI5IFNwBnAUcCnROSoSYddBLxmjPlPwF8D/ydsQZX0M2/ePIaHh1URhYAxhuHhYebNmxe3KEqCqMVCPwF4wRjzIoCI3AacCTxddcyZwIbg8RbgBhERo99cpYqFCxcyNDRUqapU5sa8efNYuHBh3GIoCaIWhX4Y8HLV8yFg8hj0yjHGmH0i8q9AJ/Db6oNEZC2wFmDRokUNiqyklfb2dpYuXRq3GIqSWWrJQ5/K4TnZ8q7lGIwxfcaYlcaYlYceemgt8imKoig1UotCHwIOr3q+EHhlumNEpA14K/C7MARUFEVRaqMWhb4DWCYiS0WkAJwLbJ10zFbg/ODx2cA29Z8riqI0F6lF74rIGuB6/LTFzcaYL4rI1fgVS1tFZB7wPeBYfMv83LEg6gzv+SqwuwYZFzDJF58gVLb6SapcoLI1isrWGI3KttgYM6XPuiaFHicistMYszJuOaZCZaufpMoFKlujqGyNEYVs2WzOpSiK0oKoQlcURckIaVDofXELMAMqW/0kVS5Q2RpFZWuM0GVLvA9dURRFqY00WOiKoihKDahCVxRFyQixKnQROV1EnhORF0Tk8in2d4jI3wX7HxGRJVX71gfbnxOR05os15+LyNMi8oSIFEVkcdW+koj8IviZXIDVDNkuEJFXq2T4b1X7zheR54Of8ye/tgmy/XWVXP8kIq9X7YvsvInIZhH5jYg8Oc1+EZGvBXI/ISIrqvZFfc5mk+28QKYnRGRARJZX7XtJRAaDc7YzBtlsEfnXqv/bVVX7ZrwWmiDbX1TJ9WRwfR0S7IvsvInI4SLygIg8IyJPicjnpzgmuuttulFGUf/gFyntAt4FFIDHgaMmHfOnwDeDx+cCfxc8Pio4vgNYGrxPvolyfQA4IHh8yZhcwfM9MZ+zC4AbpnjtIcCLwe+Dg8cHN1O2Scf/GX6RWjPO238GVgBPTrN/DXAvfk+iPwIeacY5q1G29499Jn4L60eq9r0ELIjxvNnA3XO9FqKQbdKxH8WvXo/8vAHvAFYEjw8C/mmK72hk11ucFnqlLa8xxgPG2vJWcyZwc/B4C9AtIhJsv80YM2KM+SXwQvB+TZHLGPOAMeaN4OnD+P1tmkEt52w6TgPuN8b8zhjzGnA/cHqMsn0KuDXEz58WY8w/MHNvoTOBfuPzMPA2EXkH0Z+zWWUzxgwEnw3NvdZqOW/TMZfrNArZmnmt/bMx5ufB498Dz+B3o60msustToU+VVveyX/4hLa8wFhb3lpeG6Vc1VyEf7cdY56I7BSRh0Xk4yHJVK9sZwVLuS0iMtZYLcpzVtf7By6qpcC2qs1RnrfZmE72qM9ZvUy+1gxwn4g8Kn5r6jiwRORxEblXRP4w2JaY8yYiB+ArxTurNjflvInvIj4WeGTSrsiutzhnis6lLW9N7XobpOb3FpE/BlYCq6o2LzLGvCIi7wK2icigMWZXE2X7EXCrMWZERC7GX+GsrvG1Ucs2xrnAFmNMqWpblOdtNuK4zupCRD6Ar9BPrtp8UnDO/gC4X0SeDSzXZvFz/L4ie8Tv93QXsIwEnTd8d8t2Y0y1NR/5eRORA/FvIpcZY/5t8u4pXhLK9RanhT6Xtry1vDZKuRCRDwJXAB8zxoyMbTfGvBL8fhFw8O/QYTGrbMaY4Sp5voU/57Wm10YtWxXnMmkJHPF5m43pZI/6nNWEiLwP+DZwpjFmeGx71Tn7DfBDwnM71oQx5t+MMXuCx/cA7SKygISct4CZrrVIzpuItOMr81uMMT+Y4pDorrcoAgM1Bg/a8J3+SxkPnPzhpGMuZWJQ9Pbg8R8yMSj6IuEFRWuR61j8oM+ySdsPBjqCxwuA5wkxGFSjbO+oevwJ4GEzHnD5ZSDjwcHjQ5opW3Dce/GDUtKs8xa87xKmD+59mIlBqn9sxjmrUbZF+DGi90/a/hbgoKrHA8DpTZbt7WP/R3yl+KvgHNZ0LUQpW7B/zAB8S7POW/D39wPXz3BMZNdbqCe4gT9+DX4UeBdwRbDtanyrF2AecEdwQf8j8K6q114RvO454Iwmy/VT4P8Bvwh+tgbb3w8MBhfwIHBRDOesF3gqkOEB4Iiq134mOJcvABc2W7bg+Qbgukmvi/S84Vto/wyM4ltBFwEXAxcH+wV/EPqu4PNXNvGczSbbt4HXqq61ncH2dwXn6/Hg/31FDLKtq7rWHqbqpjPVtdBM2YJjLsBPnqh+XaTnDd8lZoAnqv5na5p1vWnpv6IoSkbQSlFFUZSMoApdURQlI6hCVxRFyQiq0BVFUTKCKnRFUZSMoApdURQlI6hCVxRFyQj/H8CPXqVmyoBJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(psetRK4.lon, psetRK4.lat, 'r.', label='RK4')\n", "plt.plot(psetAA.lon, psetAA.lat, 'b.', label='Analytical')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final locations are similar, but not exactly the same. Because everything else is the same, the difference has to be due to the different kernels. Which one is more correct, however, can't be determined from this analysis alone." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bickley Jet example\n", "\n", "Let's as a second example, do a similar analysis for a Bickley Jet, as detailed in e.g. [Hadjighasem et al (2017)](https://aip.scitation.org/doi/10.1063/1.4982720)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def bickleyjet_fieldset(times, xdim=51, ydim=51):\n", " \"\"\"Bickley Jet Field as implemented in Hadjighasem et al 2017, 10.1063/1.4982720\"\"\"\n", " U0 = 0.06266\n", " L = 1770.\n", " r0 = 6371.\n", " k1 = 2 * 1 / r0\n", " k2 = 2 * 2 / r0\n", " k3 = 2 * 3 / r0\n", " eps1 = 0.075\n", " eps2 = 0.4\n", " eps3 = 0.3\n", " c3 = 0.461 * U0\n", " c2 = 0.205 * U0\n", " c1 = c3 + ((np.sqrt(5)-1)/2.) * (k2/k1) * (c2 - c3)\n", "\n", " a, b = np.pi*r0, 7000. # domain size\n", " lon = np.linspace(0, a, xdim, dtype=np.float32)\n", " lat = np.linspace(-b/2, b/2, ydim, dtype=np.float32)\n", " dx, dy = lon[2]-lon[1], lat[2]-lat[1]\n", "\n", " U = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)\n", " V = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)\n", " P = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)\n", "\n", " for i in range(lon.size):\n", " for j in range(lat.size):\n", " x1 = lon[i]-dx/2\n", " x2 = lat[j]-dy/2\n", " for t in range(len(times)):\n", " time = times[t]\n", "\n", " f1 = eps1 * np.exp(-1j * k1 * c1 * time)\n", " f2 = eps2 * np.exp(-1j * k2 * c2 * time)\n", " f3 = eps3 * np.exp(-1j * k3 * c3 * time)\n", " F1 = f1 * np.exp(1j * k1 * x1)\n", " F2 = f2 * np.exp(1j * k2 * x1)\n", " F3 = f3 * np.exp(1j * k3 * x1)\n", " G = np.real(np.sum([F1, F2, F3]))\n", " G_x = np.real(np.sum([1j * k1 * F1, 1j * k2 * F2, 1j * k3 * F3]))\n", " U[t, j, i] = U0 / (np.cosh(x2/L)**2) + 2 * U0 * np.sinh(x2/L) / (np.cosh(x2/L)**3) * G\n", " V[t, j, i] = U0 * L * (1./np.cosh(x2/L))**2 * G_x\n", "\n", " data = {'U': U, 'V': V, 'P': P}\n", " dimensions = {'lon': lon, 'lat': lat, 'time': times}\n", " allow_time_extrapolation = True if len(times) == 1 else False\n", " fieldset = FieldSet.from_data(data, dimensions, mesh='flat', allow_time_extrapolation=allow_time_extrapolation)\n", " fieldset.U.interp_method = 'cgrid_velocity'\n", " fieldset.V.interp_method = 'cgrid_velocity'\n", " return fieldset\n", "\n", "fieldsetBJ = bickleyjet_fieldset(times=np.arange(0, 1.1, 0.1)*86400)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add a zonal halo for periodic boundary conditions in the zonal direction" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "fieldsetBJ.add_constant('halo_west', fieldsetBJ.U.grid.lon[0])\n", "fieldsetBJ.add_constant('halo_east', fieldsetBJ.U.grid.lon[-1])\n", "fieldsetBJ.add_periodic_halo(zonal=True)\n", "\n", "def ZonalBC(particle, fieldset, time):\n", " if particle.lon < fieldset.halo_west:\n", " particle.lon += fieldset.halo_east - fieldset.halo_west\n", " elif particle.lon > fieldset.halo_east:\n", " particle.lon -= fieldset.halo_east - fieldset.halo_west" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And simulate a set of particles on this fieldset, using the `AdvectionAnalytical` kernel" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "X, Y = np.meshgrid(np.arange(0, 19900, 100), np.arange(-100, 100, 100))\n", "\n", "psetAA = ParticleSet(fieldsetBJ, pclass=ScipyParticle, lon=X, lat=Y, time=0)\n", "\n", "output = psetAA.ParticleFile(name='bickleyjetAA.nc', outputdt=delta(hours=1))\n", "psetAA.execute(AdvectionAnalytical+psetAA.Kernel(ZonalBC),\n", " dt=np.inf,\n", " runtime=delta(days=1),\n", " output_file=output)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then show the particle trajectories in an animation" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "output.close()\n", "plotTrajectoriesFile('bickleyjetAA.nc', mode='movie2d_notebook')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Like with the double gyre above, we can also compute these trajectories with the `AdvectionRK4` kernel" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO: Compiled JITParticleAdvectionRK4ZonalBC ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gr/T/parcels-504/2c3dc13dcb5affb33a37376ddd7e7666_0.so\n" ] } ], "source": [ "psetRK4 = ParticleSet(fieldsetBJ, pclass=JITParticle, lon=X, lat=Y)\n", "psetRK4.execute(AdvectionRK4+psetRK4.Kernel(ZonalBC),\n", " dt=delta(minutes=5), runtime=delta(days=1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally, we can again compare the end locations from the `AdvectionRK4` and `AdvectionAnalytical` simulations" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD4CAYAAADo30HgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOyde5gcVZnwf6eqL4nrskATDQsksGw+EJwNCTFSckmH4S5IXNRFWBsBMwQSJKgbmVU0azTB4GUWiclMJGyaR0WUj5sGCYypcJkSCAQJN78AEshiMDSwLEi6uqve74+q7ulJZpJJpu99fs/T09OnL3W6qrre896ViKDRaDQaTQGj1hPQaDQaTX2hBYNGo9FoBqAFg0aj0WgGoAWDRqPRaAagBYNGo9FoBhCp9QRGyn777ScHH3xwraeh0Wg0DcVjjz32uoiMGey5hhcMBx98MOvWrav1NDQajaahUEptGuo5bUrSaDQazQC0YNBoNBrNALRg0Gg0Gs0AtGDQaDQazQC0YNBoNBrNALRg0Gg0Gs0AtGAAHAcWLQruCwPOpWkWXbqpf0yj0WhahIbPYxgpTs8G2uccjutFiMUVvV0b4PJO2t1VuMSI3ejRu8bEsmo9U41GowlwHLBtSCapyLWptQWD42DPvhs3/w08FNveE9I35MA9l23EEUxc18O2wcLBSW/EZhrJ1HgtKDQaTVUpCINEAuZ+0cN1FbGYVGTh2rKCwenZgH3tsyTyW4iQx8NEgJ+s+ycURyIYgBAxhWRiA860q5ieuweXKLEbfNasNbRw0Gg0VcFxoH16IAyUEnwffAzcbA47vRnLGl/W7bWkYHB6NtB+yaG4fIgY53I6q7iDsxFMPF8R7BaFwuPCyU9irV/OpblzyRIHFNmckE5XRoXTaDSaAgUt4eVHtuBm98PDBPEw8THxiZEjyVogVdbttqRgsG/N4PIhPCK4CGN5jVFkcYlimgrl5cmHOz11cRTWD/FBjgPpNM6WQ7DH/os2MWk0mrJRqiUY7EfQhFkAA0WeL3AjqdjNWKlFZd92SwqG5DkJYqtdXCS4+M/bn9Tbvyr6D9jwLPatGZLnJLA62sBJkfrJV1mRv5AcUaJRSE16BpJJHHcy0/mhNjFVG2cQn89gYxpNA+I4MH/um2Sze+Fj4KEgFA2gEAzG/WMMK72oMqYLEWno29FHHy17Ql/3k7LwlDXS1/3kMN/QJ32zVsrCWS9JX5+ILFwoopTM4scCvoAI+DJr1h5NR7M79PVJX2yajOZdMcnJaDMrffNuGzgWzQXHqa8vOFZ9fbWetUYzLPr6REbH82KQD68t/oBrjCIvo3l3+NeuIQDWyRDX1ZbUGACsjjasjt15g4VlWfTL5iREo+CWfWqaEopheYkNsH499pbDSfzxYW51O8kSwydC1ssz/9rR/IP8Cy4xPCJkc3nmf+EV5j/fCbkctnqX5Fdew/rujFp/JY1mp6QXb2FbdgyCSaAlCKAwyBMhz0UH3kfq6vGBNaNCqEBwNC5TpkyRmvVjcBycxQ+QvPNKchIhGlVBaKs2YZSFYo5J3sQUFwXkiOBjovAQTAx8fAwMPCLkB7zGwCdCDgXkiRAhz4Uz3iQ1b38sKhwIrtHsAY4DyeNyuH7pmj0QCidxH/Oji7DWXlOWc1Yp9ZiITBnsuZbVGMqCZWHdZmEPdo0JHdMApFL64rO7hDkm2fw38DHwiYUW1iCMWIhgkOcf+BMvcgg+ETylmMlyXpRDuI/2UCwECCYeJt23f5CVd3v0SieW9yDEYtDbq4+Ppi6w05vw/ANggE9BiJJnPt/CunhiVc5VXRKjDFgWdHZuJxSSSZxlT7Bo2d44065C19bYTWybRH4LfqhOC4rSH4tBnnhc8W/zDOJRwTSE2CiT1LJjmT/rNeJxhWlC1AwCDBQeAIKB64KdOxY8D+e9o1g09zV9eDR1QeKZBzDwAT8cUSh8LuRGrPjjwSKzCmiNoRLYNo57NO3cR5YYRs5nyeL76Lit1hNrDJyeDaR7JvE400OTUYSCnVXhcfZ+DlM/dVAYeXQobTNKNbY2oI3eVGEsAhueJX1DjhXrJ+L5JrEIJP0HcXLHMJ3f4T4SJTZNR5RpaovTs4G59/8zeQyCRRBQiJzkJrjwwupptkN5pRvltqdRSRWlr08Wml8Tg1wxkiBqejowZhf09YnMmvGqRMkOiMZQ5AR8MciNKBpjQIDSrFk6okxTV8ya+pgo8sXzEUQUeZnFj0Xi8bJH1qGjkqqMZZH88fsxZoEvwUrXE+2YHoqCO+bGG3zc3AdLzEag8DhZreGc82Jktkp/bskeYFml+z8Fy58itDDtOCHtmNZUEceBFY/9U7EUT0CNtAUYucYAHASsAZ4FngauCMf3Be4FNob3+4TjCrgOeB54Ephc8lkXhK/fCFwwnO3XpcYQ0t0tEo2KGIbI6NE6lH4wCjHbCm/A6r1wi/Ge9M34bmW23f2kxMycKOVLLCb9eQ/xuIhSFVmlaTSDMWvGn3fQFsCXGdwaXEQqcB5SYY0hD3xZRB5XSv0t8JhS6l7g80CviFyjlLoKuAr4KnA6MCG8fRRYCnxUKbUv8E1gCoHIfEwpdaeIvFmGOdaEjg5oa9OLz51hpzfhZg8oxmyrMM/zWB7iCJ4lFfsF1rzyp/xDkMtib398Lk1DNguAk52EPfc1kl362Gkqh+PAirv2205bKOHKK6t+Ao5YMIjIn4E/h///r1LqWeAA4GwgGb5sJWATCIazgXQosX6vlNpbKbV/+Np7ReQNgFC4nAb8fKRzrCUDzRdoM8V2BFEYn0HIo5AwNFXxKFO5ZsYjgVCo4H7a4fiEOBxDO724j8aJteuIVk3lsG3I+4HDWYURScFCCVZxBs7b71DtU6+s4apKqYOBScDDwAdDoVEQHh8IX3YA8ErJ2zaHY0OND7adDqXUOqXUuq1bt5bzK1QWx4H2drj66uC+xWMkC1EYHkYoFIIaMEKELDHsv06t/tU4lYJYDJvpZInhiUk2G/x4NZpKkEgEIRaFsOzjeTAUEIo8Eewth1d9TmUTDEqp9wO3AnNF5O2dvXSQMdnJ+I6DIj0iMkVEpowZM2b3J1srbBtcFzwvuG/xq036P99gG3H8MFM5WCUFPxBF4GiuOpYV5FDMOC7MoQDfD368Gk0lyGSCsniFcOx9eaNoVvIxSbzx/6o+p7IIBqVUlEAo/FRE/m84/FpoIiK8/0s4vpnAYV3gQODVnYw3D8lkkGlrmjjmcSx6pB3n0nRLag5OzwZWPPPR4g9AFRN6As6a8FxFa8HsFMsiM/UMDCNYqxhG8OPVaCpBIkEYiSdFE1KQkBkIiszr1S9bNGLBoJRSwA3AsyLyg5Kn7iSIMiK8v6NkPKUCjgH+JzQ13QOcopTaRym1D3BKONY8WBb09uLMXEG73MfVt08muewzXHr8BpyeDbWeXVWxb82QLzZE8jlu9OPhM8GP4PT9/1CzuUEgw+NxMM3gPpkMn3AcWLSoJYW5pjJkMoTZzgrweYN9i8EYgslbz/25+ufbUOFKw70BxxH8mp8EnghvZwAJoJcg9LQX2Ff6w1WXAC8AG4ApJZ91EUEY6/PAhcPZfj2Hqw7FwoUipvIGltGNui0VGdk97/kBYakn7PX4wGSzg++u9RR3rNitQ1k1FaCvTyQe9Yq/hSAxtj85NkK2IiHbVDJcVUQeZHD/AED7IK8XYPYQn7UCWDHSOdU7gUVJ2Jb1QmeriesZLZUAl3k7goGHHxbDe90faMTfsm3vGs2snx0iltL9oaxks+j+rppyYFlw+scNbr898DME1YN9JPQ75DFJ37EXVk9PEANfBXQRvRpgWdC7xuSSGX8hbnpBAbi46jdXtABJ1hLHxSCPic9+739vwPNvbMmy6MibcD65uC7NNg7HsOjxU+txapoGZOzYgY/DegkUWnkul4vpufTxqv0WdD+GGtOyaQ09PfRc8iizWYKHUUztKUQCQfDDiJHnQrWS1Nn/gzXv+NrupELVXPfooPieihOLKdasabFjpyk7jgPHH+fj+f1VhBUSag2Bv8HAY+kJN9Ox9l/Lss2d9WPQGkON2aFkd6uwfj0Z9sMP8xY8zDBCqeCEI8xniNMtX6D99jn0nHATiy7dVLtVehjKmp56PVniiKiiRUmjGQmWBT9eamAon4K+EEQqGcXHPiaz7/8Mzldvr/h8tGDQ1IwEr4f/Df5DKEQoCSZZYszOd3H1sgOZfkKeSz/6eG0iuSwLJk9maLeaRrNndHTA0n97CaNoQir8BvqFRR6T9OIt0NNT0blowaCpCc6ky/giP8IvnoLb5zjmS/4PFGoPAw+TbN6k+5GJTL9kApdOe6bquSBhcjRKBfdV6p2iaQE6vnsoS+dtwiwRDuYA4WCwgs/jXHZTRc95LRg0NSG9vo0scfpPwf5VkcJnX94sJvkAxegtip3YzMDMdP9htC/7FE6ys2rCIbQo8Z3v0FKRZJrq0PHdQ5k5Y2uxLIZXrAgAhTIZ872v48z9RcXOeS0Y6pGWSaLqP9lPOOodIiqPIo9g8Cb7DBAEYGDgcQTPYpKnVEC4REnnzmXR/GzVdlnL+oY0VSE1b39GRaXYkrb/Uu3hY3If7bQ/srBiCyItGOqNFim0t7055pof/x33L3uOk401GGE7TwOPqeoxYiqHSY4oeZ7nH0PzU78JysBnhXyeq+89nvbpXu12WU8PnHpqxe2/mubHsqB3bYRLZmwlqvIU/G3BUsrHJ8I24qRz51ak5lrrCoZ6XZW3SKG9wcwxVkcb85eOJR4VTMMnHhW6lo3GfijGglmvcuGHHsYjUqwnUzA7TWI9HhE8MXGzPnZ6U/W/UE8PXHIJrF4d3GvhoBkhlgVLbxvLxWe/XjQrSXFRJAgGK+TzOIkzy77t1hQMjoOT7GTR196pqm16WISF9hzjWBapf6/IQa8XBjPHWB1t9K6NsuDbBr1ro1gdbcHrlo4ndcM0YnEVhvSBQZ5RZEmqtWGDnzwRPF5e+1L1ixPeeuvOH2s0e0hq3v6MMr0Ss1K/CdYzYtiZChSbHKpWRqPc9qRWUt+slTKad8UsNJevUOvIPaWv+0kZHXXFVF5QQ6n7yVpPqW4o1C/q7hZZOOsl6Z7xGxkdyYpBTkxcibCt/7jGplWvnlF3t4TFbYJbd3d1tqtpCfq6n5RZZrdE2Taw9a2x5zXW2EmtpJbUGGym4RLDCxvCzL/jqLqqbmpn2nC9CJ4YuDmwZ/+yvrSaGlLQMjo6Ai0iM/UMXC8S9nQwApNSodFP7tjqmeI6OqC7G045Bbq76aFDuxs0ZcPqaGPpA21c/KHfF81KCp+LWIGFdj6XhWRqPDEzMD34mNwnJ9J+2WF1c+1NJiFm5jHJESNH0v9d0/oaRkqhIKFJjgjegAYnb6m/Y9HL51XvuHZ0wD330EOHdjdoyo9lkbphGqNMr2g2nSSPa+dzubAs6P3yKk6it1jh0/VUbZyWg2BZ0Hv9cyyILKBLXYltnNjUvoaRUChIuGDWq1w89SkMFeRCGPj8wJ/L1d0H0T4tV1WNULsbNJXCsqDry69g4ONhMFd+WJFrw4jLbjckjgNdXfwD/0KEaXhIsDJnLVAfaayF7mXtcw7H9SLE5ip623Tc/GAE5bHH4zjjWTkth5vzMBDyGIgYZHM+9uxfYrW9U5UdeM45gbZQ+lijKReZvQ9FDMH3Fa5hYmfaKPdZ3ZKCwUlvpN1dhUuMCHnO4i7GmhmYdEytpzYAO9OG64Pn90euasEwNAVNy579S97K/w2LmUexb673WtV2YKFk/q23BkKhSiX0NS1CMgmxuMJ1IRarTLn+lhQMpc5nAX7NmYhnsvKL1NWqvNAiOjgBaKl+DXuK1dGG1fYOi+a+hvFIfyOgjDGmqjuwo0MLBE1lCDsEV7Rcf0v6GArOZ5McJn4YzVLD5KghKJwACxYE9/UisOoeyyI5+W3iuMWM6ZcPOhYnvVFHd2magkqXZGlJwWBZ0PvjPzLTvJGP82si5DECizSJLU/XenoD0DV59gwrNYHe2BnMVDcgwPKXTq56sb0dqNdse41mO1rSlARAWxsrI0fgemDgoZDAy/+bU2hz9IW44bEsLHsR9vws+dVRfEyygJ07FqsWzpqw+xu5HESj2mGkqWtaUmOAsCRRLjAhBV2HjSBsNSd1ZU7SjADLInFOMiy6FzqhjTdq46xJpwNnkUhwr9u+aeqYlhUMhSSywITkh+IhVxK2qmkGMhkwFAS5DR4Z9qv1lDSauqdlBYNlQdeVL2PgIygUiqN5jK7IV7BSE2o9PU2ZSCYhHsn3O6G9AwIndLVJpSAex8Fikfl1nEmXVX8OGs0waVnBAGGiiArq7OSI8igfYS5dOGVPFykz2ok5bAq5DTPNGwMnNF+gffm51a+NZVk41z3KdPN+vuZ/i+lfbNOHT1O3tLRgSCYhFukvZyuYuPn6KY0xKC3SyKecWB1tjJt5Kp6K4hXKn9SgMGF6fRtZL4KIIpvVbgZN/dLSgqGwmjxb3YkZhqzWvY+hRRr5lJtkajymIWEHXU8XJtRodkJLCwYA2tq4J3IWgsLEr38fQyEd2jR1OvRuoozgdFcAkUjV99327UxT9VGWS6PZgdbNYwixbXA9Ex+F4LP+zG+ANbbW0xqaauTDNyG2DXnfQIA8gn3GYixralXnUGhnqg+dpt5pecGQTGwg4k/AI46gWHHXfqTqPcHNsnCwggsMdT7XOqHgT3I9P2j/edcfcHpGF6vYVougEmxVN6nR7DYtb0qyMr/mQm4sdkXyvPpPcNP+593HsqD3wp/2l8jwLqR9zuF632k0g9DygoFkkknGHzDxGsP5jPY/7ylWagLjzFeL7T9dz9T7TqMZhLIIBqXUCqXUX5RST5WM7auUulcptTG83yccV0qp65RSzyulnlRKTS55zwXh6zcqpS4ox9x2hYPFXPN6fIzGcD6znf854pF8Oa3VhuFgWSS/NJlYWHU15m8jmaiDXt86L0VTZ5RLY/gv4LTtxq4CekVkAtAbPgY4HZgQ3jqApRAIEuCbwEeBqcA3C8KkkvQ7nyPkiYTO5/o2AhfLcc/cRK+0Yy2/SNuUhom197P0GqewgG/Qa5yClfl1bSdUsAt+/etwwgm6QbSmLiiLYBCR+4E3ths+G1gZ/r8SmFEynpaA3wN7K6X2B04F7hWRN0TkTeBedhQ2ZSdwPmcBCZzPvxnTENdXy4LOcT/D8h7UNqXdwEmciW2cSNJ4ACv+eO3DfW0bslnwfcjnYc4cLeA1NaeSPoYPisifAcL7D4TjBwCvlLxuczg21PgOKKU6lFLrlFLrtm7dOqJJBs7n/+p3Pucb6Pqqcxp2C8eB9i8ewdXeN2lXv8Pperj22mEyCUbJz9DzGugE1DQrtXA+q0HGZCfjOw6K9IjIFBGZMmbMmJHNJpkkFf05o8gGdmfJ1ofdeTjoFm+7hZ3ehJsVPDGDshjr96r1lIJjtmRJ0KPBMCAe1wJeU3MqmcfwmlJqfxH5c2gq+ks4vhk4qOR1BwKvhuPJ7cbtCs4PCJzP9qQr6XpkLhkSgYkh83GguvHte4wOjB82SdYS41O4SEn0WR2kH3d0BOfhrRmS5ySwrAY59zRNSyUFw53ABcA14f0dJeNzlFI3Ezia/ycUHvcAC0sczqcAnRWcX9Hv52bPJsYpgTMy/jgkr63kZsuO4+hs2uFgpSbQu+IM0rlzwTBh0jG1nhIQnodz28hmwfgdLAE6Omo9K01LIyIjvgE/B/4M5AhW/hcDCYJopI3h/b7haxXBuf8CsAGYUvI5FwHPh7cLh7Pto48+WvaUhQtFTFMEREyVl4VTbxPp69vjz6sFfX0io0cH32P06IabftXp635SRkddMQ2/bvbXwoUihhGchyASjdbHvDTNDbBOhriulkVjEJHPDvFU+yCvFWD2EJ+zAlhRjjkNh0KZhKwnGOKTeHw18Llqbb4sDJbsprWGobEzbbie4PkKNyvYtqr5/ir4n30/eFzwP9d6XprWpaUzny0Luk6/BwMfD4O5+e/VprvXCNCBSbtHMrGBmP9emOD2Xl0EGmj/s6beaPkiepmxRyIY+ERwEWym1Xv/tgHoYqu7yfr1XKAeAhFSxk/rJtCgowPa2vRx1NQHLS8YkpPeJsaY/kiVSW/Xekq7jQ5MGh6OA+03no8rQgyXVOQXdbU018dRUy+0vGCo1xWkpvzYNrh5Ew9wlcK+aCWWNb7W09Jo6o6WFgz1voLUlJdiTwY/9MektFDQaAajpZ3Ptg1uzghKMKs49kUrtS7fxFg49Ep7UEBP2rGo45pEuuKqpoa0tMYQRKgcikuUmDSmf0GzG9g2lvcglqwFz6zfmNBC5mU2G4QpLVmiM940VaWlBYOV+TVd6iVulU9yjroNK3MwDe9f0GnQQxPG9jrZydjqRJKJM+szAq204qrvBxVX29r08dRUjZYWDE7iTObKobjEeEBOoC3xQn1eKIZLscaHGxjRdWG9gVgWTtfDtM85HNeLEJur6K3H663OeNPUmNb2MWTacNWowMdgjMLONLi2oHt+7hI704brR4PM53rdRWHGm2MexyLViRM5XgdFaKpKy2oMjgMvP7IFU/YGTGJ+jmTiBRralNQoppIakkxsIGYcjisRYjFVt9dbp62D9sjFZLMKw4MlGww69MHUVImWFAxFi8u2DxDBZSY/aY4chkYxldQKx8Ga206vNzno4tb16botcW3bkM2Z+IDvaTeDprq0pCmpaHERgzwm49Tm+mjzWAYawlRSK8IDb/kP0SkLa9/veSfoxm6aWtKSgqGQ6GQqj1hUkbzksKZx1OqiejuhgXaOLqynqSUtaUoKEp06sTmWpHoIK7WoKYQC6KJ6O8PBwr7gWZKsxUpNqPudowvraWpFSwoGbBvyeRA/uG+yUEBdjG1H+iN5xxOLpehN0RCOeX0sNbWgJQWDkziTdv8KXGLEfJfeRs9f0OwS3dBIoxk+LeljaLr8haHQ9XaKDPArRbzGs9frY6mpIi2pMSQTG4hJWCOpGfIXBkNnQQ9ggF9JHsJiEY1hTEIfS13mpeq0pGCwMr+m1/gNtn88SeOBxs9fGAxtOxlIoxTQG4zSY5nNwvz5wa1R5r+HOD0bSC/azJaXtjGWfZhkpll/0FvwvveRumJfrI4m+83uBhWXlSLS0Lejjz5adpu+PpHRo0VMM7jv69v9z6h3WuE77g6NvD8KczcMEQjuG+077AZ9fSKzTnhaomwT8Ae9xXlP+mZ8t2n3wc4o16kMrJMhrqst6WMoxnQuWNC8ankrfMfdoZH3Rzh3Z8rlQe0kf2rT1sLq6YFpx/t0338YOWKAKrlR/D9LlPm3T8RJdrac38VOb8Ld5le0JJoKBEfjMmXKFFm3bl2tp1G3aPNsc+A40D7dI5sVDHyWRObScf/nmuqgOg6ccLxP3isIgp1fmww84rj0zrge67Z5VZljzXEcnGQn7e6qwEcaN+hdY+7RaaCUekxEpgz2XEv6GFqFVvdZltLoAnJA7SSEOXI9bRiN4j4fFunFW8h7YwiCJQUDj8P4I/H3RXGNOPvF/hfyOV7861j+O/8BfCJsQ5G+c28sx2nMA7u7hL6yXtqDQpkXHoZlpcq+GS0Ymhjtfw5oBgE5sEWDwhPVVMfTcWDFnQn6hULQi+I5PoT8VW33aq/4OsHgRj9Fau5VWF00zw4ZirCsi+U+ihX7A6R6K7KZ1vQxtAgNVBqoojRDm4pmr52UXryFnG8ACoWEmpFJvxu01M9Qup5V5DFJP3IYi064G6dnQzWnXX3CCsqL2u/D6Xq4YoJQ+xianEY3oZSDZtAYCjTj8XQcSB6Xw/W3N2DszM9QeE6IksNAyBMhZuTpXbqxaUNZy3kuax9DC6Nr7YRBPV0bsG/NkDwnUbc9GIZDMx7PQFsYQ3CxD9uZlvgZvsK1vM0+cPDB7PWBGD985DjyGAgmhAIBQDBxfcGe/Uustneab0dRPfOwFgya5ids0GO5LjwQg7YGVhmajIJvQUJB0G8uEkzy/FjNoeMSA1LJ4jGb0bOB+QsM7t18GEIECbWKoLuKx8ve3+OkN2I12TEudp1kPzDMinYg1D4GTfPTDE6G7WmS2kl2ehOeXxqe6hP4GXxmcgMdy46GpUsHCHKro435txyJaRTeowCP/XkVwWC5fIH25ec2lb+hEK68/Pb9UF6emeon9HZtqNj6RgsGTfPTbF74gqH561+HE04IssIaEMeBlx/PhBFIhQu8gUGeUWRJzXg7aEoxCJYFS5aaRE1B4QEmr3IAOaJBcUxPYc/+ZcMLzgLB2kbhEQm6TvovVbQDoRYMmqbGcWCRHURyNGTW82DYNs62SSzy5+HkpwQNoRvsAlhcAT86ET/UEEBh4HESvfTGzsCad/xOP6OjA9Y+YHDy1Lcx8EOfA4CHQkh4rzWHdkihOrCPSY4YOZLRhyq7wBmqVkatbsBpwB+B54GrdvX6PaqV1Ir09YksXNhStWUauTzSzujrflJG866Y5GQ070qf+lhwbBuIhbNeEpOcgIjCExNXTNzg+0y9YrcOVl+fyOh4XkyVF5OsGORFkQ8+q1nqKfX1SV9smixUndJnHifS3T3ij6RRaiUppUxgCXA6cATwWaXUEbWdVRNQMD1cfXVw32Cryz2lGV0LMLCfyDbipI3PN5x5LLHlaQrhplLiU+iNnYHV9S+7pdVZFvSuMZn5kT8ACj+MWMoSw77jf5rjnA8znjtlERYOZDIV3VxdCQZgKvC8iLwoIi5wM3B2jefU+DTrFXIXNJtroUAyCWa0JPPXuBingYpjOA7M/c0peEXTj0IwGGe8gvWj8/bI1GdZMG5yAqHfka0QXpYDcbKTG/+cr/LJXG+C4QDglZLHm8OxASilOpRS65RS67Zu3Vq1yTUsySSOeRyL1L/jmMc1zxVyFxTyFxa02xWN4Kg2lgUXXQRKBRfBvG801HXPTm8imwsczQFCBI8ka0e0Ek6mxhOPB34KAw9B0cNM2v3VOIkzyzL3WuFgseiCZ3FmrqiOn2woG1MtbsCngZ+UPP4c8KOdvUf7GHZNvw3Wk9HxfFOYXIdFszoZpHG/Wl+fyKypj4mJG/ZWEG0ZWIMAACAASURBVIG8zGJpWb7IYJ9vkJeFp6xpnJ20HZU61jSKj4FAQzio5PGBwKs1mkvTYNvg5k08MXDzZkOtLkdEE5vQiu0lZm6i94J0YHeucwqurp5HJyEUVvZ5RhdCU8uwEh7cpOTz8urnGrZ3Qy1O43oTDI8CE5RShyilYsC5wJ01nlPD06y29l3S5F/cwqFz5Yewei5siHwG24bsNsEXhY+JQuhgOb3GKVhTvbKZRwY3KX2BdncVTnpjWbZRTRKJgkvdIxbxqnIa15VgEJE8MAe4B3gWuEVEnq7trBqfRm5eNiKa/Ys3WD5DMrEBU3IUktkExTi1GSv+eFmFdiFKqWPqH8JKrSY+kSBK6fG96nofbY/jwNwvevieYOLR5V1eFe2w7molicgqYFWt59FsNGPxtWHRxF/cSZxJu1xBlljQ1S03h446bdLgOGDfmuFKfsAP+BI+BnFckmf/Hcwrv9C2LLAnJ5BH+k1KJj7Jdd+D9scbZqFQyHj2MVD4ZLx9qtJYpe4Eg0ajGR52po0sPn6Y9zuH62lL/LHuAlcLvoXstmmYfIwv8X325m2S6n6sqWdW7CKXTI0nfkOebC6Pgc+VfB/bPx6yYNWpAN2eZGIDMeMwXE+VZDwvqvh2tWBoNZqxoH+LkkyCETHw84FpxjOi2Jm2uhMMA30LUX7Al7lfTccatR6S36vYdi0Lei/+Gellf2ULH+A/uTLo2eC79CZeqLv9tANhVeBemYxtTid51t9izVtUld+tFgytRDN1rNEUu7rNmaPwPIjHK1eGeSQk3noB5GACl2bgeLZP/g7W/Hjlz79Jk1jJoWwjHpb2VriAvX6v+hcMYTiS5T+EZf4epi6o2u+1rpzPmgrTxOGb29MkVal3SUcHrL1+A98+qT6T+BwH5v5wXLFngiJPnCzJcxJVucjZmTZcY3SxwJ7CI4ZL8oZU/Z8cNYyq0xpDK1E40QoaQz0uL8tASylGdd6EKMhyPgjBCCun3sd8tQArcyZQ+U56ySTE4gp3m4eSPJNZT5I12Lljoc6b+ThY2Bc8S5K1WKkJVT2uWjC0EmEj8WZocbkzqtX+sC6wbZzsZGz/eJLvrcVKp+vnyzoOiZ/chM8SQPAxOIf/W3HfQimFiOX04r9w4+37sI4pPMJHMfCI36joTdXP7iqlpwdmzwbfH088ngrmWcXta1NSC+E40D63jat7k7TPbat7TXpPafK8tgE4iTNp91fzdRZwAjY9y1VdmEgcBxbNz7I+/08YeBR6LWROOa/qKpxlwbip+5M34vjhWtgnEiwa0puqNo/h4jhBSko+D74P2Wz1rb5aMLQQreJiaPa8tlKCkNXggpcnyhyvq+bZvcUq7/dN40Y+H/YcywV5C1XyLWxPwaRkqH5fhxKPxPLv1oUgLcW2g99oAcOo/uJGC4YWomVW0o6DZS+iM+k0tVCAMGTV7O997GFgM62mcyouQHyFS4wzWMUCvhGUvqhgO8qdUVgsfPvktZzPTSgUHiZzve/VXJBuTzIJ8WhQziNq5FnypReqfh5rwdBCtMRKusWaElkWLPmxQdQUDOUHIaup8TWdUzIJESPo4ywoVnEGSeOBspe+2F0sC5LnJPgF5w5s5lNnZTIsHHr9E/k2X2etfzwdXUdUfX7a+dxiNHGFiICW8jwHdHRAG0/XTVCBhcOF/lN0cxGCiaei2Cd9uzp5C7vAzrThh/2li2UyHr0W2tfXz2rJtrHyD2Bxf/A4p6p+HmuNoRVp5iD/lrGXleA4cPnlcO9quOyymlZZLTidJ3mPMops0Lxeqpe3sCsKZhpFHgOPz3AztpxQX13eCudwgWi06uex1hhajSYO8g+qfVgkux4ObNktUvbDSW+k3V0VFNPzfJZc+kU62pyqf/fiqZWdRoypdHEFGfYLaiJVKW9hV1gWdF35MnMWH0Qeg5/yuSB0tZ7KZFgWrFkD6XTwOFX9mFotGFqNJjW1DJR3bfT2tjXD1xoWNtPIEsMnEhTT86+jLX1z1ZO37PQmstsOwheDrBpFRn2ATq6BeLxqeQvDIbP3ofiGIH5gTvKJ4CpVX2Uyamzz1aakVqNJ+z+3SijuYCRT4zGMGkcmOQ6Jn3wXX8KLrSgSX7kIvv3tutNKk0mIRX0UeYCgAY5kSa64oObm1Xqx8mqNocVwsGhXvbgoYkroxayfVdIIaJFqH4NiWbBkqcmcSz08H+JRqhqZ5Dhgz8/ycpjM5hMJktn2PhQ6O6s2j+FiWdB1+j3Muf0k8gQtg07lt0FGWQ016P5s50DJqqU81RpDi9Gs/Z9bIhR3J3S0OVxvXM4UHuXU/F1w++1V2e7QyWxZkokNVZnDnpAZe2QYshrBI8IdnE27vxrnEbMmy/V6yHYuRWsMLUYzr6ybPhR3JzjpjVye/wEucRC4e3GWNYduwOqorMO3NJnNJ8bZ3M5UHg3yFjIfpx4czoORTI0ndkOebTkPCXMaXKLYd/wP1j3tVV9d2DZ4eZ9grS4YRm1LqGuNocVo2pV1vRhna4TNNHJECeLzVXCRuzVT8e3WazLbrrAs6Pr4agz8cEQw8HlZDqxJ6GryrduJyzYM8kTJseRf1tb0t6k1hhbEsoIkpODkTza+dGjiENzhkkyNJ7o8h+sFa70YuSB3oMLskMxmxOommW1XZMYeSUGQgoePwXJmstK/gN63VlfV92Y9sZRevotNkiQ21ta9oIalTbTG0Io0W9mIVg5JCrFw+JG6gg/xNEfwNNed8uuKm5EGTWYzPZLzk3UvFCA0J8UVpvIw8fEx8YgE2tYPHq/u7+Kcc7D4PZ1cg8Xv4ZxzqrftQdAaQyvSbLkMzew4GSYDfAzAF1f/I209lfMxFJW0bScMTGa78qNY1oyKbLPcWBb0rjFJz32cGx45Eo8IhJ3mEt5r1f1ddHQE97feGgiFwuMaoTWGViS8kDrGsUE+Q+LMWs9oZDSt42T4VNvHECSz+XhikCVGhv3oNBZj7f1sxbZZCSwLxk1O4GMS7DvwMZkrP8R560MV377jwKWXBjenrQPuuafmQgG0xtCahJ3c2uccjutFiM1V9LY1+PW0lUOSqLKPodCZTQqd2UwS6o0wwzlZmW1WkGRqPObyPJ4XFNcTjKI5yZrxwYqdV44D06cHoakAK1bUj/KuNYYWxc604frRoGZ+a5rlm4pq+RgG78zmkzn5sw2trQW6QkkTH6TfnFQhAouuFB/nclI3v0OtMbQo2izfXFTDx1BaJC/CMUTI4yHE4kbocC7bpqqKbUNeAlNSIBQUPgZz5Ye0VTA6KZnYQEwmkA2PWVSyJBMbqYfcD60xtCiWBb1dG1jQbtPbtaFRF3qakGr4GIp+BV+RI8JF3MgC8z/ove6Zhj5/+tt++hQqTvlViE6yMr9mDScyi2XMYhk2J9asw932aI2hVXEcrLntWK4LD8SgrXHNAJrQTt6TJ+8HpokI+fL6GAbxK0zicTq4ETJ/Qz2scvcUy4KuLphzqZCXINVtB3NSJX4bySRW/D+wsqHgicUg+f3yb2cP0IKhVbFtnOxkbP94ktkHsOrF67UbOGGOXou0Xdg5Gzag/P8DmOGAKttHD14kL09GfaBp7JCZDGHtpMqbk/rPWwurxn0XhkILhhbFSZxJu38FLjFi9dSkZJjoZOeB2LdmyBOhIBDymNi3ZrBGGPk4pF8hqkhefBikmmPHF0pxZ7OBxiWAEGEbivT3t2LNKE/jI8eB9ukeWVdhGnD9jy06ltbf/tM+hhbFzrThGqODTE9jNHamsUwBOtl5IIkxKqz7E1zSDHwSRx004s8d0q9w/bNYS+tnhTtSCqW4TfziHgzqPxnc6H0OZ/EDZdmOnd5ENgu+GOQ8xezL/LosPDAiwaCU+rRS6mmllK+UmrLdc51KqeeVUn9USp1aMn5aOPa8UuqqkvFDlFIPK6U2KqV+oZSKoakYBYebaQb3jWYNaMXWzkPh9Gxg7k8/gqAw8DDwEBRz//PgPb7oOA4sunQTby2/pb/5TuhX6OSaunGSlpPSUtyF0FUIBKJ91/+WxQmdZC0mHoWmSr5XPyGqpYxUY3gK+Gfg/tJBpdQRwLnAkcBpwI+VUqZSygSWAKcDRwCfDV8L8F3ghyIyAXgTuHiEc9PshGKy8MxN9F6QDorqNRA62bkf+9YMbtjaU1D4KHwie1zTv1hKq/sgfuhdgSrmK3hN5VfYnmRqPKYR6Av9PpowgW+EOQ2F4r/stRfXM5sIOQzyxKNSl7tyRD4GEXkWQKkdHF1nAzeLSBb4k1LqeWBq+NzzIvJi+L6bgbOVUs8CJwLnha9ZCcwHlo5kfpqdY+FgrQwN9SsbzFDvOFi2jaU9zySPegtzdR4/vHh7RAkuaAaJPQhMstObcLcdhCcGghGaV3JN51fYHsuCiz6Rofv2MQgmBQGh8FjPJEiYu/qIQSn4FVxXETPOoFddy/0yDVtNJ3nx4VhWqqzfoxxUysdwAPBKyePN4dhQ4wngLRHJbzc+KEqpDqXUOqXUuq1bt5Z14q2Ek97Iom1X4ngfaSxDfbNVhx0pb79dEoOkgKC0g0LI7EYqQ6n5yJB8sKLFZYnxRRbMWE/v2khT+RUGIzVvLKNMr9gPuuhn4EKcu9/ao8+005twsxJ0TfQUtjEdy3yUzlFdWKkJ5Zt8GdmlxqCUug8YO8hTXxORO4Z62yBjwuCCqFRv2358UESkB+gBmDJlypCv0wyN40D7DefhCsT4Gr3mGcHquxFotuqwIyRIbosE/RCKP7Ogcc5wNYaCrM1uOwhfvoIiKEXdxRV0qBth6niwpu76gxocy4KuL7/CnMUHkUMR7MsgYTB9595Yzu5HJyVZS4TP4KOI4JE8ay+YuqCu46x3KRhE5KQ9+NzNQGlIxIHAq+H/g42/DuytlIqEWkPp6zUVwE5vws0dENafF+wzFmM1wA/fccB++TyS5j1YPNi09u7dIbFXPqwOKlC8N1D4ZDK7zmdwHJg/H7LbBF8KQiWCkCOjxrTcPs68HcHDoH9fhlqDnyK1+Dqs24Z5MXccnPRG0s8cjVdajen006HCvTJGSqVMSXcC5yql4kqpQ4AJwCPAo8CEMAIpRuCgvlNEBFgDfCp8/wXAUNqIpgwEq5g8Ci9YxYx9rtZT2iVFC9Ly8bSrXpyZKxrLL1IhMk+8UnQQU7wfnsZQ2Kf33SsESdNBq0uDfOBTuOTwltvHiS1PlwhaKBg1XKKkb98Lenp2/SGOg5PspH3Zp+i+/3DyxEKNzsRev1cFZ18eRuR8Vkp9EvgRMAb4jVLqCRE5VUSeVkrdAjwD5IHZIuKF75kD3EMgjleIyNPhx30VuFkp9W1gPXDDSOam2QWTJg047Zk0qYaTGR4DLEiY2ONSrXS9GpLEGFXiLB2exlDIvn35kS242z5Q1BQAIrh8Qf0Xqes/htVRf47RSpMZe2Qxwxs8DAQ/LMe9nIuZdOnldLQNbVIqZornzsUNBQIICi8oh85aoL7360ijkm4Dbhviue8A3xlkfBWwapDxF+mPXNJUGDvThmcI4is8w8DOtNV95rOuCDs4mY1vovDC+Ps8BeEgKN56ejOBZbafnh6YMwe8vBCRvYng4hMNk7qCSKRxbMLKvEEj10DaU5Kp8UR6crh+wV9TyDsw8Igwx7+OtsVdg5qU+iOQjseUYzDw8UOt/GJWkIrdjJVaVOVvtPvozOcWpRET3HTuwiA4DonHVg+iMQSawrU//Xu++tX+TmGf/CRcdqlPLif4osgTZSJPcDZ3EFNu0LeZHMloX8tK3kLYqgqju/wwwquwX3OYzL39BJxpVw2IiHN6NjD/M0+TzQqemOSJ4IfHQZkGqVnvw7IXNcSJq2sltSiFstvpG3Lw93/P4IFn9UeLN2rbEdsm4+8zQGMwILwgBZezxYuF73/Px/ML60BFwQ/hY7COKcRx+ZExl8xZF5Ic+1ywqm3hHZ2aN5YVd+VxvX6tIdjHJmDwCB/l+PuP5svH/oC3x7/Flm37cPeWo8gRwccsNjEKipSYeL7XUKZPLRhaFceByztZ6a7CJcbKuz1615gNc+JqQpJJEkYa8fo1hs9yEz/j/NA0FAgAzy8IAyh0HCisgv0wMi3j70vn1F7o7KzNd6kjLAsumhlh2TKfwLxmFi/zBUOLR4TFMg9eKn2nKuZABHvXDBz5scbQygtoU1Krkk5jux/DJRaErNZzfluhnkCrJ7INhmWx/qBPhA+CC//G903iPH4ajgUO5f76PIX2lYJJnhj5EvPRQy1rPhqMVAoiJhRMSBLGFVFiVgroF7pBMZIgskuIYOBx0tT/bbhFl9YYWhHHgRUrSDKZGC4uQixm1Oc1QdfX3jXve9+Ah4/89UgeU0dwivyWFzmUj/J73uX93Mkn8EMtQjARfC7kBsZNHUty8tstbz7aHsuCL33ZYPHigiAwOZebeJf3cxefKEYqlebiGnhcqX7Ij+Ty4HcVN5jftU/D7VYtGFoQJ70RO/cVkvyOXk7CnjqPZNeM+jx5dZbzLkldsS83XpIlSwzCbF1PhNWchsLneSaE691+30IhdDIVuxmraxFYM2r7JeqUvffu32Mg/Jzz+ArXMpPlbGEsb+x3GM9mxrBV9qXQJGnvjxxG7+RfYTONZGp8Y56uItLQt6OPPlo0w6evT2R0PC8mORnNu9IXmxYM1it9fSKjR4uYZnBfz3OtIX0nfFWm4gj4AhLe/B3uFXmJ8Z7MMpZJ34zv6v25C/r6RCKR7felL4pc8f/CzSj8prqfrPW0hwWwToa4rmofQ7WoEzu5bYObNwO/gopjX7SyvlfgOkZ1WFjXnE2X8RWiuAz0JQSOUIM8MbJccsSD2LNuYemD/4R12zy9P3eBZcGSJWAYgcZQ8CVI0dhS0MJ8Ttr3CXq7X8Cq83IXw0GbkqqA07MBe/bdJP3fYcVre4ErJollhZjpk5z0dk3msVvoGNVdY1lYS1Osvewk0l5QvX4S68mQIMHrZNiPZLQP6yfXQB2Wea5nOsL2qHPmKPJ5QSQQuFJSNiNOjvmL4k0hFABUoFE0LlOmTJF169bVehpD4jjQfryL6xnEcOlVJ2Od/P6galktLnaOQ89VL3DrAx/kHH5Fx6ib9Eq8mSjUukgkgg73iQSsXx88V0fN5huRAbt2/SYSW55m/RsHw7ZtpC6ONpxQUEo9JiJTBntOawwVxk5vIusdgE+EbcB8uZr59y7AeqC9+hdkx6HnhJuYk+/Cw+ABjqUt+zRWHTp0Cz/COq5MXJ9o7apiDNy148Nbc6J9DBWmtFKjYHIvJ9Mu9+K8dxSk01Wdi5PeyOx8FzmiQetHYtjGiXUXu6778Gg0tUVrDBVm/RsHo/CL9kjBxCWKzTSsFT8IqppmMlVZGttMK8axg2AakFzyabDqSwXWEaoaTW3RGkMFcRy48cEJ2yXBBNmTSWwc92gWzdqE8/XfwPTpQZWzSi2PHSfowWAGMeym8rl+aaQu7aIFB7lp6iqqGk0t0BpDBbHTm8j5B1BYoRcaqZzBKhbzb9zFWYgo4uLSlb2CzLJ9SK7oLHsFxkJUVCK/BcVnAIOI5GjjGeqxrHIhQlX7GDSa2qAFQwVJPPMAPuczsBOUwV2chVcSB70NmM0SBIOY69Kb/hVWma6GjgPtcw7HzX8DhY8fFgTL42PfmsHqKMtmyo72oWo0tUObkiqF47D+wb8Wa7oHBHXdPUwYUJ4A/LAJSMH/UC5sG1wvgkcEHwMTv79o2jnD7BRfLeokCVCjaXW0xlAhnPRGbvRTO/gXIuTwlYkngbAw8DiOB/k9x5AHDKVIbHkanFeDJfMI4jadng08clMM5B8xgLjp03WuQ2arkDwnUV/+BV0sT6OpG7RgqBA208gTYXv/wsfV3dxtnomfl7AePjzIcZhhqd68KC6//UTafnMq1sf3xfl1Bts7nmR0mL6HUJA4b32I5OLTcIkDYJKnS82lY/bn6vOCq0ORNJq6QQuGCpFMjcdcnsfzPChqDQaMHUt+S9BXV4qaRGD9L5iXXAwW565g7O2vsYIL8Yjs6HsYRJMoLb1hyzvkOIuCGcvHIOPtU5cXXMcB++XzSJr3YPGgDkXSaGqMFgwVRJX8LXR2Gvvak8QiU9iWkwEdtpRSlFYnCRzUQXN2UGSB+fdPZ/6laaxJ23Au/xl27lgS5k1kzjRJ8Dpz7zoJ1/sGJv/OGazCJB82eqSkEUt9NSLvtyCNJxbppXfmT7FSE+pOeGk0rYQWDBXCtiHvFy7sQbVLwWSSv45UxyjSW05hxV37kRcT01RcafXxn/dPxiUWOqgV/Y3dBR+Te5/Znwee+RRdxpeY668iSww/b2Lc7mHik8dAiOBhcgdnEyPHmdzG2IPfR+q0v9RlI5YBFiTMhuqLq9E0K1owVIhkEkzDx/OgYCJSeKxnEh17PY21NEWq1BqE4tATvsTssI5Rv1Ao+CgC13UWuMH/PNuIF7OpfSLIdtUeCyGpUyNP0Pmz08E6rcp7YHgkExuIGYfjSqTh+uJqNM2KFgwVRBkGeP0Xd8HgRi4k9f2TsWY4WJZVsjq2yHzh75FlZthA3Cu+Z2DnLWE9k3eIdjIQ/JLSGwqfWBSS19dfyYsijoM1t51ebzK2cSLJrk9j1etcNZoWQucxVAjbhrxXag6CwLEcDerl2/YO70mmxhOLg0mOOFmWml/kkiMeIGbkMcgTJcdZ3FVS7yj4TIXPWdzFKLLF914y4zV619Z5KeDQjmT5D9EpC7Eyv671jDQaDVpjqBjJJJjKK0lmgwFaQ2Ij25vSLQt615jY6c0kWYuVCkJLUz0bsGf/kqT/OzBN7smdRhbwMTHIE8dlXuw65s01sJ/Yu/5yFIai2DXI1ZFIGk0doRv1VJBLP7mF7tvHhOad/nBUkxwLZr1K59LdqOdeEp7qbHg/9q0ZEkcdRObtSChEGjSSRzde0Ghqws4a9WjBUEEcB5LH53E9szimEGK4rJl1C9ZS3WJRoykll8uxefNmtm3bVuupNA2jRo3iwAMPJBqNDhjXHdxqhGXBRTMjdHcLIgUHdOhxmDSptpPTaOqQzZs387d/+7ccfPDBKKV2/QbNThERMpkMmzdv5pBDDhn2+7TzucKkUjBqlKL/HDfwjDh2pgF8AGVG18jT7Ipt27aRSCS0UCgTSikSicRua2BaY6gwhd4C6TTceKMin6cl4/V1jTzNcNFCobzsyf7UgqEKFHoLpFJB854ka7GYADvEJTUvukaeRtM4jMiUpJS6Vin1nFLqSaXUbUqpvUue61RKPa+U+qNS6tSS8dPCseeVUleVjB+ilHpYKbVRKfULpVRsJHOrRywckisuwO7+I06ys6VsKrpdp6ZRME2To446ig9/+MOcddZZvPXWWwC89NJLfPjDHy6+bvny5UyePJk333yzOPa9730PpRSvv/561eddTkbqY7gX+LCI/BPw/4BOAKXUEcC5wJHAacCPlVKmUsoElgCnA0cAnw1fC/Bd4IciMgF4E7h4hHOrO5z0RtrdVVwt/0G7uwonvbHWU6o8oWPBwqG3FxYs0GYkTZkps/Nq9OjRPPHEEzz11FPsu+++LFmyZIfX3HTTTfzoRz9i9erV7LPPPgC88sor3HvvvYwbN64s86glIxIMIrJaRPLhw98DB4b/nw3cLCJZEfkT8DwwNbw9LyIviogL3AycrQIj2InAr8L3rwRmjGRu9YjNNFxiFenUVo84PRtYdMLdOF//DbS3Y+HQ2amFgqaMFJxXV18d3JdZC7csi//+7/8eMHbLLbdwzTXXsHr1avbbb7/i+JVXXsnixYubwkdSzqiki4C7w/8PAF4peW5zODbUeAJ4q0TIFMYHRSnVoZRap5Rat3Xr1jJNv/IEJS8UpvKIxQ2Sqd1IcGswCr2mr85/g3Z/NU528qBlQDSaETGY86pMeJ5Hb28vn/jEJ4pjmzZtYs6cOaxevZqxY8cWx++8804OOOAAJk6cWLbt15JdCgal1H1KqacGuZ1d8pqvAXngp4WhQT5K9mB8UESkR0SmiMiUMWPG7Oor1A2FkhcLvmPSe90zQUe2JvUzlPaadoliGydqx4Km/FTAefXee+9x1FFHkUgkeOONNzj55JOLz40ZM4Zx48Zxyy23FMf++te/8p3vfIdvfetbI952vbBLwSAiJ4nIhwe53QGglLoAOBM4X/rTqDcDB5V8zIHAqzsZfx3YWykV2W686bAs6Ez0wGWXsehr/4tzwlfh0kubS0A4DsmX08SiPqYhQZXXJZ/WNiRN+SnEg5fReVXwMWzatAnXdQf4GN73vvdx9913s2zZMn7602Ad/MILL/CnP/2JiRMncvDBB7N582YmT57Mli1bRjyXmiEie3wjcCw/A4zZbvxI4A9AHDgEeJGgzGgk/P8QIBa+5sjwPb8Ezg3/XwZcNpw5HH300dJQ9PVJt3GJRMmKQU5G8650M1MWRq6Wvu4naz27kdPXJzJ6tIhpSl9smiyc9ZL09dV6UppG4Zlnnqn1FORv/uZviv8//vjjctBBB4nruvKnP/1JjjzySBERefHFF2XcuHHy29/+dof3jx8/XrZu3Vq1+Q6HwfYrsE6GuK6O1MdwPfC3wL1KqSeUUstCYfM0cEsoNH4LzBYRTwIfwhzgHuBZ4JbwtQBfBb6klHqewOdwwwjnVpc46Y3M9q8jRxSfCNuIMYfrA1v8nMMbX3Eosfla3oN0jvuZVhQ0DcukSZOYOHEiN99884DxQw45hDvvvJOLLrqIhx9+uEazqxwjSnATkX/cyXPfAb4zyPgqYNUg4y8SRC01NTbTSvopBO6VPBEEA9fzGz/xS5fS1jQ477zzzoDHd911V/H/p556qvj/xIkTd4hYgiDfodHRmc9VJpkaT/xGj2zWC+oniYdHFBAiyieZbKzyVTtUzS7YfHUpbY2mYdGCocoUm/HY8PIjW1h++34U+kFfeNZWLGv/Wk9x2AxZ4bYYUQAADTNJREFU/6hQA0Sj0TQkjbU8bRIsCzo7ITVvbDGvYVQcUvMaRyhARUPINRpNDdEaQw0p1R6KVpcG6GjmOEExwMSWp4lFTsXF1O4EjaaJ0IKhxgywujgOTrITO3csyWhnkABXZ8LBcaB9uoebPYAYY+iKXE5m5ldJpsbX21Q1Gs0eogVDHeEsfoDp7m9xiRJzc6xZfB2M3YjNtNpfeENNxn75PFz3IDxMXISMtw+d434GVmcNJ6fRaMqJ9jHUC45D+s6/I0scwSRLnMV3/B/al32Kq5cdQPu0fO1yHEoKlSVXXEAs4mOSI0aOZPQhbUPSNB233XYbSimee+65Pf6Mz3/+8/zqV7/a6WsWLlw44PHHPvaxPdrW/Pnz+d73vrdH7x0MLRjqBdsGf2B5qFdlbH811pxgpzfVbm4lSWu9F/+MBbNepXfWr+rS3KVpLSrRMvbnP/85xx133A6JbeVme8HQ19dX0e0NFy0Y6oVkklTsZmJkUXjETI+Ljf8ihotBHgMhseXpXX9OmRjwY9uuUJmVmkDn0vFYS1NaKGhqSiWqbr/zzjs89NBD3HDDDUXBYNs2yWSST33qUxx++OGcf/75hfI/fOtb3+IjH/kIH/7wh+no6CiOF+jt7eWTn/xk8fG9997LP//zP3PVVVcVC/adf/75ALz//e8vvm7x4sW0tbUxceJErroq6Gm2fPlyPvKRjzBx4kTOOecc/vrXv478Cw/GULUyGuXWcLWSdkZfn/TNWtlfX6i7W7qNSyRSqKsUz0tf95MDX1PezcvCWS9J94zfyOh4XkwzKHvU11d4cqHowkeaSrK7tZIWLhQxTREI7hcuHPkcbrrpJrnoootERMSyLHnsscdkzZo1stdee8krr7winufJMcccIw888ICIiGQymeJ7//Vf/1XuvPNOERG54IIL5Je//KX4vi+HHXaY/OUvfxERkc9+9rPF15TWZSp9vGrVKrEsS959990B23j99deLr/3a174m1113nYiIfPOb35Rrr712yO+0u7WStPO5nrAsLMvq7wRtdZBZvwnpjuCLgZsT0pf9npXe+WSJYfZ4XL/UpKNjhNt1nKC73A3n4f7/9s4/tsrqjOOfp79TpkJBXfkxQJItQqpQKq4RtRkJP8qqy5CFSShWJptINlzMosEsRIyJNUuM0SAuI9DJ6o9uEpfMxAqUaVZ/FFcVhwoiyaoM6+2IJa62vffZH+9pee+l99bb+957S/t8kjf33POe95zv+7znnueec973vH3TyKGUMEoE3/uZ77OH1ozRRzpWYGlsbGTLli0ArFmzhsbGRlauXMmiRYuYPt17F9n8+fM5efIkixcv5uDBg9TX1/PVV1/R1dXFvHnzqKmpGcxPRFi3bh1PP/00dXV1tLa20tDQkFDDK6+8Ql1dHcXFxQCUlJQA3pIc999/P2fOnOHs2bMsW7YsUTYjxhzDKKeqdiYFe1zFz+mHvjBfU0CEPCIR5a5NEcrKcqjEa9yHvIOp9dy+yQtmEgrFPDexZAktPXfTqxAmD6WfHCKIhCkoyLW5ZWPUEvQKLKFQiAMHDnDkyBFEhHA4jIhQXV1NYWHhYLrc3Fz6+/vp6elh06ZNtLW1MWPGDLZt20ZPT895+dbV1VFTU0NRURGrV68mLy9x06uqQ74J7rbbbmPfvn1cffXV7N69m5Y0PVVqjmGUE1XxJ38Am/byh/DtRNwCfJFwxJuU3nUfS3r/Ri8F5D3VR91NX1D7m1LPYVR5+zyHouQIFBaJt4SFm1iu0gMUsJVelAL6eDTvHkI/s+cTjNFPkCuwNDU1UVtby86dOwfjbrzxRl577bUh0w84gSlTpnD27Fmampq45ZZbzks3depUpk6dyoMPPkhzc/NgfH5+Pn19feTn50elX7p0KQ888AC33norxcXFdHV1UVJSQnd3N6WlpfT19bF3716mTYv7osuUMMdwAXCu4pcB63j8zl9yV+QxIuRQWChUcYiWvusG72AKR5Sd+y5jz0th9tcdG9wXIQ9QIir0fq20tAiVri9e2fsW+3Oraamup+rbH1BZuw4qx+6rRw1jKBobGwcnegdYtWoVO3bsYM6cOeelnzhxInfccQdlZWXMmjWLa665Jm7ea9eupbOzk7lz5w7Gbdy4kauuuory8vLBF/8ALF++nPb2dioqKigoKKC6upqHHnqI7du3c+211zJz5kzKysro7u4O4KzPR1TjvkHzgqCiokLb2tqyLSOztMYMG/l6BT0Uom5Z71wJs/3nHVTtWu/rMeSSQ5jCfGX/ofwLZhkOY3xw9OhRrrzyymzLSAubN29mwYIFbNiwIeNlD2VXETmsqhVDpTfHMFZobaW1/lUaXryEXZH1hMmloDCH/QdzPcdR/yotf/2SyZFOQrmXU/XEaio3lmVbtWFEMVYdw8KFC5kwYQLNzc1RcxWZIlnHYENJY4XKSipfqKSytZXahudiJqHP7TvXMzCnYBiZ4vDhw9mWkBTmGMYasbe8xuyzoSJjtBPvjhxjZIxkVMiefDYMY9RQVFREKBQaUWNmnI+qEgqFKCoqSuo46zEYhjFqmD59Oh0dHXR2dmZbypihqKho8MG8b4o5BsMwRg35+fnMnj072zLGPTaUZBiGYURhjsEwDMOIwhyDYRiGEcUF/4CbiHQCI32DzRTgiwDlBIXpSg7TlRymKznGqq6ZqnrpUDsueMeQCiLSFu/Jv2xiupLDdCWH6UqO8ajLhpIMwzCMKMwxGIZhGFGMd8fwVLYFxMF0JYfpSg7TlRzjTte4nmMwDMMwzme89xgMwzCMGMwxGIZhGFGMS8cgIstF5EMROS4i9w5/RMrlzRCRgyJyVETeF5FfufhtIvKpiLS7rdp3zH1O34cisixd2kXkpIi858pvc3ElItIsIsfc5yQXLyLymCv7XREp9+Wz3qU/JiLrU9T0PZ9N2kXkSxHZkg17icguEflcRI744gKzj4gsdPY/7o79RutNx9H1iIh84Mp+QUQmuvhZIvI/n92eHK78eOeYgrbArp2IzBaRN5y2Z0WkIAVdz/o0nRSR9kzaTOK3DdmtY6o6rjYgF/gYuAIoAN4B5qa5zFKg3IUvAj4C5gLbgHuGSD/X6SoEZju9uenQDpwEpsTE1QP3uvC9wMMuXA28BAjwfeANF18CnHCfk1x4UoDX6z/AzGzYC7gBKAeOpMM+wJtApTvmJWBFCrqWAnku/LBP1yx/uph8hiw/3jmmoC2wawc8B6xx4SeBO0eqK2b/74DfZtJmxG8bslrHxmOPYRFwXFVPqGov8AxwczoLVNVTqvq2C3cDR4FpCQ65GXhGVb9W1U+A4053prTfDOxx4T3Aj3zxDerxOjBRREqBZUCzqnap6n+BZmB5QFqWAB+raqKn29NmL1X9O9A1RHkp28ftu1hVW9X7BTf48kpal6q+rKr97uvrQMK1locpP945jkhbApK6du7f7g+ApmS1JdLl8v0J0Jgoj6BtlqBtyGodG4+OYRrwb9/3DhI30oEiIrOABcAbLmqz6xLu8nU942lMh3YFXhaRwyKy0cVdrqqnwKu4wGVZ0DXAGqJ/rNm2FwRnn2kuHLQ+gNvx/h0OMFtE/ikih0Tkep/eeOXHO8dUCOLaTQbO+BxgUDa7Hjitqsd8cRm1WUzbkNU6Nh4dw1Djaxm5Z1dEvgX8Gdiiql8CO4A5wHzgFF5XNpHGdGi/TlXLgRXAXSJyQ4K0mdSFGzu+CXjeRY0GeyUiWR3psttWoB/Y66JOAd9R1QXAr4E/icjF6So/DkFdu3Rp/inRf0AyarMh2oa4SeOUH6i9xqNj6ABm+L5PBz5Ld6Eiko934feq6l8AVPW0qoZVNQL8Hq/7nEhj4NpV9TP3+TnwgtNw2nVBB7rOn2dal2MF8LaqnnYas24vR1D26SB6uCdlfW7S8YfAWjd0gBumCbnwYbyx++8OU368cxwRAV67L/CGT/Ji4keMy+vHwLM+vRmz2VBtQ4K8MlPHhpuEGGsb3lvrTuBNdA1Mas1Lc5mCN7b3aEx8qS98N95YK8A8oifkTuBNxgWqHZgAXOQL/wNvbuARoie+6l14JdETX2/quYmvT/AmvSa5cEkAdnsGqMu2vYiZiAzSPsBbLu3AxGB1CrqWA/8CLo1JdymQ68JXAJ8OV368c0xBW2DXDq8H6Z983jRSXT67HcqGzYjfNmS1jqWtMRzNG97M/kd4/wK2ZqC8xXjdt3eBdrdVA38E3nPxL8b8eLY6fR/iu4sgSO2uwr/jtvcH8sMbx90PHHOfAxVMgCdc2e8BFb68bsebODyOrzFPQVsxEAIu8cVl3F54wwungD68f18bgrQPUAEcccc8jluNYIS6juONMw/UsSdd2lXu+r4DvA3UDFd+vHNMQVtg187V2zfd+T4PFI5Ul4vfDfwiJm1GbEb8tiGrdcyWxDAMwzCiGI9zDIZhGEYCzDEYhmEYUZhjMAzDMKIwx2AYhmFEYY7BMAzDiMIcg2EYhhGFOQbDMAwjiv8Dm8XU75AgkCYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(psetRK4.lon, psetRK4.lat, 'r.', label='RK4')\n", "plt.plot(psetAA.lon, psetAA.lat, 'b.', label='Analytical')\n", "plt.legend()\n", "plt.show()" ] } ], "metadata": { "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.6.10" } }, "nbformat": 4, "nbformat_minor": 4 }