Lecture.22-Molecular-Dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular Dynamics\n",
    "## Lecture 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import sympy as sp\n",
    "sp.init_printing()\n",
    "from ipywidgets import interact\n",
    "from IPython.display import display\n",
    "plt.style.use('default')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lecture we will consider the mechanics of multiparticle systems. The methods are applicable to any system with a large number of particles, such as\n",
    "* droplets in an aerosol\n",
    "* particles in a flame\n",
    "* stars in a galaxy\n",
    "* motion of molecules (hence the name, molecular dynamics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verlet method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equations of motions of each particle $i$ is\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{d v_{i,x}}{dt} &= a_{i,x}  \\\\\n",
    "\\frac{d x_{i}}{dt} &= v_{i,x} \\\\\n",
    "\\frac{d v_{i,y}}{dt} &= a_{i,y}  \\\\\n",
    "\\frac{d y_{i}}{dt} &= v_{i,y} \\\\\n",
    "\\end{align}\n",
    "\n",
    "where the $v_{i,x}$ and $v_{i,y}$ are the velocity components of the $i$-th particle located at position $(x_i, y_i)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common second-order numerical method for solving such systems is give the Verlet method:\n",
    "\n",
    "\\begin{align}\n",
    "x_i[n+1] &= 2 x_i[n] - x_i [n-1] + a_{i,x} [n] (\\Delta t)^2  \\\\\n",
    "v_{i,x} [n] &= \\frac{x_i[n+1] - x_i[n-1] }{2\\Delta t} \n",
    "\\end{align}\n",
    "\n",
    "with similar equations for $y_i$ and $v_{i,y}$.  The Verlet method can be derived starting from the central difference approximation for a second derivative:\n",
    "\n",
    "$$\n",
    "\\frac{d^2 x_i}{dt^2} \\approx \\frac{x_i[n+1] + x_i[n-1] - 2 x_i[n]}{(\\Delta t)^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lennard-Jones Potential"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In previous lectures, we considered the gravitational interaction between celestial bodies.  For molecules, the attractive forces are given by Van der Waals forces.  As atoms get close, there is a repulsive force due to the overlap of their electron clouds.  \n",
    "\n",
    "This can be modelled with the Lennard-Jones potential:\n",
    "\n",
    "$$ U(r) = 4 \\epsilon \\left[ \\left(\\frac{\\sigma}{r}\\right)^{12} -  \\left(\\frac{\\sigma}{r}\\right)^{6} \\right] $$\n",
    "\n",
    "where $\\epsilon$ and $\\sigma$ are constants that set the energy and distance scales associated with the interaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = sp.Symbol('r')\n",
    "U = 4*((1/r)**12 - (1/r)**6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAD0CAYAAABThLtwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYVPWd7/H3t7be2DdFdmxUQBGhAzqZiSsjGoOOWzCJOqMJOknMeM1MrknmkonOmEwyRmN0kseJcUsiyTWZC4OoiYpxCYqAuCCiiKAtjTT70vRW9b1/VHXT3XR3nWqqq6qbz+t56uk65/zq1Jei63z6/M45v2PujoiISDaE8l2AiIj0HgoVERHJGoWKiIhkjUJFRESyRqEiIiJZo1AREZGsUaiIiEjWKFRERCRrFCoiIpI1gUPFzJa1me5rZqdkvyQREempMtlTKQIwsx8BuPte4D+7oygREemZMgkVM7NhwBfMzFLzSrJRhJmNMrOlZrbWzNaY2T+09+ZmdpeZrTez181sWjbeW0REsieSQdtvAi8AvwbuMLN3yN4xmUbg6+6+ysz6AivN7I/u/laLNucBE1KPmcBPUz9FRKRApA0FMzvNzMzdn3D349z9RuA3QDlwbTaKcPcqd1+Ver4XWAuMaNPsQuAhT3oJGGBmw7Px/iIikh1B9lSuBu5J7Zk8ATzh7suAZZ2/rGvMbCxwCvBym0UjgA9bTFem5lV1Rx0iIpK5tKHi7tcDmNkJJLugHjCz/sBSkiHzorvHs1GMmfUBfgfc6O572i5ur7x21jEPmAcwadKk6WvWrMlGaUecRMIZ/60l3HjOBG4857h8lyMiudXe9jaQIN1fxQDu/ra73+Hus4GzSB5fuYxD9yi6xMyiJAPlV+7++3aaVAKjWkyPBDa3beTu97p7hbtXlJRk5TyCI1IoZERCRn1jIt+liEgPEqT7q8rMtgKvt3g87+5LgCXZKCJ1Ntl9wFp3/1EHzRYBXzWzBSQP0O92d3V9daNYJKRQEZGMBOn+Gmhm44ApqccXgJ+a2Z+A6919dxbq+CRwJfCGma1OzfsWMDpVw89IBtj5wHqgBvi7LLyvdCIWCVEfV6iISHCBTil29/eB94GF0Lxn8W3gTrKwcXf3F0jTh+fuDnzlcN9LgouFtaciIpnJ5DqVZqkN/L+a2dos1yMFRN1fIpKpIAfqbzKzc1JX07ecXwQUd1tlknexSIg6dX+JSAaC7KkcDZwDnJQ6Q+t14D3gEyTP1pJeSt1fIpKpIAfqv9H03MwGACcBxwO/d/c/dmNtkmdF6v4SkQxldEzF3XcBz6ce0svpmIqIZEo36ZIO6ZRiEcmUQkU6pGMqIpIphYp0SN1fIpIphYp0KBYJq/tLRDKiUJEOqftLRDKlUJEOxSIh6hQqIpIBhYp0KHmdSlZulSMiRwiFinQoFgnRED/kPmgiIh1SqEiHYmFdpyIimVGoSIdikRDxhBNPaG9FRIJRqEiHouHkr4fOABORoBQq0qFYRKEiIplRqEiHmkKlLq4zwEQkGIWKdKhI3V8ikiGFinRI3V8ikimFinSoOVR0WrGIBKRQkQ7F1P0lIhkqmFAxs1+Y2VYze7OD5WeY2W4zW516zM91jUcadX+JSKYyup1wN3sAuBt4qJM2z7v7BbkpRxQqIpKpgtlTcffngB35rkMOOnhKsUJFRIIpmFAJ6DQze83MHjezyfkuprfTMRURyVQhdX+lswoY4+77zOx84P8BE9o2MrN5wDyA0aNH57bCXqZI3V8ikqEes6fi7nvcfV/q+RIgamZD2ml3r7tXuHvF0KFDc15nb6JjKiKSqR4TKmZ2tJlZ6vkMkrVvz29VvZuuUxGRTBVM95eZPQKcAQwxs0rgO0AUwN1/BlwK/L2ZNQIHgLnurjHZu5GOqYhIpgomVNz9ijTL7yZ5yrHkSNOeSoP2VEQkoB7T/SW513xKsfZURCQghYp0SN1fIpIphYp0yMx0n3oRyYhCRToVi4S0pyIigSlUpFPRsClURCQwhYp0SnsqIpIJhYp0KhbRMRURCU6hIp2KhbWnIiLBKVSkU7FIWNepiEhgChXplLq/RCQTChXpVFE4RH1jPN9liEgPoVCRTunsLxHJhEJFOqXuLxHJhEJFOqWzv0QkEwoV6ZS6v0QkEwoV6ZRCRUQyoVCRTumYiohkQqEindIxFRHJhEJFOlWkPRURyYBCRTqlYyoikgmFinQqFg6RcGjU3oqIBKBQkU7FIqn71CtURCSAggkVM/uFmW01szc7WG5mdpeZrTez181sWq5rPBI1h4q6wEQkgIIJFeABYHYny88DJqQe84Cf5qCmI55CRUQyUTCh4u7PATs6aXIh8JAnvQQMMLPhuanuyBUNJ39FdE8VEQmiYEIlgBHAhy2mK1PzpBsV6ZiKiGSgJ4WKtTPPD2lkNs/MVpjZiurq6hyU1bvFwur+EpHgelKoVAKjWkyPBDa3beTu97p7hbtXDB06NGfF9VY6piIimehJobIIuCp1FtipwG53r8p3Ub2dTikWkUxE8l1AEzN7BDgDGGJmlcB3gCiAu/8MWAKcD6wHaoC/y0+lRxZ1f4lIJgomVNz9ijTLHfhKjsqRFHV/iUgmelL3l+RBU6jolGIRCUKhIp3SKcUikgmFinQqFg4D6v4SkWAUKtIpHVMRkUwoVKRTTaHSoO4vEQlAoSKd0p6KiGRCoSKdar5ORXsqIhKAQkU6FQ0nh1zTKcUiEoRCRTplZrpPvYgEplCRtIrCChURCUahImnFIiHq4/F8lyEiPYBCRdJS95eIBKVQkbQUKiISlEJF0oqGQzqlWEQCUahIWjEdqBeRgBQqklYsEtJ1KiISiEJF0tIxFREJSqEiaRVFdExFRIJRqEhaOqYiIkEpVCQtdX+JSFAKFUkrpu4vEQlIoSJpqftLRIIqmFAxs9lmts7M1pvZze0s/1szqzaz1anHF/NR55EoFgnpzo8iEkgk3wUAmFkYuAeYBVQCr5jZInd/q03T37j7V3Ne4BFO16mISFCFsqcyA1jv7hvcvR5YAFyY55okRQfqRSSoQgmVEcCHLaYrU/PausTMXjezR81sVG5Kk6LU2F/unu9SRKTAFUqoWDvz2m7B/gcY6+5TgKeAB9tdkdk8M1thZiuqq6uzXOaRKRYJ4Q6NCYWKiHSuUEKlEmi55zES2Nyygbtvd/e61OR/AdPbW5G73+vuFe5eMXTo0G4p9kgTiyR/TdQFJiLpFEqovAJMMLNxZhYD5gKLWjYws+EtJucAa3NY3xEtFlaoiEgwBXH2l7s3mtlXgSeBMPALd19jZrcAK9x9EfA1M5sDNAI7gL/NW8FHmFgkDKALIEUkrYIIFQB3XwIsaTNvfovn3wS+meu6RN1fIhJcoXR/SQFrChVdqyIi6ShUJC0dUxGRoBQqklYskjzjW8dURCQdhYqkFQunDtRrT0VE0lCoSFo6UC8iQSlUJK3mUInH81yJiBQ6hYqkpQP1IhKUQkXS0inFIhKUQkXSKkqFSkNcA0qKSOcUKpKWDtSLSFAKFUnr4DEVHagXkc4pVCStg2d/aU9FRDqnUJG01P0lIkEpVCStSMgwU6iISHoKFUnLzIiFQ9Sp+0tE0lCoSCCxSEh7KiKSlkJFAilSqIhIAAoVCSQWVqiISHoKFQkkFgnplGIRSUuhIoHomIqIBKFQkUAUKiIShEJFAomF1f0lIukVTKiY2WwzW2dm683s5naWF5nZb1LLXzazsbmv8sgVDYc09L2IpFUQoWJmYeAe4DxgEnCFmU1q0+xaYKe7lwN3AP+e2yqPbOr+EpEgCiJUgBnAenff4O71wALgwjZtLgQeTD1/FDjbzCyHNR7Rjh3ah2MGFOe7DBEpcJF8F5AyAviwxXQlMLOjNu7eaGa7gcHAtpxUeITbV9fIqx/syncZ0gO5O40Jp6ExQQKIJ5x4wmlMJJI/407Ck/OSPw++pmlewpPtHJLTCZrnh8xoiCdIeHKeuxMJhahtjJPw5LrcwSx5o7mmaceJhUMcaEgkX5eqNRyyZK0ODpREQ9TUx5tf46n5IYPGRHK66d9ZVhRhX21jal1goWStTa/BndJYhP31jQdfl1xCLByiNtUb0PReAEWRMHUNcVreIq9fSYTdNQ3N6wAojoWpqYu3WufA0hg79te3+v/oXxJl94EGPPXigWVRduxvaNXmlgtP7PL/d6GESnt7HG1vMxikDWY2D5gHEI1GqaioOPzqull1dTVDhw7Ndxmdqtpdy/Z9dVT8rn++S0mrJ3yekLs6mzbUTRthb7VRTm4MmzfING2YD/48UFtLWUkx9fFEcmOXWhYyozHRYiPYYqNrlgyPll/QstTGtK3iaJjahtb36imNhampbz0vZEbCD/nKt7vesliY/W1eHw3bIXcvbe99IiGjMeGdtgEIh4x4ovX62tbSXs0l0TAH2vx7DYiEQzS2czJMNBKiITXfUpvB9mpq72Sa0liYA23atfq8DUqjYWra1HPrRe++6e5dSpZCCZVKYFSL6ZHA5g7aVJpZBOgP7Gi7Ine/F7gXoKyszFesWNEtBWdTRUUFhV7nPUvX88Mn1/HCrbMpjobzXU6nesLnCcHrTCScXTX17KhpYGdNPTv217Orpp76eILNu2rZc6CB3Qca2FPbSFEkxIbqfeyra2R/XZx9dckN3MDSKLtrWv812rc4wt7agxtAAwaXxdhX20gsEqIo9dhS9RGzZ0zmg501xMIhouEQsUiIPkUR3CESNqLhENGwEQmHiIYs+TtiyQ10JBQiEjL6lUSpa4wTTk2HU4+SaJh4IrmHEAoZYTOKIiESfnBeyIyQJTfkyedGOJQc7LRp2VVXXsmvfvVLQpYcVTtshqWeh8wwnJCFMKN5HiTXadDcFppem9yIO06oZRtILUtO0zx9cBlNyzz52qbpnvK7aWa1XX1toYTKK8AEMxsHfATMBT7Xps0i4GpgGXAp8Ix7O3+2SLcYWBoDYGdNPcP7l+S5mt6jpq6Rqj21bN51gKpdtWzefYC9BxrYsG0/2/bVU723jm376igf1oe3t+xt9dqKMQN5rXIX/Uui9CuO0rckyogBxZxwdD/KisL0KYrSpyhMn+IIfYujFEdDlEQjlMbClMbCFEfDlDT9jIYpjoYoioQJh1p3ClRUVPCfdxX+hjCyt4rJxxTanvSRd9i3IEIldYzkq8CTQBj4hbuvMbNbgBXuvgi4D3jYzNaT3EOZm7+KjzwDS6MA7NzfoFDJUF1jnI3VNby3bR8bt+9n47b9bNxew46/+jp/98ArvPx+6x3uMYNK6VsSYWifIk44ui9D+xZxdP9i+pdEGVgaY1BZjIFlMQYURyktCqPzVaSQFESoALj7EmBJm3nzWzyvBS7LZJ1DhgzJTnHdbN68efkuIa0BqT2VXTX1aVrmXz4/z53763nzo928sXk3b1ft5e0te3ivej/TRw9g+cadAAzpU8S4IaWcODjMX086iitmjObo/sUc07+Eo/oXURQprO7FnvD7CT2jzp5QY8q9XX2h9eYepIqKih5xTKUneHvLHmbf+Tz3fG4an54yPN/lFAR3573q/azctIOVm3ayYtNONlTvZ8bYQSzfuIMRA0qYOLwvxx/dl8nD+zFmSBljBpfRp6hg/pYT6UiXd3/12y2BtDymciSr2n2AP62r5vl3q1mxaScf76kDkt2D08cM5NLpI5k+eiAnDO9H/5JonqsVyb1CufjxsKQb4qWuro7PfvazlJeXM3PmTDZu3JiHKuGJJ57g+OOPp7y8nO9///uHLH/ggQcYOnQoU6dOZerUqfz85z/PeY3XXHMNw4YN48QTW59N2LSB3FlTz9e+9jXKy8uZMmUKq1atynmN0HGdTZ599ln69+/f/FnecsstXXofd2fdlr3c8cd3uO7hFZz2vWe4+fdvsOqDXXx6ynB+cMkUnrrpdFb9n1n8/OpP8OUzypk5fnDz5/Xhhx9y5plnMnHiRCZPnsyPf/zjdt8jn59pkBqz9XkejtraWmbMmMHJJ5/M5MmT+c53vnNIm0L4rgepsxC+603i8TinnHIKF1xwwSHLujQ8lrv36AfJA/vvAeOBGPAaMMndmT59uru733PPPX7ddde5u/sjjzzil19+uedaY2Ojjx8/3t977z2vq6vzKVOm+Jo1a1q1uf/++/0rX/lKzmtr6U9/+pOvXLnSJ0+efMiyE/75cb/mrsU+e/ZsTyQSvmzZMp8xY0Yequy8Tnf3pUuX+qc//ekur/+jnTX+k6ff8fN+/JyP+d+LfezNi/2m37zqP3t2vb9dtccTiUSg9WzevNlXrlzp7u579uzxCRMmHPL//thjj+X1Mw1S4+F+ntmQSCR879697u5eX1/vM2bM8GXLlrVqUwjf9SB1FsJ3vcntt9/uV1xxRdv/36bt65eBn6WezwV+42m2yb1hTyXtEC8LFy7k6quvBuDSSy/l6aefbr6aNFeWL19OeXk548ePJxaLMXfuXBYuXJjTGoL41Kc+xaBBg9pdNrA0ytoNH3DVVVdhZpx66qns2rWLqqqqHFfZeZ1d1RhP8Mzaj7n6F8v5y39/hv/4wzscf1Rfbr1wMi9/62xuv3wq151+LMcf3TfwGVfDhw9n2rRpAPTt25eJEyfy0UcftWqzcOHCvH6mQWosBGZGnz59AGhoaKChoeGQ/4dC+K4HqbNQVFZW8thjj/HFL36xoyYZD4/VG0KlvSFeRrRs8NFHHzFqVPLaykgkQv/+/dm+fXvuKmxTA8DIkSPb/eL+7ne/Y8qUKVx66aV8+OGHhyzPpwGlMfbWJQL9OwrBsmXLOPnkkznvvPNYs2ZNh+321jZw3wvvc8Z/PMsPnljHOx/v5StnlvPcP53JHZ+dypWnjWVY38Mf92zjxo28+uqrzJzZegSioL8budBRjRD88+xO8XicqVOnMmzYMGbNmtXpZ5mv73qQOqEwvus33ngjP/jBDwiFOoyCVsNjAU3DY3WoN4RK2uFb2vtLJdd/OQSp4TOf+QwbN27k9ddf55xzzmn+i6tQDCyL0hg+dONaiH+FTZs2jU2bNvHaa69xww03cNFFFx3SZn9dI/csXc/f/3Ilty5+i2P6l3DjrAk8/40z+fpfH8/owaVZq2ffvn1ccskl3HnnnfTr16/VskL4/YTOawzyeeZCOBxm9erVVFZWsnz5ct58881Wywvls0xXZyF81xcvXsywYcOYPn16Z80CDY/VUm8IlbRDvIwcObL5L4HGxkZ2796d9a6TdFrWAMndzmOOOaZVm8GDB1NUVATAl770JVauXJnTGtMZUBrDivuk/XcUgn79+jV3QZx//vk0NDSwbVty7NHahkbue+F9Tv/hUn745DqKIiH++8t/wW+vP43ZJw4nEs7u16KhoYFLLrmEz3/+81x88cWHLA/yu9Hd0tXY2eeZDwMGDOCMM87giSeeaDW/EL7rLXVUZyF811988UUWLVrE2LFjmTt3Ls888wxf+MIX2jZr3r52NjxWS70hVJqHeDGzGMmDSYtaNpgzZw4PPpjsFnz00Uc566yzcv7Xyyc+8Qneffdd3n//ferr61mwYAFz5sxp1aZlP/qiRYuYOHFiTmtMZ2BpFCvqw0MPPYS789JLL9G/f3+GDy+861a2bNnS/Ffr8uXLSSQSDBo0iD+s2cJF9/yZu55+l+OO6svvv/wX3Pe3Mzhl9MBuqcPdufbaa5k4cSI33XRTu23mzJmT1880SI3tfZ6DB3faC5J11dXV7NqVHCn7wIEDPPXUU5xwwgmt2hTCdz1InYXwXf/e975HZWUlGzduZMGCBZx11ln88pe/bNusaXgsCDg8Vo+/TsU7GeLl2GOPBeDaa6/lyiuvpLy8nEGDBrFgwYKc1xmJRLj77rs599xzicfjXHPNNUyePJn58+dTUVHBnDlzuOuuu1i0aBGRSIRBgwbxwAMP5LzOK664gmeffZZt27YxcuRIvvvd79LQkByIcOD4MzkQh3Hjx1NeXk5paSn3339/zmtMV+f111/Po48+yk9/+lMikQglJSXcc/+vueGRV1m2YTtD+hRxz+dP4S/Lu3+E4BdffJGHH36Yk046ialTpwJw22238cEHHzTXev7557NkyZK8faZBamz7eS5YsCDnG+uqqiquvvpq4vE4iUSCyy+/nAsuuKDVd6gQvutB6iyE73pH5s+fz6233jrHuzg8lq6ol8Due+F9bl38Fqvnz2oetqXQuTsLV29m/sI3aYg73zr/BObOGE00y11cIr2MrqiX7tc8qGRNQ48Ild0HGvjWf7/BY69Xcdn0kdxw1oSsHnwXkUMpVCSwlkO1jKMsz9V0bs1Hu5n38Eq27KnlG7OP5/pPHdt8XwsR6T4KFQlsQGpPpdBHKv79qkpu/v0b/FX5EH7yuVOY1k0H4UXkUAoVCezg8PcNaVrmRyLh/OSZd7njqXc5dfwg/uOykxlYVvjddCK9iUJFAmt5TKXQNMQTfOPR11m1aQdXnTaG/3PBJB2MF8kDhYoE1q84SsgKr/urtiHOV3+9iqfWbuWfzj2eL59xbEFe5S9yJFCoSGChkNG/JFpQ91SpbYhzy+K3ePrtrdx60YlceeqYfJckckRT/4BkZGBprGC6v+obE3z5V6t4ZPkH/PizpyhQRLLguuuuw8w+2dXXK1QkIwNKowXR/RVPODf9djXPvL2Vf7voJOZMLbzxx0R6kng8DsDLL78M8FJX16PuL8nIwNIYVbtr81qDu3PL4jUsfr2Kb553Ap+bOTqv9Yj0VJdddhmjRo3i1Vdf5eyzz+aSSy7huOOOY/Xq1XEzOxn4CTAEOIHkVfa3uPuht7JsQaEiGRlQGmNt1Z681vBfz2/ghXe38bWzyrnu9GPzWotIT/bGG28wceJEli5dCsCPfvQjZs+ejZkVA78BrnL35WZ2K1AM/Eu6dea9+8vMBpnZH83s3dTPdq9UM7O4ma1OPRa110a638DSaF6PqTz3TjXff/xtjj+6L/9r1nF5q0Okp6utrWXHjh3Mnz+/ed6TTz7J7NmzAc4BVrn78tSi14FB6UYohgIIFeBm4Gl3nwA8nZpuzwF3n5p6zOmgjXSzgWUxDjTEqW2I5/y9N23fzw2PvMpxR/Xlh5eerNOGRQ7DmjVrmDlzJpFIssOqpqaGXbt2Nd3L50TgjRbNpwGrgqy3EEKl5T2QHwTyc0s5CeTgUC253VvZX9fIvIdWYgb3XllBWZF6bkUOxxtvvMGUKVOap5cuXcqZZ57ZNLkdmAJgZscBFwOB7iNQCKFylLtXAaR+DuugXbGZrTCzl8xMwZMnLQeVzBV35+bfvc67W/dy9xXTNNKwSBa0DZXHH3+8qesL4BGgj5m9CdwLXOHu24OsNyd/7pnZU8DR7Sz6dgarGe3um81sPPCMmb3h7u+1817zgHkAo0frrKBsG1CS+z2VX738AW9v2cu3zp/IX04YkrP3FenNbr/99lbTf/7zn7njjjsAcPd9wGe6st6chIq7n9PRMjP72MyGu3uVmQ0Htnawjs2pnxvM7FngFOCQUHH3e0kmKxUVFb33DmR5cnBQydzsqWzavp/blqxl+piBXPPJcTl5T5Ej0apVgQ6ZpFUI3V8t74F8NbCwbQMzG2hmRannQ4BPAm/lrEJpNrAsd4NKxhPO13/7GuGQ8YNLp+h+KCI9QCGEyveBWWb2LjArNY2ZVZjZz1NtJgIrzOw1YCnwfXdXqOTBgJIop44bxIH6xm5/r58/v4EVm3by3TmTGd6/pNvfT0QOX95PoUkd/Dm7nfkrgC+mnv8ZOCnHpUk7SmIR3tu2n9Ef7+3W91m3ZS+3/+Edzp18FH9zyohufS8RyZ5C2FORHmbc4DI2bqvptvXXNya46ber6Vsc4ba/OUnXo4j0IAoVydiYwaVs3L6/29b/65c38c7He7nt4pMY3Keo295HRLJPoSIZGzukjK1769hfl/3jKh9sr+G2x9/m4mkjOXdye2ehi0ghU6hIxsYOLgNg0/bsd4HdsvgtoiHjJo3rJdIjKVQkY2NSV7Rnuwts6bqtPLX2Y244ewJH9SvO6rpFJDcUKpKxsUOSeyrZDJW6xji3/M9bjB9SposcRXqwvJ9SLD1Pn6IIQ/oUsSmLZ4Dd98L7vL9tPw9eM4NYRH/riPRU+vZKl4wbUsr7WdpTqdp9gLufWc+sSUdx+nFDs7JOEckPhYp0yZjBZWzKUqjctuRtGhPO/AsmZWV9IpI/ChXpknFDyvh4Tx01hzlcy0sbtvM/r23m+tOPZdQgDWkv0tMpVKRLms4AO5zTihsa4/zr4rcYMaCEv9e95kV6BYWKdEnTtSobt3W9C2zBKx+ys6aB73xmEiWxcLZKE5E8UqhIlxy8VqVreyq7axr40R/fYdSgEmZNOiqbpYlIHilUpEv6FkcZ0ifW5YP1dz79DrsPNDD/gskaMFKkF1GoSJeNHVzG+13o/lq/dR8PL9vEZz8xmknH9OuGykQkXxQq0mXjhpSypzbzO0D+22NvURIN8/W/1vheIr2NQkW6bPIx/VlbtZctu2sDv+bZdVtZuq6ar509gSEa1l6k11GoSJdNGzMQgFc27gjUviGe4NbFbzFuSBlX/8XYbqxMRPJFoSJdNml4P0pj4cCh8suXNvFe9X6+ff5Eje8l0kvpmy1dFgmHuHjaCD7ckf604q17a/njWx9z3olHc/bEYTmoTkTyQaEih6V8aB+WrqtmQ/W+Ttv9y6I1rNi0k38893idQizSiylU5LDMmnw0/YojPLtua4dt/rBmC39Y8zH/cPYEjh3aJ4fViUiu5T1UzOwyM1tjZgkzq+ik3WwzW2dm683s5lzWKB0bMaCEk0b058Flm0gk/JDllTtr+Mf/+xrnTDyKeZ8an4cKRSSX8h4qwJvAxcBzHTUwszBwD3AeMAm4wsw0TnqB+NzMMRzVt5gn11S1mn+gPs6/P/42ITO+ef4JRMOF8OsmIt0p799yd1/r7uvSNJsBrHf3De5eDywALuz+6iSIcycfxeA+MRatrqJ6T/Kald0HGvjGo6+xtmovd86dypjUAJQi0rv1lNsJjwA+bDFdCczMUy3SRiQc4vpPjeNnz23gn/8XlVSVAAABvElEQVT7NcYOKqFyVwO7ahu44cxyzjheZ3uJHClyEipm9hRwdDuLvu3uC4Osop15h3bgJ99rHjAPYNIk9ZDlysmjB3Hd6bBk9Ufsa0gwanARXzpxPFNHD8x3aSKSQzkJFXc/5zBXUQmMajE9EtjcwXvdC9zbNHmY7ysZmDpqEFNHDcp3GSKSR3k/phLQK8AEMxtnZjFgLrAozzWJiEgbeQ8VM/sbM6sETgMeM7MnU/OPMbMlAO7eCHwVeBJYC/zW3dfkq2YREWmfuffqHqJe/Y8TEekmXR72Iu97KiIi0nsoVEREJGsUKiIikjUKFRERyRqFioiIZI1CRUREskahIiIiWaNQERGRrFGoiIhI1ihUREQkaxQqIiKSNQoVERHJGoWKiIhkjUJFRESyRqEiIiJZo1AREZGsUaiIiEjWKFRERCRrFCoiIpI1ChUREckahYqIiGSNQkVERLJGoSIiIlmjUBERkaxRqIiISNYoVEREJGsi+S6gm1m+CxAROZJoT0VERLJGoSIiIlmjUBERkaxRqIiISNYoVEREJGsUKiIikjUKFRERyRqFioiIZI1CRUREskahIiIiWfP/AbMNm+7+DLgzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x10a16c8d0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.plot(U, (r, 0, 4), ylim=(-1.5, 2), ylabel = \"$U / \\epsilon$\", xlabel=\"$r/\\sigma$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To initialize our system, we could randomly distribute the particles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random(N):\n",
    "    x = np.zeros(N)\n",
    "    y = np.zeros(N)\n",
    "    for i in range(N):\n",
    "        x[i] = (np.random.rand() - 0.5) * np.sqrt(2)\n",
    "        y[i] = (np.random.rand() - 0.5) * np.sqrt(2)\n",
    "        plt.plot(x[i],y[i], 'o')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXl4XOWZ4Pt7JZVKwovkNSpvMXK8JUQdEzUEiCYYJXZu3Iq5kBBPyECYdPsmgcE4j7mNw+Drhk6gGyfGudDpdifpwAwZQsDBaES3HYQz7THg4GCuCHjDCgTJEl4lb5Ks5bt/VB1ZVTqn1lNnq+/3PHpK9dVXdd6qOvV+33lXUUqh0Wg0msKiyG0BNBqNRuM8WvlrNBpNAaKVv0aj0RQgWvlrNBpNAaKVv0aj0RQgWvlrNBpNAaKVv0aj0RQgWvlrNBpNAaKVv0aj0RQgJW4LYMXkyZPV7Nmz3RZDo9FofMXvf//740qpKanmeVb5z549mz179rgthkaj0fgKEXkvnXna7KPRaDQFiFb+Go1GU4Bo5a/RaDQFiFb+Go1GU4Bo5a/RaDQFiGejfTQaTfo8t7edh7cd4EhXD9Mqy7l76XyuXzTdbbE0HkYrf43G5zy3t521W96kp38QgPauHtZueRNALwAaS7TZR6PxOQ9vOzCs+A16+gd5eNsBlyTS+AGt/DUan3OkqyejcY0GtPLXaHzPtMryjMY1GtDKX6PxPXcvnU95qDhurDxUzN1L57skkcYPaOXvU7obGzl0XT37Fn6UQ9fV093Y6LZIGpe4ftF0Hrzh40yvLEeA6ZXlPHjDx7WzV5MUHe3jQ7obG+m4bx2qtxeAgSNH6LhvHQAVDQ1uiqZxiesXTdfKXpMRWvn7kKMbHxlW/Aaqt5ejGx/xrfI/t/cop7e9y2BXH8WVYcYvnc2YRVMDczw7eLbzJA+2dtDe18/0cIi11RFurJrotlgan6KVv004mWQz0NGR0bjXObf3KF1bDqH6hwAY7Oqja8shgLwoZKePZwfPdp5kzYH36RlSALT19bPmwPsAegHwGy1PQ/P90N0GFTOgfh3U3OS4GIVr8295GjZeBusro7ctT6f1NDNbu5Fk097Vg+Jiks1ze9vzInpJJJLRuNc5ve3dYUVsoPqHOL3t3UAczw4ebO0YVvwGPUOKB1v9ueAXLC1PQ+Od0P0+oKK3jXemrX/spDCVf5ZfgGFrHzhyBJQatrX/3XN7HU2ymbr6LqSsLG5MysqYuvquvBwv3wx29WU07tbx3HSyt/f1ZzSu8SjN90N/Qv5Ff0903GEK0+yT7AtIcvllZWvv7AVk9Px8JdkYdv2jGx9hoKODkkiEqavv8q29v7gybKp4iyvDrhzv4O5OXtl6mLMn+xg7McxVy+fwoaOvuepknx4O0Wai6KeHQ3k/tsZGuttMh1VXG//wzZeGz7d5V1blXRRbdv4i8nkROSAi74jIPRZzbhKRt0XkLRH5hR3HzRqLL8ByPIaVTX1KT5fpeD6TbCoaGpj7UjML973N3Jeafav4AcYvnY2E4k9FCRUxfulsx493cHcnO57cz9mT0cXh7Mk+djy5nzd+st3Sye4Ea6sjlBfF7zDKi4S11f409RUsFTNMh88MTgYunm8Hd3fmXZScd/4iUgw8BnwOaANeE5HnlVJvj5gzF1gLXKOUOiUi7nrVKmbETD4m40koiUSiJh/gpemLePxj/wfHyicwdqCHULHQP3jRJquTbNLHcLI6FX2T7HivfHcXAxfi/QEDF4Y4NOE/MJUXR72WU052w6mro328jdlVY9wuvn5d1MQ8wvLQPxTm1bM3D98fuDDEK1sP5333b4fZ5wrgHaVUK4CIPAUsB94eMeevgMeUUqcAlFJHbThu9sS+gGcrr+bB6pW0h6cyve8YaycNcWOSp01dfRcd962jedJCfrToy/SVlAJwJnQJIQUTLgnRdb5fl9TNgjGLpjoaaWN1PGPHn0hf2FzJOulkv7Fqolb2Hsa4ajQ2D8YuHrioyA2zciza5/TAZF49ezOHej8T91pW56Gd2KH8pwMjt9FtwJUJc+YBiMguoBhYr5T6NxuOnR01N/Fs7zjWnJlIT1HUzttW9iHW9Ah0nrT8gRmmlcdfOjus+A36hxSXlJawd92S/MquyStjJ4ZNf3hjLlFIWVmc6cfPTnaN/byy9TCXFu3gU5OfZFzxcc4MRhX7K1tD8bv4mpuGF4Fff3cXZ3tHn29jJ+bH3zUSO2z+Jq5OVML9EmAucC3wH4GfiEjlqBcSWSkie0Rkz7Fjx2wQzZoHB2YNK36DdELnKhoaOBYeb/qYrqLof65aPoeS0vifRUlpEVevuIzIA/dTMm0aiFAybRqRB+73ta9FYy+R89tZPP7HjC85hohifMkxFo//MZHz2y2fY3W+XbV8Tr7FtWXn3wbMHHF/BnDEZM6rSql+4I8icoDoYvDayElKqc3AZoDa2trEBcRWcgmdm1ZZTruJotdVFP2PsUMzt9s2aGWvseTqil8QkvhdfKioj6srfgH8jelzkp9v+cUO5f8aMFdELgXagRXAVxPmPEd0x/9zEZlM1AzUasOxsyaX0Lm7l86P65wE2sEbJOZdWeXIj0/jP5KVBRkj5tYKq3EDt863nJW/UmpARO4AthG15/9MKfWWiNwP7FFKPR97bImIvA0MAncrpU7keuxcWFsdiUuXh/RD5wxHru6ZqtG4R0fnVloPb6C3r4OycITqOWuIVC235bXNlDyQtCyIWEQRSoooQrcQpfJqXcma2tpatWfPnqyfn04RLD8XympqbWLT65voPNdJ1ZgqVl2+imXVy9wWS+MxgnqedHRuZf/+exkaumh+LSoqZ8GC7+W8ACTWfoJoHoiEihg6PzBqfnFlmMg9V1ysHDAygTRUDg0/crR2j4j8XilVm2peIDN80y2C5dfQuabWJta/vJ7ewWjkSce5Dta/vB4gED9sjT0E+TxpPbwhTvEDDA310Hp4Q87K36r2U+KYwXC2eEIYp5tF29IhkLV9gl4Ea9Prm4Z/0Aa9g71sen2TSxJpvEiQz5PePvPfstV4JmRaUyquDEnNTbD6D7C+K3rrUcUPAVX+QS+C1XnOPPXbalxTmAT5PCkLm/vmrMYzwaqmlJQXO1qGJN8EUvlbRewEpQhW1RjzyACr8aDwbOdJal9+i8iON6h9+S2e7TzptkieY2Tl0clnzVJwgnGeVM9ZQ1FRfGh1UVE51XPW5PzaVrWfKr/4ESpvmDu8OBRXhqm8Ya5ne0CkIpA2/1wiefzAqstXxdlyAcqKy1h1+SoXpcovuplJahLbe65o7uefvlDMhRF7nqCcJ4ZdPx/RPqlqTflV2SdS0NE+fiaoURxW1L78lmlexoxwiD1Xf8wFibzHoevqhwsPGuz8qPBUfYjjY1VBnCea9KN9Aqv8NcEisuONUTVDIFpbpGPxJ5wWx5PsW/hRMPs9i7Bw34g6ix5pI6jJD+kq/0Da/DXBI+h+HDtIq72nh9oIatxFK3+NL9DNTFKTVntPD7UR1LhLIB2+muChm5mkJq32nll2sdMED638Nb7BrxnZTlLRkKLyaJZd7DT2kc+aRJmgzT4aTSFRvy5ab2YkofLouCbvGDWJevuOAIreviPs338vHZ1bHZdFK3+NppCouSlaaKxiJiDRW4cLjxUyyWoSOY02+2g0hcaINoIaZ8lnTaJM0cpf40me29tekP0SDu7udKWrk8YZSoorGBjsMh13XBbHj6jRpOC5ve1xndLau3pYu+VNgEAvAAd3d7Ljyf0MXIiWDj57so8dT+4H0AtAUBDzekuW43lE2/w1nuPhbQfiWmQC9PQP8vC2Ay5JlBlNrU0seWYJNY/XsOSZJTS1NqX1vFe2Hh5W/AYDF4Z4ZevhfIipcYGBgdG7/mTj+UTv/DWe40hXT0bjXiKXBipnT5rXkbca1/iPsnAkFukzetxptPLXeI5pleW0myj6aZXlJrO9RbIGKqmU/9iJYVNFP3aieX15jX20tLTQ3NxMd3c3FRUV1NfXU1NTY/txquesMW0/aUcp6kyxxewjIp8XkQMi8o6I3JNk3pdERIlIyqJDmuzYt3MHm2+/jR+saGDz7bexb+cOt0XKmLuXzqc8VBw3Vh4q5u6l812SKH1yaaBy1fI5lJTG/yRLSou4avkcW2TTmNPS0kJjYyPd3d0AdHd309jYSEtLi+3HilQtZ8GC71EWngYIZeFptvQdzoacd/4iUgw8BnwOaANeE5HnlVJvJ8wbB9wJ7M71mBpz9u3cwfbNjzJwIbp7PHP8GNs3PwrAwrrFboqWEYZT14/RPlVjqug4NzpsL50GKoZTV0f7RM/lnU89wZkTxxk3aTJ1K27J2znc3NxMf398ufD+/n6am5vzsvuPVC23VPZOvm87zD5XAO8opVoBROQpYDnwdsK8B4C/B5y/vikQdj71xLDiNxi40MfOp57wlfKH6AKQrrLvbmxMXs/GQXJttDPvyqqCVPYjcXoTY+z40x3PllRlHZx+33aYfaYDI4uFtMXGhhGRRcBMpdT/tOF4GgvOnDie0XgQMLpXDRw5AkoxcOQIHfeto7ux0RV5llUvY/3V64mMiSAIkTER1l+9XjdQyYBkm5h8UFFhHmNvNZ4N6ZR1cPp927HzNwtQHe4oISJFwEbg6ylfSGQlsBJg1qxZNogWHNJxSI2bNJkzx4+Neu64SZOdEtNxjm58ZLhtoYHq7eXoxkdc2/0vq16mlX0OOL2Jqa+vp7GxMc70EwqFqK+vt+0Yyco6GLt/p9+3HTv/NmDmiPszgJGxTOOAy4Dfisi7wKeA582cvkqpzUqpWqVU7ZQpU2wQLRik65CqW3ELJaXxkSElpWHqVtzimKxOM9BhnhZvNa7xPlablXxtYmpqamhoaBje6VdUVNDQ0GCrvb+3r4NdfJpV/Jib+RWr+DG7+HRcWQen37cdO//XgLkicinQDqwAvmo8qJTqBoalF5HfAmuUUrpHY5qk65Ay7IJOOYy8QEkkMqpvrTGu8Sd1K26Js31D/jcxNTU1eXHuGvyu5C/4ycAKLki02c5xpvIT9S1CJRMwri+cft85K3+l1ICI3AFsA4qBnyml3hKR+4E9Sqnncz1GoZOJQ2ph3eJAK/tEpq6+i4771sWZfuK6V8X61Z47MYfTQ7cxODSR4soyxi+dzZhFU12SWpOMIG5inpavcUHi1e0FKeNp+RprY/edft+2JHkppV4AXkgYMy0QrpS61o5jFhIVFRWmit5Oh5RfSdq9Ktav9lzPFXQN3IEiuusa7Oqja8sh3ji6l7858wM6z3VSNaaKVZevyshWX6jF55wgaJuYzgFzVZs47uT71hm+PsDMIVWsivjkudmc23u04Hewlt2rYv1qTw/cOqz4DVT/EOX/u5eOj0RtrpmUYYDCLT6nyY7p4RBtff2m426hC7v5AMMhNb58LCgYO1RGXf8Cqs9OpmvLIc7tPeq2iN4k1pd2EHOH2eT+CXH3jTIM6eD34nMaZ/n82A8QdSFurLxIWFvtnm9K7/x9Qk1NDVNe6GWwLz4OWPUPcXrbu67v/p3MTEybWL/aYo4zyOjP51jJyVFj6ZRhAH8Xn8sUbd7KjabWJra1rGdseBHnKm9iqHgSxYMnuPlDRdxY9WeuyaWVv48Y7DKv7mg17hSeLStRvw4a72T8wON0DfyXONNPX9EFfj51dN/UdMowgL+Lz2WCNm/ljlHsr+z8K5Sdf2V4/HfdEai5zjW5tNnHRxRXmld3tBp3CqczE9Mm1q92zKRWKksepbjoBKAorgzT+ZkBXp34h7jpmZRh8HPxuUxIZd4KQiHBfJNLsb98onf+PmL80tl0bTmE6r/Y8ENCRYxfOts9ofB4WYlYv9oxwJgRwxFgfesAm17flFW0j5+Lz2VCMvOWZ6/4PERLSwtfaPsC4f4w54vP84cJf6BtXNQXle5VZr4IpPJPVUDJrxh2/dPb3mWwq4/iyrAn4tX9WlYi4zIMsZwButugYgbX16/j+nuC3QjdyrxFySme/dk/UH4hvrqLXwsJ5gMjM7+sP2puHDM4hk+e+CQAxyuPp32VmS8Cp/yNAkpGHQ2jgBIQmAXAbWWfiBsZmY4TyxmgP6YIu9+P3ofo1UVAuXvp/DibPwBygdKp2yg7ZP4cT1zxeQCzzPwSVcKfdf8ZdV+oc73+U+CUfzoFlDT2EsSMzFHEcgbi6O+Jjrug/J3qPDXSvNXedR4p6aJ06jZKK97gXNl0xvaOViFev+JzCqvM/LL+MtcVPwRQ+Y8slJTOuMYegpaROYpYzkDa43nEMCcYu0qj0B9g2wJwbu/RYfPilZVhti+9jKtaPou6WLCX388/xTVvTqJk6GLcSOCu+HLA65n5gYv2sWqE7EaDZE2AqJiR2XgeSVbozw7O7T1K15ZDwyHERjmMD4XiK+3+cfp5dn38BD2XKBBh3OQpLFl5R7A3ASYc3N3J49/dxWPffInHv7uLg7ujUTz19fWEQvEZvHaXis6FwO38vdQg2Ws4ZSrIBc8662M5A3Gmn1B5dNxh8t156vS2d+MiyiCaTHjr0S+yafJ/j+tS1jFriEUrvu0JM4YbHNzdyY4n9zNwIfp5nT3Zx44n9wNQc2X0t+XV31zglL+hKDypQFzECVNBrjjtrD+4uzP9frmGXX9EtA/161yx9+fbnGCVNHht5+VM+OJHsg6PDSKvbD08rPgNBi4M8crWw8y7sirvpaJzIXDKH5I3SC5UnG5SnQ1OOuuT7diSLgAeiOzJd+ep4sqw6QJQXBl2pEuZJ0uFWHD2pPlCaTXuJQJn89eY41ST6lxw0lmfbMfmdfLdeWr80tlIKF41OJVMaCSOnTl+DJQaThzzaubw2Inm2fVW414ikDt/N8nIlOAgXo88gKhTPtrgevS43fh5xwb57TzlZjJhslIhXtz9X7V8TtwVJEBJaRFXLZ/jolTpoZW/jWRlSrBgZKidHT8+J5pUm9Hd2GjeaMUEJ531YyeGTRW9H3ZsTuBWMqGnS4WYYPyuvbjhS4VW/jaSyvmTLkaonRFxYYTaAVn/II1dopORB92NjXEtFgeOHKHjvmh0jNkC4KSz3s87tiDjx1Ih866s8oWyT0Qrfxuxy5RgFWqXa91+pyMPjm58JK63LoDq7eXoxkcsd/9OOevt2rF51cznVwqiVIhHsEX5i8jngU1EG7j/RCn1UMLj3wH+EhgAjgH/WSn1nh3HToWTkQN2mRJGRlq8d+YtWk79O+cHT3NJ8Xiu3flXnrR9mjHQYe6otRp3mlx3bHaa+TRRvFgq5NnOkzzY2kF7Xz/TwyHWVke4sWqiba9vt4k3XXJW/iJSDDwGfA5oA14TkeeVUm+PmLYXqFVKnReRbwF/D3wl12OnwomSsyO7HE0pL+VT5TC/52Kd92xMCUao3Xtn3uK1E//GoBoA4PzgaV+VzC2JRBg4MtqBWxLxXrZ1Nt2q7DLzpUMmvhM/093YSMnGR6gz3ufNN1DhsuJfc+B9eoaiZS3a+vpZc+B9AFsWgHyYeNPFjlDPK4B3lFKtSqkLwFNA3HW7UmqHUup87O6rgCM58fluMmJ0OWrv6kEBR3susO2Sfg5XRh8fOzHM4psXZKwIjFC7llP/Pqz48yF/vpm6+i6kLL5xupSVMXX1XS5JZE7i92h0q3pub3vS5zkVMWT4TgaOHAGlhn0n3bEkvaDgxff5YGvHsOI36BlSPNhqz9VrMhNvvrHD7DMdeH/E/TbgyiTzvwH8qw3HTUm+IwfMuhz1DSp2V8IPH8q+PZux4p9/6LTp416NfEjE2Jl6fcearFtVst2/UxFD2fhOEhkum9F7hOLuYsb9Gsa3z/TU92HH+7Sb9r7+jMYzxc3WrHYofzEZUyZjiMjXgFrgMxaPrwRWAsyaNStnwfIdOZDvJt6XlIzn/MDoBcDLkQ+JVDQ0eEa5WJHt9+hUxFCuvpO4shkCg5WDdH8VePJ9BpNEXzmNF31E08Mh2kwU/fRwyGR25iTLps43dph92oCZI+7PAEYZekXks8C9wBeVUqbLmlJqs1KqVilVO2XKFLMpGVG34hZKSuM/RDsjB6yadVeU535inN72LjWV/4FiiV+fiyWkIx9sxup7TNWMfd6VVSy+ecHwTj9bM18qrHwk6fpOzMpmqDCcWT44vLP2Arm+z3ywtjpCeVH8/ra8SFhbbY9MbmZT26H8XwPmisilIlIKrACeHzlBRBYB/0RU8R+14ZhpsbBuMUtW3sG4yVPyUnL27qXzCRWNvvA5d2Egpb04FYNdfXx43Mf480mf55Li8QBcUjyeP5+01BfOXj+RSzP2eVdWcev3r+H2f7yOW79/TV6ifHL1nViVxxiM+Su9En3lRR/RjVUT2TB/JjPCIQSYEQ6xYf5M26J9xiyaSuUNc4d3+sWVYSpvmOuPaB+l1ICI3AFsIxrq+TOl1Fsicj+wRyn1PPAwMBb4lYgA/Ekp9cVcj52KlpYWmve8QfeUD1PxkRo+VV/PQhvj3K9fNJ2/aXyLU+cTCqYNqpT24lQYl4MfHvcxPjzuY3HjXsVPBblG4vVm7Ln6TqzKZhSfjN56JfrKqz6iG6sm2hramYhb2dSilKl53nVqa2vVnj17sn5+YgljiJYzsLMAFsCl9zSZOjgE+OND2Vc/TAwBg+jloFO7gkxJDKuFqImtEJt7eI3EUtkA0gcVTxYz5g9jiDxwv+sKVmMfIvJ7pVRtqnmBzfB1qoTxtMpy2k0cg6nsxalws7hWNvitIFchEVc2Y2S0z5GZTH3A/Z213/BDU6R0CKzyd6qE8d1L57N2y5txoYLp2otT4dblYDb4rSBXoTGqbMYN7sniZ/zQFCldAqv8nSph7HV7sVPkM6zWrfR3r1Do799L+KEpUroEVvk7WcL4+kXTC07ZJ5Kvglxupr97AT+//yAuWn5oipQuge3kle9uR5p48hVW62b6uxfI9/vvbmzk0HX17Fv4UQ5dV29bKQVj0TISmIxF69xexyK984KV5cBLTZHSJbA7f3C+hLHT5LvaYKYsrFtsu3PXzfR3L5DP959pv4VMyFdZcldpeZrbL/yYEB/QzTiauYY3WehIU6R8ENidf9Axqg229fWjuFht8NnOk26LZitWeQ1eznewk3y+/2S1dHIlcIt2y9PQeCelPR8gQCVnaOBFrih/z7cWBa38fUq+qw16BTfT371APt9/PmvpBG7Rbr4f+uNDuksZ4Aulr/lS8UPAzT5BJt/VBr1CPvMd/FAjP5/vP5/9FsYvnW2apJi4aHnNdAkWjuruNvPJVuM+QCt/n5LvaoNeIh/5Dvm0d9tNvvI9pq6+K+4zAPtq6aSzaOW7UUo2WEVXlY+NUNQzeqGkwpHWJHlBK3+fsrY6EvfDAXurDQYdL9aOd5p819JJtWglM126pfytHNXd/bcwIbQp3vQTKof6dQ5LaB9a+fsU48fhtUtmv+DF2vFu4Ga/BS+aLq0c0ufOfpoJX50Ttf13t0V3/PXroOYmhyW0D638fUy+qw0GmUzt3Qd3d/LK1sOcPdnH2Ilhrlo+RzdpzxEvmi6TNlepucnXyj4RrfwzICgFnTSZ2bsP7u6M69h19mQfO57cT/P7J/jv73xQ0GU9ciEd06XTv7l0HdVBQCv/NAlSQadUFMIuNxN79ytbD8e1agRoURfYvuddjH2r0fQd0AtAmqQyXbrxm/NbNd1c0Mo/TYJU0CkZVrtcIJALQDr2brMm7TvLBkg0WKTT9F0TTzLTpVu/OT9V080FneSVJkEq6JQMs13uwIUhXtl62CWJ3Mfo0TuS00XmTZBSNX3PiZanYeNlsL4yetvydP6O5QEK5TfnFlr5p0mQCjolw2yXm2y8ELhq+RxKSuN/KuPV6N7NkHsTH0ti5QXofh9Q0dvGOwO9ANj1m3tubzvXPPQSl97TxDUPvZRzf+2goJV/mtTX1xMKxUch+LWgUzLMdrnJxguBeVdWsfjmBcOfwdiJYb51xYezbvqeFSblBejviY4HFDt+c8/tbWftljdp7+pBcdE3oxcAbfNPG8PGGPRon6uWz4mz+QOUlBZx1fI5LkrlPvOurBrl84jMqXSuiU8Aywukwo7f3MPbDsR12QPtmzGwRfmLyOeBTUAx8BOl1EMJj4eBJ4BPAieAryil3rXj2E4S9BLRcNGpm2m0z3N72wuum5mjTXwqZsRMPibjASbX35yVDyavvhmfkLPyF5Fi4DHgc0Ab8JqIPK+UenvEtG8Ap5RSHxGRFcDfAV/J9dia/GC2y02GcWlt7LDau3pY8z9e46Wf/xPf+pL9Nf4DS8vT1hmk9euiNv4AlRewJNnnkCHTKstpN1H0efPN+Ag7dv5XAO8opVoBROQpYDkwUvkvB9bH/n8GeFRERCllHjKh8RVml9YDRSF2lHyUOZsfBRheAILY2s8WDIeuodwNhy7EZ5YGqLyAKak+hwy5e+n8uI0J5Nk34yPsUP7TgZHXo23AlVZzlFIDItINTAKOj5wkIiuBlQCzZs2yQTSNE1hdQp8pGcvAhT52PvUEC+sW+7ofbd5J5tA1lF7AyguYks7nkAGGWS5bk2RH51ZaD2+gt6+DsnCE6jlriFQtz1gOL2KH8jeLeUvc0aczB6XUZmAzQG1trb4q8AlWl9aVF05z7dvvUd5/mEPX1XNJ3X2o/vgIGd+39rOLAnTommLxflV3G/98+22cOXGccZMmU7filrTNidn6Zjo6t7J//70MDUXP7d6+I+zffy9AIBYAO5R/GzBzxP0ZQGLFLGNOm4iUABVAYPoNxu0Oiiqo/uM5Iu915nxp7pcyC2aX1qHBfr7xZhOX9A8A0Xr5Q71FiMk2wLet/TIkaZ2aAnXojsLiczjTH+bM8WPR/48fY3uCOTEftB7eMKz4DYaGemg9vCEQyt+OOP/XgLkicqmIlAIrgOcT5jwP3Br7/0vAS0Gx9xu7g96+I4Cid6iL/TMu0DEllFMijlFmwUiuMsosHNzdafM7yJ3rF03nwRs+ztRyAaUY13+Gb7c8y2fbXo+bp3rM13vftvbLAKNOjZGdatSpaWlpiU6oXxd14I4kqA7dZJh8Dv2qmJ1HPxw3ZpgT80lkQSa2AAAgAElEQVRvn3l5b6txv5Gz8ldKDQB3ANuAfcDTSqm3ROR+EflibNpPgUki8g7wHeCeXI/rFUx3B8VC66VjoneyTMTxW5mF6xdN53f/zxf4178Yw3/p3cbS9/aMmtP31q9RA/G7fLsrJu7buYPNt9/GD1Y0sPn229i3c4dtr53Is50nqX35LSI73qD25bd4ttP6YjZZnRogenXY8COomAlI9LbhR8G38Sdi8jlsP/IR9p8ebRY8c+L46OdnQEfnVnbtqqP5pY+wa1cdHZ1b4x4vC5uX97Ya9xu2xPkrpV4AXkgYWzfi/17gy3YcKxfy4byx3B2ER6yrWdht/VpmYWFdNLTz0HX1o+rlD7T/jv4JEyiv/Wpeon327dzB9s2PMnAh+hnl0zyQaQvCtOrUFIJDNx0SPof2228Djo2aNm7S5KwPkY49v3rOmrg5AEVF5VTPWZP1cb1EwZR3GGWeiX3Ziat9pljuDvpG7NqzsNv6vczC1NV3IWVlcWNSVsakW+uJ3HMFMx6qI3LPFbY6enc+9cSw4jfIl3kgWQtCMwqlNlQ+qFtxCyWl8ed9SWmYuhW3ZP2ayez5BpGq5SxY8D3KwtMAoSw8jQULvhcIez8UkPJP58vOhuo5aygqirdRFg0qqv94LnonS7utWTExP5VZqGhoIPLA/ZRMmwYilEybRuSB+/PaMtDKDJCreWCYEVU123svmE6xakFYKLWh8sHCusUsWXkH4yZPARHGTZ7CkpV35HQ1l649P1K1nGuu2Un9de9wzTU7A6P4oYBq++TLeWOcDHHRPu+dI3KsP2q3zDLaJ9syC17C6f6w4yZNHo4ISRzPmYTko+l9H9BWNvq7sGpBWCi1ofKFYU60i7JwJGYFGD1eKBSM8rf6sqGIjs6tOa3okarl8c+/NuuXiiPTMguFTt2KW+Js/pC7eWCYhOSjta2bWTP//6an+KJpK7EFYSKFUBvKLwTdnp8OBWP2MTPPRBm0xfavcZ98mAeGSXDa33ismQ0H/p4ZvZ0IMCMcYsP8mZZdqTTeIuj2/HQQr4bb19bWqj17RocL5kJH51befvtuYHDUY2XhaVxzzU5bj6cJEBsvs0jCmgmr/+C8PBqNBSLye6VUbap5BbPzB8M+P2T6WFASNzR5wsYkrFTx5aloam1iyTNLqHm8hiXPLKGptSljGTTewo3vtGBs/gba0XORIBetsh2bqmrmWi+mqbWJ9S+vp3ewN/p65zpY//J6AJZVL0tbjkxLh+hzJX/Y9Z1mSkGZfWD0jw+ijp5Cs/fpz8Eddu2qs9h8pGd2XPLMEjrOjb5KjYyJsP1L29OSwSgdktitbfHNC0wXAH2u5Bc7vtORaLOPBdrREyVfeQ+a5OQactx5zry2k9W4GZmWDnHjXCmkput2fKfZUHBmHzAJzSxAgl60yqvkanasGlNlukusGpN+SHCmpUOcPlfMOsOt3fImQCBbg9rxnWZDwe38NVGCXrTKq5hmhGcQX77q8lWUFceXzSgrLmPV5avSliHT0iFOnyvJmq4HgYO7O3n8u7t47Jsv8fh3d/GNsu/k/J1mQ0Hu/DWFneTiZp+EURnhGTpPDQfgptc30Xmuk6oxVay6fFXajsGm1iZ+/rEfcrz/KGP7JnDln/6CuSdqk5YOcfpcCXLT9UR/y9mTffRuL2PNkgf4ae8Ps/pOs0Ur/4CQaTRGrkrIr5j9+HY8uR/A0QUgl895WfWyrBRDXFSJwNmyU/yvOU8RHlvCN677quX7d/pcCXLTdSt/S+/L49j+/cydu7mglX8AyDZ8sBB9H8mcnUEvpbHp9U3D4YQGA8X97J37b8y78jtJn2t1rnQ3NnJ04yMMdHRQEokwdfVdOddzCnLTdS+VatfK34KkLfc8RtDbzdmJl358TmN3VEl3YyMd961D9UYXlIEjR+i4L5r0lssCkGvTdS8zdmLY9Fxzo1S7Vv4mGC33jM5LRss9wJMLgI7cSR8v/ficxu6okqMbHxlW/Aaqt5ejGx/JefefbdN1r3PV8jmmORZulGrX0T4mpGy55zF05E76+L1PQi7YESk0koEO882F1bgm6ldafPOC4c3G2Ilhy+S6fKOVvwlptdzzELmGDxYSXvrxOc2y6mWsv3o9kTERBCEyJsL6q9dnHVVSEjHfXFiNa9yNNEskJ7OPiEwEfgnMBt4FblJKnUqY8wngx8B4ouU0v6eU+mUux803FRUVporejpZ7+fAlFGrkTrYUcp+EbCOFRjKswOZ9l7IPn6L68FaqjkZLsUhZGVNX32WHqIHDC5FmI8nV5n8P0KyUekhE7ond/+uEOeeBW5RSh0RkGvB7EdmmlOrK8dh5o76+Ps7mD/a03MunL8EPkTvP7W0PpBPPbzzbeZIHWzto7+tnejjE2upI2n0I4hWY0BueyP4FNwPCjJJ2W6J9gorXIs1yVf7Ludi36nHgtyQof6XUwRH/HxGRo8AUwLPKP18t95L5ErzoSLaTQkvZ9yrPdp5kzYH3h5vPt/X1s+ZAtE9BOguAmQIbKirlT5/+v1j8/WvsFzhAeC3SLFfl/yGlVAeAUqpDRKYmmywiVwClgGkFKRFZCawEmDVrVo6i5UY+Wu75zZdgJ8lS9rXyd44HWzuGFb9Bz5DiwdaOtJS/1xSYn/BapFlKh6+IvCgifzD5y8jGICIR4L8BtymlTDuqKKU2K6VqlVK1U6ZMyeTlfYGVz8AOX4LXCXLKvp9o7+vPaDyRTOsCaS7itUizlMpfKfVZpdRlJn9bgQ9iSt1Q7kfNXkNExgNNwH9VSr1q5xvwE/X19YRCobgxO3wJfsAqNT8IKft+Yno4lNF4Il5TYH7Ca5FmuZp9ngduBR6K3Y7qRycipcCvgSeUUr/K8XjepuXppJ2e8uVL8ANBTtnPBLczx9dWR+Js/gDlRcLa6vTCMw1F5ZVwRb/hpUiznDp5icgk4GlgFvAn4MtKqZMiUgt8Uyn1lyLyNeBfgLdGPPXrSqk3kr12vjp55Y2Wp6HxTugfYcYIlUPDj5K2+ssl8sJvFHq0T2K0F0Sv/BoaGhxdAArpnCtE0u3kVXBtHPPGxsug+/3R4xUzYfUfTJ+SGHkB0V3YhvkzffNjLHSFngkbN260zB9ZvXq1CxJpgohu4+g03W2ZjZM88sIPGOGb7V09KC6Gbwa55V4uFHK0l8Z7aOVvFxUzMhsn98gLtwl6xyW7KeRor2c7T1L78ltEdrxB7ctv8WznSXcFank6erW+vjJ62/K0u/K4QEEp/7yegPXrojb+kYTKo+MW5Bp5kS37du5g8+238YMVDWy+/Tb27dyR1evo8M3MKNRoL8O82dbXj+JiYplrC4Dhn+t+H1DR28Y7C24BKBjln/cTsOamqHO3YiYg0dsUzt611RHKiyRuLJPIi2zYt3MH2zc/ypnjx0Apzhw/xvbNj2a1AOjwzcyoqamhoaFheKdfUVHhuLPXDTxn3my+Pz4wA6L3m+93Rx6XKJh6/rlmNqZFzU1JlX0ixnGzibw4t/cop7e9y2BXH8WVYcYvnc2YRUkTrAHY+dQTDFyIzzIcuNDHzqeeYGHd4rRlBx2+mQ35yBz3Op4zb2bhnwsiBaP8PXcCxrixamLGi8+5vUfp2nII1R9NlB7s6qNryyGAlAvAmRPHMxpPRpA7LmnsY3o4RJvJ7yzf5k1LKmZYROZZ++eCSMEof8+dgDlwetu7w4rfQPUPcXrbuymV/7hJk6MmH5PxbAhqxyWNfeSaWGY79evMc3KS+OeCSMHY/N2wr+eLwS7zIlpW4yOpW3ELJaXxdVhKSsPUrbjFFtkKmYO7O3n8u7t47Jsv8fh3d3Fwd3a9cYPGjVUT2TB/JjPCIQSYEQ65m8uShX8uiBTMzt840f72UCudA8VMUse4WV7gahYTrUztH4orw6aKvrgydXEtw66/86knOHPiOOMmTaZuxS0Z2/s18XitUYfXyMa8mVcy9M8FkYJR/gBXs5MfDt7LkIpd7g3A/v0vAni+EcpIxi+dHWfzB5BQEeOXzk7r+QvrFmtlbzNea9ThBtkGIWjcoWDMPhBtczg0FB/iNTTUQ+vhDS5JlB1jFk2l8oa5wzv94sowlTfM9eYPrUCSaQq9zr0RhGBckRpBCOf2mhb61XiAgtr59/aZxxVbjXuZMYumelPZjySx2J2RTAOBu+T2WqMOp8klCMFNCrnIXUHt/MvC5s5dq3FNjhRQMk2h17nPJQjBLTyXeewwBaX8q+esoagoPvu0qKic6jlrXJIo4BRQMo3XGnU4jVWwQTpBCG7hucxjhykos4/h1G09vIHevg7KwhGq56zxlbPXVxRYMo2XGnU4Ta5BCG7g1cRPpygo5Q/RBUAre4fwYDJNU2sTm17fROe5TqrGVLHq8lUsq17mmjxBwbDr+ynaJ5vEzyCdPwWn/DUOYjh1k7S2dJKm1ibWv7ye3sFeADrOdbD+5fUAvv0Be4l8BCHks1lQppnHQTt/dCcvjXOk6HGcb5Y8s4SOc6PtuZExEbZ/abtjcmjSw2gWlFg48MEbPm7bApBJtI9fzp90O3npnb/GGTwQ9tl5zrzcgtW4xl2SNQuyS/lnknkctPMnp2gfEZkoIr8RkUOx2wlJ5o4XkXYReTSXY2p8igfCPqvGmDtjrcY17mLVFKi9q8eVrmBBO39yDfW8B2hWSs0FmmP3rXgA+F85Hk/jVzwQ9rnq8lWUFZfFjZUVl7Hq8lWOyaBJH8umQGXFrsTmB+38ydXssxy4Nvb/48Bvgb9OnCQinwQ+BPwbkNIWpQkgHgj7NJxy6URr5OJo1DVu7MGsWRDFwoW54+Lm2d6UyYJMzh8/kJPDV0S6lFKVI+6fUkpNSJhTBLwE/CegHqhVSt2R6rW1w9cb2Bbalmjzh2jYpwdL6ebiaExstAPReHfP1l7yOImL8LszwwxOGzNqngAdiz+R1TFaWlpobm6mu7ubiooK6uvrfd1tzTaHr4i8CJgZte5NU5ZvAy8opd4XkaQTRWQlsBJg1qxZab68Jl/YGtrmsbDPZOTiaPRrjRuvktgsqPblt7JqyrRv5w7TMuYtLS00NjbS3x99ze7ubhobGwF8vQCkQ0rlr5T6rNVjIvKBiESUUh0iEgHMSvhdBdSJyLeBsUCpiJxVSo3yDyilNgObIbrzT/dNaPLDptc3DSt+g97BXja9vim73b9PaqhbORqtxkfixxo3fsIsNr9kaJCbk2iyfTt3sH3zo8O9q88cP8b2zdG4k+Y9bwwrfoP+/n6am5sDr/xzdfg+D9wa+/9WYGviBKXUzUqpWUqp2cAa4Akzxa/xHtahbdFY5+7GRg5dV8++hR/l0HX1dMd2TH7kub3tXPPQS1x6TxNFFleo0yrLaWptYskzS6h5vIYlzyyhqbUpbo4fa9z4iRurJnJnuWJc73lQirG95/nM/tc5/5smWlpaTJ+z86knhhW/wcCFPnY+9QTd3d2mz7EaDxK5OnwfAp4WkW8AfwK+DCAitcA3lVJ/mePra1ykakyVaVJL1cAg3Y/dS8c/v4DqjV4ZDBw5Qsd90bINFQ0NjsqZK4k2/kETP1h5qJjPLupJaQY78KnJTN/eTvkIy09PEbR/ajJ+qh3r6VLHu3Zwc4Jy7gfL3fqZE8dNX+bMieNUfKTGVNFXVFTYIqqXyWnnr5Q6oZSqV0rNjd2ejI3vMVP8Sqmfp+Ps1XiDVZevoixBEZYNDbHq5CmO/suWYcVvoHp7ObrxESdFtAUzGz9AsQgCTK8s58EbPs4rZ6zNYAZrSnv424+F6SgThoCOMuFvPxZmTWlqk5FX8Hqp40x36+MmTbYcr6+vJxSK9xeEQiHq6+tzE9IH6AxfjSXLqpfBs3/FpgkVdJYUUzUwyKpTXSw7d559ZyuIxljEM9Dhv3K4Vrb8IaX440MXfRvrWqJmsMsv6ecvKgaYUKw4NSg0jQhhbe/rp21aKdumlca9lvioUmSyUsde2P1XVFRktFuvW3FLnM0foKQ0HHX6xq4UghTtky5a+WuSsqxkIsvaRsfnl4wVBs6Onl8S8ZNxI8q0ynLaTRaAxCSjqjFVRNSfWDGhH6Nvy8QSxYqJ/XR0biVStTyrSpFew+uljuvr6+MidCD5bt3oV20W7QPRqJ6amprhkM8tW7bQ3NzsyCJgFYXkBFr5a5JjUZZ56m03xNn8AaSsjKmr73JByNwwSyYqDxVz99L5cfNWXb6KnnfuIqFhFyFRtB7eQKRqecaVIr1GU2sTk49sQA0cZ6h4Eucqvkzf2GsA7yxgNVns1hfWLU6qVN0I+UwWheTEAqCVfxp0dG4t3AYwFvH5FTU3waxPcXTjIwx0dFASiTB19V2+c/YCw3HkqTJ6l1Uv48V37zR9DaMPtGEW8ayzNAlGXgeDvQhQPHiCcad+BkDR+E97agEzdut20dzc7HjIZ/PPN1tGIWnlbwPdjY05KaiOzq3s338vQ0PRnW9v3xH274/mtxXUAmASn1/R0OBLZW9GYjKRFeXhafT2HRk1PrIPdCaVIr2EWV6HqAuMP/0M6//8P/ryPaWL0yGf+3buoO/sGdPHrKKT7CbQPXy7GxvpuG8dA0eOgFLD4YiZxKO3Ht4wrPgNhoZ6aD28wW5xNT7AK32g9+3cwebbb+MHKxrYfPtt7Nu5I+fXtMrrkIETgVb8YO0szlfI586nnrB8zCo6yW4CrfyPbnwk53BE43I+3XFNsIlULWfBgu9RFp4GCGXhaSxY8D1HrwINW/GZ48dAqWFbca4LQNBKFmeC0yGfyXb3dStuycsxEwm02ccq7DCTcMSycCTlZb6msHCkD3SSrmfJMlaztRV3NzbylX89yz9cpbgQuhjC6+eSxZmQjRM5F8ZNmhxdvBMoGzdOR/vYQUkkEjX5mIwnYpXRWD1nTZzNH9y5zNcUECm6niXLWE2HxJLTJVNOcuL/XcfVvb0MnhL+x7VFnBgvTC2uZPXV9/i2ZHGm2O1EToZV7sF1t6505PgQcOU/dfVddNy3LmU4opHRaITnGRmNADfGdngFG+3jM2wrQe0mybqe1dxkuWtMx1acWHJ6sKuPgRMhiifVMND+O+reVtS9HQ15LZlWytz/ZP3ZBeKzdolUuQdOEGjlb0SipIr2SZXR6Mhlvl243CTdTWwtQe0mKbqeJctYTYVZyWkpLiX8sf+TgfbfxY0nM4/65bP2Uo2ixMjDaavvYuVj/+KKLBBw5Q/phSN6PaMxbTJokp5LpyqvYnsJardI0fUsl12jVWlpKR+tEJNla/vhs05+Re/sAmBEHnqpEGLglX86BCElH0hpLjBIrGLZ3tXD2i1vAvh6AbAuQW0+7lkssqqpXzd8N1XGqhXFlWHTBUD1noq7nypb2w+ftZdqFCWLPHRL+Qc61DNd1lZHKC+KL1Lmp5T8YdJskp6sU5WfCUyoYs1N0faWFTMBid7a1O5y/NLZSCj+Zy+hIsprwpRMmwYilEybRuSB+5MqJT981l66orcj8tBu9M4ff6fkx5Fmk/RcOlV5mVWXr4qzQ0PmoYoHd3fyytbDnD3Zx9iJYa5aPod5V7qg0EZkVZ/be5TTL7zL4C925twQ3nieWYP5Kd/4fNqvY8dnnW+8dEWfSeShU2jlH8OvKflxpGEugPSrWPoNw9acbQTKwd2d7HhyPwMXog7Rsyf72PHkfgB3FgDMo3O6thwCyGkByLWfcK6ftRN4qcheupGHTqKVPwyXcvV9Pe80m6SnW8XSjyyrXpa1Anpl6+FhxW8wcGGIV7Yedk3556shvB1hmrl81k7gpSv6dCMPnSRwyj/TQm5ulHLNK2k0SU+3imWhcfakeSSM1bgT5KMhvF/CNO3AS1f0XiuEGCjln004lRulXL1AulUsC4mxE8Omin7sROearyduXkJzGiiZ9MlR83JpCO+HME1N/skp2kdEJorIb0TkUOx2gsW8WSKyXUT2icjbIjI7l+NakU0hN6dLuWq8y1XL51CS0KmlpLSIq5bPceT4ZlVoe3f/nP4j8clXEipi/NLZWR/HD2GamvyTa6jnPUCzUmou0By7b8YTwMNKqYXAFcDRHI9rSjbhVE6XctV4l3lXVrH45gXDO/2xE8MsvnmBY/Z+083LhT76WxuHd/rFlWEqb5ibk73fD2GamvyTq9lnOXBt7P/Hgd8Cfz1ygoh8FChRSv0GQCll0vnVHrIJp8q0H6gm2My7sso1567VJmXwxFEi91xh23H8EKapyT+5Kv8PKaU6AJRSHSJith2ZB3SJyBbgUuBF4B6l1KDJ3JzIJpzK6VKuGmfxTNx+GjgVC+6HME1N/hGlVPIJIi8CZr+We4HHlVKVI+aeUkrF2f1F5EvAT4FFwJ+AXwIvKKV+anKslcBKgFmzZn3yvffey+zdkHvbRk1wSIzbh6gN38yU44UKlYkBCxDdvKTKttVoRiIiv1dK1aaal3Lnr5T6bJKDfCAikdiuP4K5Lb8N2KuUao095zngU0QXhMRjbQY2A9TW1iZflSzwWjiVxj3Sjdv3SuijF2PBnUZv3pwjV7PP88CtwEOx260mc14DJojIFKXUMeA6YE+Ox03Kvp07XK2TbTeBSUJzmHTj9r0U+ljImxcvVr4MMrlG+zwEfE5EDgGfi91HRGpF5CcAMdv+GqBZRN4EBPjnHI9rSb76m7qFkYRmhJ4aSWgtLS0uS+Z9rOLzE8ftDH1sam1iyTNLqHm8hiXPLKGptSnj1yhU7Oi5bReF8D3mpPyVUieUUvVKqbmx25Ox8T1Kqb8cMe83SqkapdTHlVJfV0pdyFVwK5L1N02ko3Mru3bV0fzSR9i1q46OTrMLF3dJloSmSU66cft2hT4a5qOOcx0o1LD5KIiKIx94pfJloXyPgSvpnG5/047Orezff2+sObuit+8I+/ff67kFQCehZU+6cfurLl9FWXFZ3Fg2oY/JzEea1FhFNTld+bJQvsdAlXcA0u5v2np4Q1xTdoChoR5aD2/wVMvGiooKU0Wvk9DSI524fbtCH4OYOeuk/8wrlS+D+D2aETjln25/094+80tJq3G30ElozmBHhcqqMVV0nBt9/vg1c9bwnxm/JcN/BuRlAfBKtFPQvkcrAqf80+1vWhaOxEw+jBr3EjoJzT8ELXM2mf8sHeWfTZSaF6KdgvY9WhE45Q/p9TetnrOG/fvvjTP9FBWVUz1nTb7Fy5iamhqt7H1A0DJn0/WfmeHlUumpsr6D9j1aEUjlnw6GXb/18AZ6+zooC0eonrPGU/Z+jf/weoOTTEjXf2aGV0ulp9utLUjfoxUFq/whugBoZe9/sk2CS+XMfG5ve0E3vEnXf2aGV6PUPNWtreXplF338klBK3+N/8nWvJDKmfnc3va4VpftXT2s3fImQMEsAOn6z8xwO0rt2c6Tpu0bPdOtreXp+H7b3e9H74NjC4BW/hpfk615IZUz8+FtB+J6HAP09A/y8LYDBaP8IT3/mRluRqk923kyrnF7W18/aw68D3ijWxsQ3fH3x4ea098THXdI+QcuyUtTWGRrXkjlzDzS1WP6uNW4Jp6amhoaGhqGd/oVFRU0NDQ4Yu9/sLVjWPEb9AwpHmztcL1b2zDdbZmN5wG989f4mmzNC6mcmdMqy2k3UfTTKsuzlLTwcCtKrb2v33J83uKoXd/1Hg8VM6KmHrNxh9A7f42vqa+vJxQKxY2lY16oW3ELJaXxl/ojnZl3L51Peag47vHyUDF3L51vg9SafDI9HEo6Pu/KKm79/jXc/o/Xcev3r3GnuU/9OgglbCRC5dFxh9DKX+NrsjUvLKxbzJKVdzBu8hQQYdzkKSxZecewffv6RdN58IaPM72yHAGmV5bz4A0fLyh7v19ZWx2hvEjixsqLhLXVHkrgrLkJGn4EFTMBid42/MjRaJ+Unbzcora2Vu3Zk9ey/xqNJqBYRfsUArZ18tJoNBq/cWPVxIJR9tmizT4ajUZTgGjlr9FoNAWIVv4ajUZTgGjlr9FoNAWIVv4ajUZTgHg21FNEjgHvuS1HApOB1MXMvYGfZAV/yatlzR9+ktersn5YKTUl1STPKn8vIiJ70omf9QJ+khX8Ja+WNX/4SV4/yWqGNvtoNBpNAaKVv0aj0RQgWvlnxma3BcgAP8kK/pJXy5o//CSvn2Qdhbb5azQaTQGid/4ajUZTgGjlnwQRmSgivxGRQ7HbCRbzZonIdhHZJyJvi8hsZyVNX9bY3PEi0i4ijzopY4IMKeUVkU+IyCsi8paItIjIVxyW8fMickBE3hGRe0weD4vIL2OP73bjex8hSypZvxM7N1tEpFlEPuyGnDFZkso6Yt6XRESJiKsRNenIKyI3xT7ft0TkF07LmBVKKf1n8Qf8PXBP7P97gL+zmPdb4HOx/8cCl3hV1tjjm4BfAI96+bMF5gFzY/9PAzqASofkKwYOA9VAKfD/AR9NmPNt4B9j/68AfunSZ5mOrIuN8xL4lpdljc0bB/w78CpQ64asGXy2c4G9wITY/aluyZvJn975J2c58Hjs/8eB6xMniMhHgRKl1G8AlFJnlVLnnRNxmJSyAojIJ4EPAdsdksuKlPIqpQ4qpQ7F/j8CHAVSJq/YxBXAO0qpVqXUBeApojKPZOR7eAaoFxHBeVLKqpTaMeK8fBVwrl9gPOl8rgAPEN0g9DopnAnpyPtXwGNKqVMASqmjDsuYFVr5J+dDSqkOgNjtVJM584AuEdkiIntF5GERKTaZl29SyioiRcAPgLsdls2MdD7bYUTkCqI7r8MOyAYwHRjZZLUtNmY6Ryk1AHQDkxyRzkKOGGayjuQbwL/mVSJrUsoqIouAmUqp/+mkYBak89nOA+aJyC4ReVVEPu+YdDlQ8M1cRORFwKyJ571pvkQJUAcsAv4E/BL4OvBTO+QbiQ2yfht4QSn1vhMbVBvkNV4nAvw34Fal1JAdsqVzWJOxxNC4dOY4QdpyiMjXgFrgM3mVyJqksq1ZNNsAAAICSURBVMY2KBuJ/oa8QDqfbQlR08+1RK+odorIZUqprjzLlhMFr/yVUp+1ekxEPhCRiFKqI6aAzC7n2oC9SqnW2HOeAz5FHpS/DbJeBdSJyLeJ+iZKReSsUsrS6eayvIjIeKAJ+K9KqVfzIacFbcDMEfdnAEcs5rSJSAlQAZx0RjxTOQzMZEVEPkt04f2MUqrPIdkSSSXrOOAy4LexDUoV8LyIfFEp5UZf13TPg1eVUv3AH0XkANHF4DVnRMwObfZJzvPArbH/bwW2msx5DZggIoYt+jrgbQdkSySlrEqpm5VSs5RSs4E1wBP5UvxpkFJeESkFfk1Uzl85KBtEv9e5InJpTI4VRGUeycj38CXgJRXz+DlMSlljppR/Ar7osk06qaxKqW6l1GSl1OzYefoqUZndauidznnwHFGHOiIymagZqNVRKbPBbY+zl/+I2m+bgUOx24mx8VrgJyPmfQ5oAd4Efg6UelXWEfO/jrvRPinlBb4G9ANvjPj7hIMyfgE4SNTPcG9s7H6iygigDPgV8A7wO6Daxc8zlawvAh+M+Byf96qsCXN/i4vRPml+tgL8kOim701ghZvypvunM3w1Go2mANFmH41GoylAtPLXaDSaAkQrf41GoylAtPLXaDSaAkQrf41GoylAtPLXaDSaAkQrf41GoylAtPLXaDSaAuT/B5x9m/6phMg8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "random(200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, atoms are not truely randomly positioned. Two atoms can never be very close to each other. One way to solve this is to consider a lattice (grid) of position initial positions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGfCAYAAAAH0zaSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAERZJREFUeJzt3F+I3Xl5x/HP0yTF2W0lFqM0WWEVJNsilLVDa7sg4rbEVtG9aEFBsVLYm9Y/pURMb7wVIqW9KELwX6FiKduwFSlNRSulUJZONkJWY1Cs1U223UhJW8pQ4/rtRSbixsnmzDlnc36P5/WCJTO/zJnfw3lO9j3nz5waYwQAuviJVQ8AAHshXAC0IlwAtCJcALQiXAC0IlwAtHLbcFXVx6vq6ap64oeO/UxVfa6qvrbz54ue3zEB4LpZ7nF9Mskbbjr2gSSfH2O8Msnndz4HgOddzfILyFV1b5LPjjFetfP5xSSvG2M8VVU/m+SLY4yjz+egAJAk++e83EvHGE8lyU68XnKrL6yqh5M8nCR33333L953331znhKAHzdnz579zhjj0F4uM2+4ZjbGOJXkVJJsbm6Ora2t5/uUADRRVf+218vM+6rC/9h5iDA7fz495/cBgD2ZN1yfSfLOnY/fmeRvljMOADy3WV4O/+kk/5zkaFU9WVW/m+RDSX69qr6W5Nd3PgeA591tn+MaY7ztFn/14JJnAYDb8s4ZALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtLJ/1QN08ei5Szl55mIuX93O4YMbOX7saB66/8iqx1ordjAN9jAN67wH4ZrBo+cu5cTp89m+9kyS5NLV7Zw4fT5J1uaGsmp2MA32MA3rvgcPFc7g5JmLP7iB3LB97ZmcPHNxRROtHzuYBnuYhnXfg3DN4PLV7T0dZ/nsYBrsYRrWfQ/CNYPDBzf2dJzls4NpsIdpWPc9CNcMjh87mo0D+551bOPAvhw/dnRFE60fO5gGe5iGdd+DF2fM4MaTnev6Cp4psINpsIdpWPc91Bjjjp1sc3NzbG1t3bHzATBtVXV2jLG5l8t4qBCAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFYWCldV/UFVfbmqnqiqT1fVC5Y1GADsZu5wVdWRJO9JsjnGeFWSfUneuqzBAGA3iz5UuD/JRlXtT3JXksuLjwQAtzZ3uMYYl5J8OMm3kjyV5L/GGH9/89dV1cNVtVVVW1euXJl/UgDIYg8VvijJW5K8PMnhJHdX1dtv/roxxqkxxuYYY/PQoUPzTwoAWeyhwl9L8q9jjCtjjGtJTif51eWMBQC7WyRc30rymqq6q6oqyYNJLixnLADY3SLPcT2W5JEkjyc5v/O9Ti1pLgDY1f5FLjzG+GCSDy5pFgC4Le+cAUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUAr+1c9QBePnruUk2cu5vLV7Rw+uJHjx47mofuPrHqstWIH02AP07DOexCuGTx67lJOnD6f7WvPJEkuXd3OidPnk2RtbiirZgfTYA/TsO578FDhDE6eufiDG8gN29eeyckzF1c00fqxg2mwh2lY9z0I1wwuX93e03GWzw6mwR6mYd33IFwzOHxwY0/HWT47mAZ7mIZ134NwzeD4saPZOLDvWcc2DuzL8WNHVzTR+rGDabCHaVj3PXhxxgxuPNm5rq/gmQI7mAZ7mIZ130ONMe7YyTY3N8fW1tYdOx8A01ZVZ8cYm3u5jIcKAWhloXBV1cGqeqSqvlpVF6rqV5Y1GADsZtHnuP40yd+NMX6rqn4yyV1LmAkAbmnucFXVC5O8NsnvJMkY47tJvrucsQBgd4s8VPiKJFeSfKKqzlXVR6vq7pu/qKoerqqtqtq6cuXKAqcDgMXCtT/Jq5N8ZIxxf5L/TfKBm79ojHFqjLE5xtg8dOjQAqcDgMXC9WSSJ8cYj+18/kiuhwwAnjdzh2uM8e9Jvl1VN35V+8EkX1nKVABwC4u+qvDdST6184rCbyR51+IjAcCtLRSuMcaXkuzpN54BYBHeOQOAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFYWDldV7auqc1X12WUMBADPZRn3uN6b5MISvg8A3NZC4aqqe5K8MclHlzMOADy3Re9x/UmS9yf5/q2+oKoerqqtqtq6cuXKgqcDYN3NHa6qelOSp8cYZ5/r68YYp8YYm2OMzUOHDs17OgBIstg9rgeSvLmqvpnkL5O8vqr+YilTAcAtzB2uMcaJMcY9Y4x7k7w1yRfGGG9f2mQAsAu/xwVAK/uX8U3GGF9M8sVlfC8AeC5LCdc6ePTcpZw8czGXr27n8MGNHD92NA/df2TVY60VO5gGe5iGdd6DcM3g0XOXcuL0+WxfeyZJcunqdk6cPp8ka3NDWTU7mAZ7mIZ134PnuGZw8szFH9xAbti+9kxOnrm4oonWjx1Mgz1Mw7rvQbhmcPnq9p6Os3x2MA32MA3rvgfhmsHhgxt7Os7y2cE02MM0rPsehGsGx48dzcaBfc86tnFgX44fO7qiidaPHUyDPUzDuu/BizNmcOPJznV9Bc8U2ME02MM0rPseaoxxx062ubk5tra27tj5AJi2qjo7xtjcy2U8VAhAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCvCBUArwgVAK8IFQCtzh6uqXlZV/1BVF6rqy1X13mUOBgC72b/AZb+X5A/HGI9X1U8nOVtVnxtjfGVJswHAj5j7HtcY46kxxuM7H/9PkgtJjixrMADYzVKe46qqe5Pcn+SxXf7u4araqqqtK1euLON0AKyxhcNVVT+V5K+TvG+M8d83//0Y49QYY3OMsXno0KFFTwfAmlsoXFV1INej9akxxunljAQAt7bIqworyceSXBhj/PHyRgKAW1vkHtcDSd6R5PVV9aWd/35zSXMBwK7mfjn8GOOfktQSZwGA2/LOGQC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0IlwAtCJcALQiXAC0sn/VA3Tx6LlLOXnmYi5f3c7hgxs5fuxoHrr/yKrHWit2MA32MA3rvAfhmsGj5y7lxOnz2b72TJLk0tXtnDh9PknW5oayanYwDfYwDeu+Bw8VzuDkmYs/uIHcsH3tmZw8c3FFE60fO5gGe5iGdd+DcM3g8tXtPR1n+exgGuxhGtZ9D8I1g8MHN/Z0nOWzg2mwh2lY9z0I1wyOHzuajQP7nnVs48C+HD92dEUTrR87mAZ7mIZ134MXZ8zgxpOd6/oKnimwg2mwh2lY9z3UGOOOnWxzc3NsbW3dsfMBMG1VdXaMsbmXy3ioEIBWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWhAuAVoQLgFaEC4BWFgpXVb2hqi5W1der6gPLGgoAbmXucFXVviR/luQ3kvx8krdV1c8vazAA2M0i97h+KcnXxxjfGGN8N8lfJnnLcsYCgN3tX+CyR5J8+4c+fzLJL9/8RVX1cJKHdz79v6p6YoFzrqsXJ/nOqodoyPU2H9fbfFxv8zm61wssEq7a5dj4kQNjnEpyKkmqamuMsbnAOdeS620+rrf5uN7m43qbT1Vt7fUyizxU+GSSl/3Q5/ckubzA9wOA21okXP+S5JVV9fKq+skkb03ymeWMBQC7m/uhwjHG96rq95OcSbIvycfHGF++zcVOzXu+Ned6m4/rbT6ut/m43uaz5+utxviRp6UAYLK8cwYArQgXAK3ckXB5a6j5VNXLquofqupCVX25qt676pm6qKp9VXWuqj676lk6qaqDVfVIVX1153b3K6ueqYOq+oOdf6NPVNWnq+oFq55piqrq41X19A//Pm9V/UxVfa6qvrbz54tu932e93B5a6iFfC/JH44xfi7Ja5L8nutuZu9NcmHVQzT0p0n+boxxX5JfiOvwtqrqSJL3JNkcY7wq11+s9tbVTjVZn0zyhpuOfSDJ58cYr0zy+Z3Pn9OduMflraHmNMZ4aozx+M7H/5Pr/xM5stqppq+q7knyxiQfXfUsnVTVC5O8NsnHkmSM8d0xxtXVTtXG/iQbVbU/yV3xO627GmP8Y5L/vOnwW5L8+c7Hf57kodt9nzsRrt3eGsr/fPeoqu5Ncn+Sx1Y7SQt/kuT9Sb6/6kGaeUWSK0k+sfMw60er6u5VDzV1Y4xLST6c5FtJnkryX2OMv1/tVK28dIzxVHL9h/UkL7ndBe5EuGZ6ayhurap+KslfJ3nfGOO/Vz3PlFXVm5I8PcY4u+pZGtqf5NVJPjLGuD/J/2aGh23W3c5zMm9J8vIkh5PcXVVvX+1UP97uRLi8NdQCqupArkfrU2OM06uep4EHkry5qr6Z6w9Lv76q/mK1I7XxZJInxxg37tU/kush47n9WpJ/HWNcGWNcS3I6ya+ueKZO/qOqfjZJdv58+nYXuBPh8tZQc6qqyvXnGy6MMf541fN0MMY4Mca4Z4xxb67f1r4wxvDT7wzGGP+e5NtVdePduh9M8pUVjtTFt5K8pqru2vk3+2C8qGUvPpPknTsfvzPJ39zuAou8O/xM5nxrKK57IMk7kpyvqi/tHPujMcbfrnAmfry9O8mndn7I/EaSd614nskbYzxWVY8keTzXXwl8Lt7+aVdV9ekkr0vy4qp6MskHk3woyV9V1e/m+g8Bv33b7+MtnwDoxDtnANCKcAHQinAB0IpwAdCKcAHQinAB0IpwAdDK/wOvQmVSQ2HfwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lattice(N, L):\n",
    "    # organize the particles in to a grid\n",
    "    sqn=np.sqrt(N)\n",
    "    if (sqn == int(sqn)):\n",
    "        grid = 0.5*L/sqn\n",
    "    else: \n",
    "        grid = 0.5*L / int(sqn + 1)\n",
    "       \n",
    "    x = np.zeros(N)\n",
    "    y = np.zeros(N)\n",
    "    x[0] = grid\n",
    "    y[0] = grid\n",
    "    for i in range(N-1):\n",
    "        \n",
    "        if x[i+1] < L:\n",
    "            x[i+1] = x[i] + 2*grid\n",
    "            y[i+1] = y[i]\n",
    "        if x[i+1] >= L:\n",
    "            x[i+1] = grid\n",
    "            y[i+1] = y[i]+2*grid\n",
    "    \n",
    "    return x, y\n",
    "   \n",
    "N = 16\n",
    "L = 10\n",
    "x, y = lattice(N, L)\n",
    "fig, axes = plt.subplots(figsize=(7,7))\n",
    "plt.plot(x,y, 'o')\n",
    "plt.xlim(0, L)\n",
    "plt.ylim(0, L)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also can perturb each particle slightly from its grid position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_lattice(n, L):\n",
    "    # organize the particles in to a grid\n",
    "    sqn=np.sqrt(n)\n",
    "    if (sqn == int(sqn)):\n",
    "        grid = 0.5*L/sqn\n",
    "    else: \n",
    "        grid = 0.5*L / int(sqn + 1)\n",
    "       \n",
    "    # random perturbation to gridded position\n",
    "    drandom = lambda: 0.25*grid*(np.random.rand() - 0.5) * np.sqrt(2)\n",
    "    \n",
    "    x = np.zeros(n)\n",
    "    y = np.zeros(n)\n",
    "    x[0] = grid + drandom()\n",
    "    y[0] = grid + drandom()\n",
    "    for i in range(n-1):\n",
    "        \n",
    "        if x[i+1] < L:\n",
    "            x[i+1] = x[i] + 2*grid + drandom()\n",
    "            y[i+1] = y[i] + drandom()\n",
    "        if x[i+1] >= L:\n",
    "            x[i+1] = grid + drandom()\n",
    "            y[i+1] = y[i]+2*grid + drandom()\n",
    "    \n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAGfCAYAAADLSgx5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+UXWV97/H3N8nEDD9MEKJCkgr2tkELQ4hTpMKFeGMJUn6klqLeaKPUlYvVS3uvUqB2RVZcrfbGFi+1lYWooCtSUg1BrN7IolYaW5UJxIQfTUEayyQUAmkiKZOaxO/945wJk8k58+vMzDmT5/1aa9Y559m/nn2evc9n9rP32ScyE0mSSjCp2RWQJGm8GHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYgwaehExJyK+HRGPRcQjEfG71fJXRMS9EfF49fG4OtMvrY7zeEQsHe0VkCRpqGKw7+lFxInAiZn5YEQcC2wAFgPvAXZm5ici4jrguMy8tt+0rwC6gE4gq9O+ITP/fdTXRJKkQQx6pJeZT2fmg9XnLwCPAbOAy4Dbq6PdTiUI+1sE3JuZO6tBdy9w4WhUXJKk4ZoynJEj4mTgTOD7wKsy82moBGNEvLLGJLOAp/q87q6W1Zr3MmAZwNFHH/2GU089dThVkyQdwTZs2PBcZs5sdD5DDr2IOAb4KvB7mfmTiBjSZDXKavanZuYtwC0AnZ2d2dXVNdSqSZKOcBHx49GYz5Cu3oyINiqBtyoz11SLn6me7+s97/dsjUm7gTl9Xs8Gto+8upIkjdxQrt4M4HPAY5n5Z30GfQ3ovRpzKXB3jcnXARdExHHVqzsvqJZJkjTuhnKkdw7wbuC/RcTG6t9FwCeAX42Ix4Ffrb4mIjoj4laAzNwJfAx4oPq3olomSdK4G/QrC83gOT1JY2Hfvn10d3ezd+/eZldFdUybNo3Zs2fT1tZ2SHlEbMjMzkbnP6yrNyVpIuvu7ubYY4/l5JNPZogX42kcZSbPP/883d3dnHLKKWOyDG9DJqkYe/fu5fjjjzfwWlREcPzxx4/pkbihJ6koBl5rG+v2MfQkScUw9CSpxdx888188YtfPKx869atnHbaaU2o0ZHDC1kkqY61D21j5botbN/Vw0kz2rlm0VwWn1nzToqjZv/+/Vx11VVjuoySGXqSVMPah7Zx/ZrN9Ow7AMC2XT1cv2YzQEPB97GPfYxVq1YxZ84cTjjhBN7whjfw9a9/nTe96U1897vf5dJLL+WFF17gmGOO4cMf/jAbNmzgyiuv5KijjuLcc88dlXUrmd2bklTDynVbDgZer559B1i5bsuI59nV1cVXv/pVHnroIdasWUPf7yPv2rWL73znO3zoQx86ZJr3vve93HTTTfzjP/7jiJerlxh6klTD9l09wyofivXr13PZZZfR3t7OscceyyWXXHJw2Nvf/vbDxt+9eze7du3i/PPPB+Dd7373iJetCkNPkmo4aUb7sMqHYqA7YB199NE1x/crFqPL0JOkGq5ZNJf2tsmHlLW3TeaaRXNHPM9zzz2Xe+65h71797Jnzx7+5m/+ZsDxZ8yYwfTp01m/fj0Aq1atGvGyVeGFLJJUQ+/FKqN59eYv//Ivc+mll3LGGWfwmte8hs7OTqZPnz7gNF/4whcOXsiyaNGiES9bFd5wWlIxHnvsMV73utc1tQ579uzhmGOO4cUXX+S8887jlltuYf78+U2tU6up1U7ecFqSJqBly5bx6KOPsnfvXpYuXWrgjTNDT5LG0Ze//OVmV6FoXsgiSSqGoSdJKoahJ0kqhqEnSSqGoSdJ4+hNb3rToOO8733v49FHHwXgj//4j4c9/THHHDOyyg1B3589uu2229i+ffvBYX3r3ar8np6kYgz7e3qbVsN9K2B3N0yfDQuXQ8cVY1fBGo455hj27Nkz5tOMxIIFC/jkJz9JZ2fDX587xFh+T88jPUmqZdNquOdq2P0UkJXHe66ulDeg9yjs7/7u71iwYAGXX345p556KkuWLDl4b84FCxbQ1dXFddddR09PD/PmzWPJkiWHTL9nzx4WLlzI/PnzOf3007n77rsHXO7WrVs59dRTWbp0KR0dHVx++eW8+OKLANx3332ceeaZnH766Vx55ZX853/+JwDXXXcdr3/96+no6ODDH/4wADfccAOf/OQn+cpXvkJXVxdLlixh3rx59PT0HKw3wB133MHpp5/OaaedxrXXXnvI+n/kIx/hjDPO4Oyzz+aZZ55p6P0cLkNPkmq5bwXs6/eLCvt6KuWj5KGHHuJTn/oUjz76KE8++STf/e53Dxn+iU98gvb2djZu3HjYfTenTZvGXXfdxYMPPsi3v/1tPvShDw14Q2uALVu2sGzZMjZt2sTLX/5y/vIv/5K9e/fynve8hzvvvJPNmzezf/9+PvOZz7Bz507uuusuHnnkETZt2sQf/uEfHjKvyy+/nM7OTlatWsXGjRtpb3/pRtzbt2/n2muv5W//9m/ZuHEjDzzwAGvXrgXgP/7jPzj77LP54Q9/yHnnncdnP/vZRt7CYTP0JKmW3d3DKx+Bs846i9mzZzNp0iTmzZvH1q1bhzxtZvIHf/AHdHR08Ja3vIVt27YNetQ0Z84czjnnHADe9a53sX79erZs2cIpp5zCL/7iLwKwdOlS7r//fl7+8pczbdo03ve+97FmzRqOOuqoIdftgQceYMGCBcycOZMpU6awZMkS7r//fgCmTp3KxRdfDMAb3vCGYa3zaDD0JKmW6bOHVz4CL3vZyw4+nzx5Mvv37x/ytKtWrWLHjh1s2LCBjRs38qpXvYq9e/cOOE3/nymKiLpHh1OmTOEHP/gBv/Ebv8HatWu58MILh1y3gY4429raDtZjuOs8Ggw9Sapl4XJo6/fbeW3tlfJx1NbWxr59+w4r3717N6985Stpa2vj29/+Nj/+8Y8Hnde//uu/HvwF9jvuuINzzz2XU089la1bt/LEE08A8KUvfYnzzz+fPXv2sHv3bi666CI+9alPsXHjxsPmd+yxx/LCCy8cVv7GN76R73znOzz33HMcOHCAO+644+AP4TaboSdJtXRcAZfcBNPnAFF5vOSmcb96c9myZXR0dBy8kKXXkiVL6OrqOnhe7dRTTx10Xq973eu4/fbb6ejoYOfOnbz//e9n2rRpfOELX+A3f/M3Of3005k0aRJXXXUVL7zwAhdffDEdHR2cf/753HjjjYfN7z3veQ9XXXXVwQtZep144ol8/OMf581vfjNnnHEG8+fP57LLLmv8zRgFfmVBUjFa4aeFmmXr1q1cfPHFPPzww82uyqD8yoIkSaPA0JOkApx88skT4ihvrBl6korSiqd09JKxbh9DT1Ixpk2bxvPPP2/wtajM5Pnnn2fatGljtgx/OV1SMWbPnk13dzc7duxodlVUx7Rp05g9e/S+C9mfoSepGG1tbZxyyinNroaayO5NSVIxDD1JUjEMPUlSMQw9SVIxDD1JUjEGvXozIj4PXAw8m5mnVcvuBOZWR5kB7MrMeTWm3Qq8ABwA9o/GfdNa2dqHtrFy3Ra27+rhpBntXLNoLovPnNXsahXHdmgNtkPz2QaHG8pXFm4DPg18sbcgM9/e+zwi/hTYPcD0b87M50ZawYli7UPbuH7NZnr2HQBg264erl+zGaD4jWw82Q6twXZoPtugtkG7NzPzfmBnrWFR+SXAK4A7RrleE87KdVsObly9evYdYOW6LU2qUZlsh9ZgOzSfbVBbo+f0/ivwTGY+Xmd4At+KiA0RsWygGUXEsojoioiuiXi3hO27eoZVrrFhO7QG26H5bIPaGg29dzLwUd45mTkfeCvwgYg4r96ImXlLZnZmZufMmTMbrNb4O2lG+7DKNTZsh9ZgOzSfbVDbiEMvIqYAbwPurDdOZm6vPj4L3AWcNdLltbprFs2lvW3yIWXtbZO5ZtHcOlNoLNgOrcF2aD7boLZG7r35FuCfMrO71sCIOBqYlJkvVJ9fAKxoYHktrffEsFdKNZft0Bpsh+azDWqLwX5iIyLuABYAJwDPAB/NzM9FxG3A9zLz5j7jngTcmpkXRcRrqRzdQSVcv5yZfzSUSnV2dmZXV9dw10Wq2LQa7lsBu7th+mxYuBw6rmh2rSQ1ICI2jMbX3gY90svMd9Ypf0+Nsu3ARdXnTwJnNFg/aXg2rYZ7roZ91ZP1u5+qvAaDT5J3ZNER5r4VLwVer309lXJJxTP0RmLTarjxNLhhRuVx0+pm10i9dtc8xVy/XFJRDL3h6u0+2/0UkC91nxl8rWF6nV9crlcuqSiG3nDZfdbaFi6Htn7fQ2prr5RLJbFHqqZGvrJQJrvPWlvvxSpevamSeUFXXYbecE2fXe3arFGu1tBxRfE7tgo3UI9U4fuG3ZvDZfeZNHR2sTWHPVJ1GXrD1XEFXHITTJ8DROXxkpuK/+9JOowXfTWPF3TVZffmSNh9Jg3OLrbmWbj80HN6YI9UlUd6ksaGXWzNY49UXR7pSRobXvTVXPZI1eSRnqSx4UVfakGGnqSxYRebWpDdm5LGjl1sajEe6UmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkoph6EmSimHoSZKKYehJkooxaOhFxOcj4tmIeLhP2Q0RsS0iNlb/Lqoz7YURsSUinoiI60az4pIkDddQjvRuAy6sUX5jZs6r/n2j/8CImAz8BfBW4PXAOyPi9Y1UVpKkRgwaepl5P7BzBPM+C3giM5/MzJ8CfwVcNoL5SJI0Kho5p/fBiNhU7f48rsbwWcBTfV53V8tqiohlEdEVEV07duxooFqSJNU20tD7DPDzwDzgaeBPa4wTNcqy3gwz85bM7MzMzpkzZ46wWpIk1Tei0MvMZzLzQGb+DPgsla7M/rqBOX1ezwa2j2R5kiSNhhGFXkSc2OflrwMP1xjtAeAXIuKUiJgKvAP42kiWJ0nSaJgy2AgRcQewADghIrqBjwILImIele7KrcD/qI57EnBrZl6Umfsj4oPAOmAy8PnMfGRM1kKSpCGIzLqn2Zqms7Mzu7q6ml0NSVKLiIgNmdnZ6Hy8I4skqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaeJKkYhp4kqRiGniSpGIaejjybVsONp8ENMyqPm1Y3u0aSWsSUwUaIiM8DFwPPZuZp1bKVwCXAT4EfAe/NzF01pt0KvAAcAPZnZufoVb31rH1oGyvXbWH7rh5OmtHONYvmsvjMWc2uVlk2rWb/3f+TKQf2Vl7vfqryGqDjimbWrDjuD81nGxxuKEd6twEX9iu7FzgtMzuAfwauH2D6N2fmvBIC7/o1m9m2q4cEtu3q4fo1m1n70LZmV60oL35z+UuBVzXlwF5e/ObyJtWoTO4PzWcb1DZo6GXm/cDOfmXfysz91ZffA2aPQd0mlJXrttCz78AhZT37DrBy3ZYm1ahM03r+bVjlGhvuD81nG9Q2Guf0rgS+WWdYAt+KiA0RsWygmUTEsojoioiuHTt2jEK1xtf2XT0Dl3ueaVxs/9nxwyrX2Bh0f9CYsw1qayj0IuIjwH5gVZ1RzsnM+cBbgQ9ExHn15pWZt2RmZ2Z2zpw5s5FqNcVJM9rrl29aDfdcDbufArLyeM/VBt8YuHXqu3gxpx5S9mJO5dap72pSjco04P6gcWEb1Dbi0IuIpVQucFmSmVlrnMzcXn18FrgLOGuky2t11yyaS3vb5EPK2tsmc82iuXDfCtjX77+rfT2Vco2qeb+2jOW5jO6fncDPMuj+2Qksz2XM+7UBOxo0ygbcHzQubIPaBr16s5aIuBC4Fjg/M1+sM87RwKTMfKH6/ALgiP2U770iquaVUnd3155od51yjVilHX6Ht69b6BVrTTTg/qBxYRvUFnUO0l4aIeIOYAFwAvAM8FEqV2u+DHi+Otr3MvOqiDgJuDUzL4qI11I5uoNKuH45M/9oKJXq7OzMrq6u4a5L67rxtGrXZj/T58D/enj86yNJE0xEbBiNbwEMeqSXme+sUfy5OuNuBy6qPn8SOKOh2h0pFi6vnMPr28XZ1l4plya6TasrXfW7u2H67Mp27Xci1aJG1L2pYer9APCDQUea3ou0ev+h671IC9y+1ZIMvfHScYUfAjryDHSRltv7+POoe1CGnqSRq3cxlhdpjT+PuofEG05LGrnpdW7GVK9cY8evRg2JoSdp5BYur1yU1ZcXaTWHR91DYuhJGrmOK+CSmypfvyEqj5fcZHdaM3jUPSSe05PUGC/Sag1+NWpIPNKTpCOBR91D4pGeJB0pPOoelEd6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIYepKkYhh6kqRiGHqSpGIMKfQi4vMR8WxEPNyn7BURcW9EPF59PK7OtEur4zweEUtHq+KSJA3XUI/0bgMu7Fd2HXBfZv4CcF/19SEi4hXAR4E3AmcBH60XjpIkjbUhhV5m3g/s7Fd8GXB79fntwOIaky4C7s3MnZn578C9HB6ekiSNi0bO6b0qM58GqD6+ssY4s4Cn+rzurpZJkjTuxvpClqhRljVHjFgWEV0R0bVjx44xrpYkqUSNhN4zEXEiQPXx2RrjdANz+ryeDWyvNbPMvCUzOzOzc+bMmQ1US5Kk2hoJva8BvVdjLgXurjHOOuCCiDiuegHLBdUySZLG3VC/snAH8I/A3IjojojfBj4B/GpEPA78avU1EdEZEbcCZOZO4GPAA9W/FdUySZLGXWTWPMXWVJ2dndnV1dXsakiSWkREbMjMzkbn4x1ZJEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxZgy0gkjYi5wZ5+i1wLLM/NTfcZZANwN/Eu1aE1mrhjpMlvd2oe2sXLdFrbv6uGkGe1cs2gui8+c1exqFcd2aA22Q/PZBocbcehl5hZgHkBETAa2AXfVGPXvM/PikS5nolj70DauX7OZnn0HANi2q4fr12wGKH4jG0+2Q2uwHZrPNqhttLo3FwI/yswfj9L8JpyV67Yc3Lh69ew7wMp1W2DTarjxNLhhRuVx0+om1fLI17cdLp20nvVTr+aRSW/n7LvP930fRwPuDxoXQ2qDAj+bRnyk1887gDvqDPuViPghsB34cGY+UmukiFgGLAP4uZ/7uVGq1vjZvqunZnnnT+6Fe74A+6rDdz8F91xded5xxTjVrhy97XDppPV8ou1WjoqfAvBqdvi+j6N6+0O9co2+Qdtg0+rKPlHYZ1PDR3oRMRW4FPjrGoMfBF6TmWcAfw6srTefzLwlMzszs3PmzJmNVmvcnTSjvWb59VP/+qWNqte+HrjviD212VS97fD7U1YfDLyDfN/HTb394WB5gUcY423QNrhvRZGfTaPRvflW4MHMfKb/gMz8SWbuqT7/BtAWESeMwjJbzjWL5tLeNvmQsva2ybyK52pPsLt7HGpVnt52OCl835up3v5wzaK5Lx1h7H4KyJeOMAy+UTVgG0D9feEI30dGI/TeSZ2uzYh4dURE9flZ1eU9PwrLbDmLz5zFx992OrNmtBPArBntfPxtpxPTZ9eeoF65GtLbDs9Gnd4C3/dxUW9/WHzmrGKPMMbbgG0A9feFI3wficwc+cQRRwFPAa/NzN3VsqsAMvPmiPgg8H5gP9AD/O/M/IfB5tvZ2ZldXV0jrldL6d9vDtDWDpfcdET3mzed73vrumEGUOtzJ+CGXeNdm3JNsH0kIjZkZmej82noQpbMfBE4vl/ZzX2efxr4dCPLmPB6N577VlS6DabPhoXLW3KjOqL4vreu6bOrXZs1yjV+Ct1HGjrSGytH1JGepENNsCMMtYbROtLzNmSSxlfHFZWAmz4HiMqjgadxMlrf05Okoeu4wpBTU3ikJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKkbDoRcRWyNic0RsjIiuGsMjIm6KiCciYlNEzG90mZIkjcSUUZrPmzPzuTrD3gr8QvXvjcBnqo+SJI2r8ejevAz4YlZ8D5gRESeOw3IlSTrEaIReAt+KiA0RsazG8FnAU31ed1fLDhERyyKiKyK6duzYMQrVkiTpUKMReudk5nwq3ZgfiIjz+g2PGtPkYQWZt2RmZ2Z2zpw5cxSqJUnSoRoOvczcXn18FrgLOKvfKN3AnD6vZwPbG12uJEnD1VDoRcTREXFs73PgAuDhfqN9Dfit6lWcZwO7M/PpRpYrSdJINHr15quAuyKid15fzsz/FxFXAWTmzcA3gIuAJ4AXgfc2uExJkkakodDLzCeBM2qU39zneQIfaGQ5kiSNBu/IIkkqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSqGoSdJKoahJ0kqhqEnSSrGlJFOGBFzgC8CrwZ+BtySmf+33zgLgLuBf6kWrcnMFSNdZqtb+9A2Vq7bwvZdPZw0o51rFs1l8Zmzml2t4tgOrcF2aD7b4HAjDj1gP/ChzHwwIo4FNkTEvZn5aL/x/j4zL25gORPC2oe2cf2azfTsOwDAtl09XL9mM0DxG9l4sh1ag+3QfLZBbSPu3szMpzPzwerzF4DHgGLfyZXrthzcuHr17DvAynVbmlSjMtkOrcF2aD7boLZROacXEScDZwLfrzH4VyLihxHxzYj4pQHmsSwiuiKia8eOHaNRrXG1fVfPsMo1NmyH1mA7NJ9tUFvDoRcRxwBfBX4vM3/Sb/CDwGsy8wzgz4G19eaTmbdkZmdmds6cObPRao27k2a0D6tcY6P3/b500nrWT72aJ1/231k/9WqWHvODJtesLO4Pzdf/ve7dJ340bQnceBpsWt2kmjVXQ6EXEW1UAm9VZq7pPzwzf5KZe6rPvwG0RcQJjSyzVV2zaC7tbZMPKWtvm8w1i+Y2qUZlumbRXC6f+g98ou1WZk96jkkBsyc9xx/mzcXu5M3g/tB8fdvg0knrX9onSNj9FNxzdZH7xIhDLyIC+BzwWGb+WZ1xXl0dj4g4q7q850e6zFa2+MxZfPxtpzNrRjsBzJrRzsffdnrRJ4ybYfGZs1hx9Fc5Kn56SPmUA3vhviP2wuGW4/7QfH3b4PenrD5sn2BfT5H7RCNXb54DvBvYHBEbq2V/APwcQGbeDFwOvD8i9gM9wDsyMxtYZktbPPm7LH7ZCpjWDS+bDZOXA1c0u1rFOarn32oP2N09vhUpzabVlQ/R3d0wfTaLFy5n8XVu/820+MxZlX80bqhzrFHgPjHi0MvM9UAMMs6ngU+PdBkTyqbVle6CfdWTxL3dBwAd7vjjavrsyvtfq1xjw+2/tblPHOQdWUbLfSte2uF7Fdp90HQLl0Nbvwsm2tor5Robbv+tzX3iIENvtNTrJiiw+6DpOq6AS26C6XOAqDxecpNHHGPJ7b+1uU8c1Mg5PfVl90Fr6biiyB26adz+W5/7BOCR3uix+0Alc/vXBGHojRa7D1Qyt39NENGK3yDo7OzMrq6uZldDktQiImJDZnY2Oh+P9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScUw9CRJxTD0JEnFMPQkScVoKPQi4sKI2BIRT0TEdTWGvywi7qwO/35EnNzI8iRJasSIQy8iJgN/AbwVeD3wzoh4fb/Rfhv498z8L8CNwJ+MdHmSJDWqkSO9s4AnMvPJzPwp8FfAZf3GuQy4vfr8K8DCiIgGlilJ0ohNaWDaWcBTfV53A2+sN05m7o+I3cDxwHP9ZxYRy4Bl1Zf/GREPN1C3ZjuBGus4wbgOrcF1aA0TfR0mev0B5o7GTBoJvVpHbDmCcSqFmbcAtwBERFdmdjZQt6aa6PUH16FVuA6tYaKvw0SvP1TWYTTm00j3Zjcwp8/r2cD2euNExBRgOrCzgWVKkjRijYTeA8AvRMQpETEVeAfwtX7jfA1YWn1+OfC3mVnzSE+SpLE24u7N6jm6DwLrgMnA5zPzkYhYAXQE8fmPAAAE60lEQVRl5teAzwFfiognqBzhvWOIs79lpPVqERO9/uA6tArXoTVM9HWY6PWHUVqH8MBLklQK78giSSqGoSdJKkbTQm+i38IsIuZExLcj4rGIeCQifrfGOAsiYndEbKz+LW9GXQcSEVsjYnO1foddEhwVN1XbYVNEzG9GPeuJiLl93t+NEfGTiPi9fuO0XDtExOcj4tm+30eNiFdExL0R8Xj18bg60y6tjvN4RCytNc54qLMOKyPin6rbyl0RMaPOtANud+OlzjrcEBHb+mwvF9WZdsDPsPFQp/539qn71ojYWGfaVmmDmp+lY7Y/ZOa4/1G58OVHwGuBqcAPgdf3G+d3gJurz98B3NmMug6wDicC86vPjwX+ucY6LAC+3uy6DrIeW4ETBhh+EfBNKt+5PBv4frPrPMh29W/Aa1q9HYDzgPnAw33K/g9wXfX5dcCf1JjuFcCT1cfjqs+Pa6F1uACYUn3+J7XWYSjbXZPX4Qbgw0PY1gb8DGtW/fsN/1NgeYu3Qc3P0rHaH5p1pDfhb2GWmU9n5oPV5y8Aj1G5A82R5jLgi1nxPWBGRJzY7ErVsRD4UWb+uNkVGUxm3s/h31ntu83fDiyuMeki4N7M3JmZ/w7cC1w4ZhUdQK11yMxvZeb+6svvUfn+bsuq0w5DMZTPsDE3UP2rn5dXAHeMa6WGaYDP0jHZH5oVerVuYdY/MA65hRnQewuzllPtej0T+H6Nwb8SET+MiG9GxC+Na8WGJoFvRcSGqNwKrr+htFWreAf1d/BWbweAV2Xm01D5IABeWWOcidQeV1LpJahlsO2u2T5Y7aL9fJ1utYnQDv8VeCYzH68zvOXaoN9n6ZjsD80KvVG9hVkzRcQxwFeB38vMn/Qb/CCVrrYzgD8H1o53/YbgnMycT+XXMj4QEef1Gz5R2mEqcCnw1zUGT4R2GKqJ0h4fAfYDq+qMMth210yfAX4emAc8TaWLsL+J0A7vZOCjvJZqg0E+S+tOVqNswHZoVugdEbcwi4g2Ko20KjPX9B+emT/JzD3V598A2iLihHGu5oAyc3v18VngLirdNn0Npa1awVuBBzPzmf4DJkI7VD3T23VcfXy2xjgt3x7ViwkuBpZk9cRLf0PY7pomM5/JzAOZ+TPgs9SuW0u3Q/Uz823AnfXGaaU2qPNZOib7Q7NCb8LfwqzaX/454LHM/LM647y69zxkRJxF5f1+fvxqObCIODoiju19TuUihP6/bvE14Lei4mxgd2+XQ4up+19tq7dDH323+aXA3TXGWQdcEBHHVbvdLqiWtYSIuBC4Frg0M1+sM85Qtrum6XfO+tepXbehfIY101uAf8rM7loDW6kNBvgsHZv9oYlX7FxE5SqdHwEfqZatoLKzAEyj0lX1BPAD4LXNqmud+p9L5TB6E7Cx+ncRcBVwVXWcDwKPULmy63vAm5pd737r8Npq3X5YrWdvO/Rdh6DyY8E/AjYDnc2ud431OIpKiE3vU9bS7UAloJ8G9lH5b/W3qZyzvg94vPr4iuq4ncCtfaa9srpfPAG8t8XW4Qkq51h694neK7BPAr4x0HbXQuvwpeq2vonKB++J/deh+vqwz7BWqH+1/Lbe7b/PuK3aBvU+S8dkf/A2ZJKkYnhHFklSMQw9SVIxDD1JUjEMPUlSMQw9SVIxDD1JUjEMPUlSMf4/P215qQxFkzoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(figsize=(7,7))\n",
    "\n",
    "N = 16\n",
    "L = 20\n",
    "\n",
    "x_grid, y_grid = lattice(N, L)\n",
    "plt.plot(x_grid, y_grid, 'o', label = 'grid')\n",
    "\n",
    "x, y = random_lattice(N, L)\n",
    "plt.plot(x, y, 'o', label = 'initial position')\n",
    "\n",
    "plt.xlim(0, L)\n",
    "plt.ylim(0, L)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to initialize the initial velocities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "v0 = 1.0 # maximum velocity\n",
    "\n",
    "vx = 2*(np.random.rand(N) - 0.5)\n",
    "vy = 2*(np.random.rand(N) - 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, for the Verlet method, we the need the (fictious) previous position.  We can calculate this as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.02\n",
    "\n",
    "x_prev = x - vx * dt\n",
    "y_prev = y - vy * dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Verlet method can calculate the subsequent position and velocities of all the particles given all the accelerations.\n",
    "\n",
    "These accelerations are given by the forces due to all of the forces acting on the particles.\n",
    "\n",
    "\\begin{align}\n",
    "a_{j, x} &= \\frac{1}{m} \\sum_{k \\neq j}  f_{k,j} \\cos \\theta_{k, j}\\\\\n",
    "a_{j, y} &= \\frac{1}{m} \\sum_{k \\neq j}  f_{k,j} \\sin \\theta_{k, j}\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $f_{k, j}$ and $\\theta_{k,j}$ is the magnitude and direction, respectively, of the pairwise interaction force between two particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pairwise force is given by\n",
    "\n",
    "$$ f_{k,j} = - \\frac{ \\partial U} {\\partial r_{k,j} } = 24 \\left( \\frac{2}{r_{k,j}^{13}} - \\frac{1}{r_{k,j}^7} \\right) $$\n",
    "\n",
    "where $r_{k, j}$ is the distance between particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "#@numba.jit\n",
    "def F(x1, y1, x2, y2, L):\n",
    "    # calculate pairwise force between (x1, y1) and (x2, y2)\n",
    "    \n",
    "    time.sleep(100.0)\n",
    "    \n",
    "    dx = abs(x2-x1)\n",
    "    dx = min(L-dx, dx)\n",
    "    dy = abs(y2-y1)\n",
    "    dy = min(L-dy, dy)\n",
    "    \n",
    "    r = np.sqrt(dx**2 + dy**2)\n",
    "    \n",
    "    if r > 3:\n",
    "        return 0, 0\n",
    "    \n",
    "    f = 24 * (2 / r**13 - 1/r**7)\n",
    "    \n",
    "    # fx = f cos(θ)\n",
    "    fx = f * (x2-x1) / r\n",
    "    # fy = f sin(θ)\n",
    "    fy = f * (y2-y1) / r\n",
    "    \n",
    "    return fx, fy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANcAAAAVCAYAAADVXXsKAAAABHNCSVQICAgIfAhkiAAABLFJREFUeJzt2n/s1VMcx/FHCbVQZLT5FakYhjG0Jl+hmNkM/Wd+TH+xYtj8nvgjzQzJJs1osZmNxWz5tRaSaWaamYh0GX2tYlKt/MjXH+/z1e1+7+d77/187/3eb+3z3D47d+e8zznvndc9vw8FBQX9ymJsxPB2O1JQMMA5C124qR7js/Evbm+BI9dgPlbgj+TUS3XkOxrPYwP+RAlP4tAW1z8KM7AE32EHtuAj0ZiDq+QppXKrfb806O/eTl69m1FWHu3Ip98SdOKg8sghVQznJOef6cXxvNyP07ENP+GkOvKMxcc4Am/ga5yDW3EpJuHXFtU/XbRDJ5bjRxyJq/AcLks2XRX5tojOX8m2Ov3cV8ijd7PKyqsdjev3CFZhlug/VRkvZq2FNRzPy4UYh0HoUN9I9k6ym1kR/3iKX9DC+qfgCj1HudFCrC5cXZFWSl9BPr2bVVYe7civ3xr8gP2yDOamSi/KUXijdKjdQCckm/V6NtLBYiTZLt/esJ76e+PelH9+RXxJ0bmq0aFv7d3MsrK0I79+D6Yyp3VHVC4LL8YufJKj8FYwJYXvihm1nK1Yiak4D8v60S/4O4X/VEk7ENfiWNH5v8CHom0L2k9v2pFPv5UpvESstvboXMNxhpjetudyuflMSOHajPRvRecar3871xBcl36/XSV9NF6siFuPG/FBC/0qqE0t7cin36cpnNwdUb7UOkqsFzsb8bTFjEjhloz07viR/eBLOXNxKpZKo1QZL4hl9WgxYJ2GZzEGb4lNeUH76E078uu3BTvFbNeDiWLN+EpG5pLsI8pqX631cEcddguTzYyM9Dkp/e4adeWtvxqzUr41OKyBfI+lfEsarK8/KWmuxuV05MjT7LLyakd9+v2sbKlZvizckcKhGRnXiZ5ZLxsasM2ie2YakZF+SIVdq7kF8/CVGN1+ayDvAtyhbNkwAGmHxv1FX7SjPv2G2d2P9uhcG1M4KiNjf5wgVvJNCsdnpI9LYdaerJnchifwpWiLjb2b96DbfiC/emmHxv1BX7Wjtn6DxfZkfXlEN53YZPchwkBgeQqnqn4UP0mMFK0+3bxLiLNa3LfkEWdiCr9vllMFddEM7ait3wRxB7e6O6L8D9sljhsPx4k5HWg268Qx/BgxrZfzkBhFFut5ujlW3ODv3wQfHhCb4M/EqLe5F9tTVF/LH4en0+9q+4RFov1vyOvkPkS7tCO/fsR1ELsnhB73XK+JW+tp4j1Ws7kyfcRpDDEiLEq/N+POijw3i+dPT4kGWoNzxSi0FvdVqWeZaJDj7Xkh2Gj91+Nhcb+xQmyIKymV5Z8uDleWi+XBVvFnuVzsZZeKjXEl3YNc1r3L3koevdulHfn1I1ZXu8QTvaocIB4nrsoy6COz9X76VMrId4w4Iu3EX+KZyTzZJz6lVN6YPtZfy74L75fZX4CXxfvH38Vl5Sa8J+5WBmX4+7l4z9noQ+SBzmyN613SHu3Ir98IsT15PSP9f+5JFZ9Zy7CgKYwUI96j7XakIDczRZ85v5bhUDEzvNlqjwoQj0t32r3UKdi7GCauJF6tN8Nk8RBxIB8bFxQMBE4WS9Ax7XWjoKCgoKCgr/wHZ7mxaXclV7IAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left ( -10.125, \\quad -1.125\\right )$$"
      ],
      "text/plain": [
       "(-10.125, -1.125)"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F(0, 0, 9, 1, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 s ± 1.82 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "F(0, 0, 1, 1, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The net force is calculated by summing up all of the components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "#@numba.jit\n",
    "def net_force(x, y, L):\n",
    "    n = len(x)\n",
    "    fx_net = np.zeros(n)\n",
    "    fy_net = np.zeros(n)\n",
    "    for k in range(n):\n",
    "        for j in range(n):\n",
    "            if k != j:\n",
    "                fx, fy = F(x[k], y[k], x[j], y[j], L)\n",
    "                fx_net[j] += fx\n",
    "                fy_net[j] += fy\n",
    "                \n",
    "    return fx_net, fy_net"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now apply the Verlet method to calculate the evolution of the system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "#@numba.jit\n",
    "def Verlet(n=20, L=10, v0=1, tmax=100, dt=0.02):\n",
    "    \n",
    "    # initial position\n",
    "    x0, y0 = random_lattice(n, L)\n",
    "\n",
    "    # initial velocity\n",
    "    vx0 = 2*(np.random.rand(n) - 0.5)\n",
    "    vy0 = 2*(np.random.rand(n) - 0.5)\n",
    "    \n",
    "    # previous position\n",
    "    x_prev0 = x0 - vx0 * dt\n",
    "    y_prev0 = y0 - vy0 * dt\n",
    "    \n",
    "    N = round(tmax/dt)\n",
    "    \n",
    "    # allocate memory for n particles over N timesteps\n",
    "    x = np.zeros((N, n))\n",
    "    y = np.zeros((N, n))\n",
    "    vx = np.zeros((N, n))\n",
    "    vy = np.zeros((N, n))\n",
    "    \n",
    "    x[0, :] = x_prev0  # previous position\n",
    "    y[0, :] = y_prev0\n",
    "    x[1, :] = x0  # current position\n",
    "    y[1, :] = y0\n",
    "    vx[1, :] = vx0 \n",
    "    vy[1, :] = vy0 \n",
    "    \n",
    "    for i in range(1, N-1):\n",
    "\n",
    "        fx, fy = net_force(x[i,:], y[i,:], L)\n",
    "        \n",
    "        x[i+1,:] = 2 * x[i,:] - x[i-1,:] + fx * dt**2\n",
    "        y[i+1,:] = 2 * y[i,:] - y[i-1,:] + fy * dt**2\n",
    "        \n",
    "        vx[i+1,:] = (x[i+1,:] - x[i,:])/(dt)\n",
    "        vy[i+1,:] = (y[i+1,:] - y[i,:])/(dt)\n",
    "        \n",
    "        # periodic boundary conditions\n",
    "        for j in range(n):\n",
    "            if x[i+1, j] < 0:\n",
    "                x[i+1, j] += L\n",
    "            if x[i+1, j] > L:\n",
    "                x[i+1, j] -= L\n",
    "            if y[i+1, j] < 0:\n",
    "                y[i+1, j] += L\n",
    "            if y[i+1, j] > L:\n",
    "                y[i+1, j] -= L\n",
    "        \n",
    "   # vx[N-1, :] = (x[N-1,:] - x[N-2,:])/(dt)\n",
    "   # vy[N-1, :] = (y[N-1,:] - y[N-2,:])/(dt)\n",
    "    \n",
    "    return x, y, vx, vy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 214 ms, sys: 5.92 ms, total: 220 ms\n",
      "Wall time: 221 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "x, y, vx, vy = Verlet(n=6, tmax=20, v0=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 58.5 s, sys: 315 ms, total: 58.8 s\n",
      "Wall time: 59.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "x, y, vx, vy = Verlet(n=50)  # old, not jit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 427 ms, sys: 7.26 ms, total: 434 ms\n",
      "Wall time: 440 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "x, y, vx, vy = Verlet(n=25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " "
     ]
    }
   ],
   "source": [
    "%%prun\n",
    "x, y, vx, vy = Verlet(n=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discuss different programming languages\n",
    "\n",
    "* Interpreted\n",
    "\n",
    "* Compiled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXl8VPW9///8nDMzCQkQsrAFCCHsRoQIAi4I4m5VXGpd6L3a1mLvUu/jtrduWNvaUm/b29tvb3+9bdWqtRdtRQWsbV0QUVxYEkAxQlgCJCRhT0LINjPnfH5/nDmTmcmZyewzCfPqw5KZOXPmZGbyOu/z+rzer7eQUpJBBhlkkMHAgpLqA8gggwwyyCD+yJB7BhlkkMEARIbcM8gggwwGIDLknkEGGWQwAJEh9wwyyCCDAYgMuWeQQQYZDED0Se5CiGeEEMeEEJ/53FcghHhbCLHX829+Yg8zgwwyyCCDSBBO5f4ccE3AfQ8B70gpJwPveG5nkEEGGWSQJhDhNDEJIUqB16WU53pu1wCLpJRNQojRwAYp5dREHmgGGWSQQQbhwxbl80ZKKZsAPAQ/ItiGQohlwDKA3Nzc2dOmTYvyJRMMZzs4z4BjMDhy/W8DnNwHUoIQUDjJb5vjTjtHOlVy6CKXLtRBQxheUNB7n3GEq6uL5qYGpJQIIcgfPQZ7dnZcXwMI/3dwtlu/RxlEBVdXF51nTtPZdhrM+ktAweixifmcM0hrVFVVnZBSDo/kOdGSe9iQUj4JPAkwZ84cWVlZmeiX7Bv1W+DgRihdAOPmGrf/cCNoTlBdcM2D8MZDPbdn3QFVR0FqIFRYfK/xXM9znNh4rPvLfM/+R+y40ZUzOL7wmP8+7v6D8dqBrxvqdghsXv0SH/z5j+Ah94tv/wfm3fyl+L9Pvu/L3X8Iflwbfw7rVxjvEQLmXA3X/7/4Hs9ZgMY9u6h+bz3VG95Gc7v9HxSCSxLxOWeQ9hBCHIr0OdGS+1EhxGgfWeZYlPuJHiYRDiqEzpNhEaL3eV7CcsDdrxn70ZwGMWlO2LXW/zbC2NZ8TukCv+eoSK5Vt2DHjU3oaNLVex+fvAg7XuzZxzX/6UP+Frfvfs043kCy9/zeg1w5RpUMSCkZNGRI/N/jwPfl4Mbg73HpAlBU0DRAwvYXYOZd4X0mZzlMQm9vaebgjsrepO6BqqqMK5+R5KPLoL8iWnJ/Dbgb+E/Pv2vjdkSh4EvobzwE7m5AB6GAmmUQpBXR+1bEVoRVusCfvKcvgUMf99yeeafxXyDRep6jYePv2lzmKjUg3eiKHTVwH8jQJ5C+TgYm2XtOTEebJgE9atiBj9Zx3hXXRlT99wnv+9JtSC2DCoNvO24uVHwZKp81flfdHfpkkIGX1D979y10Teu9gRCoqo3SWbPJHZZP+cLFFE+ZnvwDzaBfok9yF0K8CCwCioQQh4HvYZD6S0KIrwF1wG2JPEjAv+IWAqQO6MZjUjeI/m/fNqpZkwx7SS6eCjmwCh83t6eCN2+PPKc3SfoSlc9zfrB9GH9qHMUe5zjmK7toHTmPFXPu8d8HwI4/BT+B9HUyOLjR2Id5H75kINlXvYvGtT+neOfPQlf/kZD/uLnG+/W3b4OuGyfUkecEf97MO3tOSIoKrYeN18sQPCtXrmT58uXU1dUx95yp/ON1V+NsOGhdpXtIvXzRFRlCzyBq9EnuUso7gzx0eZyPJTjqt8CGJ4wKUuogFVAUz0KTp3IXiucxD9FveAIWPdy7Uu882ZvIwfg3kLzDIb9xc9m8ZQPQzjY5hW3aFCa5cq330dcJJNTJwLzPc2Iqz2/h0+bRGG+CAGDrOxtYUhiBFBQO+Xee9Mg/et/SjHnC++QFQ5ap+oPx+uaJ9izFypUrWbZsGcMHObi54hwuKB1H+8G9KEL02lZRVc697KoMqWcQMxK+oBozzMo7mARjau6BUs3+d+HgBzD/n41KH8W/Uo8j2RTkOuB4u/9tK/R1AunrZOBzX3HpAiY++2f2V+/2br7vqKQxN5/i7ObwpKBwyH9QoedKSfRIM6Gq/3Fzjcd0t/E6vifas5Tgly9fzvBBDu5bOB+bqiAAIYThdFKUjPSSQUKQ3uTuW7GjAwqULQpOFCPPMcipoQovsX30K6PyVBSDvBJAMMNyHCFvRw2rk5DPfXPvGML+x77juYLxVO+DbmLJZcXhSUHhkr/uWSTVNfj7d4zX0t3Bq3/fEwI61G4wXvcsreDr6uq4bGoZqiJQPKSu6zqalMy+8roMoWeQEKQvuQer2AOJPbCKHD3TQ+4eSB2DmDDIK5Rm3M9QPGU6E2fPZ3/lJu99+z7fQ+Md91I8zkMWoaQgCI/8vUZr87bo+TlY9S91D8Hjkcq6DLlmgLz3kaCkpIT9x0+i6Ya8pUvJ1gP1NHRrPPTyv6T68DIYoEhfcje18lAVu5WtceadsH2lZ1HPZhCM5savggzmqokSLR3OkLcTiblLbmV/1WavLRIpqX5vfU8lGKn0Y0n+nrUOL6TnZBuq+vddCzm77ZErVqxg2bJl/O69TUwcXsj+4yc53unkySefTPWhZTCAkZ7kXr8FWusNctYxSMRKirGyNS74Ntzzuj9BbXjCIPZQrpoYcKrdGfJ2ImFU7/P8qvf2lubwdxAu+TdUwe6/YZTiPidbCF39T7ocdv+Vs9keuXTpUsDQ3t+tqaWkpIQnf7nCe38GGSQC6UfuvtW4osLsu41qPKCJh9IFvf3pJpkHEtSih3sIx7RRSo/745MXYvaFh72gmiDMXXIrB7Zv9XqlD2zfSuOeXfHRcc33sn4L7F3nea/t/ifbvqr/fe/0/ozOMixdujRD5hkkFelH7r7VuA7kjfUn9kAZxspREghfH7vpqjFPHttf8F8cjILgE7agGiaKp0xnQsUF3upd1zS2vvYKS/7j0Ti/kgz414NoXD8ZDFhUHWpmU+1J5pcVMnt8fkT3ZRA/pB+5B6vGIbgMEw5h+BKOWVm2Hja82APAspc7zP+PY1/l5vhV7+CxN/q4ZiKRV+JsPc0gfoiEdMPZtupQM0uf3oTTreOwKay8dz5AWPclgvQ7tm+nY8tWcuZeQE5Fhd9twO+xgYb0IndTcvH1sJudmePmhib+SOArNex4sceRE+WCayoXVE2UL1zMp++84bewGtfqPZIoggwShliq4miJOJJtN9WexOnW0SW43Dqbak8ChHVfsNcx7w+H8APJu+4rX0U6nQiHg5EPP8TRJ/7TuG2zGdefbjfC4aDk2WeM5w8gsk8fcreKCQhsqLGKCYgF5v5iXHBN5YKqCauF1bhW75FGEWRgiVjJOdqq2Oq+cIl49vj8sLedX1aIw6bgcuvYbQrzy4wiINz7gr1OOITfsX27H5nnLVmCdDpB15EuF21vve13GwApkS4XrWvW0rp2rd+JQGtp7ddEnz7kfnBjTwXt7u7dUGPKAPG+xB83N/iCa5hSTaoXVE1Y2SLjWr37RhG4uwyPe4bcgeRIFrFUxVb3RULE4W47e3y+98Th+16Ee5/V64Qi/LKjtTScqkX52hLG1e3uRd7C4UC6XAi7nSFXXUlHVZVxW1WNyl3TEHY7QM9znU6O/PBHoOv9uqpPD3Kv3+JpPDK91DqMOs/fUpdIl4XVgmsEUk2qF1RNJLx6L11g2FO9sb7/5+9kGoCIJ2nHKlnEWhXHQsSRbhsonURyn9U+rQi/7GgtKz74HXbdDQ+sQ330EYTDge50oisqR+YvZupNS/xIOWvKFEvNHTAqd5fL0xujBa3qS559pl8QfMrI3UzJK9YP884/5pJl8xnoKhTIHppcl0XggmuCvfGJQp9NTbFg3FyouAsqnyOqhdUUIBkySLIki1ir4liIONJtY0HgPoP93g2narHrblQkUnOjtbTS+cQvWfX7tewoKKP2o3ZW3juD2ff1EHFORYUfMfv+XPLsM3Rs2Yo6LM/Q5j0VP+B3RdC6Zm2/qOJTQu5mSt55Bd08ttCBXZEoCE++oTBiBhIQ8BU2Qkk1Ft74dFhQNRFzU1NfmHlXT9dqmi2spkoGSZZkAbFXxf3Vcmgee8f27Zx4YxXT517AbV9bAg+sQ2puFIeDnLkX8FHLEP40eTG6BNXnqicc+BJ/YIVvVvVCVWlZvdq7EJvO2nxKyH358uWcV9DNO/+Yg0MFVYBbk2hSkDX/K+lxqR+BNz4dFlR9kfCmpjRcWLUi8mTJIMmSLM52BC6YTn32GXj+Ob8qev6hZsvPLFJ7ZWCFb1b1rsZGWlat6hfafErIva6ujjsvtuNQwaYINF2y7oCbx9938tGPfpGKQ7JGmN74efbr2cco79NStaBqwqqpKW7SDESW8Z4kWBF5MmWQDEEnBr7Wxo4tW/3kkY4tWym6b5kfiVp9ZsHslZHAJPuO7dv7jTafEnIvKSnhREcjqjDmfyoCXtnlplEZm4rDCQ8hvPHf40M05csUiDNs0qczLCf9WuxPNtTFb2fehVXd+DcNIgVMIne6dIQQ5Oc4kiqDZBB/BFbqIx9+yM/9YkomgQj8fKxO/NF+fjkVFf1Gm08Jua9YsYI9v78PTRqVu1uTFOc5WPGtFak4nMhg4Y1XcfJD+3MIJC5s/LLrv4A5KT3MwI7Vht2fx7djNVgUQYowe3w+j11fzmNrP0OXksdfr2bqqCEZGaQfI7BS11pavcQaCWEGu4KD6Lph+4s2nxJyX7p0KZtaPkAefQG3LnFJwUV3PcjV/SVYKWDBVUpQ0FGFBOlmfttbsPFESrNUrDpW4ybNxBJFkEA0dzjRpYxLhZZBauAXFzD3gl6VeqAWHg6CXcHFU64xEa42nwyCT40Vsn4L85tXG9ORFAXbdT/n6jn3pORQoobPguvTW1q4+/RvsUs3GgoXtb0J6/9uLL5WfDklC8RWrpm4IU2jCKykmQz6DwJlmJJnn4mqUreC1ZVZPOUaE+Fq88mQa1JD7r6DOKQwFuj6Izw6/MtbNvCWs4D5yi6KxQnutL3rsU1qUPlsyoZEl1XM9iP3kRPK4rPjNHXMhJJmMkh/hLNgGk8Ek2viEVwWVJsPkGsS6bBJDbnHKwAsTVCQ62DL8Sls06ZwvtjDbbYPUDHH00lLb3wycPRAbcjbMSFNowgy0kz/Ql8yTCKRKGeNCStt3k+uSbDDJvnkHpj8OAAyvn3jBrbJKfyi+Gc8VLy9xw8fp9z4SBHYvBTXZqY0jSIItXiWQXohkTJMuEiks8YXVnJNoh02ySV3Z3vvYRv9nNitUJtdDtffbXRzpjA3PtAxE3g7JqRpFIEpzfz9syauPXd0pmpPYyRbhgkH8XbWBMJXronEYRMNkkzuZ6yTHgcY6k91GD/EOTc+UgRq7HHT3E3MvMu4IjGHkaeBvFZ1qJnHX6/G6dbZevBURnNPM6RShgkHiXTWmIjGYZOrKLmRvk5yyd0xmI5TOh1NgpzRkpw0IIN4YPiQLL/bu4608cLmOu6aV2LcEafc+EhxYEeV3+24au5epJffPVGX1RnEjnSQYcJBspw1JsJx2AxWlCGR7lfpe5P4QXcL6jYUcXznUOo2FNH80QFO/O5JOrZvT+ZhxB23nN+7s/aZDwKI1PTGq1kgVMMGKvWeq5hPXoCNPzeq/Digcc8u9lfFZ19BYeV3TzHMy2oFMnbINIOVDJNTUZFwKWblypWUlpaiKAqlpaWsXLky4n2Y3ytV0MtZ8+t391F1KPb1LFOyGX7//Yz67qOIrCxQVYRNJTtbRvxFTmrlrre3I3MHe85G7rQO3YkEs8fnM7c0ny0Hez7g/cfb/at3CC+MLE7e+PrqnUjpU00LQfnCxVHvzxJp6HfP2CHTC6mWYcwEWjlaUnhdIcd2H2PZsmWA0UwZLhLtrDHRy2Hz9qvkNDxD0XE5PNJ9JZXcldxc74eb7qE7keLBa6dz228/QjcbQoFH1+wE6E3wocLI4uSNH1c+A1W1obmNiTSKonCi7iD11TsZVz5jQCdEZuyQ6YF0kGGWL1+OHC2Z8MAEhE0g3ZIDPz3A8uXLIyJ3SJ6zxkRORQU5Z9bD+s6onp9ccs/JoeSpp/pF6E6kmD0+nx/dNINH1+z0ErwugxC8iV4Lrl3EyxtfPGU65Yuu8EYQSF3nnWd+i5QS1Wbjsru/TmdbW+xEn4YJkZlO1dQh0hTHRKOuro7C6woRNoFQBQC503Kp+1vsQXqJdtZQvwVa60GxIaNY1Eq6z72/hO5EA5PAAwn+sbWfhZYGTLnmkxfi6o0vX7iYz99/B83tRgiBrusgJZrLFT+iz0gzGXgQbYpjIlFSUsKx3ceQugQBUpe0726npMSi2IoQCXXW1G/psY0rKs2d8kSkx5fSGarpHLoTLe6aV0LdyXZ++37Pgqqmy74v2cwqPo7e+OIp07ntuyuor97JoCFDePcPT4UkekVRKV90BeULF4dP8hlpJgMP4pXiGE+sWLGCbz7xTT9TV3Z2Nisejk8CbcKcNWZEi9RAB6dGxBOA0mNAtgfhWII6tmwF0nPh1bwUO93tRtDzfVIVEX6nZJy98cVTpnuJuqikNDTR6zqfvvMG1RvWRUbyaRhFML+sEJsicGkysvc/g4iQiBTHeGLp0qV8rH/MBtcGhCJAhbsevCtivT0SxKVLOiCi5XS3bIt0F8LPUZFgzJkzR1ZWVoa1rfmlCdTmRz78kHHb59Iv1ZJN1aFmXt12mFWV9bh1iU1VQErcukQRgseXnGutufeF+i3+3ngUj4Uydm98455dfkTvdrl64oEBhMBmt4cn2dRvgeeuN6QZMI7tnr+mlOCrDjVz51ObvH9gL349didDBv6wWjCF9Cu8dhzbwdfe/Bou3YVdsfP7q3/PrBGzEvqa4Q5ht4QZ0TKo0FvIiZJ5VVLKiIZEpFXl7otg2rzfpV+KJRtfUndp0lupa5rOHXNLKB42KLYFlQgHdUeCwIq++r31VG94G81zlRSRZJOGUQSbak/i1nQkxueRkWXij3RYME1XBMo1Yevwvlp7jAVcTOQuhPh34F4MBWIn8BUpZVcs+7RC4KVdKDtlMqqGYKQOIDCaHG45f2x8yCQJ3niT6MsXLvaSvK7rkUk2aRZFkHHMxB9+EkxFRVrGB1ih8mglbt2NRKJJjcqjlQmv3AMRtg7vq7XH6D6LmtyFEGOA+4FzpJSdQoiXgDuA56LdZzhIdU6yeQbudumWpP7F2WO5NV7EbiJJ3nhfkreUbKREc7tCkHz6RBFkHDPxhZUEExiCla4V+5yRc7ApNly6C1WozBmZ/BGYYdsm4xiHHqssYwMGCSFcQA7QGOP+wkIqc5Jf2XbYj9gTSupWSKA33kRYkk0gyecdpTjNRu9lHDPxQ7DogFQvmPYXhGObnGvbx8/ntjEmToGCUZO7lLJBCPFfQB3QCbwlpXwrcDshxDJgGRAXb2kgkpWTbEoxL1XWe4ndpgq+NGdcckg9EH154+MUY2Al2fQi+XV/5zNF4fJRxZw3rCFt/O4ZaSY2pDo6IF5IB1kGQtsmZ7GHZ5Ufk1XlBltWXIIEo3bLCCHygVeA24EWYBXwspTy/4I9JxK3TDTw/TICxmWkR7KREJVkE0xfFxie9hU3z0jY7xM2zNV1X6kGAAG27LgmTjbu2dWb5M1XE4IZwxopH3qc4qHOtMjrf2FznVeaiVf+x9mA/uKECQepcMuEC7Ny/5pczb+rq7AJHYRKw/nfYs3gO7xVvhAiqW6ZK4ADUsrjAEKIV4GLgKDknmiE1RQVgWRjEoNb9z8BCiDLrlimQaYESZBqTARW8jvXv4nUdQCklHzaPIrqlpGU5x2lfN0zFH8lI830R2ScMMnB7PH5rLnRjmubG+WoDaSGrtj59pYhbHHXeAuSaBALudcB84UQORiyzOVA4sryKBCLZPPC5jq/GAFIgb4eKRIUY2AFk+RHTijjnWd+i671XC1oEj5tGUX1m3WUt3yP8i/cEZ+gsiiQaWYKHwNFhglEusgylqjfwrQ3v9zjgDv/bl51L2DL5iy/giQaxKK5bxZCvAxsA9zAduDJcJ4bl1CdCBDJaKvOJ37JhprjHH3/Q6YUTWR3QSlgdJnefkGK9PVIECrGIAGDus+74lr/RVe32/OIh+Q3VVFd+WnksQbxhPD0CwuR/NfuJ0iHBMdEIR3cMkEREDNA3lgmlFyGo2pTzHOAY3LLSCm/B3wvkuckIgM5HIQj2ehOJ+v/5zkWH6rErru5XbHx8CX3saewlB9G22WaKgRKNQnKjYcAqeavf6J60xY0nzkwoe2TiUWmmSk8ZGSYFMHC+jh7nLWzJlIkvUM1Xcag+Uo2zWvWIJ0uXEJFlxK77kZFInU3V9ZX8Z0R7VzoKAb6Ebmb8G2CSoA33hfFU6ZTPOgLlB97nupTBVS3jEBDNR7s0yOfGGQcM8ExUGWYQPjKMm7dnT6yjGmEsLA+WjlrIkXSyT0uoTpxxBpnISvnL6P8+D4+LZoIwBV1VUjdjVBtXHu4EnFoM3Wvv9h/HQNJXHDl4EaKs1spHn2K8mHHqc69huqaY6E98gkk+UwzkzUGsgwTiDxHHjrGor+OTp4jL8VHRFxjBoIh6eQeysyfTB3etDj+aUsdWv54qvPHA8ai6aMLvsG3hp/hHLUjKY1RSUMyvPE++e7Fue0UX3cR5XfMC68RKkEkn3HMGEi3QRrJQquzFYFAIlFQaHW2pvqQ4hozEAwpCQ6LOlQnTghmcexZNL2I2ePz+3TZpHP8cFCEXHCNg1Rjke9efPdrFH/9X/puhHr3LS7/6jc474pr4/orZ6SZ9BykkSzkOfJQhYoudRyqI/ULqj4TltCJOWYgGNIiFTKZOnwwi6Oq9I7mDemysdtRh+X1+oNJdfxw2EikVBNk9F7IbldA1zTW/f43HD1QG9cqPiPNpOcgjWRgx7Ed/HTrT9GkhipUHrjggdTq7QETlph9d1xMDVZIC3JP+CxC/GUYX2Lvy+IYzGWTbvHDUSMR3vg+Ru+FbISKdmBIHzgbpZl0H6SRDFQercSpOZGe/6VckrGwPiaqizstyD2hswgJLsMogogtjuHGD/crySbe3vgwR+8FbYRKgB5/tjUznU0LpqGQdh73OKY+9oW0IHdI3CzCSGSYSBE0fthCsukXTpt4euMjGL3XqxEqUYuuZ1Ez09m0YNrvMOsOQCRMjjGRNuRuhVjkmmhlmEgR1sSo/ua0iYc3vnSBsWCkeWKAt/9fyC9zWOmTMZD82dDMdLb41iNBWkUPBNofZ96Z0JdLa3KPVq6JpwwTCYJJNv3SaRPrgmuUo/cSRfID3TGTkWGskVayTBLsj75Ia3KHyOQaP+96AmSYSBCp0yZtZZtYFlxjGL0Xdo58mPbJge6Yycgw/QBJ1NuhH5C7Fazkmhc21/HdNTv9SB1SG/gVttMm3WWbUAuu7m7Y8IQxyNuyColt9F487ZMDzTGTkWH6RtrIMiGiBhKFfknugXJNzZE2lq/e2Ys+bEmu1vtCv5dtekk13YAOtRvg0Me9v7gHNxpyTBxG78XDPjmQpJmMDBMe0kKWSULUgBX6JbmbaGzp5Hfv7eedXUd7EfvMsXk8dkN52lZm/bpBypRqNjxhELvUDaL/27eNitr8Avfhd48GsdgnB5I0k5Fh+hGSrLWbSCm5r1y5kuXLl1NXV0dJSQkrVqxg6dKlIZ9j6up/rqzHHajBeGBXhR+xJzu3Jlz06wapcXMNKebQx8YXVgiD5KXes+CaNw7m3Qcf/38h/e7RIFr7ZH+WZjIyTORIuSyTpKgBK6SM3FeuXMmyZcsYOaSUK2fewd6mT1i2bBlAUIIP5oIBY9FUEXD59JHct3Bi3BuhkoF+1yDla5kcVGiQd6A33iT9gCiCeCAaZ838slH9spkpI8NEh5TKMkmMGrBCysh9+fLljBxSyje/8DNU1Y6mufjVX7/D8uXLLcndqhnJhE0VfGmO9aJpuuTHR4pIG6RSSvDml3XkOb0XXKXiaRryNA7FQZoJRCQkP2L2AkZ25XPYPqJfNTNlZJh+iCRGDVghZeReV1fHlTPvQFXtqIoKSCaPnsnbn/zJb7tQzUiXTxtB0ZCskE6YdMuPjwThNkilVSVv1eFqLqpKDf7+QNykmUCEQ/JNm9dzAwobCi+hJq88rU/2GRkmdqRUlkmy9TEQKSP3kpIS9jZ9gqa5jDsknOlqpaSkx9kSj2akdMmPjxWhnDZp55kP7HCtfLbnsSQsKPVln1TQWXRyIyOdJ8htyQUmJexYokVGhokP8hx5KEIBCXbFnjxZJgXWx0AIKaPzH0eDOXPmyMrKSqBHc59VsogvXXI/QihoupuC2W3c/U9fSmgmTH/S4YMhcPDC8V/+0li0VFWGffGL6eOZr98Cz33BIHUANQvueT2pX/TGPbt62yc9/+lCZdzchSy6/trUDO/2QajPdPj991N037KUHl9/w45jO/j6W1+nW+tGFSqPzHuE26belvgXToD1UQhRJaWM6MyUssrd1NVf+/2HCKGgKiqqopA/bAbLV+9MaCZMf9XhfRGJZ751zdrUVvH3/NVwzyQhLMkKweyTCiCkRtPm9bxUtTHpw7t9cTYP00gUUhb3myLrYyBSaoVcunQpl194PWt/sR23WwcEz2w9xA6H5rddvDNhkpEfn0yE9MyrKi2rV4Pbnboq3nfRNYUw7ZMb/vp36je9h0rP9ywVw7t9cbYO00gkUuKUSaH1MRApb2IaVZbHJV+azIYXatClZFGHnWOKpNGmJywTJtH58alAMM+8q7HRbw5s2iy+pgjFU6Zz15TprH//Ej5+4w2yDlR5FnxJ+vDuzDCNAYYUWx8DkXJyrzrUzFtbDjNYShTPENtxboWjDpnQTJhE5cenC0xiCJwDm3aLrylC3vjJ/FE9Sf6oMZzTvofyMzVIPXnDuzMLpomHr1PGrbsT75RJsfUxECkld9MNM8Ip+BIOzHCpLkUmNJo3GILJNf1RqjERKNn0q8CyBMI8kTdljeJY9ihmXXYF087sSdrw7oxvPfHIc+ShYyyg6+jkOfIS+4Iptj4GImXk7uuGabRJ1g9ycUWnHQFc1Z3F4uEJ/iAsYCXX9GepxkS/DyxLAALH7l100Rxmj78yocO7M7715KLV2YrwqAEKSmIXVNPA+hiIpJN7sKYiz6nLAAAgAElEQVSkHAQKAgEIKWnY08yostQQvC95DySpBqLLmR+oBG81di9Rw7szMkzykefIQxUqutRxqI7ELaimKPWxLySV3E+1O7n9dx9bNiVdt7iUzreacLt1BILsXHsyDy0oBpqzBiLLmR+olXxfY/fiPbw7I8MkFzuO7eCnW3+KJjVUofLABQ8kTm9PE+tjIJJK7g0tnYz2IfZAN0z1sMG8/+IedF3ywUt7KRwzOCXVuy8GorMmEJF2vw4EQgo32z3a9Ek/CaaiIiPDJBlJ9binmdZuImWau1VTUle7C7NjVtP0lEkzgRjozhpf9LUA27Fl64Ag90iy3SNNn5w4eizOR7/f64SYkWGSh6R53BOotfsqA9EgJeQerClpzJR8VJuSdtKMFQais8ZEX5X8id89OSAIKtJs93BJ/jMEY4qGMPbUafK7e06IGd/6AEMCtXbTSajpkiy7gnAMyo10H0kn91Cj78yGpnSTZqwwUJ01gbCMHh4gHvlox+5ZkrzbbTwoJTqS+oIhHC4YwrlHWhifkWCSjqSkQSZIaw/M1XK6dBRHzpBI95NUch85NJs/33dhSMJLV2nGCgPdWWPCrDhP/O7JsDzygXpz4G3orUmHs028EevYPZPkJ44eS9X//IL6vBykEIb7Rggk8FlxAaLyI8pzs1MeTHY2IeGyTAJiBoI5CRVFoDs72iLdX1LJfcSQrD7/cExpRtN0FEXQdqqLI7WtaUvwvujP2fHhIHBREKw98oEBWFbVfqTbWF0lxIP84zF2b2jDUc49fJwxJ+w0FAylrnCo9zEpZUozazJIABIQMxAq3vzxJeey9InO9kj3GRO5CyGGAU8D52K0l35VSvlxLPscVZbHkn+vYPemJnZ/2MTnGxup+fgIS/69Iu0JPpizBgaOFh/KI2+1ANv21tuWJ4BItwm8SrA6IYR71eCLaKUZq4ak/G4XBSfaKF36ZT54+68x2ycziB4JlWXiGDPgW637joQOdBKGnixtjVgr918Cb0gpvyiEcAA5Me4PMAi+YU8zui6RMv3lGV9YOWsGkhYfyiNv3u9b3Q+56ko6qqp6WQD72ubYETeHSq9h2KkahnU1cNI5hPqRC8lv3kNee13QE0I4VwS+xzv91EF+LXfy382D2VUwPixpJpyGpDEXL4jYPplB/JDQIR0xWh9NQj/e1s2GmmO4NIlvvR6vePOoyV0IMRS4FLgHQErpBJxRH0kA+pNzpi8MVC0eepO9leXPd0SguW2obU7nlfHmM9vRxn8BZfy1zLtkMJs/6kArPQ+l5BrO3/VbRlmcEPq6aghW/Y/qdvIjReXhi+9DBZr+dwcdd1wTfH0gjIakaIZ3Z0g+PkhoA1MM1keT1P9cWY9bk70ej3cKbiyVexlwHHhWCDETqAL+TUrppw0JIZYBywC/EXp9oT85Z/rCQLZNWsGK8AMlkVDb7HrjIJpbR0rQFYX6k1noegcIFakqyK8+QP7tF9E6tIy6LQcpmVsa1lUDWMtBitSxaZIr6iq5or4Kh9Soe/Mlf+2/uxtUlVHffRR1WB4oivF6fdhDMySffCSsgSlK66NJ6qsq63tV6WCQut2m8MXZY+OaghsLuduA84FvSik3CyF+CTwEfNd3Iynlk8CTYIzZi+QF+pNzJhTOFttkvOC7qK6qChMrRtC0t9VzW6V0UTlHalt582/taO5Cdv2tnSVTWxnVx1UD+K8RmOSvO51oQgXAprmNqClf7b+72yBjt5sjj/8QVBU0DVSV/C8v7VP66di+HceWrVw896LQJJ+A9MmzEQlzykRofeyL1AFsquBLcxITbR4LuR8GDkspN3tuv4xB7nHDQJJmzhbbZLww9cJRAEybP5pRZXkUjhlMw55mxkzJZ1RZHlU+1b154gdoOJrPmCtvJ8dTBPS1RpA1ZQrr//QG/9s8FB24or4KITUUVcXV2Ej2OdMNMjd97Lpu/Och5e5duyNe+L14zkWMaW7jkNZNzWfbezzyxCd9MoMEIQKtPZj7xazSF00ZTtGQrITNq4AYyF1KeUQIUS+EmCqlrAEuBz6P36ENLGkmEAMxkCweOFLbytpfbEdz66g2hWnzRwPGd8H3cw+s7rNz7X7PM91VR2pb/U4KVnLQ6IJSap/ehNOls/ySb/CgUsuID9+mZdUqhMNBwT13c+q5P4CuI2w2owLTtIikn0DyF04nExwOypZ9jU/XvEzdME8DohBG+uS6v1P97luUX3YVIyeU0dnWxrjyGRmyDwMJdcrMuoNQs4BDuV8SIb2EQqxumW8CKz1OmVrgK7Efkj8GijQTiP4cSGYSZnauna52F9m5do7XGz0Ww8cN8fvZ6vGudpeXbAPRsKe5V0VutZ1pmTWJ2+p5gCXhB8K3mWlXwXje2LufL7vdCA8hq0OGMv6Pz/vJO5FKP8HIP6dqB+cePs7QMx1UjykymqCkEUOsud18uu7v3uNUVZXyy65i4uixDG042i+7gpOBhMgygXr7zDv9Hg4lwcTL/RIpYiJ3KeUOIKFTZweSNBOI/hRIZhJ6d4ebT9bVo+sRLZ/0gqIIZl4xDmeXIUmYpH+qsd37h6EogjFTgv/ugdW8byUfivB9K3kT2mefcuvu9XxSNJFPCsq4S7WhCi3oPNNQC8bR6P4lrR0M1U9ybOY57D911F+fDSD7z3TJmOY2xj37e2b+5rdA/Ju7MghAEL09FKknagZ0uEj5DNW+MJClGSuko1xTvbHB+/7HC7ou2f5WXchtNE2y7a1D5Ax1WF4R+BJ0YCVv3h+OdHPw71spevldbjxZwx271/Hogm+w76GfcMmZuqgIsi+3UDDyNzNoRvzTNzg8OIv6/MFGJQ/eiVAAuoD6giE05A/m2G9/yYhPdjHs9Jmwvf2Q+GiHVCLuskyQqIG+dPVkSjBWSHtyh4ErzVgh3eSa6o0NvPdCjS+3JA8SDuw4Yf2YMKqikhmFvcj/ZMMZL8mHI9288VorWsk1KGOv5LxPfsW5x/fx0O4JrLz3toS8x32R/8zf/JbJW7ZyesxIGusO0PGnVbQ6VH+yFwId2Ndygv0lRZQflpS0dvTp7Q832qE/I64NTD5yjC5UqkfdRNP4m3ivchB/2rIz5bp6KPQLch/I0owV0kWuCUXspqySlWOLSHN3ZNviIusgjco+KPkHHOPJhjMcPXja+AvEyOzI/fx9DjZOQkcBYVTEp/Kn8knRRIq6YNPrtYy5dlLSCwmT7Iswmkk6LlxIh4fs9zcdpvrdtwwbpYfoJVA9pojTgzs4b/YsRAiNP5xoh0Tl+CQDcW9gOrgRqXUjpI6uS96ot/HrAyrgf9WZKl09FPoFuZ9t0owVkt0IZUXsQsCsK0vIyrEFXRANB2WzhnsXZI/Xt3G8ro1jB3tC73KHORg6fBBHa1vRteh/h1DSz+BT+9hV34Rdr4FxlwMgFMEHEy7CMWgwX2yzw2etrN29PeW5RoFk32u+qxBIIajLH0zDh+8w9fYljFezGX/l1UDv/B8IPSQ9HDsnpCf5+zYw6VKPqYGp6lAzWw8Vc7duw44bFzY+1vzdSqnW1UOhX5A7nF3SjBWS2QgVjNgX3jWV8gVjYt5/4EJo9cYGTtTtQZcSm03hmmUzvDbG3ZuagN5XBMfr29j1YWOU5C9pHTqR1qETPTclIJAojCWfMZ2gGJvhdutsfGkPucOyyBnq8PruU4mg810Bze3i80+qqFFtlOcPoXzhYsv8n1ABcBCenTMS8k8W8hx56BgDzXV08hzRfVYvbK5j9WuvcgGf8wP5DxSIM2zSp7NdTgGM74ctjSQYK/Qbcu/PUcDxQjIaoY7UtiaU2K1e74OX9qJLiSIEl3xpst8iaajPd9r80ZbkH1z6MUgchFeeMYndhGL0p/o9xbiqMPb9+cZGSmcWMb68kOP1bXS0OsnJSw3pRzTf9bLFFPl45EMt8kJkds5wdf5kVPqtzlaE5zNUUCKu3E33S83WdfzR/mNvxb7U+QifqdO4c/ZYzi3Oo7nDmfZ9KP2G3PtzFHCikIj8+I9e3ddLY08UsUOPrx0JUki62l1hPzcU+ftKP217D9H+4vOcHjQaZ1YeJ4vO8/yO0qtbm9CRHvoXVrtFehZ5A7X+zzc2xkWyihTxyq6Jh50zHWSePEceqlDRpY5DdYS9mBpoafwndRd23NiEDtLNfeObKLr2nrQm80D0G3KH3lHAbpfO7k1NZy25xzs//khtK037/Cud/FE5CSP2I7WttJ3qMibNSOn1p8cDvsR/4vO/crxhI8W6DqqK+NoDHBlldBcGyj0fVTbSUXMa1VO9ByP5QEiJV98XAkZNyqNgdG7SKvp4B5RFa+dMlsxjRf7RLKZa+dTPF3soFifQMBw3qA6u/sIXYVz/IXboZ+QOhjyjqALNbfzx7f6wKS100FQhnvnxpsThi5mXj4vbsfrCN2ZAUQXnXFwct88x8A8/cIJUyWUzmFYxzfK5G9ydvNR4gundKgow87zhlOcM4lRTO0f2tYZlCZUSmva20rS3lV0fNjL94uI+O3PjhbBIPg4BZZGQP8RX5glG/k1/+V/Gq13UjCVkGmSoPPXzxR5WOgw5RgqVU1PvZMQl98Rt8HUy0e/IfVRZHtMuGk31+40A6FKedYurfSFaLf5Uo/8kr8KxuQmXY6Q0XC1DCrLjRuyBHm6rjPlgmF9WyK+y97JedaEogkUzClnkcUFYLfCeamznyP7gpK9reL+rYNgzL71zSsLeVxMhSZ7EB5QlUuYJRv7jnU4eVXR+eKeKTYF5Jw7T4d7urfT3vvU+67PG8uvjg4LmqV9s241DuFExXEgjxk7sl8QO/ZDcwVhIq/n4yFnje48U0XS5HqltpWm/f6VzOGcfP/z4LW6YeEPcJ8dn59oRwlj4iqccYzVIwySacLTbUEOzg2n8vtEMO9bVIfXg+9d1yXsv1HC8vi0pV5yBJO+1T4IRUJakLPl4yjzByF/oOjYJl+7UWfhZN6q+irrn19J13/2o//v/UF0uFik23rzkPgDOO7GfT4smUlNQyrmtddzpOE7FeXNQa9dGPWUpndAvyT3jew+NaLpcd29qwrw2NUYc6Lzh+BNH9xzklT2vsHz+cm6beltcji+UQyYaWM0zDfR1R4JIh2b7kn7ZrOHs3tREx2knLUc7aG7q6LW9lEY1n0xDQFD7ZJoMDIkH+WvObtyK8SW2aSCkjtbtpPrFNcx0uVCRSN3N5XWVXFFXhV1341ZtbLvhK8z/6A8IlwuxXtDxyDehu5WOY1nknLCTkxhlMuHol+QOGd97X+iry1WzHeB/t3/KNV1l7Dq1i+zPp5JDkXcBsWlILUeHHARAQ+NHm34EEBeCj8UhE4hw5plGivllhdgUgUuTqIqIyIXkS/RHaltZ89/b0NwSoRhpj77yTSq+txHZJ9MoTz4c8t//7lqe6F6DS3ez8DMNqQncQuX94hmcc+IAUnfjVmwoisAhNRQkitS44ng17S6XcbUnJa0vPE1r3VCky4144fW0btoKhX5L7hnfe+QoKmwie/gGNNmBveADtrTqbNkEI9tKufHEhX7btuQc9buto/PjzT9mcv7kmCWaeEoy4cwzjQpCYGWVjASjyvK46Vvne7NtwLhC2v1RE7oeXzkqEgzE0X85FRVoYwQ1f1+LWwoev8PO5E/n8Mmw2ewuKKVu6Ghmnapl+CUXsXjqcNSHtyNdLhRT5tmyCSklQpGga0inyzgZ9+OO3X5L7hnfe3jYcWwHr+1/jZOdJ/mg4QNshS6vzc/ElOMXoKB4mz8kOvtHbGP2iNlsO7bN29QTj4S9WCWZvpww0cgwgdhUexKXW0Mi6O52cus3HmTFlxeydOnSiPcVqNOPKstj2vzRlrHDycZAIvmqQ8385OPXcUsNIWDPWPgsqwDXyVIcNoWKay/l1vPv8l7NdowKkHmGaXQ8+wg5wztBsdF62I50a0m1csYb/ZbcobfvPSPP9MAk9dV7V+OW7pDbDusY6XfbPaKNn9z2fWaNmMV/V/43z1Y/CxhafLTt3CZikWRidcKEi9N7K9FdDlCM0+DJxjqWLVsGEBXBB6KvzttkI1n2yUTA16cuhzjJGiUM6UvakJ0TuXNeiWU8gJ+sU7+FnJwmcv7j+9B5EkoXUHKPPelWzniTfb8mdzj7EiPDwaqaVfx4848tSd2szgUCVahcn/1VhrdN9Hu8Yto5zBoxjapDzVQfdmGYxIzqfdepXVEfV6RNS4HVTaxOmHDxux8/xMm8aRRe9U8IISi4YhlHTxxi+fLlcSH3dEWq7ZPhwHR75ec4qG5s9TYfiUGHyBn5OqCDVHAdvYEfXP2FvsO8Aics3f0ajJtLzjji7uaJhPy1lla/qV+5ipIb6XuVdHKP9+VIxjljwFd+ea/+PTR6J2rZhI1bJt/CtIJptDpbmTNyDi3rsqmmx4eNMKymprPGbRtCdomKUIwTxdp9a7lx4o0RSzORNi1ZVukJkGCsUFdXx5B5c0EIhKKClGSXzKBu88sJeb10Q7rYJwNhDsfQdEmgS92eUwvCbSyRCLjroqLwUhqDTFgKRDKsnOg60unkyA9/5D+v1+1mvN0+JcK3K7nkrnd0+P3BWp2hoiH9s9k505f8ogqVhWMXUjSoyNKv/mpjld/t0ZMMyeCVd/cZzhrXeNwts7EXbAEkLt3Fa/tfi5jcI21aCrZQGm8JxgolJSU01e1Eam5jUU3X6KrbSUlJekW6JhrpYp8MNnTahABE9yRU8Q46Gg7FxpJpYfjTg0xYigbxIn+EAM/VknR5JEuD2yJe2U8uube3G7FMIc5QoTSo5j//mVPP/xGttRXhcKAMHQpOJzK3BJF/E1KooElO//637P9tDTidqAUFAGinTqEWFKDm5WErKiLvpiVpa2EKF33JL6pQeWTeI0Hti1aNSwWjjas/30Yo0TEHW+E23NKFRPLqnldxtZ7PTdMvDju3pi+HTLgLpfGWYKywYsUKli1bxtE/LSe7ZAZddTuxtdSx4sknE/q66YpU2SdDzScF/9jdmROH8tNPFFzhZkD7yjGKCrPvNoZeJ6gbNVzyV4flGXq8y4VQPdYH46Qa8XSbpJK7kpvr/YMNdoYKpkGdfOZZXIcO+e+w0ZATcqll8mgneybfjhSC3cMWMeiTT8k7fQBqa3u29/m55aWXcEydankC6Otnx4QJFN77tZSeHHYc2xGU2E35pa/OUt/GJcAryUBgI9RFPL9/F+vr1wPglhqrdq9h1YdqWLk1fTlkkrVQGi5MXX358uXUbX6ZkpISVjz55IDW2/tCMp014QydvveSCQwZZPc26T298+nI5qb6yjE6kDc2pTEDvmQfWNV3bNnKoX/+pz2R7jO55J6TQ8lTT4U8Q1lqUN//gd+AYCu47J75kkJFF9A8bLJB7sEgJc7du42fg5wAgv3srK3lzHvvMf6Pz6eM4P+y/y9osqdK6Ut+sUJglowpyZjwbYT6W2OR/5PVtrBza/pyyCRroTQSLF269Kwm82BIJMn3Reqh5pPOGTkHm2LDpbtQhdp31G/pAkOGScOYAasqv/0b97WHeIolkr6g2tcZylKDcoe28gHkt+xF0TV0AQKdrqx8WodOCE3wscDt9pJQsrGqZhWv7H3F6z9XUVk+L7J4gFCSjBVumHgDr+591XulYBtcgzJ6NV3ZC3h65waGyKl81tiKLecAS6Yt8J5crBwyhfI4J373VkK96mc7Eu2jjrd90lwsdQcMWEnI0On6LUblfs1/eq2P/TUcLBSE7KMijifmzJkjKysrw9rW/HKqw/I48vgPTd3JC5GVhX3CBHA6weEAp5PTRVM5bJtIvW0yUigoUmN+15vk68fRTp0ChwNnTU2fVwFhwWZLSeW+qmYVP9r0I+8oMYDbptzGYxc+FtF+Nryw2y+tEAG3fmd2yEXOH378Q1btWUVv9VMgdRWjo1PHrti4ZfLNTOw6l+ZXhoCmoKqCaReNpjS/ja6Hl/WSYNKho2+goGP7dg7dfQ+4XGC3M/4Pzxn3J/D9bdyzy9I+CSAUhRmLr7as4oMtlkZK6k/vfJpfbfsVOjqqUPnXin/l3hn39t4wiPUx3SGEqJJShjd5xIO09bkHVvhNP3i8R0ax2Sh57lnLL6n+xkHq19YaEoCqoN91P6XXlHof79i+ndY1a3GfOIHW2hq2zp4Omrups/sSu03YuHHijRHvqy9Jxgo3TLyB1/a/RrfWHUDwEoRR0QsBbuli1Z5VzGo4xQXu61AAza1xUNtL59bXGO9J8EsXCWagoXXNWqPoAXA6aXr0uzjr6/0MC/F+r8O1T46YvYCTo2dQWDbVz6fu+21SFcHtF4yLqFLPc+ShCGO4hl2xB5dlwrQ+DgSkLbn7IqeigolrVodV3fXV1NSfSSRQZ1dQeGTeIxHbEqs3NvSauBRKkjExa8QsnrrqKV7b/xpr963FrbvRMd5nKY3KXWK0f0skQ9rOoOrGfYrU2LD3d7RlHeS7io4d0FVomDyMoj5eN4PY4Ny/3/uz7wk1EVdLfdknGzavR+c9KodMpWbwVJzZo7zPNRdLH19ybngedQ8imsCUxlp7vNEvyN1EOMQ8UJuaLHX2KGJ4j9S28t6LNX73CR+XTF+YNWIWs0bM4saJN1J5tJI8Rx6tzlZG1tg48/Gn1E7p4NWcrRS0jmH6qVuRQkFInYn7VtGoHGDNRQqP36lSXiepLoEDh5/glo9r/Bqr4pEdb2aspzq/JRXIu2kJratXG4vUvhDCu6YRzKEUL/jaJ3e++za6ZlzZKYBA49y2z5l+Zg9rRt3A0exRMenqlUcrcWpOTy5SkAlMps5eusCQYsyfB2jVDv2M3MPFQGtqMuUY36r9lim3RBW/u+3NQ/7DJIQxADvS98ckefBYGX9kEEXFOw6u/++H+ajFgSZsgIKUEmfWYKpLjD6MvWMFe8d6ejKk26vjh+PNDwfVGxu8J3ebXTnrwuRyKioo+cNznHz695x55x3v/YNmz2bEt79FTkUFJ373pKVDKZ7VfFPWKDYUXcr7owYx8fRuprfVoKAhMEgeqVHibGLxwrkxLZb26ZSx0tkXfDum360/YECS+0DLmwmUY6LV2as3NnDgkxN+902YWRTxyLe+Ml+ydnQhaod7MyalkKwqH8TFi7/DLA716qY1r0YkErd0s2LzCnaf2h1VNV+9sYH3Xqjxrue53f3/5B4Ncioq6DjvPM6sX+9d3Oz0rDcBqMPyQFEAEHY76rA8mr73fVpWrzYWYlWVUd99lKwpU7zPyT5nOlpLK+qwPG9nudUJoJel0TGS+qKR7B48lWlnapjeVoMqdBTFxjfuuIrFl85I7JtxFunsvhiQ5N6fpZkdx3ZQebTSW32Y0QK+cky0Ovt7L/SWY86/anxE+wkn86Vl2CR0rct4DQRi0jDu/9KD3spsau5l/GDD80hxGnVwDYpixAyb0KTWq5r/h3P+gaFZQ0MSfSCxAyhCpCQzPR2QM/cCUNUeK7Gm0fLSS7SsXm30sus6KAr5X15q9Jx0d+NzVjRcaoFWZGEMHUFR/CJETMKvL5nG43+pZtqRvZQVTWR3QanxNODU4GKKF87lgsJOck8eZFz5jLh0s762/zVcusu6gSmOEQP9DQOS3CF+0ozvUOREzrw0M2LMhUqbYnw05pcWDKKMRo45UtvK+y/u8SM9EaYcE04yo2/my5mSCrqah6AoTV5f+5JbpzJqfM/rnDg5mq6mm9Al2Acd4tKZbZSPtfPHz//oregDq/lnq58NSfRWxC4EXHrnlH5xUk8EcioqGPXdR42YDx8POi6X8e5KYxhJ967dxmcaaBHW9d73mbd1Hdnd7b9vRcEmFB6XoOoadwqFdeMv4L0Jc7ls6ggu7z7M5NJR5FTMBeJTOe84toO1+9b2FD++skySIwbSDQOW3OMxqelIbSuv/KzK26K/+8Mmbvr2+QBxW6zzJXVfInfpRhenL7FnqVlRyTG7NxmTf0yYxN6XHBNJMmNORQWn88p48xfb0dxtIZMffXNrVPcE/rnCiDBYXLK4lxPHF1ZE/8i8Rzjn2EWWxB7O7zjQkX/77V5ppXX1aqSmgaoaIc6ejnDflEJvZS4l2GwGwQdrIjSGKBixIYDQdQQ6xkoLKFLjuoObuK6hCuVDgXS7ObTyGYbdfLNX4olV2688Wolb99hwEdw06aaeqj3NIgaSjQFL7vGY1PTR6n1+2SuaJtm9qYmaj4+guXVUm8IlX5pMV7vLe+kfLukHI3UwvqR2xVgnMO1dN026KexYAV9Ub2zg843+zUqTLxjJ8fo2Nrywm+HjhtDV7qJQHmdw3Xa/P7ZIkxnDTX4MNsB71ohZaJ3jcbWejy3nAENyXPzx8z+iSa1X45RJ9C+vXceC2nzPu2Ys0maI3R+myyzvpiWWHeFWKYV+3eI+mnvbW29z5sMPexL+FcU4SWDwp1uoCCFQdbfHGQO43T1XCk4nLX/+s3FgQaSdSAjfdzHVrti5YeINPQ+eRbZHKwxYcofYJzUdO3C613111Sdxu4yK0u3WPXKHRFEFeEitL9IPluZokrpJ5IBXf4/GHmjaHv2urCXs2XK017ZCd1F48gyse4/CS9oZff4E2uznIgomM7R5X1jJjJHMRrUa4G1myDvdAodtEivvnc/iaxZTebSS092nexH99KMXckntbZjEbtwvaZ9Xi2tqAZAhd19YZZaE+5g5JGPM4jxKN21B1dxoqo2/TLyEJXvfQ5E6ulD4zcybqBs6misPV3FNfaVR+ZtXCi6Xv8zjmw7rI+1gszHs5ptjS249SyIGQmFAkzvE3znTdrK754YE3dT13T1f2lCkP+jS07yz92MKh44FoPj0JJqG7udUXoNldR6L53v3piZ/22MISGHjRNFMAE7s6qZml7H4Kmbcz+ghZ6i4bHTIP7RYZ6OCMbvU6dbRJd5Qsn+5rMdyubikh+i3rN/DxbVf9NC610+BzAsAACAASURBVJfDxrKX2CU+5sU3bFw69tKIgtQy6Jl0NL+sEMA79ejx16txunUUIZh88X3MOLGfnUUTOe/EfhQpUQEpJQXuTgqvvZTLsmeQv8loRMq7aQlAjzTkdnsXc1EUf21f143q3rPw25eEY8oyfoup3c5+GTEQb8RM7kIIFagEGqSU18d+SPFFLM6ZieePsKxy+0RQ0tdoWZ/NHHkt54urAVCkglAEo8/P5uKRM+AMVG07mNjmm0Bx2upnn00bTw+mcW0boz+v4qKbJ1keVyyzUU34avF2m+IlGDBJZzDzy24lu66DrAM1XjeNH7GP/BgAt3R7I4pf3ftqWBHIZxt8R9Y1dzj9SNymKiAlbt04WfdMP5LsLZrAnsJSVFXBpgjcNTZjsIndzpe+toSpo9zUfeXfaPGs1Zjkbi8uZuTyR/wkGG86rEf+89P8AyQc34oeDOlo1Ki23rEDNe+cldbHQMSjcv83YBcwtK8Nm9qb+Lf1/waQ1IoqWufMlV8tp+VoB8cOtcXhKCS6lAgUFBSEFEiMCAF0OFLpZM2ObZbSTnau3SvxREL40+aPZtcHjei6RFEEl945lZPrP8K54U3aBo8FoVA8bwqdE2bx+cZGz3sUfOBL095WXvlpFRVXlXDRLZO890c6GzUYgmnxPXKNziyXjcvb7UacjfdYJcdmf8LerK2WIw3cnkapV/e+GpalciDDl9BNItdlT+u/SeIutycXxvP/qiKQUmK3KTx2fTnNHU7PyXc+Oy8o4bwT+5l81aWWDVKB8xkCu2FNvV8dlkfX57uM6t5CwvGt6AWgu92MEzr3zBB8eJ6d2271xA50O2PS2n3tyLNGzOp1O17bBNqeA3/Oc+Sx69QuBAIlO/IZqjGlQgohxgJ/AFYA3+qrch80YZCc9P0eUlBQqBhRQV5WD2ElgvTN+Z1uz2XlpXdOiWixrXpjA59/2IhqU0Bg5LJE9Lb1zlEEk0aDkKkwPNq6lN7XstmtCR+sF3I7tm9n5zd/wKncUgraDzLjV98DMBwwHreL+YdmWj47TjvJGerAkW1jx7q6oLLOlLkjufKr5b1mo067aHTcLaO/fncfP3+rhnO7VK7stKP4vGe+i6e+c2TfP/y+5SAT6JlSNdCJPlRl7l+NG1A83znpOUEjJZrem9D76iT1uqw837G8JUtoWbXKq78Pv/9+IwIhSCesGe7nJ+GYFT14Yr1BeP40JOBWoe7CMVwyoRTNPpqGiUPYaq9mTtnVzDr3LqA32Vrdt+PYDr7+1tdxak4cqoMHLniAn279qff2U1c9BRDzNr63TduzrwU60C2273v7ZOfBTiWSzz/Wyv3/AQ8AQ4JtIIRYBiwDyC7N9ntMR6fqWFWv57y852UWjVvEJWMuYdepXZzsPAlET/yxNjWVLxjjdzLwJUKA0yc6OXk4WJa+UQ1bCB6EHIvoI+2YcLt0w/JnPh1QbT2aPsCwQV2MnZyHnldI8456KKhg9JHNDG2vC+l2GVXWOxWybNZwPnp1X6+QMTAWZYsnD6Or3RXRbNRokJ/jYJRL4YpOu987FuiK8Y1ECDUwPFzvfH9BXzq5VWVuVuO6lOjSsC46elXl9LqKCgeBU7T2Nu9FXa2guEHxdMMevOcebyRx6XPPkVNR0UO0Y+Yw6wff59jCc2j84G3GFE9jeH0bza++gu52oxk1FjYd8By7XYOyDxpo+uCwcRCqoGOGwk9m7eTBEecA/mRrRcBPXfWUN6dGR8elu1hXt87vduVRI7I81m0CbxufiOxlgfbCikL6QNTkLoS4HjgmpawSQiwKtp2U8kngSTAq93D2raOzvn69VzP1xct7Xo6q2o9n3owVEQYSPkD9sUacjSoC/xNuyIq9DwReaPlq+khJc3sWzTu6QDQCg2H0JTSOuojxjRtotZ+Lu7YV8spoGJ/PmLx8ckK81qiyPG75j9lUb2yg6u8HaTvV7ff4J+/UM/PycWE7ZCKBL2E1dzg5x2k4LsK1O/oSvelOCmWpDPTOx5Jtk0gEErnZ5u/WZQid3KgFdGktr5iVvRWJh0PqlnKDrYo5Nxs699erf8b422HGYRs33fYwh3d8jnA6USVoTif7312L3ryXtatWsHOszu/GZxuVbePPcI534lB38NS/PEXNDIXqt1/isxJQhcLddeMZt3EfNg1vXo1ZStk1uHyHzsKdnZxseILjC8v7JGDzd3CoDq+18oqSK9h2dJv3timbxLqN721VqECP7Rl6V+5BLv5DImpZRgjxBPAPgBvIxtDcX5VSfjnYc8ZMHyO/+Osv0trdyvZj23s1qcQCFZV/LP9HzrjO+FX6Zj7JlO6Z7PpDu7epKRHygS92HNvBV974CoWnxzLl+AUMcg4hy53LOErR2+Kx1BFY+XuuUc3LVyEst/Wt9H1lFAju0T9S28q2tw5xYId/Lk3OEDvOLi1quSsYfDV2h03hwcljaPvouPcUKYRgUYQ+dpOArCyVgQgcWQj0qa8Gvk48HjOJvKiwid0tOzhytJh3dgzyI3LNcRA1pxZ3Rxmy04iSkBgSi21QHWLQflwdZdA5Hoe9h9CLCptoEzVBr1TCPVYILT/cOPFGXt7zst8QjSE1jUz9/ovYNENOqf/qlUx47l1wuXGrsOIuO+VF5ehVO/msBGrH2fjXin9lzsg5fP2tr/t52j9d/xILduos2ilRNVC8n6no8eIDqCqvzxW0Z0n2lDp48KvPeI/d3N9TVz3VpzaeKs39+4u/v1vr1CLKaojLJCZP5f4ffWnuvpOYfC+bTSSC9MH4mBWhcHX2TUw5PpfOarsh3KmSgi+eYV/2TgQirlr/0zuf5n+2/Y+XQBQUHp3/KLdNvc2yyu9qd4XQ8i0knF4ETnAXTB8IJPyScwu9L9vV4aJpr0WEKpA/OoeWIx3GoSgw78YyZvsMRokEvhXpptqT/PytGnQJM50qV3U4/LYtv7SYRXdNi+p1IEKiR0UKiS51HIqDpdOXep9jV+wsmbSEGyfeaKnXmmRhvmY4j9kUO9cN/z7T88/j8dercdsOkF3ytDEMRdroqLsXvXO8UakOOsQgn8dcDcuQnSVousQxuJ7scU+hSTeqsHH9iB9w0/SLmT0+P+SxRHKsVuQ9b/Q8NjVu8t6+ZfIt/GX/X/wIFOAnz3yVKQed7Cl18IjjRtSnVyF0HU3Ahgobl30G0kP2P/nyIB786jO9CfHY53yt8sc4gWkNkn9unsSs/OF0NTtofnMj0u1GyJ6/HAkePcdG/i23knfTEvaMETH1kiQL/WoSk+9lsy9M0hcIphVM82rusRC/6YH9W+crNLW2coF+HQoKuqbxxsfvs33MOsCwzF069lKgp+o3z5yREn+eI8+PMO4uv9t7qW8l60BvaUdrbqFj2zYc3a3YdSeHxl3usYmZv5j5zTVWmIQQXs398DGV5qYO/xcQoHp895rPTDNfaUdzy14VejBkD7Z7/3JikWQCK/XHri/HYVOY2i64osO/L0FRRNjZ88Hg+93z9c77ZtuY0NC877dTd/Js9bPex5y6k5f3vMxf9v+FBy54gHV167xTqrq1bv6y/y9eQvrNJ7/p9ZjWOZ41uz5kn/NVurRuQNKtOXlp5wa0ZsVYzBxaC8KNEBKJG1tOLa7O8dhtCmruAe9jCI27LnVx/fgL2VR7khNqLa8ecCPRkWiMH3PEK7EE6sp+QVt9PB74mESGlB9unHijN/vfl0Af/OozVB6t5OaRc5jYIKl7fi260wk2hYXjFqLseNcztVHwyKmLGL16Cx1zJbMqKnpONJ+8CFIihKC2GPLnTiP/yp8AkLfUWJRtefll74hOAXisQF7Xzcibb+aum5aQk8bEHi3iQu5Syg3AhnjsKxjpQ/Bqf8exHejoYclSjUP3oSka6AIpoMvWsxDq640OhOmV9iV8U/LxvXwyTwIfNHzg9/x2V9/Dy4e21nJu21Zy5l3gtZQdr37K6zKYtHgqh9Uyzry7gcGnD3FmaAmDL7uM0edP6GWVPFLbyqv/VeXndhk9KY+LbjbcSmYsgy4Nm2Qg4YcDb0UvYMZlRlNW1Rvhe/TNar2xpdOveam5w8lvr5lB9f/t8ds+EUFggURvpnAGc9oEQiJxak5j/KHs+Q5KJGv2rWFawTR+uvWnfqMJJZJX9q7mhQ801KLXPCMKPe+9VAwZxdTHO8pA2pAY1bnonsSd80q49fyx7GsdzBM71iOlG4fq8Awmz/dU5wt4ve75Xjow0EtXDsw/D/V44GNW5D05f3IvMg/8m/b7Ox+B3wIsQN2bHyJdLlRVRXljI8fd7/ayUVZmZ+MWAikEmue2+Spmx232OdP9O2ADffQ+zVIxdcSmIdJ2QHYkMC/VTJL11dxz7bk8//nzfnno049eyCUHvoiQAk3ReP2cX3N0yME+X6enxT34fTZhdEa+W/+u3/3nFp3LLZNuCZpPbhXSBda2xXAHKqz6z60cO+jv0b/1gdl+JwBTZwfY+NKeXtuHC6F4rJthavm+1bpNESAEmmY0L628dz7samXTmtqe/ScxLyZcS6X32DyXL9Irn/VIccVZ59HQ/alRQXu+Dga/CLT2Sai5+4yq3Hs1ptB9dAlq+4V++vjulh24OyZ4pRXfY423/h/rc+MB8zvuamzsZaMsum8ZHW++wOY3/8hvCg9TMxYcip2nrn7G8njMfYX00QsBdjuDFyzAVlSUdkQfjSwzIMj9/2/vzOOjqs/9//6emQRICCHsAgFEQEIUiUK0UpfihhVFxd36U+vy6l36a63357Var7e1ePtrf231dW9vb93QWrSigOAuiyCILAlBYggECEsghDWEEAjJnPP9/XHmnDln9kkms2TO+1VKZjLLYTzzzDOf53k+TyT8pZ5dK5vJ2ThC72IRkpaSPewf801Ub+TOYnRkPJf/Ay5pyDMD9JG/vMThF18MPIk7sRmnatV+Vsy1e7gXnJXDtPuKQspCC35Xbjvn8/r1YEBhHgjYW3lUH8QSIIO4wfpj1fKNoayd+05Q19LKiaY2Pmw4xn6XhkvAnaUjGNa3l9mxsWbBDio+32s+lv/gVKKwBvrV+1fTprWFvb2UAqQCQgPppu3gDLIHf+jVxb1WWkLVNfKDN5I1+AMQevubaZUuXDxd8p/cOTGzjK6C4d8zP2LOa3Comt0/+xWaJvC44Pl73Nw+69+i7nCy9dH7B3kDt5u+s2alTJBPK809kfhLPQ1ZTSyqNIaaFG4ouo7iy35oeyNbNXfjjR3MhjZYNh8OieTsunaG//YVDmkKShRWuh09uYovG8a2dQ22PvXGA6dY8LvyoFnwkNH55PXvyYkjreZ1A0bk8f0fTQR8NYGtaw6gIX1zJVFo+R6Pxsq3t6FqEgH0Bm4nmy9y2uktFK4dVMCVU4ebz/PN0jrbsfXISc6pajhVrq09yqUX3MZXDZ+y4sBHSFS/TNy4hwSh4Wkuov3YFainR6KdGYI7dxeelrPRgKycWpS2MTx79fepbpxCk7KWlQc+0JeN61s0qGlZATjB3b9nPqekhCM//wNo4JKACkV7NHNvajTJkNUlM2SQ93jSXrLJiODuT6ihpkmjw+v9/u1J/pp7sOxu7D5J8V7JlhGCGu/e0OK9Uu/NldFZ6XaGS28dE6C9Swkr5m6jcsU+hpyTb7aELnmtyhbYAXL6+DpVTJdN1VvUVaD4sqHk9etpSjtWLR8haFc1s5ClSIHi/TBUvGLGtaf1x9/29x0M1AStLe0c3H3C5j+vKInbphTZb2UqWvYwsvI3Il3NuPOqkVLaOlCllLh7b9V/bs9DnJrMM5f9OERP+QhgBu9uu4DZa2ejono1+fmM7zc+ZfvtE0mAY+WV18PiDaia3kq5daSbWwZPjnnpt3+QPz5/vt273uJvc3z+fIY88wsK7ryzK/+pcSUjgzvEPtQUrtBrZdOhTTy75llqm2oZu0/yb2+ruFWQWS5+dadg63BJ1Qj966SiKShxytJDMWR0PlfcfW6APANwdH8LR/e3UPVlPT1yXJw5pQbcxr8zxd/W139WYMjofBjZm6pNBzmcI1iy5SDD2xVaFcnVrdlIb+auGX4l3g8dw0nTarcAXbdNqSNTnabfyumRaK0j9RpDy1qyBy/CuiZQCEBouPO2AODqV8ZO7TA3FofuuLr93NupPlbNuzXvAvpAy/PrnmdswdiUbtFLBjnX3UP7r/ax8PO/UjkCdhfqYSzY/gHr+ylUVu8f5D1HjtCyapU9m/euHWzdUp02WXzGBveuWKJ9qqKC4es3cE3OSP5CrZmhuyRoHpWivZKtwxV2DFf46l+nceep8+KapYfCkF9WvLUtpCdOsMB+1lh7y2Y0tr7lexr50aeVvkJpD0GD20OWW+GxWeM4sbeZupZWRuT25JzhfVg9bzuqqv83sGbrBhMuG9rpImq8pjqD+62Mo0VMZdvxTQzoo/H5vnlB2ykN47Jw0683nXOTrVNHlarZTulgoW49G7KqWHipoicJXqvfcSGkTQixVczvfWdNrgxd3tpKiaqmVRafscG9s34zAbtFLSfPtCw3K+8UZoaOCh4Fvh2hyzJZShZTr/shAxL4pi2+bBj9h/VmzcIdIYeSbAjMtknQA/v6D2vNRSVWW1//ASQj81U1yZ2lhbZCKVPtT9N/WG/21zTSMzeL1fO24/HaBiPA7VZi6mmPlI0bgbxd9QXwYO6HsfutjAF0C+e7D10ftJ3SsDmYvW42q/evDmqZMWnQJJ66+ClHngmHdy/qZDe4hwygXQhzb2rOoEkhpc1Ys/qQrZSQNll8xgZ36LjfTLAswHryCI/Kva0X8sw5FfzqbhfFe3UpZvtwEbjnMYEMGZ3PrY/rXjHfLKujseFU0NudNSafS28dY2uZNFw1AW8BVR9aCjWAZHiyz7pweFh/EutAlzXQR7I47kg2bg/kuuyS5Q7ufthRvxVDvrvpnJuCGpepUg3rMx9Mnvn12l+bv8t4jL2oLlfQKe1Q0maohgUIn9Vbd9AGZPHz5tG0aFFEbT9ZZHRwj0aaCabTBc0C/E6eMxPHIFo2sX04bB/uOwFdwmXf85gEDJdLo/vl2IEWWpvb6Ts4hwuvGxkQUK2LOBBQWFRAn4sGMn/P4aADSME82aMh2ORuvLJxf1nF5VK47aLhzLpwuPm4sR5vOIIZl/nLNYbP/OKdi20j/v7yjIbG7LWzqT5WbVodZCzevahlvbL1ASbwbWAK87oE67oxiJTVg75opN+DD3Ds9Td8WbyUyDNnaHp/kRPcU41I0kyoT/RgWYD/ydM6TOD6NHDS8b4J96XMmzOUDYI/1iKq263Q56KBNl3d7VLMASQjQHYkSHZkM1As2Xg4G9t4BfVg3H7u7YwtGBtSrjmjnmHxzsW2ac6nLn6KX6/9tdl6a+j2C7cvTGnXyoQw6S4mtzfhbt5Eu9dJ0X/KNhjxyOqHPPMLWrdU+zprpOT4e+/Rc0JRymnwGR3cwU+a8WjUvPk5fW4ao+voIT7RQ2UB1pNnEgS8QQWCPj0iLqxKKcwiqrd9pPCqYVS1tYbX1SPQqU4V7/93JhvvykAeimByjTE0J5Es3L7QZl9hBG9/q2KP9GRuFu/V21HboFcODBkQl4eNJatXjzdx1i//HcC3AlBVU1KDz/jgbkgzqkdDeNpoW/EZX65by/gfexgU5hM9mrZF6xtUkxrZruyoMoxUYn9No6m1q1Iyd9UubrhjfNS6esf2dEbbqZK8bLyjWOWa575+jndr3jWLrf4dNUbG/8HOD5i/fb5poaGi8l7NeyzascjmStntMfR2qVKWreCRatSyTCRizep7TvBz37Vo8IN//mTIhd6JJCPsB0Jh6OknR5Swd/1u2lZ8xvYxt6EpLlyK4Ob/M4U+TbWdHixKhBdHV1C+p5G1S3ajlTcipd7J/V5eG3d/f6zZGWME1I5k41YrVmPFm61TJSs+m4FSFcM+12oqBro/0Zzpc2znSriFI27hzgypxpK5b+qVw0NDBtLutV5+9bpXu9znxhoDjvzlJQ7/8Y+BN1YUvcirabpdwhuv64/RyRji2A9EgdVEyNi6LrKzOf/nT7JuXT6a4gLhQkPPWhkX3daicEQ7AJVKlO9p5PE/r+Pm41m4vNd9kdPOkZ762jsjeK+tPcq2hua4ZuOd3QyULkwaNImXr32ZxTsXs6BmgdlR45EemwYP2LL493e8b9rtGrfPGKlm0l2AgFFToOL3IAPnM+JNsKw+p3QKomdPXbIRwu5B4e2okW1tHH3lVVq++sqMM4nM6jMquFuLIyiKWfU2tLTxP76L3Yub0ACXW6Fnbpa5/NnlVpj5WEmXbW5KBfz71Ye0ClzgtQyA8QW5fOeiATFr412xp7O7YHzwF/UrYva62ab04q/BW2974zk3Zp5UY9XbXdmU5ffBo3nMXQ2dlWVixV+nBz07P715MyeXLTNv5zl0yKfZt7XR8KvnzKx+8NNPdWmgz4jgbrUPNV5owOyVNbS0ASUl3HyuzwbXaAGU0tcHD6HX0aUr5Xsamb9xH++V78Oj6ivzbhs5kDxNoO8g1iWZ1SdOcnB1s5ONdwG3n3s7W49tDanB+3ukTxo0ifH9xtukGok0F4h0uyBv0dtR28g/cQBFKCAJ6kmfCAI8b7zDjCdXrTKXf/e9bRYHt2/XrQzAltU3/Oo5Pdt3u7vEnKzba+62Via3Ww9Kqqp/ckb4imQM76iqhsulW9aunre9W2XyxhDSmXbNDNhDPQp3nMzGBWhAZbbKlmyVereWkdp4ogilwbuEC5dw4dE85p5SYy/A9sbtLNyxkC1Ht9h2FoB3UbTfKsC0xU9vf+SsQZzRPCm50DzY9HqwrD4At5t+D9yPK6+P7dtATukUci+8MLX93CeNHCnXvP9+QirIocz++952G1lDh0b9Vci60GJ/TSPrFtXadoYa16drJv+nL3aY+0oNLm51891WNwoCDcnXOSplvTwxZeMOHcOwnbZq8OCzllZQUISCJjXciv7FO5gVtRVrkC/qVxRyYUzKUrdez9x79YfTR3mFE/xn7fu2pdsPn/9wso8yIqcqKthz/wN6Vu9yIRQluJ+8ooDbra998XgQ2dlMrty89aQa24LshMoynkOH2PvgD7u0qGAa/ixcCB4Pwu0Gt9vM1mP96uM/6ONyK2Ym3x00+UtG9zfbGoUikJqkVfjyRsWl8Oisc3l8eK6TjScAqwZvtNAaQVyVKgJhyjDtmv5VP9I+AUOuMSwNFBTcijs9Mno/rZ37FzO5RzbZez4OuSYwVckpKWHkG6/bdPoAWwPQE9H2dv2/qncKtq/L1T/W50to5n5ez17y3dGj9U8mTQNFoWdREX1vmxWX6S5Tgjlzxvdp2IFsPRzhMvkJ3/V5m6dTkLcWUpv3tVD11nZQ9d2ql989LiGr7RwCsbbQAuY+gd9u+C3tWjsuofcxRcrcw5GtZKd2kF/1e1g+W9fahYtNlz5C2ZCx5Gfnp983kBA0vvOO3ZzMm7mjaaa//G27d8mq1lYllsdNbEFV4Avs3j+tlZU0VFbSVldnak0dCcCnKio48l9/sgd2ITqUrYcjVCavKEJfUuFdKZdOWbzVLqC8ugnhLZhKfM6PDonHv4XW+NlaXAVsS2QW7VgUU7A3MnqjAJtyso3XS8bU2huW0Fb/KdmubJsfTzpjmJMZLdqGqtH0/iKOz5sXeZ9lCBIa3N2DBjHkmV/o/eWt9o0/x16bA4DIzqbgB/dypnoreddeEzGj95dhjE8+4XaT38XrsYaMzmfmYyXsr2mk+VgrW1bVp3VnTUNtE83HWlG87YuG82PaYmi1oy6DwtJkH03cCBX0QTcd60iwT2nZxtvbXpbfhzav1t6utSe8/bErCTUh27RokdFpE3OET0q3zKmKCo6+8qq9cmzfU2Ze3fuqq+j/8ENB/+FBZRhFIfc732HAP/9TQkd/072zxjx+j4biEoy/9KyALUtpRd16eH0GqGdAuOCGP8DkB5J9VAnHfz2kv3FZJNzCzS1jb0lOkPfrkPmg5Gbeb/ga1TuV2l0y93AYjSGj/vEfYi6oJrUVsvGdd2j+fAk9isbT+Le5vsqxZs8yRM+e+tZz7GO8B579d/vXFiEQPXokzV85nfX48k93B3QCXTR9VLIPq+N8+BiUvWa5QsCMFzIywFuxLoH33/cbDkObT6hs49XbN2W7eGTIINoUBXeq1wi6iLSzHyi4805Tdsm76ipOrd+A2nyCY6+8arudbGvTt5QvWmQb421auNAX2N1u+s6alVRXtnTW4/13o6a1HAMEfouV8OFPYMcSmPqTbiXTxIJV0jECfTTZfFJkG4t3e5vQB+pUqTK099CMCuwdJWUmVK2aU3ZhIcfe+CtttbX6LzUNz5EjNuvN5s+XII1N5ULQd9Ys04ozFeisHl9fU01dVSWFxeczdFxM38ZiJprdqGnHBfdA2RvgrzVv/RC2fgzjv5/RQR6C2xCv3r86okavoZmTsB/s/MAcquqSDhavd3v2yW9ol1patT4mm5QJ7lYK7rwT9XgTh194wdTh3QMG2Kw38669hlPl5ebl/JtnJvuwAzAy+YbaJrZ93RB1f/zmpZ+w7LX/QdM03FlZ3P7M7C4N8NZNS9bdqGlNYSnM+CN8+FMCs3hND/LbPoYb/pjxUo1/Nl92sIwTZ07w5pY3w2b0Ekmr2srsdbPRpGYOWVmzeqBjjqhevX2TG8p69eKJi/+Bptx+qdPFkwakZHAHcPXN90kuUtJzQhH5N8+0ae5G+1CyfZMjYc3iQ3nWbF+3nO3rv2LImDFs++oTNO9Qg6e9nbqqyi4L7t2uQ8bK5AegYROUzQn+e6k5Uo0f1kA/bcS0qGQbq+2BNatftGMRoPfhxyzj7F7FJjc8Mrg/bUKQXfMWL09/3QnsMZCywV093mTroFGPNwU16knloG4lYKWdbKD91LcIIdi96Vt2V3wOQGP9Vtv9FEWh2EOcEwAAGyRJREFUsPj8Ljkm/w6ZCVOHpneHTDAuuAc2vglaGE3ZyeKDEmnZdzj8J2iDGZpBmKx+1GWUbfwvXWsXgnYpu1XrYyJI2eDun7m7+qZ3wLFq6Ef27qb1+N+RUtc1d1cIv1sL7/8Ek298oMuy9q1rD+Bp149B0yR5/Xp2r8AOejb+/d/DR4/pmXoojCy+8h24+pdOFm8h2LLvYEtDrAgECgpCEahacNdKCJHVe+cT8otuQTmwFICsNNxilmxSNri3bqkOezldqK+ppmrlcjYv/cS8TiiKGdh17G+SUSXXUF9zGpRhbFlTwNiLm+IedBtqm9i65oB5WRGi+8gx/hjZ+Ec/i7zcYc8aePVaXaa55pddfmjphrE0xNDl36h6I2jxVSJRUXFpLi4cdCGVRyrNBSOhsvoPdn7AE2Nup2ndn8n3tPPbfn1RFYFLuHliyhNO1h4jKRvc053NSz+hcvkSDu3eaernBtKvj9/lzmLcpdfRsGMHY0unklNQwsE9tTZNPt7BfX9NI5rq+1AZP7WbyTH+mPr760Qe9pPw1QuwcxkMnwIX3O1k8hasmXxhXmHYTF5FZdPhTUwaOInGM43sbd6r+/57fXGsAb9NbeP5bX9Dy89FkRJNYNpgNLU1JfKf2C1I2eDuv4A2YCFtCmJIL2dOtbBh8fywtxXe1VyjLyyldOYsm/TSUNsU4D5Z/unuuA5AWfva3W6F8ZecFZfHTWmC6u8C+p0Nx2oDb99Qqf8pex1Gfgd6FUDvQU6wt+Cfyb+55c2AQK9KlfJD5QC4cDFhwARuGXOLbW2g4XipSRVN6DKlAggUp/2xg6Tsso4jf3mJwy++qE+rCkHfO+5IqT52K/U11WxYPJ+d5euRUuo+zBFe14lXX0/xFdNC6unGtGvP3Ky42xgYhVSPR9+6lFHOj2WvB+rvQoHzboO9a6Fpb+THULJg3HV6oB9yAZw+qnuNN2wCREYH/2haKQUCl3Bx34T76NOjj9kfn5+dz2/X/4du5YvCE+fe67Q/eknohKoQohD4KzAEfVLkJSnlix19PH9ySqfotpdtbSAlTQsXJnX61J/NSz+hZt0aBo0azcaPF6N6fL3hET8uvZlJuEKp0V1T/unuuK/665Z97dFi6u+WAC81+PY9vVumcRd89SJh/ytq7XqHjYmw375iLjzwYUYG+GhaKY01gnOq5uiFV6FwxfAreHDQxbx84BBl2QqT2zQmTbsyI1/DeNEZWcYDPC6l3CiEyAPKhRBLpJRb4nFgOSUl9L3lFo6/8w6gbyQ5tX5DUoO7URw9un8v+6urANizuSLgdorLxUU33Myh3bW0NB3nyJ5d9htIybdffB42czcYNq4grgtCunVfe7SECvAfPaYH+Ic+h2/egsM1sOdrAqZcA/D7IFDP6G6UGR6YommlNBZcL69bzpd1X3BrryxuPHmSSW2q8xp2kg4HdynlAeCA9+dmIUQ1MAyIS3AHP51d05LaDrl56ScsffXPAcVQwOZkKRSFq374IyZefT31NdW8/cwTQR9PU1U2LJ7PzH/5RdjnjWYACqLL5DOirz1aIgX4GS/o19Wt1wP9ycO++9Z8pmfvJn6ZO0DriS468PQj2lZKD5J383JZkJfLfc2n6MMJJh/alPGSTEeJS0FVCDEKKAHWBfndo8CjACNGjIjpcZPdDmmVXso/ej94YAem3HgrZ06dArBl41UrlxPu6/3O8vXU11RHzN7jterP+sHQbfvaYyFcgDd+X1gamD0aAR/h09z3b7RLNWv+E8bf4GSefvgXYP+65a+2CVcpBB5gTp9cRO1CXLsWp9wC7HSh08FdCNEbmA/8VEoZkK5IKV8CXgK9oNqZ5/IcOdKZu0fEGswbD+xnx4a1QHDpRSguBp99DudPu4aJV18f9PFajjcGXmnJ8qWmxWwt0JlMvvs5P8aBkAH+ZzB4QvDgHCrgb/vE10cvVb2d8q63uuzQ05VguvzRY9tZebhCF2y8NSlDm5+9bjar969mQK8B3HjOjU4mHyWdCu5CiCz0wD5XSrkgPofkI//mmRyfP9/cI9iyahWnKiriqrsHa18MGsy9J1yo9sVgNB0+aLs8rKiYoWPH29oke+XlxXzMHcnku6XzY7wwArzVZEyqemH1rrnRPUZhKZx7vT173/qJHvSd7D0kZqBf9Xve3foFz/fLR0XP4A0MTR5gwfYFXD78cifQR0FnumUE8CpQLaX8Q/wOyUdOSQl9Z83qkqKqURytWrEETdMi7imcfOOt9MjJjdqCd/PSTwIKqf2HjTDlG4Ndm8pDZv7REG0mv/7DWtNqIOM6ZKJh8gOw43PY+pHvuq0f6a2T0frNTP2JbidsFmA1WPosPPhJuHs5AIy6jNtX/o6xB49S1qsXJybexpt1SwNaKT3SYwb6+TXzefqSpx3JJgSdydynAvcBlUKITd7rnpJSftz5w/IRz6KqkaX3ysvjizdextNm2UIj/P1dQHG5GTRqdFjpJRRrF74bcF3xFdO8OryPHWXrotLdwxFNJu/xeAOOwJFkQjH1p7DtU4tFgbTr75EoLNV94q3Z+541MOd6x68mHMau2+m/YdLpo0zy7rydFmGZiIrqSDZh6Ey3zGr0NoEuJVJRde7cuTz99NPs3buXESNGMHv2bO69917bbfyzdCGEnq0bCIHL5UbTVKSmIRSF86ddF1WrYjA2L/2E5iOHbNcNKyo2H2vzsk9tpmhVK5fHzRwsVCaPBAQUFhVQOmO0I8kEo7BU37UarsAaiak/8Wb/lm+Ce9bAa9c5rpPBsOxJxZUN9y82PwSDtVJ+ue9LW6D3l2xuHXurE+S9pKz9QCisRdW5c+fy6KOPMlbTuLaggPUHD/Loo48CcO+999qCuuqxfPILgaIoSClRFBfFV14NQOXyz8yb9BkwMKaAW76nkbW1R7lkdH/Kg2Ttl9/zAKAPLp1z0cXsLFtr/i5o4bUThMvkncAegVAF1o8fD11gtVJYqgf4r16wXy81r6aPE+Ct7F6lB3ap6n8H6W33XyYy59s5Qa2HPdLDuzXvsmD7AnP6NZOnW1M+uPsXVU+uXGkWVZ9++mnGahqvF44gy+v5/C8nj7PgT3/k/MH9AqUX0LP0rCy+d/8jnG5uNjX0+ppqtny5DNXjweV2x+ShXr6nkbtf+pp2VTK87SA3+2XtfQYOsn1QlM6cxa6KDaah2K6KDZ2WZkLhn8k7gT0KggV4TYt+qMZwk/QP8MYe18ZdjuOkgXdPqpm5j7os7M0nDZrEi9NeDLn/1X/6tYerBy9f+3JGBviUD+45JSX0vuIKTi5bpl/h8dD0/iJySkrYu3cvPzznbOoG9yNbVWnq2YPp/UajCWGuqTPxSi/FV14dVG4ZOq6I25+Z3aG9pfM37qPN67A4+cgaQwExufhme8Fn6Lgizi6ZYmbv0Q40dZSARSEOkTEC/MeP64Hd3SNi4LFxzS+h4Ozga/6+egGqFsB3H3eyeIBJdxGrJ0+006/tWnvGLvlI+eAO4B4wwHb58OmT7F04j+kXTSR35FBqFLv07xICqWkB0kskDd34XV1Vpe1yJI40nwFgSGsDw8402H6X07cgaDE2t6+9oBmPwqpDnJn8gC7F7F6lB/ZYC6Jmi+VjBFgYHN+rZ/EbXoEZf8jMYqu/3n7B3TE/RLjpVyXDHSXTIrgbHTONOT3YX5DHvroaZF0NV40ebq7hs/0NQaWXSNTXVPPuc0+b0ky0i6mPn9KlnyuOrATsWfvQsecGvU/xFdMCCqtdmb07dJBgA0uxYHxALH1WL6z6c7Ayc5eDRKG3x4J1+tVwmnQ09xTnwJ5dVA0fSF1Bb324wZBbpERIn0OFlBLhcjGxg50udVWVqB4PUtNQPZ6op0ePtbRxydGvGdh+LOB3U26aFfQ+wQqrTvbeTSks1XvdlzwbRIcHcznIjmWZlcXHqLdHgzWTz3RSNrjbetK/WYunn3eS05KlK1IyYf8R2twusjVJ7j13Mfrq6zocHAuLz8fldsdcVB3hOUTRCb3V37CQEuiBPdyxlM6cxc7ydU72nilc80vdb+ajn+lLQPw5WAmvXgMjL+3+ffGW3nZOH+2Y7OUQlpQL7iF70o2gLiWKlIxQFc7aVU/BKV3v7vfwQwz+x5926rk7WlQdtWslArscM3Dk2Vx+74MRn88/e4/WTMwhTSkshR+t1idfV/0++HKQPWv0ID9+hi7XdLegF6a3PRzWJfPO+yMyKRHcI02OKoqCVFWElAw71szwxmYzqBu3ceX1icuxxFpUra+pJqdxj3lZev9c/fA/RvV8pTNnUbtxvek42REzMYc0ZPID+p/5j0DlvOC32fqhPhDV3fT4GLV2/4TP5XbHXE/LRJIa3KOaHPUWRg+9/z45a9bZg7qBouibm+J0TLEUVb+c+7p+qPga3g4UXhL1CTd0XBGTZ9xiMxM7c6qlYwfvkH7MehlGTg2dxRt6fHdqnYxSaw81hKi2t7Pstf/RF23H0PiQaSQ0uLccb6S+RrcPiHZy1CiMnuo/hD3LVwd93N5XXhk3p8hYiqr1NdXs31pl/zcqOWw+67sxPWePnFzb5bIPFzJmSvQfEA5pjpHFhyy40v1aJyP0tm9e+ok+q6Lap1ARAqEoptlfLI0PmUZCg/vJY0eZ98unAGkP6hByctQgp6SE3t/7nm+YycDlov/DD8XtGGMpqi555b/1Q7dc19BzEP1ys2N+TqEoNmnGKaxmIEbB9asX7O6UVtK9dTJCb7uRrVcu/8y+HMcyhDj47NF88cbL5nu0V14e6xbO63YSjbXG0BESLsuoapCgHmZy1Er/hx/i5MqVphUBLhdD/u2ZuPq7G0VVf/dGf76cO8dm6Wto7RX5JUzMiS24Dx1XxOgLS522SAc9i73rLT0IhuqNT+fWyRB6e0gfKAhq5DdgxChbnS7W2ZRUxb/+aPy7erjduZHvbSfhwd3lcgMSTdOinhw1yCkpYeSbf6Xp/UWA7jvTVQuzDZ+ZLV8uCzhh6muqbRq5obWX95lEQ88hTOzA8zltkQ42jN74uvXdq3XSorfXtxZQVXaClrW/ZvemssBv8+jL5o2dxFaGjiti6Lgi1i2c16HZlFQjZP3RKz31zHLHvNUnocG9d7/+3PHs8wAdbmnKKSnpsoBuEEl3N+QYK4ey+rO2/3cA38RqLDhDTQ5B6Yatk/XDZlK1s5lv95xEqy0PvEEM3+atMqqiuDhx5HBavWeica51ud20tnuaY33shAb33L4F5oueyi9+ON3dX44BPXP/csDl5uVjLbEHdwievcfT690hjekGrZP1q+ZTNff/UtXYHxWFgHUQMQR1A6uMWrViCZXLPwv6bTuVCNv6DUHrjz/928KYW+hSos891Qilu/vLMQYHCi+hwT3EvBxrQdX6vP7Zu4ODjTRrnTQy05bjjeyuWIeqDvT+xh7YFZeL8753bYdsQ4aOK6KuqhJN01Janol2aVAsH27hcIJ7GPx192ByzMCRZ7PtvOmwxbcMu2+MBVUro0susgX3Hjk5HX6smDFGwp1R8NQmDVonjUD27Ref+7UzGhMhutwwatIUcvsWdDqY+X/bTqUOmliWBsUjqBs4wT0EVt3d097O0lf+O0COAX0SdfXa+A0dnW62S2sJ63mvWw+vz/C1qD3woRPgU52oWycTV3QN1/Vi4HIpFE+5kOIb7orbeW21DklmB41Vcjnd3ByT9BLvY3SCewgKi89HUVyo3or14SCB3TAGO77c3q7WkYKq9XmT0vP+zdugeqd/1TPw6ZO6qZMT4FObqFon8RZdr4UZL8RVqrEGs4O7asMEdYmC5LyCIxTf+68MvSy4W2pnCNVBU7VyeZd60gRo6O3tPnNDY+DKIM7SSzic4B6CoeOKGDXpopD6t9UYzL+A2tGCqvG8yel599sWtH+jPmwSpamTQ5KJpnXSWPO3Y0lcumpCTpEaeAPZqMJ8co9XUtznIENzTwG7O/W8kfDvoLF60sQziw/XvgiAlB1aGhQvnOAeBv9tSVasxmD9crPhcIv9cidISs/7BfdAxVxdltGfFDxnYMV/wJU/dwJ8uhBN62Qnu2pCTpEa+GenJ9bBxx/o+2hdMa4r7ABWiebEkcPmccYri4+ooWsa0pu5d7X0Eg4nuIeh+IppbF76ScD1/j7t/gXUzhRUIUk974WlcP3v9KxPGpmYBrUrYM/XTgafbkQsukbfVRO19BJMcqhbD/Oe1BfsKErCpD5DorEuvrdm8bFm0bG2LxqaezILuk5wD8PQcUUUffdKqlevMK8bObEkok97PEhK9n76qN8VQs+24rACzSFJmIu6g+xxhYhdNRGlF/Q2xotuuJkeObmBwcywG0ADKYKcY11LyCxe09i87FO2fLksaEAGAj7QEtG+GE+c4B6B7//4X+jdrz/b169hbOmlQQO7fwG1MwVVg6Qs8jBHw88AAoQCmqovSunVv2ue06HrMfa4xtBVU386LzbpJdQ52QWr9GLFP4s3C55SmvbBplbu/XeBRFVVX3Jl0MXti/FESP+D70ImT54sy8rKEvZ8ieLq369gh0VzHzMwl6WPX9npx62vqebvz/6r7c313bv+FxffckenHzskZa/Dx4/7vkZLryWaq4cjzXQHQnTV1J/Ko+5UPr1c7Rxs7U3ViaGoQWJ6zJlq3Xr45i3C2fsmkohFUNCTGQga2N1J0tCFEOVSysmx3MfJ3ONAvAuqBklZ5HH6qPek1kAzV487xdXuQpCums2Ng1nWMAbNNjVqbALWCSu9hCKCvW8yMLL44ium2Xvi29t9RVDD3FBVbdelcpYeDCe4x4F4F1StJHyRRyhpximudi8KS6mf9jJV7/2Zyoad3kZY63p3/QPeJaB4GBTf/GDsvekxrtNLJEaQB7t9cDDNPdmF0Y7iBPcUJ+FDTYWlekeDKc0IGFYC9RV6cdXJ4NOewFY+I0M3vql5g3p+A8X5hxia0wzL1kDr5uhbJ+vWQ1MdKG69jpskvT0arIHe//p0xgnucaArCqoGSRlqOn1UD+RooHngrIlwcIse2J0MPu2ItpVRUVxcdI6LHsdrKMxp0oO6SQyGZFY5RnHBRfenhN6eaTjBPQ7Ec0I1GAlvi+zV3xvc0f8eMgnuv0fP2GtXeDP4Vr1Q5rxhU5pop0gDetNDWRlEY0j2zdv6+YHUs/b84c55kgSUZB9Ad8C/gBqvgqqB0RZpxcjeu4TTR/GdGop+ubBUl2IUIx+QUOH1NHFIOeprqlny8p9Y+uqfgwd2IXC5s5h41XTuePZ5rnnkn3zfBI2i64wXCfBcNzBaJ+dcbz8H6tZDxd8wJR7FnbJyTHfHydzjQFcWVA0Smr2PugzcPfSiqlB8Pe6FpVDyAyibA0hQ2x39PUXo1BRpKIz++IiGZNfAwPEwbrq3iNpuPBmU3OOcG0nCCe5pQkKHmvyLqp8+qb/JC0t17XTT247+nkJ0eoo0HFEZkgGHt+p/rLiydc8ih6TgBPc40JUFVSulM2dRu3G9rXOmyzbOWPvdrW1shaV6ILfq7ynW5pYpxGzg1ZnzJBpDMn+crD2pdCq4CyGmAy8CLuAVKeVv4nJUaUZXF1QNEjrUZO1397cfMPT3PV/7OiKa9unZnfNm7jK6RHqJlWi2QAEIl5O1J5kOB3chhAv4E3ANsA/YIIRYLKXcEq+DSxe6akI1GAkbagonzRi/v3+x3jFT8RaUv6HLNY480yXU11Tz7nNP+3xRgtBh6aUjGFugvnkL9pXZ5RqhwA2JX+3nYKczmXspsENKWQsghPg7MBPIuOA+dnAe63c32i53FcGGmqpWLk+sNGMeTKl+neZJySnE7oSx9jGY30nSRuMNmQ5SzkPGoRPGYUKI24DpUsqHvZfvAy6WUv6z3+0eBR71XjwP+Lbjh5uaiOxeuVkFQ8/V9Qsp2xvrt8m205H0kgHAkY4834Deuef0zHL3NS6fams7fKzldBQiaGzkZZM7tp8yDoFAIrcf02qa22iJ9TZR0OHXohsS9LXo4XbnDuidM04IYbYvS6Q83dZ+pOVM+9EzHk8Xmw4lBee88HGulDKmrLEzmXuwBtiATwop5UvASwBCiLJYnc26K85r4cN5LXw4r4UP57XwIYSI2U63M0NM+4BCy+XhQH0nHs/BwcHBIU50JrhvAMYKIc4WQmQDdwGL43NYDg4ODg6docOyjJTSI4T4Z+Az9FbI16SUVRHu9lJHn68b4rwWPpzXwofzWvhwXgsfMb8WCd3E5ODg4OCQGBzjMAcHB4duiBPcHRwcHLohCQnuQojpQohtQogdQognE/GcqYgQolAI8YUQoloIUSWE+EmyjynZCCFcQogKIcSHyT6WZCKE6CuEeE8IsdV7fnwn2ceULIQQj3nfH98KId4WQvRM9jElCiHEa0KIQ0KIby3X9RNCLBFCbPf+XRDNY3V5cLfYFFwPTADuFkJM6OrnTVE8wONSyiLgEuCfMvi1MPgJ0EXG9GnFi8CnUsrxwAVk6GsihBgG/G9gspTyPPRmjbuSe1QJ5XVgut91TwLLpJRjgWXeyxFJROZu2hRIKdsAw6Yg45BSHpBSbvT+3Iz+Bh6W3KNKHkKI4cANwCvJPpZkIoToA1wOvAogpWyTUh5P7lElFTfQSwjhBnLIoPkZKeWXwDG/q2cCb3h/fgO4OZrHSkRwHwbUWS7vI4MDmoEQYhRQAqxL7pEklReAJ9CXsWUyo4HDwByvRPWKECI30p26I1LK/cD/A/YCB4AmKeXnyT2qpDNYSnkA9AQRGBTNnRIR3KOyKcgkhBC9gfnAT6WUJ5J9PMlACDEDOCSlLE/2saQAbuBC4M9SyhKghSi/enc3vHryTOBsYCiQK4T4QXKPKj1JRHB3bAosCCGy0AP7XCnlgmQfTxKZCtwkhNiNLtVNE0L8LbmHlDT2AfuklMa3uPfQg30mcjWwS0p5WErZDiwALk3yMSWbg0KIswC8fx+K5k6JCO6OTYEXIYRA11WrpZR/SPbxJBMp5c+llMOllKPQz4nlUsqMzNCklA1AnRDiXO9VV5GB1tle9gKXCCFyvO+Xq8jQ4rKFxcD93p/vBxZFc6cuX7PXQZuC7spU4D6gUgixyXvdU1LKj5N4TA6pwY+Bud4EqBZ4MMnHkxSklOuEEO8BG9G7yyrIIBsCIcTbwJXAACHEPuBZ4DfAPCHEQ+gffrdH9ViO/YCDg4ND98OZUHVwcHDohjjB3cHBwaEb4gR3BwcHh26IE9wdHBwcuiFOcHdwcHDohjjB3cHBwaEb4gR3BwcHh27I/wfTcM0E6NxAwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots()\n",
    "N, n = x.shape\n",
    "\n",
    "plt.plot(x[-1, :], y[-1, :], 'ko')\n",
    "for i in range(n):\n",
    "    plt.plot( x[::3, i], \n",
    "              y[::3, i], '.')\n",
    "    \n",
    "plt.ylim(0, 10)\n",
    "plt.xlim(0, 10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x619540a90>]"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADopJREFUeJzt3X9o3Pd9x/HXS5bUhCRbhCxW4x/SxML+aBhJJbwrhhGybKRZSAbLaIqXraPCUFqW0MLo+kdK+9/+CWPELHhOWcJEu5JkJcuSbR51SA09JzpPduq4G0ZExCRgVVHiiHWRL/feH7qtRpF9J+m++ureej7g8J3u4/u+D+Onv3z0PdkRIQBALj1lDwAA6DziDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgod6yDrxz584YGRkp6/AA0JVqtdrPImKo1brS4j4yMqKpqamyDg8AXcn2bDvr2JYBgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB3YImqzCzp8/Lxqswtlj4IESrvOHcAv1GYXdPBoVUv1hvp7ezQ5UdHY8EDZY6GLceYObAHVmXkt1RtqhHS53lB1Zr7skdDliDuwBVRGB9Xf26Mdlvp6e1QZHSx7JHQ5tmWALWBseECTExVVZ+ZVGR1kSwYbRtyBLWJseICoo2PYlgGAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEioZdxtX2f7VdunbZ+1/a1V1nzB9pzt6eZtophxAQDtaOdH/n4o6c6IWLTdJ+mE7Zciorpi3T9ExFc6PyIAYK1axj0iQtJi82Ff8xZFDgUA2Ji29txt77A9LemipGMRcXKVZX9g+4ztZ2zv7eiUAIA1aSvuEfFRRNwmaY+k/bZvXbHknySNRMRvSPp3SU+t9jq2D9mesj01Nze3kbkBANewpqtlIuI9SS9LunvF1+cj4sPmw7+VNHaV338kIsYjYnxoaGgd4wIA2tHO1TJDtm9u3r9e0l2Sfrpiza4rHt4n6VwnhwQArE07V8vskvSU7R1a/sfg+xHxgu1vS5qKiOcl/Znt+yTVJb0r6QtFDQwAaM3LF8NsvvHx8Ziamirl2ADQrWzXImK81To+oQoACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEio6+Jem13Q4ePnVZtdKHsUANiyesseYC1qsws6eLSqpXpD/b09mpyoaGx4oOyxAGDL6aoz9+rMvJbqDTVCulxvqDozX/ZIALAldVXcK6OD6u/t0Q5Lfb09qowOlj0SAGxJLbdlbF8n6RVJn2iufyYivrlizSckPS1pTNK8pM9FxJudHnZseECTExVVZ+ZVGR1kSwYArqKdPfcPJd0ZEYu2+ySdsP1SRFSvWPNFSQsR8Wu2H5T0l5I+V8C8GhseIOoA0ELLbZlYtth82Ne8xYpl90t6qnn/GUm/bdsdmxIAsCZt7bnb3mF7WtJFScci4uSKJbslvSVJEVGX9L6kj22I2z5ke8r21Nzc3MYmBwBcVVtxj4iPIuI2SXsk7bd964olq52lrzy7V0QciYjxiBgfGhpa+7QAgLas6WqZiHhP0suS7l7x1AVJeyXJdq+kX5b0bgfmg/jgFoC1a+dqmSFJlyPiPdvXS7pLy98wvdLzkv5E0o8lPSDphxHxsTN3rB0f3AKwHu2cue+SdNz2GUmvaXnP/QXb37Z9X3PNk5IGbZ+X9FVJXy9m3O2HD24BWI+WZ+4RcUbS7at8/dEr7v+PpD/s7GiQfvHBrcv1Bh/cAtC2rvrZMtsRH9wCsB7EvQvwwS0Aa9VVP1sGANAe4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AsIlqsws6fPy8arMLhR6nt9BXBwD8v9rsgg4erWqp3lB/b48mJyoaGx4o5FicuQPAJqnOzGup3lAjpMv1hqoz84Udi7gDwCapjA6qv7dHOyz19faoMjpY2LHYlgGATTI2PKDJiYqqM/OqjA4WtiUjEXcA2FRjwwOFRv3/tNyWsb3X9nHb52yftf3wKmvusP2+7enm7dFixgUAtKOdM/e6pK9FxCnbN0mq2T4WEW+sWPejiLi38yMCANaq5Zl7RLwTEaea9z+QdE7S7qIHAwCs35qulrE9Iul2SSdXefoztk/bfsn2pzowGwBgndr+hqrtGyU9K+mRiLi04ulTkoYjYtH2PZJ+IOmWVV7jkKRDkrRv3751Dw0AuLa2ztxt92k57JMR8dzK5yPiUkQsNu+/KKnP9s5V1h2JiPGIGB8aGtrg6ACAq2nnahlLelLSuYh47CprPtlcJ9v7m69b3EevAADX1M62zAFJD0l63fZ082vfkLRPkiLiCUkPSPqS7bqkn0t6MCKigHkBAG1oGfeIOCHJLdY8LunxTg0FANgYfrYMACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASahl323ttH7d9zvZZ2w+vssa2/9r2edtnbH+6mHEBAO3obWNNXdLXIuKU7Zsk1Wwfi4g3rljzWUm3NG+/Kelvmr8CAErQ8sw9It6JiFPN+x9IOidp94pl90t6OpZVJd1se1fHpwUAtGVNe+62RyTdLunkiqd2S3rriscX9PF/AGT7kO0p21Nzc3NrmxQA0La24277RknPSnokIi6tfHqV3xIf+0LEkYgYj4jxoaGhtU0KAGhbW3G33aflsE9GxHOrLLkgae8Vj/dIenvj4wEA1qOdq2Us6UlJ5yLisasse17SHzevmqlIej8i3ungnACANWjnapkDkh6S9Lrt6ebXviFpnyRFxBOSXpR0j6Tzkv5b0p92flQAQLtaxj0iTmj1PfUr14SkL3dqKADAxvAJVQBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3bEm12QUdPn5etdmFskcBulI7P34A2FS12QUdPFrVUr2h/t4eTU5UNDY8UPZYQFfhzB1bTnVmXkv1hhohXa43VJ2ZL3skoOsQd2w5ldFB9ff2aIelvt4eVUYHyx4J6Dpsy2DLGRse0ORERdWZeVVGB9mSAdaBuGNLGhseIOrABrAtAwAJEXcASIi4A0BCxB0AEiLuAJAQcQeAhIg7ACRE3AEgIeIOAAkRdwBIiLgDQELEHQASIu4AkBBxB4CEiDsAJNQy7ra/Y/ui7Z9c5fk7bL9ve7p5e7TzYwIA1qKd/6zj7yQ9Lunpa6z5UUTc25GJAAAb1vLMPSJekfTuJswCAOiQTu25f8b2adsv2f5Uh14TALBOnfg/VE9JGo6IRdv3SPqBpFtWW2j7kKRDkrRv374OHBoAsJoNn7lHxKWIWGzef1FSn+2dV1l7JCLGI2J8aGhoo4cGAFzFhuNu+5O23by/v/ma8xt9XQDA+rXclrH9XUl3SNpp+4Kkb0rqk6SIeELSA5K+ZLsu6eeSHoyIKGxiAEBLLeMeEZ9v8fzjWr5UEgCwRfAJVQBIiLgDQELEHQASIu4AkBBxB4CEiDsAJETcASAh4g4ACRF3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBICHiDgAJEXcApanNLujw8fOqzS6UPUo6Lf8nJgAoQm12QQePVrVUb6i/t0eTExWNDQ+UPVYanLkDKEV1Zl5L9YYaIV2uN1SdmS97pFSIO4BSVEYH1d/box2W+np7VBkdLHukVNiWAVCKseEBTU5UVJ2ZV2V0kC2ZDiPuAEozNjxA1AvCtgwAJETcASAh4g4ACRF3AEiIuANAQsQdABJyRJRzYHtO0uw6f/tOST/r4DjdgPe8PfCet4eNvOfhiBhqtai0uG+E7amIGC97js3Ee94eeM/bw2a8Z7ZlACAh4g4ACXVr3I+UPUAJeM/bA+95eyj8PXflnjsA4Nq69cwdAHANXRd323fb/k/b521/vex5imb7O7Yv2v5J2bNsFtt7bR+3fc72WdsPlz1T0WxfZ/tV26eb7/lbZc+0GWzvsP0ftl8oe5bNYPtN26/bnrY9VeixumlbxvYOSf8l6XckXZD0mqTPR8QbpQ5WINu/JWlR0tMRcWvZ82wG27sk7YqIU7ZvklST9PvJ/5wt6YaIWLTdJ+mEpIcjolryaIWy/VVJ45J+KSLuLXueotl+U9J4RBR+XX+3nbnvl3Q+ImYiYknS9yTdX/JMhYqIVyS9W/Ycmyki3omIU837H0g6J2l3uVMVK5YtNh/2NW/dc+a1Drb3SPo9SUfLniWjbov7bklvXfH4gpL/pd/ubI9Iul3SyXInKV5zi2Ja0kVJxyIi+3v+K0l/LqlR9iCbKCT9m+2a7UNFHqjb4u5Vvpb67GY7s32jpGclPRIRl8qep2gR8VFE3CZpj6T9ttNuw9m+V9LFiKiVPcsmOxARn5b0WUlfbm67FqLb4n5B0t4rHu+R9HZJs6BAzX3nZyVNRsRzZc+zmSLiPUkvS7q75FGKdEDSfc096O9JutP235c7UvEi4u3mrxcl/aOWt5oL0W1xf03SLbZ/1Xa/pAclPV/yTOiw5jcXn5R0LiIeK3uezWB7yPbNzfvXS7pL0k/Lnao4EfEXEbEnIka0/Pf4hxHxRyWPVSjbNzQvEJDtGyT9rqTCroLrqrhHRF3SVyT9q5a/yfb9iDhb7lTFsv1dST+W9Ou2L9j+YtkzbYIDkh7S8tncdPN2T9lDFWyXpOO2z2j5JOZYRGyLywO3kV+RdML2aUmvSvrniPiXog7WVZdCAgDa01Vn7gCA9hB3AEiIuANAQsQdABIi7gCQEHEHgISIOwAkRNwBIKH/BVxg3Nn+vunIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v = np.sqrt(vx**2 + vy**2)\n",
    "v_avg = np.mean(v, 0)\n",
    "plt.plot(v_avg, '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAG4AAAAVCAYAAACnvtv5AAAABHNCSVQICAgIfAhkiAAAA3FJREFUaIHt2V2IVVUUB/DfGImT2UxmEpRRmYbgQ9EHFThNpQRFUJBQ0EPRi1AS+BJFUA9REhFGL+FDiQg9qGQUfUwUpBiF0MdDmEVqEE4q5ZiZfU8Pax/mzuWcOfeeuZczNfOHYXHXWv+11pl199l7r8sM/lfYjCOYW3ch0xxXYhQPtOJ8Ff7Buib9XXgRu/BzCrilhXgX4GUcwu84iA04u4OcKjnqxApsx7CodxhDuDXH97VkP7Ms6BBG0Nuk/1w06wT2aq1xi3E4+e7AenyQPn+FczrAqZKjTjwuajuKV/A0NmIPns3xvyb5PzZR0KVitW3Msd2IJejBoNYa927yW9ukfz7pX+oAp0qOurBa1PQe5uXYTy/g7cV3OK0o8PoU+OaSAgaVN+6S5HMAs5ps8/ALThq/j7bLqZKjLszCflHPuW1ynxDPeUtjsEasxN/4eBIFZrgpySGxihtxArtxBq6dBKdKjrpwPS7GWziG2/AIHsZ1JdzdSa7KFI2Nm4vLxbI82YFCL0vy6wL7N0kunQSnSo66cHWSh/Ep3hRvuA34CB8qXol7khzIFI2NO1+8Q4c7VGhfkscL7Jm+fxKcKjnqwsIk14iD30rxOl8u9ukBbC3gHsdvuDBTNDYuO30d62CxE6EnydEucqrk6Bayg0WPuFq9L/bgL3EnvscNil+bP2FB9qGxcaeSnNOhQrNve1+B/awmvyqcKjnqQrYg9uOLJtspseqI438eeo31aFzjjiTZqXvPviSL9pclSTbuT+1yquSoC1mtIwX2rLHN92eiT/3GejQOPclwtIUiBpVfBxYrP6r/avxRvV1OlRx1YQH+FI2bnWN/WzzL3Tm2Zcm2vSj4tuRwaUkRg/7bF/BNyXbfRMV3AVtS3qea9KvEdWZE/kHq/sR7qCjwPcnhwRzbHeKBN+Gd5Pdtg+65HE7zOOoZY+OofVobeZVxquTYnOz35ti6iYXiijKKneJ/thV/idW4uoD3avJZVBR4Nn7AJzm2J1PCor+DBTEXiZncMP4Qo5sXML+oiAqcdv0/E4PyOobQ88Xb4ICo9Ue8rnhI0CcOJTvKAj8qGnFFR8qceugX06G8ge5UxFrRjxVljnPEN/aNbldUE24Xl9nz6i6kBfSKn6q2tUoYEIPNqXAam85YJraoi+otYwYzmO74Fz8gLD7NLQlCAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left ( 1000, \\quad 6\\right )$$"
      ],
      "text/plain": [
       "(1000, 6)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v.shape"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}