{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Band Diagram\n",
    "\n",
    "In this tutorial, we will build off of the `holey-wvg-cavity.ipynb` tutorial by considering a smaller, more abstract calculation that we really should have done first. In particular, we compute the band diagram of the infinite periodic waveguide by itself with no defects. This is very similar to the types of calculations that MPB performs, but with a different method that has its own strengths and weaknesses. By analyzing what solutions can propagate in the periodic structure, one gains fundamental insight into the aperiodic structures from `holey-wvg-cavity.ipynb`.\n",
    "\n",
    "Let us briefly review the problem. In a periodic system of this sort, the eigen-solutions can be expressed in the form of Bloch modes: a periodic Bloch envelope multiplied by a planewave $exp[i(k⋅x−\\omega t)]$, where $k$ is the Bloch wavevector. We wish to find the bands $\\omega(k)$. In this case, there is only one direction of periodicity, so we only have one wavevector component $k_x$. Moreover, the solutions are periodic functions of this wavevector: for a unit-period structure, $k_x$ and $k_x+2\\pi$ are redundant. Also, $k_x$ and $−k_x$ are redundant by time-reversal symmetry, so we only need to look for solutions in the irreducible Brillouin zone from $k_x=0$ to $k_x=\\pi$.\n",
    "\n",
    "Solving for these eigenmodes is very similar to solving for the resonant modes of a cavity. We put in a pulse and analyze the response via Harminv except that our computational cell and boundary conditions are different. In particular, our computational cell is simply the unit cell of the periodicity. The $\\epsilon$ function then obeys periodic boundary conditions, but the fields obey Bloch-periodic boundary conditions: the fields at the right side are $exp(ik_x⋅1)$ times the fields at the left side. For each $k_x$, we will do a separate computation to get the frequencies at that $k_x$.\n",
    "\n",
    "We'll first load our modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using MPI version 3.1, 1 processes\n"
     ]
    }
   ],
   "source": [
    "import meep as mp\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from IPython.display import Video\n",
    "\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And define our computational cell as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some parameters to describe the geometry:\n",
    "eps = 13  # dielectric constant of waveguide\n",
    "w = 1.2  # width of waveguide\n",
    "r = 0.36  # radius of holes\n",
    "\n",
    "# The cell dimensions\n",
    "sy = 12  # size of cell in y direction (perpendicular to wvg.)\n",
    "dpml = 1  # PML thickness (y direction only!)\n",
    "cell = mp.Vector3(1, sy)\n",
    "\n",
    "b = mp.Block(size=mp.Vector3(1e20, w, 1e20), material=mp.Medium(epsilon=eps))\n",
    "c = mp.Cylinder(radius=r)\n",
    "\n",
    "geometry = [b, c]\n",
    "\n",
    "resolution = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that our cell is now size `1` in the x direction, and there is no need for any loops to duplicate the geometry. We just have a single air hole in the unit cell. The PML absorbing boundaries have something new:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pml_layers = [mp.PML(dpml, direction=mp.Y)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our structure is periodic, we don't want any absorbing layers in the $x$ direction: adding `direction=mp.Y` just specifies PML layers on the $y$ boundaries.\n",
    "\n",
    "As before, our source will be a Gaussian pulse from an $H_z$ point source:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcen = 0.25  # pulse center frequency\n",
    "df = 1.5  # pulse freq. width: large df = short impulse\n",
    "\n",
    "s = [\n",
    "    mp.Source(\n",
    "        src=mp.GaussianSource(fcen, fwidth=df),\n",
    "        component=mp.Hz,\n",
    "        center=mp.Vector3(0.1234, 0),\n",
    "    )\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we put our source at $(0.1234,0)$. The $x$ coordinate is random, to help ensure that the source will couple to an arbitrary mode, but the $y$ coordinate is 0. This means that we will only be looking at $H_z$-polarized odd-symmetry modes (recalling the pseudovector subtlety discussed above). As usual, we will exploit this via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sym = [mp.Mirror(direction=mp.Y, phase=-1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, regardless of the source, we don't have an $X$ symmetry plane because this symmetry is broken by our boundary condition for $0<k_x<\\pi$.\n",
    "\n",
    "Let's set up our simulation and visualize the domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,1.2,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.36, height 1e+20, axis (0, 0, 1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ4AAAIpCAYAAABT1YGfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAfvUlEQVR4nO3de3hcdb3v8fc36S1JKU0vJtCGltbH3pSLXBUELLC5bC/gBvXoOeDt7O3d5xH0PB5v6Gbr9rLde3vER0VFVPQonI3CRqtCRVrcXIWKRW4tLUlpAzRpS5M2bZLv+WOtCZOQaZPOdL4zsz6v55lnJbPWrPklfXetmVkzK+buiJRbXfQAJJsUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeICZzTazr5rZo2a2y8y6zOxPZvaV6LHVKsv6+/HM7DjgN8BMYC3wF2AasBSY6+4TAodXszL9SzWz2cAKoAF4o7vfNGL+iSEDy4BMhwd8DpgFfGBkdADufk/5h5QNmd3VmlkD0EnyOHe2u+8KHlKmZHmLdzxwCLDa3XeZ2XnA2cAU4DHg5+7+dOQAa1mWw1uaTp8xs18Abxwx/wtm9m53/+lYVmZmawvMehmwC2g/sGEeNG1Ar7u3Rtx5lsNrTqdvAAaADwDXA43AB4HLgWvN7K/u/mAR91NHPYcwYyj0ytBF8lMHyXJ4udcwJwCfdPdv5s37mJnNAy4GPga8fX8rc/dlo11vZmuZwdLOJztpmthU7JiL1rGjg3N+fA4bv7ARno3bCmc5vJ15X18zyvxrSMI7vRR31jSxiaZJseG1b2/n/J+cz8btG0PHAdk+cpH77fe6+7OjzN+QTl9SnuEcXO3b2znj2jNY372e+dPnRw8n0+E9kE4bzGzyKPNnpNOdo8yrKvnRLWhewIq3r4geUnbDc/engDWAMfruNHfdA6PMqxojo7v90tuZO21u9LCyG17qy+n0q2Z2WO5KMzsGuCz99ltlH1WJjBZd26Ft0cMCsv3kAnf/iZn9DXAp8LCZ/ZHkuO2rgcnA1e5+feQYD1QlRwcZDy/1TuBO4B+AMwAH/gR8292vDRzXAav06EDh4cnB6qvTS9WrhuhAj/FqSrVEBwqvZlRTdKDwakK1RQcKr+pVY3Sg8KpatUYHCq9qVXN0oPCqUrVHBwqv6tRCdKDwqkqtRAcKr2rUUnSg8KpCrUUHCq/i1WJ0oPAqWq1GBwqvYtVydKDwKlKtRwcKr+JkITpQeBUlK9GBwqsYWYoOFF5FyFp0oPDCZTE6UHihshodKLwwWY4OFF6IrEcHCq/sFF1C4ZWRonuBwiuTjh0dii5P5k9hUS7n/PgcNm7fyPzp8/nV237FjIYZ9OzpCRlLz96Y+82n8Mokd/rXDds2sPiqxcGjiaddrYTQFq8MFs9azH2fuC96GMMc/7PjeeTZR8LuX+GVQZ3VhZ/xfaQ6i93ZaVcrIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFl8fMZprZM2bmZvZE9HhqmcIb7l+AWdGDyAKFlzKzM4FLgaujx5IFCg8wswbg28DDwFeDh5MJOjFj4rPAAuB0YG/wWDIh81s8MzsKuAy4xt1XRY8nKzK9xTOzOuC7wDbg40Wua22BWQuLWW+tynR4wIeAE4B3uvvW6MFkSWbDM7MjgCuBP7j7D4pdn7svK3A/a4Glxa6/1mT5Md5VwCTgvdEDyaLMbvGA15E8tvuWmeVfPyWdzjGz29Ov3+ruW8o4tpqX5fAAppO8hDKaKXnzphRYRg5QZne17m6jXYAj00XW5V2/IXCoNSmz4UkshSchFJ6EyPqTixdJH8/Z/paT4miLJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehMhseGbWaGYXmNn3zOxRM9ttZj1mtsbMPmNmU6PHWMsyGx7wNuBG4F3AAHATsAo4EvgccK+ZvSRueLUty+HtBb4DLHX3pe7+Znc/F1gEPAAsBv4tcoC1LLPhufu17v4P7v7XEddvBj6QfvsmM5tU/tHVvsyGtx9r0ulkYGbkQGqVwhvdgnS6F+iKHEitmhA9gAr1kXS6wt37xnIDM1tbYNbC0gyptmiLN4KZnQ+8m2Rr9+ng4dQsbfHymNli4MeAAR9z9zX7uckQd19WYJ1rgaWlGWHt0BYvZWZzgBVAM/A1d//34CHVNIUHmNkM4LfAPOAa4PLYEdW+zIeXHhr7Ncnu8D+A/+nuHjuq2pfp8MxsMvBL4ETgN8B/c/eB2FFlQ2bDM7N64KfAcpJjtG9y9z2xo8qOLD+r/SBwYfr1c8A3zWy05S539+fKNqqMyHJ4zXlfX1hwKbiCJEwpoczuat39Cne3MVw2RI+1FmU2PIml8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CRElv/ORVnpz6MNp/DKROENp/DKpK5Oj2ry6bchIbTFK4Pe3l5Wr14dPYxhent7Q+9f4ZVBV1cXV111Fe6OmQ17vJf/ff5fjxy57Mi/LFnsurq6ukr8U46PwiuDnp4ebr311uhhDNPT0xN6/wqvDAYGBnjuOf3l0Xx6ciEhFJ6E0K62gpgZ9fX1AAwODjI4ODim29XV1Q29TjgwMFAVL1YrvAri7vT39486b7RntTnjibRSKLwKMPJlkZxJkybR2tpKS0sL06dPp6GhAYBdu3axbds2Ojs72bJlC3v27BnzOiuFwqsA+YFMnDiRxsZGmpubaWtrY+HChSxYsIDW1lYOOeQQAJ5//nm2bNnC+vXreeKJJ+jo6KC7u5ve3l727t37onVWIoVXQebMmcNxxx3Hsccey9y5c5kxYwYzZ85k+vTpNDU1MWnSJAD27NlDT08P27ZtY+vWrWzdupWOjg4efPBB7r//fjZt2jS0zkrd8im8CjBlyhTa2to48cQTOffcczn11FM5/PDDh0KDZAuWf+Qh/zHfnj172LRpE6tXr6apqYn77ruPp556ir6+voqMDhReuOnTp3PUUUdx2mmnceqpp7JkyRJaW1uHRQcvji3fpEmTmDNnDqeffjqzZ89mwYIF3HHHHTz00ENs27atHD/GuCm8YIsWLeLcc8/l7LPPZunSpTQ2NtLf309/f//QSyS54HLT3FYsNx0cHKSuro4jjjiCWbNmMXPmTBobG+nr6+Oee+4J+Kn2T+EFMTPa2to46aSTWL58OUuWLBl61lpfXz90YD+37Mjb5qurqxu6rqGhgaVLlzIwMMCWLVvYvHkz7e3tZfiJxkdHLoK0tLRwyimn8JrXvIZFixbR1NTE4OAgAwMDmNlQTIV2r/DC7je37MDAAIODgzQ1NbFo0SJOO+00TjnlFFpbW8v4k42NtnhlMnI3uXDhQs477zxOOOGEoS1dse9Szr99Y2MjJ5xwAj09PWzcuJEtW7aMOo4omd/imVmDmX3ezB4zs91m9rSZfd/M5hzM+21ra+Poo49m3rx51NfXD23p9rWF25/c7QcGBqivr2fevHkcc8wxtLW1lXDkpXHA4ZnZv5pZYykHU25mNgVYCXwamAr8EmgH3gk8YGYLSnVf+S+HTJs2jZaWFmbMmAHAhAml3/Hk1tnc3ExLSwvTpk170TgiFbPF+wjwkJmdVarBBPgUcDLwX8DL3P0t7n4ScBkwG/h+qe9w8uTJzJ8/nzlz5jBlypSh68e1pXOH3/0O/v7v4fWvT6a/+11y/Yh1TZkyhTlz5jBv3jwmT55csp+jaLn/AeO9AD8BBoEB4Bqg+UDXFXEBJgHbAAeOHWX+mnTecUXez9p0PQ747Nmz/cILL/TrrrvOu7u7PWdwcNDHpKvLffly9ySz4Zfly927unxgYGBo8e7ubv/Rj37kF1xwgc+ePdvzxwKsjfr9H/AWz93fBrwe2ARcCjxsZm8+0PUFOAU4FFjn7g+MMv+GdPr6Ut5pU1MTbW1ttLS0DL0Faszc4aKLYOXK0eevXAkXXUT+trOuro7W1lba2tpoamo64HGXWlFPLtz9FmAJcBXJrumnZvbLg/3AvESOTqd/KjA/d/1RpbzTKVOmMGPGDKZNmzb+Z7G33lo4upyVK+G224a+ra+vZ9q0acyYMaOidrVFP6p19x7gQ2Z2HfA9ki3E6Wb2baDgJ0rc/fPF3neRjkinHQXm566fN5aVmdnaArMW5n9TX1/PpEmTmDBhwvifwV5//diW+/nP4eyzh76dOHEikyZNqqgPlZfs6ZS732VmxwJ3ACcClxdY1EgeX0SHNzWdFvqAae4/zSGlvNOBgQH27NlDf3//+J9dbt48tuU6O4d9u3fvXvr6+irqzaIlCy996eFq4ASSJxw3so8tXq1x92WjXZ9uCZfmvt+9ezddXV1s3759/CEcdtjYlmtpGfpyYGCA7du309XVRV9f3/ju7yAqOjxL9heXAVcAjSTPBt/j7vcXu+6DbGc6LfRaZO6R+POlvNOenh6eeuopnnnmGQYGBsZ344svhquv3v9yb37hOd7g4CCdnZ20t7ezc+fOfdyovIra6ZvZK4C7gS8B9cAngeOrIDqAp9Lp3ALzc9dvLOWdbt++nXXr1tHR0TH0bmEY4yGss86C5cv3vcyZZ+J5y+Teq7du3Tp27NhxoMMuuWKOXFwJ3AccD6wGjnb3L7r7OP8bh1mTTl9ZYH7u+j+X8k737NnDxo0befrpp9m9e/fQ9WMKzwxuuKFwfGeeCddfT/6adu/ezaZNm9i4ceOon82IUsyu9n8DO4APu/u3SzSecroT2A4sNLNj3P3BEfMvSqc3l+LO8g/O79ixg87OTrq7u2lra6O/v3/sz3Cbm5OXVW67LXn2umULtLYmu9czz0ziHBigv7+fCRMm0NXVRWdnJ88///yLxhGpmPBuBt7n7k+XajDl5O57zOwbJA8PrjKzv0lfGsLMPkry+t0fDtbDhvb2dh588EEOPfRQWltbmTx58lAM+43QLNntnjX8aKWnxzDq6+vp6+sb+hxGJb4f74DDc/c3lnIgQa4EzgJeDTxuZqtIXrc7CXgWeFep7mjkFmbdunX8+te/prGxkbPOOovJkycPPcsd9xGNVP7te3t7uffee1mxYgVPPvlkwXFEqZxXFAO4+27gtcA/kryedwFJeD8AXunu6w/WfW/ZsoU777yTO+64g0cffZSenh7q6uqG3n08ODiYf7x3tLEPXXLL1tfXU1dXR09PD4888girVq3izjvvHHovXiXJ/BtB3X0X8Jn0Ulbt7e3cfffdtLa2UldXx7Jly2hsbBx6mSX3zuKRu+D8GHPhQfJWqF27dvHwww9z2223cdddd9HRUejATKzMhxftscceY8WKFfT29rJ161aWLFnCYYcd9qJPmeXLfwyYe+t77tnyww8/zKpVq1i1ahWPPfZYOX6EA6Lwgm3bto177rmHzZs3s379es477zxOPfVU5syZMxTfaLvc/Hcr9/X1DX2udsWKFdx77720t7dX1JGKkRReBejr6+OJJ55g165d7Ny5k8cff3zoTAKzZs1i+vTpTJ06dSjEvr4+enp66O7uHjqTwKZNm3jggQe4//772Zx3TLdSzyRglTioWjLyWO3+TJw4kYaGBpqbm5k7dy4vfelLh507xd3ZuXMnmzdv5sknnxx27pRdu3YNOxoyBg8XOsZ8sCm8g2ws4e3rbFEtLS1DZ4vKvVV+9+7dQ2eL6uzsLOZsUQqvVo13i1dmYeHpMV4F0RlBJYR74TOC7ks1nhE000cuJI7CkxDa1ZZBfX09zc3N0cMYpru7e/zvgC4hhVcGTU1NnJ33qa9KcMstt4S+I1nhlcHMmTP50Ic+FD2MYe666y6FV+saGhp41ateFT2MYXKnRouiJxcSQlu8Mqm219kONoVXJsWccLEWKbwyUXjD6TGehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEyG56ZLTaz/2Vmvzez58xsr5ltMbP/MLPXRI+v1mX5L/vcCswBdgJ3AV3AUuBC4AIz+6i7/1vg+GpaZrd4wCPAJcBsdz/b3d/i7q8A3gsY8FUzWxo6whqW2fDc/Sx3/5G77x5x/beB3wL1wMUhg8uAzIa3H2vS6eGho6hhCm90C9LpltBR1LAsP7kYlZktBF6XfnvTOG63tsCshUUPqgZpi5fHzCYAPwAmAz9z9/tjR1S7qnaLZ2Y3AkvGebNL3P2efcz/OnAqsB54/3hW7O7LRrs+3RLq2fEIVRsecCSwaJy3aSw0w8w+CbwP6ATOcfeuIsYm+1G14bn7MaVal5m9F7gS2A6c6+5PlGrdMrrMP8Yzs7cCVwG9wN+6+4PBQ8qETIdnZucDPwT6gQvd/c7gIWVG1e5qi2VmpwA3kBwee7O7/zZ4SJmS2fCA/wQagCdJ3hRwwSjLrHb375Z3WNmQ5fCmp9Mj00shCu8gyGx47m7RY8iyTD+5kDgKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHw8pjZp83M08t/jx5PLVN4KTNbBHwS8OixZIHCA8zMgO8A24CbgoeTCQov8R7gNOAykvjkIMt8eGbWCnwZuM3dr4seT1ZkPjzg60AD8L7ogWTJhOgBRDKz1wEXA59198eLXNfaArMWFrPeWpXZLZ6ZTQW+CTwGfCl4OJlTtVs8M7sRWDLOm13i7vekX38BaAPOdPe+Ysfj7stGuz7dEi4tdv21pmrDA44EFo3zNo0AZnYi8AHgR+6+stQDk/2r2vDc/Zgibn4+ycOMV5jZ7SPmLU6nnzSz9wAr3P2fi7gvGUXVhlci+4p3cXrZUJ6hZEsmn1y4+xXubqNdgGvTxf5Het07AodaszIZnsRTeBJC4UmIrD+5eJH0Md07godR87TFkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwJITCkxA6k0AZDPogPXt6oocxzKAPht6/wiuDR557hKlfnBo9jOGei7177WolhLZ4ZdJ5eSdNE5tKsq6OHR2ce925bNi2gfnT57Pi7SuYO23umG/fs7eHlqtaSjKWA6XwyqRpYhNNk4oPr317O+f/5Hw2bNvAguYF3H7p7bQd2laCEZaXdrVVpH17O2dcewbru9dXdXSg8KpGLUUHCq8q1Fp0oPAqXi1GBwqvotVqdKDwKlYtRwcKryLVenSg8CpOFqIDhVdRshIdKLyKkaXoQOFVhKxFBwovXBajA4UXKqvRgcILk+XoQOGFyHp0oPDKTtElFF4ZKboXKLwy6djRoejy6K3vZZL7jISiS2iLVyaKbjiFVybzp89XdHkUXpmsePsKRZdH4ZXJeD73mgUKT0IoPAmh8CSEwpMQCk9CKDzAzC4wsxVm9qyZ7TazdjO70cxOjR5brcr0ITMzqwOuBt4F9ACrgW3AEcD5wP3pdVJimQ4P+AxJdDcD73D3rtwMM2sGZkUNrNZlNjwzmwt8AngKeIu778qf7+7dQHfE2LIgy4/xLgUmAd8dGV2t69jRET2E7G7xgOXp9I9mdhjwduClwHbg98Bv3N1LdWc9eyvjrO8dOzo458fnRA8DK+HvtqqY2WagFfgw8I/AoSMWuR240N23jXF9awvMWkw9dcw40JEeJF3AAM+7+7SIu89yeLuByUA/8F/AR4AngBNJnukeCdzg7hePcX2FwlsKDAKPFDtmYGE6XVeCdS0GBt19YgnWNW5VG56Z3QgsGefNLnH3e9Lb7wEmAs8CR7r70L7QzF4O/BkwYJG7P1bEONcCuPuyA11Hpa/rQFTzY7wjgUXjvE1j3tc7gWbg+vzoANz9L2Z2L8nW7zTggMOT0VVteO5+TJGr2EgS3oYC8zeQhPeSIu9HRpHll1MeSKfNBebnng7sLMNYMifL4d2UTk8fOcPMpgKvTL99YOR8KV6Ww7sZ+CvwajN7f+5KM6sHvkayxfsLOlZ7UFTts9pSMLNjgD8A04A1JC+nHAssALYCr3X3h+JGWLuyvMXD3R8EjgF+CLQAbyA9jAYcp+gOnkxv8SROprd4EkfhSQiFJyEUnoRQeBJC4UkIhVdiZnaKmf3KzLrMbKeZ3WNml4yyXIOZfd7MHks/Uvm0mX3fzOak899hZr6Py//NW9dhZna3me1N5/Wb2VozG9fbxsxsw37uc3Hxv6FE1b47pRKZ2d8BPyP5D30H8BxwJnCtmR3l7peny00BVgInA5uBXwLzgXcCrzOzk/NWuwZ4cJS7uztd1xySIy5TgAGgneTTcUuBh8zsBHcf7/Hmawtcv32c6ynM3XUpwYXk2O52wIE35V3fAjyeXn9Get2V6fd/BKbmLfvR9PrbgXekX1+xn/t9KF1uK3Bo3vV3pNdvGcfPsCFJogy/r+h/sFq5AB9P/6F/Mcq8C9N5N5McktuWfn/sKMuuSed9dn/hkXzw3NPLmSPmTQb2pvPOGuPPULbw9BivdP42nd4wyrxbgN3AWcAZJB8sWuej7wJztz92DPf5gXS6y91vy5/h7n0kb98H+OAY1lVWeoxXOken0z+NnOHue8zsL8DxJI/5Rl1uxPW5U4geZ2ZfIXkHzRZgpbv/IZ2Xeyz4VIF13U3yvsJxfa7CzD5G8sGiPmAtcKO7PzuedexX9C6qFi4kUeR2edMKLHNjOv8X6fRrBZY7Op2/IW+dIy+3M/yx48oC6/pIOr97jD9HofvsAd5Vyt+ZdrWlMTXv694Cy+Q+UHTIGJebAFxBsss9lOQzwG8g+Zjk6cB/Ag37WVfuFByTCswf6SbgTcA8kg9GvZzkTbGTge+a2RvHuJ790q42VezHJQ+CXnf/XN73O4Cbzez3JGexOp7kY9kl4+4fHnHVWuAyM3sE+A7wJZKXfoqm8F5QzMcld464bscoyzal0+dH3HZ/yw3j7jvN7OvAN0g+F7yvdeU+yLSnwPyx+h7JS0CLzGy+u28ocn3a1ea4+zHubuO83J7edgcvvLha6O8K5K5/dIzLbdzHcB9Pp/3p9PACy+X+Iz23j3Xtl7sP8sLZCw4rZl05Cq901qTTV46cYWYTSR4v7QZuK7TciOv/XGA+vLAlyz3TPKLAciel00Kn1xiP3H2W5OxDCq90bkmnF40y73Ukh7RuJXlGuh1YmH7YaKTc7W/ex339XTrNPd5qMLPX5i9gZpOBo9Jvv7HPke+HmS0j2Xr2UppzwOjllFJdKHzI7CUUPmR2J8mu9xFgDsMPmX0CmJUu/8F0mS/xwhGN3vQ2uUNmz5H3Ug4FDpnlreuLI64/H1g+ys91FPBwuq5/L9nvK/ofrJYuJFuiAZKzQ60Erid5ScOBf8lbbgpwF8NfK3sgnT5D8vFKJ9k1r86Lqyed7srFnca3O72+n+TF5N68748dMcYr0nk/KHD9BpIt6U9JXoDOHXb7PdCg8Cr0ApwC/DoNrge4F7h0lOUagM/nhfcMcA0wN53/OeC3JE8ycv/4O4BvkZzBKn9dh4+IZIDkcd3SUe63UHivInn2+ud067mX5I0HvwfeA9SX8vekjzdKCD25kBAKT0IoPAmh8CSEwpMQCk9CKDwJofAkhMKTEApPQig8CaHwKpSZXZKer+Sh9I2koy1zspkNmNlzZja73GMshsKrUO7+Q5I3jr6c5CwFw6QxXk3yb3iZl/pzrweZ3p1SwcxsAcnf2jDgKHd/PG/ep0j+3Omt7n520BAPmMKrcOmn+r8M3O7ur02vW0TyGY9B4BXuXoo/I1pW2tVWvn8leXfyGWb2bjMzks+4TiY5oU/VRQfa4lUFMzuO5B3GO0hC/DzJOfNOcPf+fd22Uim8KmFmXwUuS78dAE529/sCh1QUhVclzOxwoIPkicb33f3dwUMqih7jVY/PkUQHcI6ZHbKvhSudwqsCZnYa8G6S8yX/guQjjf8UOqgiaVdb4dIzAqwh+ST/RSSfs/0ryanLTnb3ewOHd8C0xat8nyKJ7iZ3/3/u3klyJKMOuNrMqvKMX9riVTAzeznJqWl3k3w4uyO93kj+wPNrgI+7+1fiRnlgFF6FMrM6knOrnAx82N3/z4j5S0hey+sHlnkJzllXTtrVVq73k0R3N3DVyJnu/lfgn0lOyvjN8g6teNriVSAzm0tyhqYG4JVe4E/Up088/gy8DHiru/+sfKMsjsKTENrVSgiFJyEUnoRQeBJC4UkIhSchFJ6EUHgSQuFJCIUnIRSehFB4EkLhSQiFJyEUnoRQeBJC4UmI/w/RMSUYnRcg7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim = mp.Simulation(\n",
    "    cell_size=cell,\n",
    "    geometry=geometry,\n",
    "    boundary_layers=pml_layers,\n",
    "    sources=s,\n",
    "    symmetries=sym,\n",
    "    resolution=resolution,\n",
    ")\n",
    "f = plt.figure(dpi=150)\n",
    "sim.plot2D(ax=f.gca())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, there are two ways to proceed. First, we could set the value of $k$ via the `k_point` variable, and then use `until_after_sources` with Harminv just as we did to calculate a resonant mode:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error\n",
      "harminv0:, 0.3050132495988499, 3.480691775115153e-06, -43815.03294539187, 0.011878313829683136, -0.01179628264240186-0.0013935764266843182i, 5.295036558380241e-07+0.0i\n",
      "harminv0:, 0.43960436531181446, -2.7363657356159734e-06, 80326.31741985629, 0.01745119452862526, 0.017126214985327293-0.0033521561348289074i, 9.312403773661967e-08+0.0i\n",
      "harminv0:, 0.5005873843994147, -0.0009565526575684324, 261.6622202858719, 0.00018410949058675333, -0.00011025456283318074-0.00014744570491736274i, 1.242711680938466e-05+0.0i\n",
      "harminv0:, 0.6404729374549927, -0.005167317620050075, 61.97344391699944, 0.0031051734641799815, 0.0007710721023072853-0.0030079145692141806i, 3.514904299909929e-05+0.0i\n",
      "harminv0:, 0.700226243241067, -0.0037510897788062153, 93.3363748313049, 0.004369347983945716, -5.4869994286284585e-05+0.0043690034434110386i, 1.3205854700252928e-05+0.0i\n",
      "harminv0:, 0.762852889565266, -0.006722784235381431, 56.73638055721298, 0.0316955396820757, 0.02965548904826291-0.011187457488026218i, 4.404983351452081e-06+0.0i\n",
      "harminv0:, 0.8234265979478261, -0.006957003127114023, 59.179691521096714, 0.000615608025293912, 0.0006152489526665892-2.102301236585398e-05i, 0.000399758756541917+0.0i\n",
      "harminv0:, 0.8243479880387249, -0.0005179204557115729, 795.8248983486759, 0.030110588821435868, -0.005890629837883119-0.029528766301466433i, 5.666365518267213e-08+0.0i\n",
      "harminv0:, 0.8878738611231001, -0.0002512055373818404, 1767.2258947331713, 0.0032343510475491863, -0.003225648132430831-0.00023710930923945473i, 8.018541351827511e-07+0.0i\n",
      "harminv0:, 0.9495092084271758, -0.005235806245342398, 90.67459374301981, 0.005382664822832947, -0.0016788910104550295+0.005114137812962859i, 1.7544167778335838e-05+0.0i\n",
      "harminv0:, 0.9758414410168742, 0.0009253350005759116, -527.2908948702515, 9.453933579227094e-06, 3.226576270516566e-06+8.886285258249198e-06i, 0.00012202208460776741+0.0i\n",
      "harminv0:, 0.9858676174323542, -0.0037464233316342127, 131.5745085596497, 0.0012188508341391812, -0.0004042226659787743+0.0011498701631883445i, 2.708082943685142e-05+0.0i\n",
      "run 0 finished at t = 306.675 (12267 timesteps)\n"
     ]
    }
   ],
   "source": [
    "kx = 0.4\n",
    "sim.k_point = mp.Vector3(kx)\n",
    "\n",
    "sim.run(\n",
    "    mp.after_sources(mp.Harminv(mp.Hz, mp.Vector3(0.1234), fcen, df)),\n",
    "    until_after_sources=300,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which gives us the frequencies at a single $k=0.4⋅2\\pi x̂$. Note that, in Meep, $k$ is specified as a vector in Cartesian coordinates, with units of $2\\pi$/distance. This is different from MPB, which uses the basis of the _reciprocal lattice_ vectors. However, this only gives us one $k$. Instead, there is a built-in function which takes as input a time to run after the sources finish, like the 300 above, and a list of $k$ points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "sim.restart_fields()\n",
    "k_interp = 19\n",
    "kpts = mp.interpolate(k_interp, [mp.Vector3(0), mp.Vector3(0.5)])\n",
    "all_freqs = sim.run_k_points(300, kpts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have used Meep's built-in `interpolate` function to interpolate a set of 19 $k$ points between $k=0$ and $k=0.5⋅2\\pi x̂$ , to cover the irreducible Brillouin zone. This function automatically runs Harminv, using the frequency range and location taken from the Gaussian source in the sources list. It returns the calculated modes as a list of lists.\n",
    "\n",
    "Plotting the real parts of $\\omega$, where the light cone $\\omega > ck$ is shaded gray, we find:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHqCAYAAADLbQ06AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde5hc1Xnn++/qVl/U3aW+XwS6tozUkuDgBAMRR1ack/HEhtjBBI9t4gnYkTKOTQLxYBtChiQiHOyEwcHBGT+R7Ih4jOPBI3ROxthO8hx7wLEGJszEAVoy5iLAAkmouyX1/VK1zh+7t1Rd6uquXbVW995dv8/z8BRVvWtp122/e639rncZay0iIiISTxWLvQMiIiKSnwK1iIhIjClQi4iIxJgCtYiISIwpUIuIiMSYArWIiEiMKVCLiIjEmAK1iIhIjClQi4iIxJgCtYiISIzFIlAbY3YYY/7WGPO6McYaY64t4DnvMMb8L2PMuDHmBWPMTQuwqyIiIgsqFoEaqAd+BHyikI2NMeuBbwHfA94K/Bmw1xjzS972UEREZBGYuC3KYYyxwPustQfm2OZzwDXW2ouzHvsboMla+64F2E0REZEFEZcedVTbgH/Ieey704+LiIgsGcsWeweK1AUcz3nsOLDCGLPcWjua+wRjTA1Qk/NwC9DvZxdFRKQMpYDXrcPh6qQG6mLcAfzBYu+EiIgseauAo64aS2qgPgZ05jzWCZyZrTc97V7g/qz7KeCnV175/zA42HX2wY99DK68MvoOZTLw/PNw+jQ0NsLGjVAR8cLC4cPwuc/Nv91nPgM9PQu/f6Gnn4aHH4bh4TSf+tQR/vRP11FfX8kNN8Bll0Xfr099CvrzjGsYA83N8Kd/Wvz+lsLHZ/Lkk/ClL527X1Nz7n0cH688+3ih38Xc9vJZrPZ8vIe5beZ7D6O0+fTT8MUvQm4/yJjg9hOfKPz7vRDfm3yK/VxKfQ99/JZdfia57T788Mx9bWmhqGNYtlOnTvFLv/RLAIPFt3K+pAbqg8DVOY+9c/rxWVlrx4Hx8L6Z/qQzmQamphrPbrdyZRDIIu3MQdizB06ePPdYWxvs2gXbIlw1v/xyaGqa2U6utrZguyhBy9X+hW098EDww6mpSVNXV8fUVCNvvlnJAw/A7bdHa/OZZ+DEibm3OXECfvpTuOSSaPvqwvg4TE0Vtl2h35uVK2e2WVl57n2cmqqcsV0hbea2N9d2i9Fe+L3u6zv/gAvBQbe1Ndr3OrfN3PcwapuZDPzn/wyTk7P/3Zjg7+94R2Ht+XjNrj+X3O92vu9hod9t179l159JKPsYlu3NNynqGJYtnU4X98R5xCKZzBjTYIx5qzHmrdMPrZ++v2b67/caY/466ylfArqNMX9ijOkxxnwc+DfA50vZj1QKtm6N9pyDB+Gznz0/uPb1BY8fzHvqcL6KCtixY+5tduyI/qV0tX+ZTBDwZzvwhI/t2RNsV6iBAbfbudbc7HY7CL5jbW3negW5jAn+Xuh3Me7tVVQEJ4Xhc3PbguDvUb7Xrtt87rm5T5CtDf7+3HOLs3/g/nNx/d12/Vt2/ZmAn2NYdts//nH05xUiFoEaeBvwv6f/g2CI+n8Du6fvrwTWhBtba18GriHoRf8I+PfATmvtdxdqh8H9h57JwOOPz73N449Ha8/l/vn44fgIhC65PjiC+4N43NuDoIdy++1BLzJba2vxPRiXbfo4YXT9ml1/LuUW+MHPMQyCDs/OnXD//fNvW4xYDH1ba78P5Pm6gLX2pjzP+RmX+zE4GHxAhQ6xRvnQC2lzvvbAbXtR98/HDyc8WMw3RBh1pMOV8OD42c8G+5K9j8UGLTh3EN+zJ/jehVpbi7skkd1e9mcel/bCNq+8Mvi+DQwEB+ytW0vLPchuE+Cuu4pr09cJo+vX7PJzyf1uZysl8Lv6Lfv4THwcw8JRy+ByYOHPiyIWgTpOonxArj/0uG/n44fjKxC65CNohe26CDKztecqKLgOrBUV7nMNKipgyxb4yU+C22L2z+cJo+vX7PJzcXnC6Pq37OMzcX0Mm2vU0iUF6hxRgozrDz3u2/k6mPkKhC75CFrgJsjktucyKPgIrHGUhBPGbC4/F5cnjD57/C4+k/AYNl/CbqHHsEJGQV1QoJ5WTJDZujVIQBucIxE/SoKa60Douj3XQ2XZfAVCl8olaJWrJJww+uLyhNF1j/9974MDB84P1NdeW9x1/h07YP/+/NtESdhdqCRXBWpKCzL5pg4U+vdsrs8gfZyRur62mru/CoSymJJwwuhDJgO9vVBVFdzG5TUfPAiPPnp+RyOTCR7ftCnaMafQhN1f//XCXv9CJbkqUFN8kHnmGRgbm3ubsbFgu0svLazNpCUFQenXVkXipNxOGMM6C4ODcPfdsHt3MBJY7PHBVd2GQq7/7tkTHItcTcODaAm2hQylu1D2gXrXLti+vbgg8+yzhW9XaKAO5X45S0lW8JUU5PLaqogsvHwZy2GdhajTyLLby1ZMe65nrYD7BNtChtJdKPvD64YN8Qoy4Re9r2/m4/390QuUZAt7CTt2BLdxes0isvB81IFw2Z6PqVQ+sr7nG0p3QYfrElx88fzbRNnOZ9UcEZFsrot/uG4vCVUBFyrrW4G6BJdcElzLmUsq5aeAiohIKeJetyEJVQEXKutbgboEFRVw881zb3Pzze4/9MWqey0i0WUyQULp448Ht3EZEYt73QYfpWzBbWlXZX0nxLZtcMcd8Jd/OfO6cjFZjnGvey0i0bhcuc61uNdtAL9VAS+/HB57DI4dg64uuPpqWBYxIirrO0FcZVVv3hw8Z64z7oqKYDsRiTeXGdC5MpnSjzeuCxiF7d177+x/t7b4gkguZ8HA7CdQBw4UVz1tIbK+FagdcTH38tCh+YfFMplgu3Ka5ymSNPMlhhoTfQ5wyGUvPbvyV7ZiK3+5lu9kJ5wFU8zJjssTKGV9L5AXX3RzzcjFdaj+frfbiUh0Ln7LPpdTdLW+fNjeo4+e/xrDyl/FrFefT3hysljTvXy0qVrfC2TPHti3r7TrHa7OcE+fdrudiETj6rfsIzHUdS/ddeWvhV72t5iCJ4u19G+pyr5HDcWfjYLbM9zGRrfbiZSDsE41BLfFjpDF/bcc93nPcZ/u5aPNhUrsVaAmPsMoLS1ut8sW1ykiUl5cfw8PHoSdO4P61BDc7twZ/aQ7CcWG4h4I4z7dy0eb8831dqXsh75DcRhG8ZX17WOKSFxX2xG3XGQXh1x/D13WqXb9W/ZxGSvugTAJ070Wculfl3RozbGYwyhRsr4L5Tr5JGzTRS9G4i38nO+8E+67L7gt9nN2/T2Me13pJJS/dN2e6wIlPgqe+GgzXwEVlxSocyzmMIrrg8V8BzNrow/n+Qj82fur4fnSuHoPXX7OPoaV415XOgnlLxcyaBVT9ctHez7b3LsXPvnJ6M8thIa+p8VhGMX1wcL12qtJmRuavb8ul/Z03V7YpqtLCAuxDnAxn3MSliv0PSSa3aaL8pc+1qsfHCy9vbBNl8vq+lim19fSv5s2Ff/8uShQU/owiqsfo+uDhet52T4OuOCngpOP66GuTyTCNgcH4e67g0sIqVRxbcZ5HeAkLFfoI7D6LH/pKxAC3HWXm6DlsiiT6/Z8temLhr6JzzCK66Eo1wktizE3FBZ/eN7Xdf64Di2X4/Va8Dskes89cNttwe3evaVX/HK9vnxFBWzZEvz/li1KDI2bsu9R79oF27fHZxjF5Vm467mcPg64rntvC1kUotjh/rgPLcc9Gxjc16kO+RoSTUrPTeKp7AP1hg1uzh5d/hhdHSyamtxu5+OA67r3pupIpW+XxOu1rq6vhvurwCpxUvaBOq5cHCyi1NQthI9eTNwz55dSdaSlfr0W3FxfFYkbBeolLCyrWMh2P/MzhW3ruhcT98z5JFVHivs6wD6GleHc9dWf/ETXV2VpUqCWorhaHzbumfNJrI7kcmhZ12tFFp/OPZewiy92ux2cy1bu65v5eLg+bDEZ0HHOnE96daRSMpbDfXWZXSwi0alHvYRdckkwLzd7iDpXKlV4D8dnwZO4Zs77aC+3zbgWmhCReNBPeAmrqICbb557m5tvLvxA7rpso2+uhud9teejTfWARZYe/YyXuG3b4I47zh8SbWsLHo/Sc/ORAR3ysQCEq+F5H8P9PtoE1UsXWYo09F0GXA2J+siABrflL8ux4EnIR5lTEVl86lGXCRdDouF62fP9O1HWy3Zd/tL18LyP4X4fbfpa1Uw9dJHFpx51mXCx8tNzzxW2XvZzz8GllxbeZpyrdCWh4Il66CJLW9kH6hdfhJUr3SxXGNclFV0dcJ99tvDtCg3Uca/SlYSCJz7KnPpY0Qz8LBUqstSVfaDeswf27XOzXGEcl1T0dcB1Je5VupJQ8EQ9dJGlTeeyxGe5Qtftub7+Gy6D52o7CIJRKjX3NqlU9CpdUD4FTxazh14oX9fQIfj+huVye3t1HV2WHgVqil/32HUgjHtiVZTgtphcV+nytU6xqzZdr8280D10iP7bC4XT+nbvDu7v3l38tD6RuCr7oe9QHJYrjHti1enTbreD4LXMVTkNgr9H+VzAfZUuH1W/XK385LrWdxKuocPMyzo1Necej8tlHRFXFKhzLGb2bjkmVvksouJ6AQgfC0q4WvnJZZnTuF9DB7/lbJXwJnGjQJ1jMbN3yzGxylcRlXLkqtcf9x46+O2lK+FN4kbnidOiXscD99cGXbeXhMQqH0VUypmrWt9xvoYOfnrpPhPeREqhQE18sneTsPyh6/YOHSqsiMqhQ9HaldJt2wZ798I998BttwW3e/dG/4x9fK9d99J9JrypupuUSkPfuFuuMM5LKvpKrILik6AgWITC5Xbilqvr8q6/164vw2goXeKs7AP1rl2wfbu77N24Zhj7SKxykQTlI5Nc4snl9zr3Onq2YnrpPofSfRQbUsJbeSn7QB1XPjKM46ix0e12Em8uv9fZvfTsKX7F9NIXeii9lKx09dLLT9kH6jiWEIV41w53qaXF7XZSXlxdhknSUHqcSwKLH2UfqKG0L7mPH06ca4eHwrKNVVXBbbGBPzxAznVQi5oRLOXFxWUY11PSkjh33MXvWfzQR0F8SohCvGuHZ7fpqmxjeICca+pO1IxgkWK4nNGw2HPHo1AZ1vjT4W9aMV9y1z+cuNcOBz+Bf9s2eN/7zg/GFRXB46X0+l1Oi9E0m6XP1ZQ0zR0XlzT0nWMxS4jGvXa4z+UPH330/HYzmeDxTZuW5uUDiScXCW+uh9IhWQlvcc2JSSq9dTkWs4Ro3GuH+xh6m+tgEVrsXr96HVIM18WBXPfSfQ+l33kn3HdfcKuh9NKoRz2tmDrVrjNF41473MfQ20L3+iFaL8Fnr0OWPp9zx+OY8KasdD90aCE+JUTjXjs8CatnzRf4IVovwVevQ8qHq/rrEO+EN59lWMudetTEp4So6zNm1+0lYfUs1yVJfS7DKVIMV730pMwdFwXq2JUQjXPtcNdlG8H9wcJ1SVItwylx5CPhLVtchtJD5Z6cVvaBesMGNx+469KIca0d7rJsI7jv9bsuSepjFEEkLuJchjWkGRcK1LHlYxENHycSUNrqWWF7rnr9rkuS+phmIxIncS3DCkpOCylQS1FcrZ4Vcn3dzWVJUh9Lj4rESRzLsPqecZGk4XQFaokN19fd8p3VF9MD9rH0qMhS4/Kk1mdyWtKG0xWoZcnJd7Ao9YdYLkuPipTC1Umtr+S0JA6nK1DLkqQesMjicXFS6yM5zXfZ1B//ONpzCqVALUuWesAiyeUjOc3nOuG5mfMuqX8hRQnXr4XgVtWGRMQl15UfYWFXIHOp7AP1iy+6CTLltKSir/Vr4/yaIf77J7LUuF7YZCHLprpU9kPfe/bAvn2lJRklYUlFV1MRshMxamrOPV5qIkacX7Ov/ROR+bnMN9m6FVKpuYeoUyl3ZVNdKftADaUFGdcZhD4yEl0FGZ/rUcf1NfvaPxFJvoWq8V/2Q99Q/MourleL8bH6jMu1lBd6Peo4vOb59s/a0lYE0rV+kbm5XN/6uefmT/gaHCz8GLZQNf4VqKcVE2RcBy7X7bkOgou9HnUhXL9m18tmZvNxrV/X0WUpcXnSDe6PYfMtJeyKhr5zRAkyrj/0hV6bOepUBNcLXkD8X7PrZTNDPq71x/06v5QPF98bH5faXCeTzbUCmUsK1DmiDGW4/tBdb5eEtZRdB3/Xr9n1spng5wAU9+v8Uj5cfW98zHn2MTc73wpkLunceJox0RdrmG/YI2qbrttzHfh9BC3XXL9mH6MIcR/uB/dDjtn7qqH5pcvl98ZHR8PH3GwIgvXevfDJT0Z7XqFiE6iNMZ8wxhwxxowZY540xlwxz/a3GmN+bIwZNca8Zoz5vDGmtrh/O7iN+gG5/tBdtxf3wA/ug//mzfO/PxUVwXaFcL1sJizucH8hfAR+cJsUJG65SGp0/b3xtb6167nZoYoK2LSpuOfO27afZqMxxnwAuB/4I+BngR8B3zXGdOTZ/gbgs9PbbwZ+A/gA8H8X8++X8gG5/tBdthf3wA/uf4yHDs1/IMhkgu0KEb7muSz2a4574Ad/PXRQ5nypXCU1uv7e+DjehMIe8D33wG23Bbd798b3kk5crlF/Ethjrf0rAGPMx4BrgI8SBORcVwH/aK19ePr+EWPM14Ero/7Du3bB9u2lJci4XgDCZXsul52bK3Gi2FEJ19eMXAet8DXfe2/+baK+5rDXP1dAidLrj3vg97kQQnaN5bvvDoJMKqXr6IVymdTo67fnan3r2dpPyloAix6ojTHVwGXA2UOhtTZjjPkHIN9X5IfAh40xV1hrnzLGdANXA1+d49+pAbK+iqQAurst1qZJp0t8IQQLroespeQ2XbV3xRXwtrfB4cNw6hQ0NUFPT/BFjdrmFVcEP96HHoKhoeDJNTVpWlvhxhuDv0dtc9cu+Pzng//P92Ms9PU3N8884My1XZT9nK/NKG0dOgRVVdltp2fcZm+X/R3Ip6cHLrggyDzPd7LT0hJstxjvYW9vEEjnajOcu1rI6w099VTwvQmCzLn3cGjo3PfpijkvoM0uk5n9txInLvYxkwl+x9XVwf3c76Exwd/f9rbC2vbx28s+3vT1nXu8lOONTxlPwznG+i5SOt8OGHMBcBS4ylp7MOvxPwF+3lo7ay/ZGPM7wH2AITjh+JK19rfm+Hf+EPiD3Mcffvhh6urqSnoNIiIiIyMj3HDDDQCN1tozrtpd9B51MYwx7wB+D/g48CTwFuABY8x/sNbenedp9xJcBw+lgJ92dHSwcuVKn7u7ZKXTaV566SW6u7uprKwsuT1XPZmwpwWz99B/93cL72n19p67djeXu+4qvDeY22ZNTZrf//2X+OM/7mZ8/Nz7GKVNCF73XD2PKMrxPczuoWcr5jVnS6fh7/4Ojh+Hzk741/8aiv25uNzHH/4QvvCFc/fzvYe/8ztw1VXR9g9K/94k0YCnua5xCNQngTTQmfN4J3Asz3PuBr5qrd07ff8ZY0w98JfGmHusteeNP1hrx4Hx8L6Z/uYYY5wEmXJWWVnp5D2srHRzzSi8puZiLufAAIyPF7ZdoW9BuDBAbuLN+Hjl2QNkmCQT5UTFdW4DuHkPw9dbSB5Cofua73PJfg/D7Qr5XMLr6GNjs/+92Ovo+/bBgQMz8xH27oVrr4Wbbiq8HR/72NRU2HvY1FT4d9vl9yaJKjxdI1n0QG2tnTDGPA38InAAwBhTMX3/wTxPqwNyg3F4pcJzMTdJAldBy8cUkYoK2LED9u/Pv82OHcUFWJcJMq7eQx9JQa4/Fx/FNfbtm/0zzmTOPR4lWLveRx/FP8B9cq3EZHoWwZD0LmPMjcaYzcB/AuqBMAv8r40x2Xm3fwv8ljHmg8aY9caYdxL0sv/WWhuj1AJZTGHQ2rEjuC12WTzXU0QymaDgx1wefzwe04xcvIfgfhqj68/Fdcby1FTQk57LgQPBdoXylVUNbot/hG27+N5IYNF71ADW2m8YY9qB3UAX8M/Au6y1x6c3WcPMHvQfA3b69kLgTYLgfeeC7bSUBR+9wSgLfUTtHce5NrfLnpbrqYKue+iPPVbYfP7HHoP3vtftvx1ldCdf+ctipm+KP7EI1ADW2gfJM9RtrX1Hzv0pgmInf+R/z6TcuZyLDv5qsCehNrfroXlXQSa8jj5XreZUqvAe+rF82TVFbgfu599nyx36XuTJQJIjNoFaJM5c9gZ99Ix8LMqRBNu2weWXw3e+E9y/8UZ417tg2SIf2bq63G4H0aruFXoylK/gSX9/ad+bOI/sJJECtUiBXPUGXSfx+Kz8FXe5lckeeihI1Cpmpab5Vj4Ki7IU8h24+mr4ylfm7/1efXXh+5iUinFJGNlJmiX2sxWJP9dJPD5qcydBnFdqWrYsmII1l2uvjdbzX8xM90Jp1TU/1KMWWQQur68mYd3xkKshUde9QR+XI8KpV7nzqCsqiptHHfe6+EnroSdpeF6BWmSRZF/3hqCKVjEHCx9rZvvg8oCblDnFN90EH/5wkN197FhwTfrqq4u7hu56BkIS5qL7yr3wEfwzGfjxj4t77nxiev4gUh4qKs6VuNyyJX5n9K6GHF0PiSZpTvGyZcEUrN/8zeC2lEQ3l/PR4z4X3ee66K6H58OlQu+/f/5ti6EetUjCnT7tdruQq16HjyFR33OKXUzD88VXxbhsS7WH7uO7mC9z3iUFapGEi/t0Lx8H3HIvf+lqBoLruehxvoYO7r+LcwV+lxSoRRIu7tO9fBxwXfcGc9t2VZQlCVzlSsT9Gjq4/y4WUmXQhZidJ4pIVOEBMt9ZvbWLO93LxwEX3NcPL2euciXifA0d3H8XF2omhXrUIjKD616Hr2FqcNcbFHfivOqa6+9i1JPLYunrLLKIMhno7Q3+v7e3uKzqcKg6n3CoutC2Xfc6fGZUh+3HOXO+HMV11TXX38X5ev2uqEctskhyy1/u3h0s/FBM+UuXCTI+Fn9ISka1xI/rBD+X38W5ciVcUqAWKZDLSkb5pnQUk1XtI0GmkMUfnnsOLr20sDYhORnVEj+uE/zCxVxcFKLJlznvkgK1SAFcVjKKe/nLZ58tfLsogRrKL6Na4mm23/OBA6WN7vicoqVzWZF5uK5k5Dqr2kd2rEicuFyUw/XvOWyvr6/4fZqPetSyZLkYqvZRychX+UtX2bEXXwzf+EZh24n4FufRLBU8ESmBqx+3j6pacS9/ecklQVLbXNfbUikNYUt+rvI5XC/K4fr3vFAFTxSoJTbi+OP2UVUr7uUvKyrg5pvh3nvzb3PzzcUvT6lksqUtzjXiXf+eVfBEykpcf9w+er9JKH+5bRvccQf85V/OvPZWylKAvtYVltKF8/mrqoLbOJwkJ2E0SwVPpGy4TO5ISqJWEspfbtsGX/4y3HMP3HZbcLt3b/FB2vXSguJGuETj7t3B/d27g/tRPxPXy1L6HM1y9XteqIInCtRSFBcVtcJ24vzj9llVa9u2IPDddVdw/667ig+EvrioMDXfZ2xtcesKS+nifJLsczQL3Pye52rPJQVqiczVGTgk48fts/dbDuUvC0m4ifIZi5vpSnE/SU7KaFa+9lzSNeoy4SNRq9SKWpCcBSBUVat4/f1utyt3cZ3R4KtGvMtFOUI+ypJeeSU8+SR85CPFtTEXBeoyENdELUjWj1tVtYpz+rTb7ZLGV+nZbHGY0eDjJNlnjXjXv+eKCti0yV172RSol7i4Z2Em7cct0TU2ut0uSeJcrCMpJ8kazVKgjq04VtXykYWpH/fS19LidrukiHuxjiSdJJf7aJYCdQyVyzWokH7cS1sYEOb6LsapFnkcT5LBf+nZbDpJjhcFakfKoaqWr0Qt0I97KcsOCPm+N8XmDbgW15Nk8F96NrtkrE6S40WB2oG4Jmv5vgaVrdRErbB9/biXpnyjJnGqTBbnk2RYmBkNEMzn10lyvJT9R/Hii/FZMi0JVbWSUFFL4iks8OKi0lm2qSn49reD///2t4P7UbmeU5yEYh25bS/1+fxJVvYfx549xRfriHvBAF8/7CRU1JJ4clHpLNu+fXD99fDQQ8H9hx4K7u/bF62dJJwkg06Uy5WGvolPFmbclz/MFp6B/+QnOgOXwrmcU7xvH+zfP/u/ET5+002FtRX3NcKzKZ+j/ChQE58sTFXVkqXM5ZziqSk4cGDubQ4cgA9/GJYVcJRL0kkyKJ+j3ChQT4tDFqaqakncxHE2A8Bjj81/SSmTCbZ773vnb08nyRJnCtQ5omRhbt4c/ODmOmBUVATbFWrbNnjf+4LeQG6gvvZaXYOShRPX2QwAx4653U4nyRJnOq/LEWVo69Chws7qDx0qvM2DB+HRR89vN5MJHi923V4Xq+1I+YjzbAaAri6324EStSS+1KOeVszQlutr1HP1PEJRex7g9tpg9r729kJVVXCr4bx4iGNVLR9ziq++Gr7ylflHs66+uvA2QUPVEk8K1BQ/tOX6GrWPakaurw2GbYaVjO6+O1iPOpUqPUHGZUZwUrg84YlrVS0fiVrLlgWXgmbL+g5de21hiWQicaevMcVnYbpOQFnIHnqx1wZdr0ed3a6PXr/LwO+6PZcnPHGuquUrUSucepWb/V1REQTpQqdmZfPxPRQpVdkH6l27YPv24pNEXCagxL2H7iPwg99ev6sDro/2XJ3wJK30rMtErZtuCqZgfec7wf0bb4R3vau4nrSP76GIC0t8YHF+GzaU1itymYDiupqR656Rj6Qg19XdwG0ilI/2XL/mJFTV8pmotWwZvPvdwf+/+93FBWkf30MRV8q+R+2CqwSUuPfQfSQFxb3X72MUwfVrTkpVrTgnavnIDxFxJQY/kaXBVQ3jOPfQfSQFxb3X72MUwfVr9llVy3UP2HWtb1d8nISKuKIedQzFtYfuIyko7r1+Hwdw169ZVbVK5+NkR8SVJfiTWxri2EP3sRpX3Hv9Pg7grl+z7+UP49gDds3XalciLizRn51kc7kOsOshUddBxvUB18cB3EdgVVWt0vg82REplYa+y4TLesPZQ6IQrEddypCoy1WGXA/3+0ysCl/z4OC5x03C10gAACAASURBVEtZWWnbNrj88mAhimPHgvKZV1+toh+F8rnalUgp9BOWorhej9rl9VDXB1yfB/Dca8pzlY+dz2xzvQ8cUJCJyuVnIuKCArXEhutev+vepeugOlvBk/7+4oprqFhH6fK9h8V+JiKu6IqLxIbLFb4OHoTf/M3gWvx/+2/B7W/+ZnGrj4UH8L6+mY+HB/DFLniiYh2l03socaYetcSCyxKdLnuXSSh4omIdpdN7KHGmHrUUJVz1CYLbUnu/rkp0xr08JyRjrne50XsocaZAXSZcDyvv3Bms9gTB7c6dxQ0rxz2wJqHgiYp1lE7vocSZhr7LgK9hZRfLXMa97rXPgieuKon5qkxWTvQeSpypR+2Iyx6ry/biPKwM8Q+sSSh4omIdpdN7KHGmHrUDcV372HUilI+Em7jXvU5KwRMV6yid3kOJq7IP1C++CCtXFn+m7Hr+6sGDcO+95z9+8mTw+B13FN5e3IeVIRmB1dcB3EeFt3JZRMMXvYcSR2UfqPfsgX37ijvgFjoUHGXt4wcfnHubBx8svL24DytDcgKrrwO46wpvLovGlCu9hxI3ZR+owV8iFETrsT7zzMxh0NkMDgbbXXrp/O35GFZOpebex1SquOUU3/e+oNxlbqC+9tr41L3WAVxEFoMCNcUXrujvd7vds88Wvl0hgXrz5uC1zJXcVVERbLeYDh6ERx89f2Qikwke37SpuOx013WvMxkNiYrIwtNhZloxhStOn3a7nWuHDs2fgZ3JBNsV4rnnCuvxR3kP57p8EIqaSe4y0z27zZ074c474b77gtti545nc1k4RkSWJgXqHFESoRob3W538cVut0tCBSzXBUp8TCHzEfjDdl0Vjgm5niYoIotPQ985oiRCtbS43e6SSwq7BrxYU598JJO5Dv6uM9191PoG94VjwjZdThMUkXhQj3paMYUrwqlFc4nSZkUF3Hzz3NvcfHPhAcF1sQ4fxT9cB//FDPyFSlKvXz10kcWnQE3x04DCqUVzBa5iphbdcUcwlShbW1u0OdTZ+xfuS+6+weJXwHId/OMe+CEZw/3g77q8iESjQE0QFItdFD6cs5vbs25rK63NL38Z7rkHbrstuN27t7T9yw38xb5m1+25Dv5xD/yQjF6/eugi8VH216h37YLt20ubZuOjGIbLObuu989HRS1XBUpcF1DxsVhD3Hv9Pq/L6xq6SHRlH6g3bHAzFzbuxTBc75/rilouTybiHPjBffB3Hfh91HR3XWo3WzjFraoquNX89uj0HsZb2QdqiY+4jiK4LkmaG/yzxaHXn5QeOpzrpQ8Owt13B1PcUin10qPQexh/CtSyZMU18IftuVo9y3WvPwk9dPAzxQ3KqwKdr/dQ3FKgFimQ68sHLq/1u+z1x72HDrqO7oLPkQ5xKzZvvzHmE8aYI8aYMWPMk8aYK+bZvskY80VjzBvGmHFjzPPGmKsXan9FXAiv9UPp1/q3bQtmB5Q6W8B1Jr6PzPkkZbrHlY/3UPyIRaA2xnwAuB/4I+BngR8B3zXGdOTZvhr4e2AdcD2wCdgFHF2I/RWJq7DXv2NHcFvq0LyLaXg+CuUs9HV0KG4ueth2HKek+RjpED/iMvT9SWCPtfavAIwxHwOuAT4KfHaW7T8KtABXWWsnpx87sgD7KVI2XF2X95E5n6Tr6HEdSvcx0iF+LHqPerp3fBnwD+Fj1trM9P18X+X3AgeBLxpjjhtjnjXG/J4xptL7DouUkTj20MF9L91H79LnULqLXrqPkQ7xIw496jagEjie8/hxoCfPc7qB/wv4GnA18BbgL4AqguHz8xhjaoCsvEZSANZa0ul0sfte1sL3Te9facrlfbziCnjb2+DwYTh1CpqaoKcnCP7FvPRdu+Dznw/+v7o6aKCmJj2jl25tYW03N8/Mep5ru0Lay2TgoYegunr2vxsT/P1tb4t+8vPUU8Fz+/rOPdbaCjfeGLzHUbh8DwUynq5rGDvXQsALwBhzAcG15austQezHv8T4OettVfO8pzngVpgvbU2Pf3YJ4FPWWtX5vl3/hD4g9zHH374Yerq6ly8FBERKWMjIyPccMMNAI3W2jOu2o1Dj/okkAY6cx7vBI7lec4bwGQYpKcdArqMMdXW2olZnnMvQcJaKAX8tKOjg5UrZ43tMo90Os1LL71Ed3c3lZW66lAsvY+lyWTg8OE0VVUvMTnZTU9PZVFD9E89da53Odt19N/93cJ7rD/8IXzhC/Nv9zu/A1ddVVibmQz89m/P7ElnMyZYUvfP/zx6L93Ve1juBjxl3i16oLbWThhjngZ+ETgAYIypmL7/YJ6n/SNwgzGmYvp6NsBG4I08QRpr7TgwHt43078+Y4wOjiWqrKzUe+iA3sfiVFZml7Mt/j0Mr5O7SP5qaoLx8cK2K3R3e3vh9dfn3ub114NLC1Hn+8+8Th28hwrU0VV4etMWPVBPux94yBjzT8BTwK1APRBmgf81cNRae8f09v8JuBl4wBjz58BFwO8BBZzDiojMzlWmu4/FXHxNp1IJ0fiLxTmTtfYbwG3AbuCfgbcC77LWhglma4CVWdu/BvwScDnwLwQB+gFmn8olIlIwF5nuPtZu9zGdqtyKvCRVXHrUWGsfJM9Qt7X2HbM8dhD4Oc+7JSJSFNeLubjupauEaHLEJlCLiCw1LhdzcV04xleRF3FPgVpExCPXq7i56qX7LCFaTiuQLQQFahGRBHHVS/dVQjTOZVOTSuc4IiIJ4yLhzUcJUSWn+aFALSJShlxnpvtcgazcKVCLiJQpl4ulaH1rf3SNWkSkjGVf8wa4667irnlrfWt/1KMWESlzFRVBGVYIbou55u1zfWsXy3ommXrUIiJSMh9lU0FZ5KAetYiIOOCjbKqyyAMK1CIi4oTL5LT5ssitLZ8scg19i4iIM64KssyXRQ7lU+JUgVpERJxyUTa1v9/tdklWUqA2xlQBXUAd8Ka1tgzeMhER8e30abfbJVnka9TGmJQx5reMMf8dOAMcAQ4BbxpjXjHG7DHGXO54P0VEpIw0NrrdLskiBWpjzCcJAvNHgH8ArgXeCmwEtgF/RNBL/ztjzHeMMRc53VsRESkLLS1ut0uyqEPflwM7rLX5isA9BXzFGPMxgmD+duAnJeyfiIiUoXBe9lwJZVEXDUmqSD1qa+2HwiBtjLnKGHNxnu3GrbVfstZ+xcVOiohIeQnnZc+1ulfUedlJVcpL/CJwZe6DxpgNxphUCe2KiIicnZfd1jbz8ba26POyk6yUrO9NwPdnefxfAe8BfrmEtkVERJzNy06yUgL1GWC28upPAPeU0K6IiMhZLuZl58pkkhP8SwnU3wFuAz6Y83gGqC6hXREREW+SttBHKecP/wH4eWPMfzXGXAJgjKkFPgP8i4udExERcSmJC30UHaitta8BPwcsB35kjBkFBgmuT3/Kze6JiIi4Md9CHxDPhT5KKiFqrX0FuNoYs4ag8Mkk8KRKiYqISNzMt9CHtfFc6KPoQG2M+Tiwx1o7aa19FXjV3W6JiIi4NTDgdruFUso16j8HftkY0577B2NMDC/Hi4hIOWuebZ5SCdstlFKGvg3wCGCMMSeAZwiSyI4Q1Pxuzf9UERGRhRWWJe3rm/06tTHQ2hq/sqSlzhrrBq4Afg/oBS4DPg38jxLbFRERcSosSwrnlyYN78exLGlJyWTAmLX2aeBpFzsjIiLiU1iWNHcedWtrfOdRlxqoNxtjBqy1k072RkRExDMfZUkzGfjxj93tY7ZSA/X/B0wZY57n3DXqfwH+xVr701J3TkRExAeXZUnDSmeDg27ay1VKoD4M/CrQBVwM/B/A+4DfJyiCUlny3onIonJdD9lHfeVMBnp7oaoquI3jPsrSFVY6sxZqavz8G0UHamvtlun/PQR8L3zcGGOADSXul4gUwWWQcV0P2Ud95eyezN13w+7dkErFax9l6Zqr0plLRf2EjTHLjTHbjTFbZvlzDXBVabslUh7C3iAEt6WULjx4EHbuhDvvhPvuC2537iyudrHresg+6isnYR9laZuv0pkrkQO1MWYjQS/6ceAZY8x/N8ZckLVJI/BXjvZPJDYyGXjmGXj88eC21HrAYWDdvTu4v3t3PAKr63rIPuorJ2EfZembWcHMUl097uXfKaZH/TngWaAD2ESwEMcPput9iyxJLnurYXtxDaxR6iEvRntJ2UdZ+sIKZrW1Y7S29mGMnzHwYgL1VcAd1tqT1toXCFbL+i7whDGm2+neiZTAVQ/Y9ZBo3AOr63rIPuorJ2EfZenbsGGM7u4+qqomOXp0FS++uNHLv1NMMtlyYCq8Y621wG8ZYx4E/jtwg6N9kzITx0So+YKqMcHfr7yy8H11vYKP6yDjuh6yj/rKSdhHWbrGx8cZGhqiqqqKD31oJR//eAcjI/XU1PR5+feKCdSHgbcRXKc+y1p7c5Dwzf/rYL+kzLjMts2eLpEt7AHffnvhbfpYFi/ugdV1PWQf9ZWTsI+y9ExMTDA0NMSyZcvo6uqio6ODyy9voKEBbrnFX2JZMf2VR4EPzfYHa+3NwNcJFuwQKUicr9f6GBL1FVhzaxeHjAn+XmiQcV0P2Ud95STsoywdk5OT9Pf3Mzo6SkdHBz09Paxfv56GhgYArrsOjhyBRx7x8+9H/tpZa++11l49x98/bq3V13mJm5qCb387+P9vfzu4X4y4X6/1MSQa98AK5+oht+asgdfaGm1Ewld7SdlHia9CclgmJycZGBhgaGiItrY2enp66O7uJpVKYXJ+bJWV/r4jpZYQlTK0bx8cOBBUgrr7bnjoIdi7F669Fm66KVpbcb9e62NINAysn/2s+8DqcqEB1/WQfdRXzm4T4K674rePEj/zXWqbmppicLoeaHNzM11dXaxYseK84LxQIgVqY8waa+2rEba/0Fp7NPpuiWuuErX27YP9+2dvP3w8SrCO+/Xa3KCaHaxLGRLNDqzZ9YHjFFjBbT1kH+2FbW7ZAj/5SXBbalD1sY8SH3PlsPzpn05xyy2DXHKJpaWlhc7OThobGxctQIei9qj/pzHmALDXWvs/Z9vAGNMI/BvgFuAvgS+UtotSKleJWlNTQU96LgcOwIc/DMsK/GbFPREK/C2L57o3CAoyInPJd6mtoiJNKjVIZWWGr361iUcf7aKlpZGKmAylRA3UW4A7gb83xowRrEP9OjAGNE//fSvwv4BPW2sfc7ivsRbXxQtcZkA/9tj814ozmWC79763sDZdB1afPWAfQ6Kue4Mikl/upbaKijQNDUNUVU1x6lQTx493cupUM88+W8E73rFou3meSIHaWtsHfNIYcydwDbAdWEswt/ok8DXgu9baZ13vqC8vvggrV5Z2gIzr4gWu5wAfO1bYv1voduAnsPrqAau3KpJs4SU0YzKkUkNUVU1y5kwjL7/cycBAM5lMsOjjG28s4k7OoqhkMmvtKPDN6f8Sbc+e4LprKUHVVY/VdXuuE7W6ugr7dwvdLpSERCgRWTyuRhibmjKkUsNUV08wOJjiyJF1DAy0kE7PXJV55UpHO+5I0VnfxphPA28lWI96FHgOeNRam7g1ZooNqq57rK7bc52odfXV8JWvzD38XVERbBdVEhKhRGThuRhhzGQyDA8Ps3LlBA0NDfT2rqG/v4V0+vwQuHo1vP3t0fcznfa3wlop/YvfBtqAE9P3P0SwOMd3phPKEqPY1XHivjCA60StZcuCKVhzufbawhPJcoWBdceO4Fa9X5HyVmoxJGstQ0ND9Pf3U1VVxUUXbeCKK7bw5psdswZpgA9+MJgTHcX+/bBuHbz//dGeV6iiD4XW2tXW2n9trf2gtfYaa+1q4P8EOoEvOtvDBVLM6jhxXxjAdWENCKZeXXfd+UG0oiJ4POo8ahGR2ZRSDMlay/DwMH19fVRWVtLd3c2WLVtoa+vk61+fuyfxN38T9I4LtX8/XH89/PSnhT8nKqd9Fmvt/wA+AhSY8xs/i1kK0tccYHBbFvGmm+Cb34Qbbwzu33hjcF9BWkRcKWaE0VrLyMgIfX19GGNYv349W7ZsYeXKlVRVVfHEE/MH1NdegyeeKGwf0+mgxvdsJxMuOalMZoz5CMG61GPAtYCfJUQWQDGlIOO8MICvDOhly+Dd7w6mFb373dGHikRE5hJlhNFay9jYGMPDwyxfvpy1a9fS3t5OTU3NjG0LzeYudLtCAr8LrkqIXgm8H2gCvkUCe9QuSkGWOrUoaXOARUR8aSww06m+fpT+/hFqampYs2YN7e3t1NbWzrptodnchW63UNO4nByqrbUfI0gs+2WgG/hZF+0uFBelIOO+MIAStURkKampGaO1tY9MZopVq1axZcsWVq9enTdIQ5DNvWrV3Hk7UbK+F2oaVynTsx4HPmWtfRLAWmuBbxtj3gAeAx5ys4v+uSwFGdfFC0REkuT06dkfr64eJ5UaYmKiitdfvwDoYM2auoLarKyEBx4Ikr/yjVr+2Z8VfikvDPxHj/q9Tl3K0PdzwD8aY54C/ivwDDBEME1ruYN9WxC7dsH27fEr5K85wCJSznLzhaqqJlixYpDJySqOHevixIlOhofrWbUqWrvXXRckv95yy8zry6tWBUH6uusKbys38PtSdKC21v6WMeZB4FPAXUAq/BPwew72bUFs2KCeqohI3ITJtWfOTNDQMEQ6Xcnx452cONHJ0FBD5GHqbNddB7/8y/AXfxGUkd6wAT7+caiuLq6t226D+++P/txClZRMZq19DrjJGPMbwAaCZLJXrLXHXeyciIiUp3R6kl//9UG+8pUKTp5s5/jxDgYHU4Apapg62/795/eo/+N/DHrHUXrUYVv33RcMfVdVRd+XQrhKJktba5+31j6lIC0iIsWamppiYGCAoaEhfuEXWvnd3+1hbGwDg4MrgCBCr1oVDF9HDaqQv0DJ0aPB4/v3F95WouZRi4iIlGJqaoqhoSGstTQ3N9PZ2UljYyMbNxre//5gzvIbbwSZ1m9/e3E96bkCa7iewq23wq/8SmHtJ20etYiISGTpdJrBwUEymczZAN3U1ITxkJ01X2C19lxlskLWo16oedQK1CIi4lQhy1Km02mGhoaYmpqiqamJrq4umpqaqMjZcLbryatWFXc92XVlstjPoxYREck137KUmUyGwcFBpqamaGxspKuri+bm5vMCNJy7npw7VB1eT456ndp1ZbKFmketiUkiIuLEXMtSfu5zGR5//AwDAwPU1dWxceNGenp6aG1tnTVIz3c9GYLryVFWunJdmSycRx0+1xcFahERKVm+ZSmNyVBfP0hLSz9f+1oN3d1voaenh7a2NirnyNiKcj25UHMF1mKnfIUFVC68sPDnRKVALSIiJTt/WUpLff0QLS39TE1V88ILb+Ef/3ELhw93sGzZ/FddXV9PDuULrKVM+bruOjhyBB55JPpzC6Fr1CIiUrJzy1Ja6uuHqa0dY2Sknpde2kBfXytTU0E1ENeJWsUkdF13XTAFy8WUr1BlZfHrRcxHgVpERErW1GSpqxth+fJRRkfrOHJkPSdPtjE5ObMup6tELWOCvxdTQhSCwFrIFKw40NC3iIgUzVrLyMgIXV19NDTAK6+so7d3C2+8ccF5Qbq11U2iVqklRJNGgVpERIoyOjpKX18fmUyGNWvW8uKLW3j99QuZmKhx0r6P68lJpKFvERGJZGxsjKGhIWpra1m9ejXt7e08+eTyea8/9/UVXvUr5ON6ctIoUIuISEHGx8cZGhqiurqaVatW0d7eTl1dHeAvSxuSdT3Zh1gNfRtjPmGMOWKMGTPGPGmMuaLA533QGGONMQd876OISLkZHx+nr6+P8fFxVq5cyebNm1m7du3ZIA1+s7TLXWwCtTHmA8D9wB8BPwv8CPiuMaZjnuetA+4DIkx7FxGRUCYDvb3B//f2BvcBJiYm6OvrY2xsjK6uLjZv3sz69eupr68/rw3XVb/knNgEauCTwB5r7V9Za3uBjwEjwEfzPcEYUwl8DfgD4KUF2UsRkSXk4EHYuRN27w7u794N/+7fTfL97/czMjJCR0fH2QDd0NCQtx1lafsTi2vUxphq4DLg3vAxa23GGPMPwFxTyO8CTlhrv2yMmfM8zRhTA2SnIqam/x3SUYrFylnh+6b3rzR6H0un97A4Tz0Fn/98ME+5piZ47zo6+pmYqODBB5v5/d/v5Gd+pgFjDFNTU/O29573BNnYn/lMMP85tGpVUAP8Pe+ByUlfr2bxTXp6ccb6XPKj0J0w5gLgKHCVtfZg1uN/Avy8tfbKWZ6zHfgb4K3W2pPGmH1Ak7X22jz/xh8S9LxnePjhh2dcZxERESnGyMgIN9xwA0CjtfaMq3Zj0aOOyhiTAr4K7LLWnpxv+2n3ElwDD6WAn3Z0dLBS2Q1FSafTvPTSS3R3d89ZXF/mpvexdHoPo+vthXvuSdPQMERFRYbR0RXceedL/MZv/CtGRs4VKvnWt2D79kXc0QTp6+vz0m5cAvVJIA105jzeCRybZfsNwDrgb825iyEVAMaYKWCTtfbF7CdYa8eB8fB++DxjjH7YJaqsrNR76IDex9LpPSxMOp3mxIlBamsznDzZxPHjnYyP1wMvMTJSzeho1dltjx2Dqqr8bck5VZ7eqFgEamvthDHmaeAXgQMAxpiK6fsPzvKUw8AlOY/9MUEv+RbgNX97KyKSTOl0mqGhIdLpNG1tK3j++S5OnWomk6lg+fLZr69qwHHxxSJQT7sfeMgY80/AU8CtQD3wVwDGmL8Gjlpr77DWjgHPZj/ZGHMKwFo743ERkXKXyWQYGhpicnKSFStW0NXVxWWXNfPv/31l1qpXM5W66IW4E5tAba39hjGmHdgNdAH/DLzLWnt8epM1QGax9k9EJGkymQzDw8OMj4+TSqVYt24dLS0tZy8PPPAAXH+9plPFXWwCNYC19kFmH+rGWvuOeZ57k4ddEhGJnUwGnnsuWAO6uRm2boWKiuy/zwzQa9asoaWlhWXLZh7yw0UvbrklqMMdWrUqCNLlsuhF3MUqUIuIyNwOHoQ9e+Bk1nyXtjbYtQt+7ucsw8PDjI2N0dDQwKpVq2htbZ0zySlc9OLxx+HMmSDLe8cO9aTjRIFaRCQhDh4MCofklr/o67N84QsjjI+PcsUVdaxfv562tjaqq6tnbyhHZWUwBeuxx4JbBel4iVMJURERySOTCXrSM4O0ZfnyEVpb+wDDl760jk2btnDBBRcUHKQl/tSjFhFJgOeeyx7utixfPkZ9/TCjo7W88spaTp5sZ3y8hqeeKu8lIZciBWoRkQQIp1HV1o5RXz/E+Hgtr766mpMnOxgbqz27XTHrPUu8KVCLiCRAKjVGa+swExPVHD26ijff7GB0dPl526lAydKjQC0iEmPj4+MMDQ2xfn0Vk5MrOXy4g+Hh89eDVoGSpUvJZCIiMTQxMUF/fz/j4+N0dXWxdetm7rhjPSMj9SpQUmYUqEVEYmRycpL+/n5GR0fp6Oigp6eH9evX09DQcLZAyYUXznzOqlXB4ypQsjRp6FtExKP5qoiFJicnGRwcxBhDW1sbnZ2dpFIpTE73OSxQ8sQTQeLYypXBcLd60kuXArWIiCdzVRHbti24PzU1xeDgIAAtLS10dXWxYsWK8wJ0tspKTcEqJwrUIiIe5K8iFjz+6U9PsWXLINZaWlpa6OzspLGxcc4ALeVJ16hFRBybvYpYwJg0jY2neOSRM6RSjfT09HDRRRfR1NSkIC2zUo9aRMSxmVXEAhUVaRoahqiqmuLUqSZeeKGTEyea2bpV/SWZmwK1iIhjYRUxAGMypFJDVFVNcuZMIy+/3MnAQDOZTCXHji3ePkpyKFCLiDjW3BwE6IaGIaqrJxkcTHHkyDoGBlpIp8+lZ6uKmBRCgVpExKFMJsPatcN0d09w9GgDr766lv7+FtLpc4dbVRGTKHRxRETEAWstQ0ND9Pf3U11dxa/92gYOHdrMyZMd5wVpUBUxKZx61CIiJbDWMjIywujoKPX19XR3d9PW1sall1axfDnccgv89Kfntl+1KgjSqiImhVKgFhHJUmglMWsto6OjjIyMUFdXx/r162lra6O6uvrsNqoiJi4oUIuITCukkpi1lrGxMYaHh1m+fDlr166lvb2dmpqaWdtUFTEplQK1iAjzVxK7/XZ461tHGR4epra2ljVr1tDe3k5tbe2c7abT6lFLaRSoRaTszVVJzFqorR3jv/yXYbZurWb16tW0t7ezfPnyedvdv3/2a9QPPKBr1FI4BWoRKXuzVRIDqK4eJ5UaYmKiit7eCzh1qoPLL68rqM39++H6688P/kePBo+XsiyleunlRdOzRKTsZVcSA6iqmqC1tY/a2nGOHevi8OEtvPLKOk6eLCxIp9NBTzpfDx3g1luD7aLavx/WrYNf+AW44Ybgdt264PFipdPwgx8E//+DHxS3X+KPArWIlL3GxuC2qmqClpZ+6upGOX68k8OHN/Pyy90MD9cD0NFRWHtPPDFzuDuXtfDaa8F2UYS99Ny2w156McE6DPzXXBPcv+aa0gO/uKVALSJlL52epLm5n/r6EU6ebOfw4R5eeqmboaGGotp74w232wX76L6X7iPwh/v6/e/D178e3KqHXhoFahEpW1NTUwwMDNDfP0R/fyuHD/fwwgsbGBxcAZy/5OSJE4W1W2gN7yi1vl330n0Nz/sami/nwK9ALSJlJwzQg4ODNDc3s2FDDy+8cBFnzjQyW4AOFRpY3/72ILs73/LSxsDq1dFqfbvupfsYnvc5NO8y8IP74J9OB1P8fFCgFpGykU6nOXXqFGfOnKGpqYlNmzaxceNG3vnOJi680DgLrJWVwRSs8Lm5bUH0Wt+ue+muA3+ShuZdB/+wvfe/v7jnz0eBWkQSLZOB3t7g/3t7g/u50uk0p0+f5tSpU6RSqbMBurm5GWOMl8B63XXBFKwLL5z5+KpVxU3Nct1Ldx34RtSPpgAAIABJREFUkzQ07zL452vPJQVqEUmsgwdh507YvTu4v3t3cD8cgsxkMmcDdH19PZs2baKnp4eWlhYqcgp4uw6sYZtHjsD3vgcPPxzcvvxycW25PplwHfiTMDTvOvjP1Z5LCtQikkhhyc/cQiV9ffC5z2V4/PEzDAwMUFdXx8aNG+np6aG1tfW8AJ3NZWANhbW+P/Sh4LaUwiQuTyZcB/64D82D++A/X3uuqDKZiCROvpKfxmSorx+htnacr30txd/8zVra2lpYtqzwQ13cF9FwuSJXGPhvuSU4wQkVsxRn2EM/enT2HqYxwd8Xa2ge3Af/KCcJpVCgFpHEOb/kZxAZmpsHGBho5IUXVk9Pt1oW66BbLJcnE2Hgf/xxOHMGvvUt2LEjeuAPe+jXXx8E5exgXcrQvKvAD+6Df5SThFJo6FtEEudcyU9Lff0QLS39ABw5sp7e3s2cONHJ1NSyBevxJF1lJWzfHvz/9u3FD8/HeWge3F+Xn689VxSoRSRxmposdXXDtLb2YW0Fr766DmA6QFed3W6hejxyjsvr/K4T/FwH/7nac0lD3yKSGNZaRkdH6eoaoaVlOYcPr+PkyTYqK2f2OYoZFhV3fAzNu1otLPu6fO7yo1Gvy+e2N9sKbC4oUItIIoyOjjI8PExtbS3r16/lYx9r4wMfqAWgsnLy7HbFDotKfLlO8PMR/H/lV+Db34b3vMfdfoYUqEUk1sbGxhgaGqK2tpbVq1fT3t7O8uXLef/7gwOri4xlKT+ug39lJWzb5q69bArUIrJgMpkgY3tgAJqbYetWyDeteXx8nKGhIaqrq1m1ahXt7e3U1c1cD9pVxrJInClQi8iCOHgwmPucfR2vrQ127ZrZEwkDdFVVFStXrqSjo4P6+vq87YYZy489VlrGskhcKVCLiHdhFbHc+bB9fcHjt98Ol102weDgIMuWLaOrq4uOjg4aGopbD1pkKVGgFhGv8lURg+CxqqpJvvGNQXp6Kujo6KCzs5OGhgaM78mpIgmhQC0iXp1fRSywbNkkqdQQ1sILL7Ry+nQXV16ZUoAWyaFALSJenasiFqisnCKVGsIYS39/MydOdHHmzAr6+/OvBy1SzhSoRcSr5ubgNgzQFRUZBgaaOX68k9Onm7A2iM6qIiYyOwVqEfGqpyfNunWDnDmT4dSpJo4f7+TUqSasDeZlqYqYyNwUqEXEi3Q6zdDQEOl0mg9+cAWf+lQXAwPNZDLnJk6ripjI/BSoRcSpTCbD0NAQk5OTrFixgq6uLq64opnGxkpn9ZVFyokCtYg4kclkGB4eZnx8nBUrVrBu3TpaWlqonO4qu66vLFIuFKhFJK9CSn5mB+hUKsWaNWtoaWlh2bLzDy+u6yuLlAMFahGZ1XwlP621DA8PMzY2RkNDA6tWraK1tZWqqqr8jYpIZArUInKeuUt+Wm67bYTNm0epr69n/fr1tLW1UV1dvTg7K7LEKVCLyAz5S35aamtHqa8f4etfr2PfvnV0drYrQIt4pkAtIjOcX/LTUls7Rn39MGNjtbzyylpOnmznxRdrWL16sfZSpHwoUIvIDNklP4MAPcT4eC2vvbaakyc7GBurBYLMbRHxT4FaRGZoboaamjEaGoaZmKjm6NFVvPlmB6Ojy2dsp5KfIgtDgVpEzhofH2flyiE6O6t44YWVnDjRychI3YxtVPJTZGFVzL+JiCx1ExMT9Pf3TwfqLnbu3Mwrr6xndPT8IA0q+SmykNSjFiljk5OTDA4OUllZSUdHBx0dHTQ0NNDdbaiqQiU/RWJAgVpkCSmkkhicC9DGGNra2ujs7CSVSmGyFoRWyU+ReFCgFlki5qskBjA1NcXg4CAALS0tdHV1sWLFihkBOptKfoosPgVqkSVg7kpi8OlPT7FlyyDWWlpaWujs7KSxsTFvgBaR+FAymUjC5a8kBsakaWw8xSOPnCGVaqSnp4eLLrqIpqYmBWmRhFCPWiThzq8kBhUVaRoahqiqmuLUqSZeeKGTEyea2bpV5+YiSaNALZJw2ZXEjMmQSg1RVTXJmTONvPxyJwMDzWQylRw7tnj7KCLFU6AWSbjm5iBANzQMUV09yeDgCo4cWcfAQAvp9LkUbVUSE0kmBWqRBMtkMqxdO0x39wRHjzbw6qtr6e9vIZ0+99NWJTGRZNMFK5EEstYyNDREf38/1dVV/NqvbaC3dwsnT3acF6RBlcREkkw9apEEsdYyMjLC6Ogo9fX1dHd309bWxqWXVrF8uSqJiSxFCtQiCZAdoOvq6li/fj1tbW1UV1ef3UaVxESWplgFamPMJ4BPAV3Aj4DfttY+lWfbXcCvAxdPP/Q08Hv5theJo0wGenuhqiq4zS35aa1lbGyM4eFhli9fztq1a2lvb6empmbW9lRJTKQw6bTbk9p0Oig85ENsArUx5gPA/cDHgCeBW4HvGmM2WWtPzPKUdwBfB34IjAGfAf7OGLPVWnt0YfZapHhhyc/BQbj7bti9G1KpcyU/R0dHGRkZoaamhjVr1tDe3k5tbe1i77bIgnMdVPfvn/0y0QMPFHeZKGwvt56BK7EJ1MAngT3W2r8CMMZ8DLgG+Cjw2dyNrbW/ln3fGLMT+FXgF4G/9r63IiXILvmZ3Tnu64PPf36M8fFh3va2GlatWkV7ezvLly9fvJ0VichlYPURVK+//vxKfkePBo9/85vR2s1uz9d5dCyyvo0x1cBlwD+Ej1lrM9P3txXYTB1QBfQ730ERh/KV/KyunqClpY/q6gn27r2QjRs3s2bNGgVpSZT9+2HdOviFX4Abbghu160LHi+mreuvnxmk4VxQjdpmOh0E/dnK7YaP3XprsF2p7bkUlx51G1AJHM95/DjQU2AbnwNeJyvYZzPG1ADZF/ZSEFwDTBf6qcgM4fum9y+a3t5guDvsSdfXjwOQSo3yxhvt9PV1MDxcxz/9E2zfPrmIe5ock5OTM24lmnQafvjD4L174olJrrqquB7w3/4t/Nt/GwSu7PPL/v7gcYD3vKfwffrMZ/L3Uo2B22+Hq68ufF9/8INg1Gquc9+TJ+Hxx2H79ujt1dZOMjZW2L5EYazvU4FCdsKYC4CjwFXW2oNZj/8J8PPW2ivnef7twKeBd1hr/yXPNn8I/EHu4w8//DB1dXUl7L2IiAiMjIxwww03ADRaa8+4ajcuPeqTQBrozHm8E5izQrEx5jbgduBf5QvS0+4lSFYLpYCfdnR0sFK1FYuSTqd56aWX6O7uplJzgAr2zDNTfPGLg2QyFQwMtDA42MKf//kP+ehH38noaNXZ7b71rcLO6iXoSf/93/8973znO6mqqpr/CQLk9oAn+cpX/p6PfvSdjI0F7+FXv1p4D/gHP4Brrpl/u0K/19/8JvzGb8y/3Ze/HAyDF8L1Pua2V1vbV9iORBSLQG2tnTDGPE2QCHYAwBhTMX3/wXzPM8Z8GrgT+CVr7T/N82+MA+NZzz17qyBTmsrKSr2HBZiammJwcJDVq2Fysp1Dhzo5fXoFy5dPATA6WsXoaNXZkp87dmgOdFRVVVVlEahdJGuF11dHRmY+nv09vPXWYG5+IW0fOwajo4VtV8hHtHJlYe2tXFlYexD8plpbg2vcsy8LG+23l9uetX6+e7FIJpt2P7DLGHOjMWYz8J+AeiDMAv9rY8y94cbGmM8AdxNkhR8xxnRN/9ewCPsuktfU1BSnTp1icHCQ5uZmtmzp4bbbLuLMmcbz1oRWyU+Zj6tkrSeeOD9JK5u18NprwXaFKHRgstDt3v72IGjmWzbdGFi9OloN+8rKIFs8fH5uexDttzdXey7FJlBba78B3AbsBv4ZeCvwLmttmGC2Bsj+iH8LqAa+CbyR9d9tC7XPInNJp9NnA3RjYyObNm1i48aNNDU18au/avjmN+HCC2c+Z9Wq6NNDJN7Safj+9+HrXw9uS8m9dJkF/cYbbrdzHVhdB9XQddfh9LeXrz2XYhOoAay1D1pr11pra6y1V1prn8z62zustTdl3V9nrTWz/PeHi7HvsvRlMvDMM0FG6DPPBPdnk06nOX36NKdOnSKVSrFx40Y2btxIc3PzjB70ddfBkSPB9TAIbl9+WUF6KXE5Vcn11CLXPWAfgdV1UM1u98gR+N734OGHg9tSfnthe488Utzz5xOLa9QicRdWEcuuPNTWdq6KGAQBemhoiKmpKRobG+nq6qK5uZmKivznw5WVQdLKY48FtxruXjpcF9aIMlRdSBnZsAc83/XaKEPLYWB1uTiMrxr2rsvtVlaeOxa4pkAtMo/sKmLZ+vqCxz/zmQwXXzzE5OQkK1asOBuglWCXTC4TtfL1fqMmaoH7oeqwB3z99e6Hll0H1nKvYR+roW+RuMlXRWz6rzQ0DPLII/1UVdXwlre8hZ6eHtra2hSkEyquiVrgfqga/A0th4H1Qx8KbvVzKI0Ctcgcnnvu/EL7xmSorx+ipaWfyclqnn76LfT3b6Gjo4NlyzRIlVRxTtQCP1nQoFyJJFCgFpnDwED2PXs2QGcyy3jppQ309m7mxIlOjh9XgF4s6XRQeAKC22KyquOeqAX+sqDDtsMCH8qViB8FapE5NDcDWOrqhmlr68PaCl5+eT29vZs5fryLqamgwIGK2y2OcKg6rA51zTXxGKr22fvVtL7yo26ASB7WWtavH2XDhhGOHVvOkSPrOHmyjYmJc2u7FJMZK25kZ1VnL7JQTFa170St7J56qb1fX1nQEl/qUYvMYnR0lL6+PiDDhz60lt7erbzxxoXnBWlQFbEoXBX/SMJQtc/er5K1yot61CJZxsbGGBoaora2ltWrV9Pe3s5lly2nvt7t3NBytH//7O/hAw9Efw+TMKcY1PsVNxSoRYDx8XGGhoaorq5m1apVtLe3z1j+VAfc0rgu/pGkoepynwMspVOgliUrkwmmVw0MBElhW7dCbpGwMEBXVVWxcuVKOjo6qK+vn7W9cjzgxrX4h8+hao2cSNwoUMuSNF/Jz4mJCYaGhli2bBldXV10dHTQ0KCF17K5Gqp2PUwNGqqW8qJALUvOXCU/77tvkt/+7UHe+tYK2tvb6ezspKGh4bzlJsudy6FqH8U/fJW/DNsut5ETiTdlfcuSkq/k57JlkzQ1DdDQMMS+fW1s3LiZDRs2kEqlllSQdpFVnYSMatCcYikfCtSypOSW/KysnKK5OQjQ/f3NHD68mSeffAs/+tGKJRWgIb51qn0V/wCVv5TyoEAtS0pY8rOycoqmpgFWrDjDqVNNPP98Dy+8sJHTpxux1kQaZk2CONep9ln6Mmxf5S9lKVOgliWlsTFNY+MpGhvPcOZMI88/v4nnn9/IqVNNWHsuSsSl5Ge51KnWMLVI8RSoZUlIp9OcPn2aCy88TU1Nip/8ZCPPP7+JgYEWrD33NS9lmNW1cqxTfeQIfO978PDDwa2GqUXmp0AtiZbJZDhz5gynTp2irq6OTZsu4tZbN9Hf3zojQEO8Sn6W61C1Sl+KRKdALYmUyWQYHBxkYGCA2tpaLrroIjZv3kxbWxu/+quVsR5m1VC1iEShedSSKJlMhuHhYcbHx0mlUqxZs4aWlhaWLZv5VfZRuMJFlS5QnWoRiUaBWmJjrpKf1lqGh4cZGxujoaGB1atX09LSQlVVVd72XBaucLmghOpUi0gUCtQSC/lKfu7cabn00hFGR0epr6+nu7ubtra2OQO0a64XlFCdahGJQoFaFt3sJT8tIyOj7N07wkc+Use7372O9vZ2qqurC243rgtKaKhaRKJQMpksqvNLflqWLx+ltbUPYyyvvLKWe+/dQlfXhZGCdFyrdIGyqkUkGgVqWVTZJT9ra8dobe2jsjLNa6+t4dChLRw9uoqXXqqJFAjjPPUppKxqESmUhr5lUQ0MQE3NGA0Nw0xMVHP06CrefLOD0dHlM7YrNBC6Hqr2taAEnBuqfvxxOHMmqFO9Y4d6wSIyk3rUsmjGx8epqemjunqC119fyaFDW3j11bXnBWkoPBAmpUpXSHWqRWQ+CtSy4CYmJujv72d8fJyrruri9OktvPrqekZG6s7bNmogTFKVLhGRQihQy4IJA/To6CgdHR309PTwlrd087nP1c86VA1BDzhKIFSVLhFZanSNWrybnJxkcHCQiooK2tvb6ejoIJVKeVkPWlOfRGSpUaCWomQy0NsLVVXBbXYVsdDU1BSDg4MAtLa20tnZyYoVK2YE6DD5K5//v717D46rPO84/n0s64Kk9eqy9spFxrIdsC1ubtKGmjppUtqUCZ2BFFKYDE2ZBNrQwpSGtIXJkDYkrRNI0nbGdAoZLk0zU5jWGaeUTqDMlLtyoVwiLg7FQMPN2F5JllbX1e7TP86uvZa8sna1Kx3t/j4zZ6xz9j3nvOfxSs++Z9/zvsV2/tIoXSJSbXTrW4rW1wdXXgk33xys33xzsN7XF6xPT08zODjI8PAw7e3tbNmyhVNPPZVoNDqrFV2J55R1q1pEqola1FKU/FHEGhuPbk8k4JZb0lx33Qhnnpmhvb2deDxONBplxcymdp5KPqesW9UiUg2UqGXeZo8iFlixIs2qVUlWrpzmu99tY/fuLjo72+ZM0DmVfE5Zt6pFpBro1rfMW/4oYgBmGQDa2g4zNtbCK69s5sknN/Piix3zStJQ+eeURUSWOyVqmbfBweBfswyRyDAdHUMA7Nu3ib17t5BIdJLJ1BV1m1rPKYuIzE2JWuZt1aoMra0jdHYOMDXVyL59GwFIJGKk00cz6Zo1xR1Xnb9ERArTd9RyQu7O6Ogo4+MTTE9HePXVdQwMdFJfX2CUkhKo85eIyPEpUUtB7s7Y2Bjj4+O0tLTgvpGXX46RStUDUF+fOu5+Bw6Udj51/hIRmU2JWmbJT9DNzc1s2LCBWCzG0FADqePn5mOU0kNbRESOT4m6RmQyQa/twUFobz/+SGLuzvj4OGNjY5x00kn09PQQi8VozD4wXanhOUVEpDAl6hrQ1xc8/5z/aFUsBlddBdu3B+vj4+OMjo7S1NTE+vXricViNDU1HXOcmcNz5lMPbRGRylCv7yqXG0ksP0lDMJLY174GTzwxwaFDh0in06xbt47e3l66u7tnJekc9dAWEVlcalFXsUIjiQE0NEzS2prk3nsbuP32brq6VtPcPHs+6OPJ9dB+7DEYHoYHHoAPf1gtaRGRSlCLuorNHEkMggTd2ZmgsXGSd99dS1/fVn7+8/XzTtI5dXWwY0fw844dStIiIpWiFnUVGxg4+nN9/RSRyAjp9Er2749z8GCcZLIVCDqHiYhIOClRV7HDh2HlyhSRyAiZzAoOHlzDgQNrGBmJAEd7gx08uHR1FBGRuSlRV6np6Wmam0eIRGBgoJP33utiePjYBJ2zevXi109EROZHibrKTE9Pk0wmcXfi8Xb27o0zPBzFvcD0VMzuwS0iIuGhzmRVIp1OMzQ0xPDwMNFolM2bN3PRRacRibTNmaQ1haSISLipRb3MpdNpkskk6XSaaDRKV1cXbW1tR+aDzg1QUmgkMQ1QIiISbmpRh1QmA/39wbPK/f3Ber50Os3hw4cZGhqipaWF0047jS1bttDR0XEkScPRAUq6u4/df906DVAiIrIcqEUdQnMN+XnOORmSySSpVIpVq1bR1dVFe3s7dXM0izWFpIjI8qVEHTJ9fbBz5+ztiUSGXbtGmZycZPv2CD09PXR0dMyZoPNpCkkRkeVJiTpEMhnYtevYbWYZmpvHaGqaYHQ0wje/eQr9/R00Nuq/TkSkFuivfYj098PISG7NaWkZpalpgrGxVl57bROJRCfT0/U88QScd95S1lRERBaLEnWIvPACgNPcPEZz8zhjY828/voGEokYqVTDkXKPPKJELSJSK5SoQ8LdgXFisTHGx5t5/fUeEokYU1ONS101ERFZQkrUS8zdmZiYYHR0lM2bm7jrrvUcOrSaycnCCVqdwkREaoeeo15CExMTJBIJ0uk069at48ILe5mY6J4zSXd2KlGLiNQStaiXQK4F3dDQQHd3N6tXrz4yH/Qdd8DFFxfe94479PyziEgtUYt6EU1OTpJIJJiammLt2rVs3bqV9evXH0nSEAxOsnv37IkyuruD7RpJTESktqhFvQimpqZIJpOsXLmSrq4u1qxZQ2tra8HyGklMRERylKgrKJVKMTIyQl1dHatXryYej9Pa2opZ4dmscjSSmIiIgBJ1RaRSKZLJJACxWIx4PE4kEplXghYREcmnRF1G09PTjGSHFmtvb6erq4tVq1YpQYuISMmUqMsgl6DdnY6ODuLxONFoVAlaREQWTIl6AdLpNCMjI2QyGdra2ujq6iIajR4zH7SIiMhCKFGXIJ1Ok0wmSafTRKNR4vE47e3tStAiIlJ2StRFyGQyJJNJUqnUMQl6vnNCi4iIFEuJeh4ymQyjo6NMTU0RiUTo6emho6NDCVpERCpOiXoO+Qm6tbWVU045hY6ODlauVNhERGRxhOpLVTP7YzN7w8wmzOxHZvbBE5T/pJntzZbvN7OPl6Me7k4ymWRgYID6+no2bdrE1q1bWbNmjZK0iIgsqtAkajO7FPgW8GXg/cDzwINmtqZA+XOBfwHuBH4R2APsMbMzSq2DuzM6OkoikaCuro6NGzfS29tLPB6nvr6+1MOKiIiULDSJGvg88G13v9vdXwI+B4wBnylQ/k+AH7j7re7+srvfBDwDXFPsid2dsbExEokEZsaGDRvo7e1l7dq1StAiIrKkQnEf18wagA8AO3Pb3D1jZg8D2wvstp2gBZ7vQeCiAudoBPIneo4ADA4OMjU1RVNTE7FYjM7OThoaGo6MMCaFpVKpIx9w9IGmdIrjwimGC6cYLtzAwEBFjhuKRA3EgDrgvRnb3wO2FNinq0D5rgLlbwT+cubGSy65ZP61FBERObEOYLhcBwtLol4MOzm2BR4B3gK6ATWfS6MYlofiuHCK4cIphguXi2FZm9ZhSdSHgDQQn7E9DuwvsM/+Ysq7+yQwmVvPG4d7xN3L9smnliiG5aE4LpxiuHCK4cJVan6HUHQmc/cp4H+A83LbzGxFdr2vwG59+eWzfnOO8iIiIstOWFrUENyW/iczexr4MXAd0ALcDWBm3wHedvcbs+X/HnjUzK4HHgAuA34J+IPFrriIiEilhCZRu/t9ZrYauJmgQ9hzwPnunuswdgqQySv/lJl9Cvgq8DfA/wIXufsL8zzlJMEz25MnKigFKYbloTgunGK4cIrhwlUkhubu5TyeiIiIlFEovqMWERGR41OiFhERCTElahERkRBTohYREQmxqk7UYZk2czkrJoZmdrqZ7c6WdzO7bjHrGmZFxvEqM3vczAazy8Mneu/WgiJj+Dtm9rSZDZnZqJk9Z2a/t5j1DaNi/ybm7XdZ9nd6T6XrGHZFvg+vyMYtf5ko9pxVm6jDMG3mcldsDIFm4DXgBgqPKFdzSojjRwjeix8lmHzmTeAhMzu58rUNpxJiOAD8NUH8ziIYj+FuM/utRahuKJUQw9x+PcA3gMcrXMXQKzGGw8DavGV90Sd296pcgB8Bu/LWVwBvAzcUKH8f8B8ztv0Q+MelvpblEsMZ+74BXLfU1xCGZSFxzJavy/6yf3qpr2W5xjC7zzPAV5b6WpZTDLPvvSeBzwL3AHuW+jqWUwyBK4ChhZ63KlvUedNmPpzb5u6Z7Ppc02Y+PGPbg3OUr2olxlBmKFMcm4F6yjzQ/3Kx0Bha4DxgM/BYpeoZZguI4ZeAA+5+Z2VrGH4LiGGrmf2fmb1pZt83s9OLPXdVJmrmnjaz0DSYxU6bWe1KiaHMVo44fh14h9kfJGtFSTE0s6iZJYEpgmGGr3X3/6pYLcOt6Bia2Q6ClvRVla3aslHK+/BnwGeAC4HLCXLuU2bWXcyJQzOEqIjMZmY3EIxj/xF3L7oTSo0bAbYBrQQT+HzLzF5z90eWtFbLgJlFgH8GrnL3Q0tdn+XK3fvImyjKzJ4CXgb+ELhpvsep1kRd8Wkza0ApMZTZSo6jmX2BoGPeb7j7TytTvWWhpBhmb0u+ml19zsy2AjcCj1SgjmFXbAw3AT3A/XlTN64AMLNpYLO776tITcNrwX8T3T1lZs8C7yvmxFV569s1beaClRhDmaHUOJrZnxN84j7f3Z+udD3DrIzvxRVAY3lrtzyUEMO9wJkEdyRyy78D/539+c0KVzl0yvE+NLM6gri+W+zJq3IBLgUmgN8HtgK3A4NAPPv6d4CdeeXPBVLA9cAW4K8Ivts6Y6mvZRnFsIGjv9TvALdmf37fUl/LMovjXxDMvnMxwXdfuaV1qa9lGcXwRoIP2huz5a/P/n5fudTXslxieJz970G9vot9H34J+Fj2ffh+gscux4HeYs5brbe+8cWfNrPqFBtD4BeAZ/PWv5BdHiV4NrgmlRDHqwk+9PzbjEN9meADZM0pIYYtwD8A3QR/GPcCl7v7fYtX63ApIYYyQwkxbAe+nS07SNAiP9fdXyrmvJrmUkREJMSq8jtqERGRaqFELSIiEmJK1CIiIiGmRC0iIhJiStQiIiIhpkQtIiISYkrUIiIiIaZELSIiEmJK1CIiIiGmRC1SxczsG2a2p4T9Os3sgJn1lLEu95rZ9eU6nkitUKIWqW7bCMYjLtYXge+7+xsAZnajmf3EzEayCXyPmW0u8phfBb5oZtES6iNSs5SoRarb2RSZqM2sGfgscGfe5l8DbgN+hWBWqnrgITNrme9xsxPc7AMuL6Y+IrVOiVqkSplZNxAjm6jNrM3M7jezJ8ysa45dPw5MuvsPcxvc/Xx3v8fdX3T354ErCGYK+kDe+V42My+wXJMtdj9wWXmvVKS6KVGLVK9twJC7v2FmZwI/Ad4GPuru++fY70ME0/HNJXf7eiBv28XZf88D1gI9BFP+fZJgqj+AHwMfNLPG+V6ESK1TohapXtuA57PzrD8K3OLun3P31An2Ww+8U+hFM1sB/B3w5Iz52uPAdHb7foLW/AqGqllnAAACF0lEQVTgcXefzJZ5h2Cu7bla9CKSZ+VSV0BEKmYbcBawC7jA3fvmud9JwMQcr98GnAHsmLH9TOCVvKR8NnDA3d/LKzOe/bd5nnURqXlqUYtUr23A94AmoG3mi2b2pJmdk/35TjP70+xLh4D24x3QzHYBv01w+/ytGS+fBfTnrZ89Yx2gI/vvwSKuQ6SmKVGLVCEziwAbCVq/1wD3mtnpM4p9BbjBzD4PZNz9b7PbnwV6ZxzPskn6E8Cvu/vrxzntWcBP89bPnrEOQUv8LXc/VMJlidQkJWqR6nQ2kAZecve7CB61ut/MYrkC7v4Dgp7bFwB/lLfvg8DpZpbfqr6N4LGqTwEjZtaVXU6CI99bn86xiXkT8MaMen0IeGjhlydSO5SoRarTNmBv3vfFfwb8DPiemTUAmNkvE9yKPpzfwczd+4FngN/NO97VBD29HwHezVsuzb6+ieB75/xE3Q982cx+NXu+JuAijvYAF5F5MHdf6jqIyCIzs5OB/yRInLuBT+f34DazC4BbgTPcPVOmc14NfMLdP1aO44nUCrWoRWpM9nb1vwLXZr9r3gnclF/G3R8A7gBOLuOpU8C1ZTyeSE1Qi1pERCTE1KIWEREJMSVqERGREFOiFhERCTElahERkRBTohYREQkxJWoREZEQU6IWEREJMSVqERGREFOiFhERCTElahERkRBTohYREQkxJWoREZEQ+3+CA0VseFvYZwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "kx = [k.x for k in kpts]\n",
    "fig = plt.figure(dpi=100, figsize=(5, 5))\n",
    "ax = plt.subplot(111)\n",
    "for i in range(len(all_freqs)):\n",
    "    for ii in range(len(all_freqs[i])):\n",
    "        plt.scatter(kx[i], np.real(all_freqs[i][ii]), color=\"b\")\n",
    "\n",
    "ax.fill_between(kx, kx, 1.0, interpolate=True, color=\"gray\", alpha=0.3)\n",
    "plt.xlim(0, 0.5)\n",
    "plt.ylim(0, 1)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"$k_x(2\\pi)$\")\n",
    "plt.ylabel(\"$\\omega(2\\pi c)$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gray shaded region is the light cone, $\\omega>ck_x$, which is the region corresponding to modes that are extended in the air surrounding the waveguide. Below the light cone, we see several discrete guided bands, which must have field patterns localized to the vicinity of the waveguide. The imaginary part of $\\omega$ for bands below the light cone is very small, due to either numerical error or the finite computational cell size. Some tiny portion of the guided mode overlaps the PML. Note the band gap between the first and second guided mode, from about 0.2 to 0.3.\n",
    "\n",
    "Inside the light cone, we also see several discrete bands. These are leaky modes, or resonances, which have some intrinsic lifetime/loss because they couple with radiating states inside the light cone, which is reflected in the imaginary parts of their $\\omega$. Twice the imaginary part of $\\omega$ is the energy loss rate per unit time; for a waveguide, it is more conventional to report loss per unit distance; to get this you would divide the loss per unit time by the group velocity $|d\\omega/dk_x| = |slope|$. Harminv only identifies leaky modes that have a substantial lifetime. The default threshold is a lifetime, or $Q$, of 50 periods.\n",
    "\n",
    "Computing band diagrams, especially for leaky modes, with a time-domain program like Meep involves several subtleties. For example, the accuracy of Harminv will go down if we specify too large a `df` (too narrow a source), because the increased number of modes makes the signal-processing more ill-conditioned. Sometimes, Harminv will report a spurious mode, which will appear as an isolated dot on the plot. Second, we sometimes have to be careful with modes and especially the imaginary parts to make sure they aren't an artifact of the cell being too small, or the signal-processing error being too large (either because the run is too short or because the bandwidth being searched is too large). Third, currently Meep doesn't attempt to \"connect the dots\" for the bands — the frequencies are printed in increasing order, but since modes disappear when their losses become too large this means that a single band may be split across several columns.\n",
    "\n",
    "For example, there seem to be some bands that run right along the edge of the light cone. These are not leaky modes, but are artifacts of the fact that PML boundaries do not absorb well for light that is travelling parallel to the boundary, corresponding to extended modes at the boundary of the light cone. Below, we will see that these modes are not localized to the waveguide."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is usually a good idea to examine the field patterns for any modes that you are particularly interested in. That is, re-run the simulation with a narrow-bandwidth source, at the particular ω and k you want, and output the field patterns just as we did for the resonant-cavity modes. We'll do that for several modes below:\n",
    "\n",
    "* $k_x=0.4, \\omega=0.1896$ guided mode\n",
    "* $k_x=0.4, \\omega=0.3175$ guided mode\n",
    "* $k_x=0.1, \\omega=0.4811−0.0017i$ leaky mode\n",
    "* $k_x=0.3, \\omega=0.8838−0.0018i$ leaky mode\n",
    "* $k_x=0.25, \\omega=0.2506$ light-cone (extended) mode\n",
    "\n",
    "First, we'll wrap the routine in a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_sim(kx, omega, filename):\n",
    "    s = [\n",
    "        mp.Source(\n",
    "            src=mp.GaussianSource(omega, fwidth=0.01),\n",
    "            component=mp.Hz,\n",
    "            center=mp.Vector3(0.1234, 0),\n",
    "        )\n",
    "    ]\n",
    "    sim = mp.Simulation(\n",
    "        cell_size=cell,\n",
    "        geometry=geometry,\n",
    "        boundary_layers=pml_layers,\n",
    "        sources=s,\n",
    "        symmetries=sym,\n",
    "        k_point=mp.Vector3(kx),\n",
    "        resolution=resolution,\n",
    "    )\n",
    "    f = plt.figure(dpi=100)\n",
    "    animate = mp.Animate2D(fields=mp.Hz, f=f, normalize=True, realtime=False)\n",
    "    sim.run(mp.at_every(5, animate), until_after_sources=1)\n",
    "    animate.to_mp4(10, filename)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,1.2,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.36, height 1e+20, axis (0, 0, 1)\n",
      "Meep: using complex fields.\n",
      "Normalizing field data...\n",
      "run 0 finished at t = 1001.0 (40040 timesteps)\n",
      "Generating MP4...\n",
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,1.2,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.36, height 1e+20, axis (0, 0, 1)\n",
      "Meep: using complex fields.\n",
      "Normalizing field data...\n",
      "run 0 finished at t = 1001.0 (40040 timesteps)\n",
      "Generating MP4...\n",
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,1.2,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.36, height 1e+20, axis (0, 0, 1)\n",
      "Meep: using complex fields.\n",
      "Meep progress: 929.575/1001.0 = 92.9% done in 4.0s, 0.3s to go\n",
      "Normalizing field data...\n",
      "run 0 finished at t = 1001.0 (40040 timesteps)\n",
      "Generating MP4...\n",
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,1.2,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.36, height 1e+20, axis (0, 0, 1)\n",
      "Meep: using complex fields.\n",
      "Normalizing field data...\n",
      "run 0 finished at t = 1001.0 (40040 timesteps)\n",
      "Generating MP4...\n",
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,1.2,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     cylinder, center = (0,0,0)\n",
      "          radius 0.36, height 1e+20, axis (0, 0, 1)\n",
      "Meep: using complex fields.\n",
      "Normalizing field data...\n",
      "run 0 finished at t = 1001.0 (40040 timesteps)\n",
      "Generating MP4...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAAFyCAYAAAD25PTWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZPUlEQVR4nO2de7AcVZ2Av1/P6859kpuYYMJCgMgjibwDC4goLG4spXyVhbpr+dgSpESL2tISWF3j1i5BqsQXiizlmjXi6lLlYpVosHioxWNZeYWQgJrEBCOQkHtvbu5rnv3bP7pn0rczc+/MzfTMPenzVfW9M919zpzpb36nz5w+fUZUFcv8xul0ASyzYyUZgJVkAFaSAVhJBmAlGYCVZABWkgFYSQZgJRmAUZJEZJmI/FBEhkRkSkS2iMh5nS5X1CQ7XYBGEZEFwKPAw8DbgdeANwAjnSxXOxBTOlhF5BbgYlW9pNNlaTcmSdoG3A8cB1wK/AX4jqreNUOaDJAJrR4EhqMqp08f8LK26uCqqhELkPOXm4GzgauBKeAjM6RZB2iHlmWteu8mRVIBeFJVLwqs+yawRlUvrJMmHEl9wJ7bfnEbZx5/ZkvLt3HLRjZs3sDa49ay6bpNAAOqerAVeRvTcABeAbaF1r0AvK9eAlXNA/nKcxEB4Mzjz+SSN7ypZQW7+dH1bPj9Bta9bR2XLL6ETWxqWd5gVhP8UeDU0LpTgN2tfJFm67SbH13Put+uY92b13HjxTe2sihVTIqkrwGPichNwH8D5+Odl66eS2atqOTXBwTdFJEgMCiSVPV3wHuADwLPA18ErlfVu+eSnxzh0i5BYFYkoao/B37e6XKEq7iom15GSYqCZg/w+jacg8LEVtJcPv2dEAQxltQs7TwHhbGSfGSGbTd3UBDEWNJMUoJ0WhAY1ATvBM0KiqqVF9tImo3ZBLWzx9NKqkEtQZ3shraSQrT7i2oj2HNSgHZ0ls6F2EdSJVI69UW1EWIrKViNzWdBYKu7eS8IYi7JBEEQY0mmCIKYnpM2btnojUkwQBDENJI2bDZHEMQ0kq58w5WsPXktz7z6TMvz/sPQH1qepzHj7lqBiPQDo9wAdEX0IjngFqCF4+5iKenB5x+kp68nkteYGJvg8tWXQ0wHR7aMM5acQX9/fyR5H+xuiZdpxLLhYBpWkgFYSQZgJRmAlWQAVpIBWEkGYCUZgJVkAFaSAVhJBmAlGYCVZABWkgFYSQZgJRmAlWQAVpIBGCtJRG4QERWRr3e6LFFjpCQRWQNcAzzX6bK0A+MkiUgvcDfwCWIwtScYKAn4NnCfqj4w244ikhGR/sqCN9+dcRg1pEtEPgCcA6xpMMmNwJeiK1F7MCaSROSvgG8Af6equQaTrQcGAstxERUvUkyKpHOBxcDTlRkggQTwZhG5DsioajmYoN7MkaZhkqQHgTeG1n0feBH4SljQ0YQxklR1DG8ywioiMgEMqerztVMdHRhzToozxkRSLVT1LZ0uQzuwkWQAVpIBWEkGYCUZgJVkAFaSAVhJBmAlGYCVZABWkgFYSQZgJRmAlWQAVpIBWEkGYCUZgJVkAFaSAVhJBmAlGYCVZABWkgFYSQZgJRmAlWQAVpIBWEkGYCUZgJVkAFaSAVhJBmAlGYCVZABWkgFYSQZgJRmAlWQAVpIBWEkGYCUZgDGSRORGEfmdiIyJyD4RuVdETu10udqBMZKAS/EmJPxr4AogBfxKRKL5wdh5hDHT1qjq2uBzEfkosA9virXf1kojIhkgE1hl5MyRJkVSmAH///AM+9wIjAaWPVEXKgqMlCQiDvB14NFZplGzM0d2kG8Dq4E3zbSTnTmyQ4jI7cA7gTerqpHVV7MYI0m8MPgW8B7gLar6p6hfU8NliPoF62CMJLwq7kPAu4AxETnWXz+qqlOtfKGwnPD6dssyqeFwLd7J/9fAK4Hlqla9gFJfUHi/dmJMJKlqpB/gZg+80r6IMimSImOukdGuiLKSjpB2iIq9pFYc5KhFxVpSuxsAc8WYhkO7qSewXmMhyoaElRRitujqxHel2FZ3tWQ0U/21s6qMraQwczno4TRRibPVHTN0A4U21OpEb8eX2lhK0jqP4XAxtba1+4qHre4CzCRopv3s96QImRZR9VoSwWW2/Q/frSXEWlKFww54ve7w0PpguiijyUoKModrFY1WkUdCLBsOcOg4Vw9y4GC7M6SrfqoDzTrVaBsTNpKgYUGV7dV92tSAiLWkcFU1m6Ca+9rqrg34B7mWIA1ZDA8Jczn0KY+yyoulJGV6FAUFhcVQY1tQlgs4EXc7xFJSFa0taNYaTLWtAy3jLSlAWNBM15PUS1AVVYmmRlvwzRJrSZUoCgo6rGkeJhhAbYqo+Eqq0XwOnqvqn5tk2r+KqGAjotXEUtK07hzVaYJmPS9VE8thvqJqjsdSEgSqOqYLmnZeqTEYvNqQU++RiJcuymovlpK8aFFc15Pi+id9N2BImR5R4j8RBBVwfDmq3mMHcAVct/XhFEtJ3qfejwL8Fptq9cuoF2US+HsIxRNSSScigS+xEklExVJSLpcjlUp7/XAKrqofTVqNplrxUDn8jucYQbwoEqlGUz6fr5HyyIilpEcefYx0V5ayq5RcpayerJLr+o+9dbW6hRICCUdIiPj/8R97/wu5lt6FA8RU0q6XXkKSGQquMlUskyu5TBVdCmWXYtmlUFbKrksxdH5JOULCccgkHdIJIZ3wHmeTCbIpb51btJHUEvKlMlBmNFfkYL7MwVyJsVyRsVyJqUKZfLFMqeSirk7rr3MSQiLhkE0nyKYT9HWlGOhO0ZtOMJBJ0ZtJ4JSa6UtvjFhKcl2YKpQZzZV4baLA8Li3FHJFivkypUKZUqGE65ZRtwSAOEmcZJJkMsF4Jkkq4zDclWKwN83r+rypIpKOkIngy1IsJZVVmSx5kvaN5hgZyzM1XiA/UaAwdZDi1Djl/BSUi4c6FgASKVJdPSS6ukln+8n0pCn5kZNKOnSnEiQTrS9vLCW5ClPFMgcmChwYLzAxmmNydJTC+AiF8VFEYOHChSxevJje3l4AxsbGee21fQwNDcHYMMXeYyiXFgKwP+HQk0l6VV4EfUOxlFR2lVzJZXSqSG6ywNTBMXIj+0hSZtXKlSxdupTunm6y2SypVAqAYrHI1NQUExMTvPLyK+zYsYOp4SIiDql0gtHJFPl+FzeCu0ZjKalYVnLFMhO5EvnJEvmxIfq606xYsYKTTz6ZgYGBmukq6xctWkR3Tzfbt29nYvQ10tksY70lcqUyrtv6QxpLSYoyWShTyJfITxwk6RZYvfpcVqxYQTqdxnXdaq/EYSlVGRgYYPXq1WQyGZ5++mnykxMU8l3kim61O6mVxFJS2YViyaVULENhihNOWM6JJ55IKpVCVXGceicWr9tHVUml0px00kkMDQ2xZ/8YxbzXdHe19Yc0lqOFXFXyJZdSoUw2qaxauZJ0OuP3w81+Tqn016XTGVauXEk24VIqlSmU3EgGSzYsSUSWtv7lm0dEPiUiu0QkJyJPiMj5c8mnUHIp5kr09fYwODiI40hTB1gVHEdYOLiQnmzGi6SSSwSd4E1Vd1tF5FOq+qPWF6MxROQq4Dbgk8ATwPXA/SJyqqruazQfRSmUXMQtsGDRIOLMHD1nbX2Sdzx8LwtH9jO0YBH3vfXdPLvqPCoX/gYXDDI+kvciKYJzUjPV3T8Bd4rIPSIy2PKSNMY/Anep6vdVdRuerEng481lI5RcJZtKMLhgwYx7nrX1Sa7+8e0s27uHrkKOZXv3cPWPb+fMbU9V91kwuIBMQii70YztaliSqn4HOANYCGwTkStbXpoZEJE03lSeDwTK5PrPL6yTJiMi/ZUFf3pP70qskkomyHZnZ6zm3vHwvTXXv/Oh/6mUge7ublJJJ5ILftBk686fvuwyEbkO+KmIvACUQvuc08LyBVkEJIC9ofV7gdPqpLkR+NKRvOjCkf1NrY+CptuLInIC8F5gBPgZIUnzjPV457AKfcCeyhXVYqnM1OTUjMODhxYsYtnew+c+HFqwCPx8JicnKZTK9M5ybpsrTUkSkU8AX8WrYlap6muRlKo2+4EysCS0fgnwaq0E9af3VJKOMFUsMzwyzEmcVPdF73vru7n6x7cftv7nl72n+nhkeIRC2bsYGMWQoWaa4JuArwDXqep72ywIVS0ATwGXB8rk+M8fbyYvQUgnHTSZYWRkpHrdqNa56dlV53HnBz/NX5YcRy7dxV+WHMedH/w0z55+rpfGVYZHhnGdBOmkg0TQcGgmkhLAGR2e9/Q24D9F5Eng//Ca4D3A95vNKJ30OkbHRiYYGh5iwYJBEonaB3jzynPZvPLcaesEKJeVkZFhJqbypAYTZJIOUdR4zbTuruj0xLSq+hPgs8C/AM8CZwFrVTXcmJgRR4RM0iGZTjBVErZt20ahkK8bTTXKgSoU8nm2bd3GZElIJv1I6qSk+YKq3q6qJ6hqRlUvUNUnms0j4XgX6ZKpBKSz7N69m507d1IsFhEB13WrIqYvWu18LRYL7Ni5k927d5Po6iWVSZBJJXDskK7WIAjd6QTpTJJMTz9j4yNs3bqVfC7PySvqX6qodLCOjo6yfft2dmzfgZvqJp3Nks4k6Uo5fuOhtcRSUiohdKUS9HQlmexOUuhbyNjwqzy35TkmJyd5/dLX09PTQ7YrSzLlHaJSscRUzrvo9/LLL7Njxw5cJ033wsV09aTo60rSFcW1c2IqKeEIXUmHgWyK8Z405WI/qi6F8RG2btvGthe2MTg4OO3y+fj4OPv27WN4eBgQUr3HkO1fSLa/h0x3moHuFJlENA2HWEpyBLKpBMf0pJkslHFdRZwF5DM9pLr7KU6Nc2BskgMH/si00eCJFJm+hSSzvaS6+8h2p+nqTbOoz5OUTTmI2CFdLSEhQncqwUBXkkJ/BscRhpMJ0l1JivkMpcIgpVIZt1RC3RIiDohTHdKV8od0ZQJDuvozSbLJBImm7mFvjFhKchzoSSYouUlEhK5kgr6uJBP5EhP+FdZy2cUtNzc4sivl4JRsw6ElZJIJJJnA6XboSpXp70qRL2XIlaYPM66MCwfvu1VCvKZ7yjk0zDib8oYZV4YeuxF8q4mlpOXHH19zwH5Zvee1BuyLCA7eXA12wH4buOTii+nt759+y4tL4Iay2rdk+nc1IVJ57FWdwVtgxscOtry8sZSU6cqQyWS8HgS83gTXvxUzPLdDLRyo3jJbuTepevmjkJk58RyIpaTK7ZjVm5j9kBG8SPIuONS4HZPA7ZiBvJTKZRCdVkW2ilhKqgx8dBzBVfXvf8Ufc1e51VJq13ccquokFEUO3giiVhNLSeDfiExg2Ij4f1QPzXpS43hXIypwr2zUE27EUlLwmIrI9PouONNgzbRSzeOwgch2lq4WE7rSLfUfTE9Wo8qDaK/5xFcSgSrPH9/dSFRI4H+7ZuqKtaQgQVGzDXEMC6o2ySWaGi/eksSbAq3yvahmRNVK1uYp9mMpyYuEQ22DSrUH0wXMNr1nJW3UP1YRS0nT8Ou3oKjqplkiJthYsFNOR0T4wDZzMKr7tqHmi7WkKoEDPdsBcagvKCpfsa3uqq24yrkp8L2poU/utC/ErS7ddGwkBal1L3O9/SoPbXXXHg470PVkhdZP615qfbGqxFpS8MDWjAgJLbPtf/huLSHWksI0WnWF94u6xotlw0FCj4NfWYMCGvnhxXYQS0lh6t361YiUdniz1Z3PXA52rfZGFMRWUr12wpGkjwpb3YWoHPxGfqqiXVhJdWhWhv2eFBEdaqw1TawlQWtERS079pKOFNsEbxNzPdDtqi6tJJ/51FAIYyUFmMOVirZghCQRWS4i3xORP4nIlIjsEJEv+9Ortf71aOhKRdsw5XvSaXgfqGuA7cBq4C68KWs+G9WLzpcmuhGSVHUTsCmwaqeInApcS4SS5gtGSKrDADA80w4ikgGCd3X1RVqiiDDinBRGRFYAnwbunGXXG4HRwNLRCazmSkclicgtIqKzLKeF0izDq/ruUdW7ZnmJ9XgRV1mOi+SNREynq7uvAhtm2Wdn5YE/N/nDwGPA1bNlXn/mSLPoqCR/9smGZqD0I+hhvNkjP+bPZBwLOh1JDeEL+jWwG68197pKVKhqzflXjyaMkARcAazwl/DJ38w6rAmMaN2p6gZVlVpLp8vWDoyQFHesJAOwkgzASjIAK8kArCQDsJIMwEoyACvJAKwkA7CSDMBKMgAryQCsJAOwkgzASjIAK8kArCQDsJIMwEoyACvJAKwkA7CSDMBKMgAryQCsJAOwkgzASjIAK8kArCQDsJIMwEoyACvJAKwkA7CSDMBKMgAryQCsJAOwkgzASjIA4ySJSEZEnvVn8Dqr0+VpB8ZJAm4FXu50IdqJKXMLASAibwfeBrwPeHsD+9uZI9uJiCzBmxz3w8Bkg8nszJHtQrx50zYA31XVJ5tIameOPFJE5Bbg87PsdjpeFdeHd9Abxs4c2Roand7zMuBCIB860E+KyN2q+pFoijc/MGJ6TxH5DPCFwKqlwP3AVcAT0ZRu/tDpSGoIVX0p+FxExv2HO1TVyMZAMxjRcIg7RkRSGFXdRQzmXq1gI8kArCQDsJIMwEoyACvJAKwkA7CSDMBKMgAryQCsJAOwkgzASjIAK8kArCQDsJIMwEoyACvJAIy8MnukPLf3OXomeyLJe2JsouV5xlLS5T+8HLoiyjzX+ixFVVuf6zxFRPqB0Tt+cwdnLz971v0nChNc/6vr2XlgJ9/822+y6nWrZk3zzK5nuPbSawEGVPXgEReamEbSKQtP4exjZ5Y0lh/jyp9cye6Du7n/Q/ezZumahvKOorqzDYcaVARt3b+VX3zgFw0LigorKcR8EwRW0jTmoyCwkqrMV0EQY0nCoSGw81kQxFhShfF5Lghi2gSvMJYf452+oF+GBM2nb4+xlTSTIDj8boBOSoulpMnCxIyCalHrFo52iYulpM8/dAN7insaFlSPdkVbLBsOu0Z3HbGgWkR1w1QsI+nWy2/lvKVrpn3y5/MdabGUdPqi0w9b12hV1QmZsZR0JHSilRfLcxJM73GY78Q+ko5ElG2CG0C7IjG2klodBVEKM+qcJCLvEJEnRGRKREZE5N455zXDMheU6Ko/YyJJRN6HN9/dTcBDeGVfHclrNbl/UM7GLRtbWRTAEEkikgS+AXxOVb8X2LRtrnm24lMvof83P7qeDZs3tCDn6RghCTgHWAa4IvIMcCzwLJ605+slqje95+aXNre8gBu3bGTD5g2sPW4tm9jU2sxVdd4vwAfwPvy78eZfPRf4EbAfGJwh3ToOnS7avSxv2fvv8MG/pYE3exrwIf/x1YG0Gby58q6ZIf8M0B9Ylvn5LAutr7c0u38wTX+rjlOnq7tGZ458vf+4eg5S1byI7ASOr5dQ60/vOdbI6NJm9w+laRmmzBz5FN7BPhV4xF+XApbjVYFHNZ2OpIZQ1YMi8l3gyyLyZzwxn/M339O5krUHIyT5fA4oARuBLN7cq5ep6kgTeeSBLxOoAlu8/1zTzEis7qowFaO6heKKlWQAVpIBWEkGcFRLEpFBEblbRA6KyAER+Z6I9Ib2+ZSI7BKRnH8Z5Cn/B7SCy3cD+79fRF7117si8qKInD9DGT5aI7/m7qztdL9cxN1Ov8TriL0AeBPwR+BHge1X4TWVPwasBP4dKAI/wOvErSz9/v4XAWW8rwI3AXf4z0eBxXXK8FF/ezC/JU29j04fyAgFnY7Xh3ZeYN1awAWW+s+fAG4PbHd8aY/UyfMnwEgozf8CE8ANM0g6cCTv5Wiu7i7EOzjB31t6AE/SBSKSxutNf6CyUVVdPAlrRGS/iDwvIutFpDuQZ38wDd4Pm+T9bfXoFZHdIvJnEfmZiMx+G3uAo1nSscC+4ApVLQHD/rZFQALYG0q3GfgT8Fa832v6MPDDQJ5OKM1evN72Y+uU4/fAx4F3AX/vp39MRBr+wS2TuoWApn4Ya65sBo5R1S3AFhF5BXhQRE6eS2aq+jjweOW5iDwGvABcA3yxkTyMk0TjlzdeBRYHV/qX4Qf9bfvxTvpLQmmX+NsrVH6faYW/flkozRK86i6Ypi6qWvSvLq9oZH8wsLpT1ddU9cVZlgLep/cYETk3kPwyvPf8hL/PU8DllY0i4vjPHw+kqfyW7Sv++oPBNMAVQDqUpi4ikgDe6OfX8Js+ahe8JvjTwPnAxcAfmN4E/yReQ2IdXhX5X8CUf+CX441KGgF+4+9/EV7zuwjcANzOoSb4En+fHwDrA6/xz3i/SXgS3liNymusbPh9dPpARixpEG8sxJh/IP8D6A1sX47XTH8Vr8p6xpc6hDeV0ySwhcClcOD9eI0F9QW/CFwQ2P5rYEPg+dfwrn9VqsT7gLObeR/2UoUBGHdOiiNWkgFYSQZgJRmAlWQAVpIBWEkGYCUZgJVkAFZSCBFJiMhjIvLT0PoB/6Ldv7W9TLZb6HBE5BS8sRGfUNW7/XU/AM4E1qjXg96+8lhJtRGRz+D1jq/C60W/B09Q628TnK0sVlJtxLvR6CG8SxFvBL6lqv/akbJYSfURkdPwLnVvAc5Rb4xE27ENh5n5ON41pROBhgeOtBobSXUQkYuA3+BdVf2Cv/pvtAMHzEZSDfxxdhuAO1T1YeAf8BoPn+xEeayk2qzHm0PjBgBV3QV8FrhVRJa3uzC2ugshIpcCDwJvUdVHQtvuxxsG19Zqz0oyAFvdGYCVZABWkgFYSQZgJRmAlWQAVpIBWEkGYCUZgJVkAFaSAfw/Fgv+1U5wAuUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAAFyCAYAAAD25PTWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAVzklEQVR4nO2de5BcVZ3HP9/ueWTymIQJJAhsHhB5hYegAUFeC4sbS9FSy4JldX1sKVKyFrWlK2F1DVu7BtkSXyi6lGuWgFsWVZasisECQZfHsoIYEgJKEhJEIJBkkgwzk5nM9G//uLcndzrdM91Jd9/+5Z5P1a1Mn3vP6dP3k9+5p8+597TMjEBrk0u7AoHJCZIcECQ5IEhyQJDkgCDJAUGSA4IkBwRJDgiSHOBKkqSjJd0habukQUlrJb0l7Xo1mra0K1Atkg4DHgYeAN4BvAa8EehNs17NQF4GWCXdCLzNzM5Puy7NxpOk9cC9wDHAhcCfgG+b2W0T5OkEOkuSe4AdjapnzAzgJavXyTUzFxuwJ96+BJwBfAIYBD48QZ7lgKW0HV2vz+4pkoaBx83s3ETaN4AlZnZOhTylkTQDePHme27m9Hmn17V+q9auYuWalSw9Zimrr1kNMNPMdtejbDcdB+BlYH1J2jPA+ytlMLMhYKj4WhIAp887nfPeeF7dKrbi4RWs/P1Klr99OefPOZ/VrK5b2eCrC/4wcEJJ2vHAlnq9gQ5gW/HwCpb/ejnLL1jO9W9bVq+qjMOTpK8Cb5V0vaRFkq4kui5960AKK3fCa+VLTRAEjiSZ2W+A9wJ/BawDvgBca2Z3plGfZgkCX9ckzOynwE/TrkczBYGjSGoVmi0InEVS2pQKataXlyBpEooikr24ZU0UBEHSOCqd+FJBzSZIorIcSF8QZFzSZE1WKwiCDPfuvAiCDEuaiFYSBEHSfrSaIAiSxtGKgiBIGqNVBUGQBLS2IAiSWl4QZFySB0GQYUleBEFGRxxWrV0V3ZPgQBBkNJJWrvEjCDIaSZe98TKWHreUJ195su5l/2H7H+peppv77uqBpG5gF9cBUxr0JnuAG4E63neXSUn3r7ufaTOmNeQ9+vv6ueSUSyCjN0fWjdPmnkZ3d3dDyu6bWhcv48hkx8EbQZIDgiQHBEkOCJIcECQ5IEhyQJDkgCDJAUGSA4IkBwRJDgiSHBAkOSBIckCQ5IAgyQFBkgPcSpJ0nSST9LW069JoXEqStAS4Cngq7bo0A3eSJE0H7gQ+TgaW9gSHkogWfPqZmd032YGSOiV1Fzei9e7c4eqWLklXAGcCS6rMsgz4YuNq1BzcRJKkPwO+Dvy1me2pMtsKYGZiO6ZB1WsoniLpzcAc4LfFFSCBPHCBpGuATjMbTWaotHKkNzxJuh84tSTt+8CzwJdLBR1KuJFkZn1EixGOIakf2G5m68rnOjRwc03KMm4iqRxmdlHadWgGIZIcECQ5IEhyQJDkgCDJAUGSA4IkBwRJDgiSHBAkOSBIckCQ5IAgyQFBkgOCJAcESQ4IkhwQJDkgSHJAkOSAIMkBQZIDgiQHBEkOCJIcECQ5IEhyQJDkgCDJAUGSA4IkBwRJDgiSHBAk1ZFGPdseJDkgSHJAkOSAIMkBQZIDgiQHuJEkaZmk30jqk/SqpB9LOiHtejUDN5KAC4kWJHwrcCnQDvxCUmN+MLaFcLNsjZktTb6W9BHgVaIl1n5dLo+kTqAzkeRy5UhPkVTKzPjfHRMcswzYldhebHSlGoFLSZJywNeAhydZRi2sHJki3wJOAc6b6KCwcmRKSLoFeBdwgZm5bL5qxY0kRWHwTeC9wEVm9nzKVWoabiQRNXFXAu8B+iQdGafvMrPB9KrVeDx1HK4muvg/CLyc2C5PsU5NwU0kmZnPq34d8BRJDcfSrkAF3ERSoygVU3zdSmGbWUmTRU0rycpkc1dLs9YKTWAmJXkjSKqCtKMpSHJAkFQlaUZTkOSAIKkG0oqmIMkBQZIDwohDCQc6wtDIkYlMRtJE15a0vxOVI5OSJqPVRAVJFagkKg2BQZIDgqQJaJVmL5O9O7NoS3Iwt+Q1es4pRFJMOXHlSGMSMEiahFZo8oKkEqqJpiTNiKxMXpMwxodIK9zIMAEhkmA/adVGU7PcZlJSId5qJa2Ay2ZzF4dOUZSSyaVplO88qEK62YHon5iMShJJDUaiSdH4f5OyKv29D2vIM1CZlDS0Zw9DHR37EkqESCp1NfGX3cRFbM+eoQkOPDAyKenB/3mIadOih9alfXElKYoojU9n7N/JexcDAwN1r28mJW3a8gJTuroQkJPIKflv/DeJ6DEDCmCGinLKSZIYCJFUHwaHC9A2Sl5CgrYc5CVyEm2xHQmsEMlRoRBLsbGOgUocmQDlsNGRutc3k5KGRkfJjxRoy+Voy0PBks2cEhE0GkWOFcAKsZj4S1VJJIlcdFwDBpIyKWl4tEB7wZAK5CwXnde4qcOMnAQ2Sh5QTsiEFUShMMKYhP2aOwPlY1H1JZOSCgajBcNy+7psxY5CXkIYs2fNYs4RhzNj+nSwAn19u3nt1a1s3749MWRe+oBM/QVBZiXti4KoVxd1FqZ2TeGE4xZy1Nw5TO2awtSuKbS3Rado795hBgf66e8f4OWX/sTGjRsY2rMnLiUORfa/VtWDTErKxZ2EfE7kJdpyYk7PYZx4/HEsWjCf7hnTAdC4a08XM7u7ATj88NlMm9rFhueeo3fnzmi3WWMMkVlJ0JaLJOVyMGPaVN502iksWjiPjvZ2rDAajRwUz3l8nTEzzIyZM2ex+JRT6eycwm+feJyBwcG4vy4a8fh1JgdY23I52vNRFHV1tLPo2IUsnB8JwgwpFzVeyeEGKfqym8thVqCjo4OFxx3H/AULaG9vb2h9MympPSfacjnyOZh92CxOPfkEOjvax2RMRvGIzo4OTlp8Ct3d8YJhEqj+p7TqEiUdVfd3PwAkfUrSZkl7JD0m6axay+hoiwS1SfTMmsnsw2aRy+Vq/oqjXI7ZPT3M6J451tw1glq0Py3pyobUokokXQ7cDNwAnAmsAe6VNKeWcorDP1O7Ojlidk90/ZngHOvun5B/63nkj5pH/pzz0X//NEqP6kRPz2F0dBTXPqy/qFok/SPwXUl3Seqpe02q4++B28zs+2a2HvgkMAB8rJZCFI/TdU+fTs+sWRMfe/dPyH/ww+jp9ej1fvT0evIf+ij6yc/GjjmspycesFVDgqlqSWb2beA0YDawXtJl9a9OZSR1EC3leV+iToX49TkV8nRK6i5uxMt7RkEjOjs76JraFY/0lG/rciturJB+U1QHYGpXFx3JqY86U1MXPF6+7GJJ1wA/kvQMMFJyzJl1rF+Sw4E8sLUkfStwYoU8y4AvHtS7bnmhQvqWgyq2Fmr+niRpPvA+oBe4mxJJLcYKomtYkRnAi/F4NkNDwwwODDLW3y4XTfPnwdPry6TPJ87FwOAgw8P1n6IoUpMkSR8HvkLUxCw2s9caUqvybANGgbkl6XOBV8plqLS8p5lRMNj9+uts7+3l2IXzKr5pYdl15D/4N2XS/2Hs794dO+jvHyAaHa/241RPLV3w1cCXgWvM7H1NFoSZDQNPAJck6pSLXz9aS1kFi7aBwSG27ejFzPa/F6/4vu+5jNE7b8cWn4xNn4YtPpnRO1Zil70zvpQZO3b0JiIp3amKPHBayuue3gz8p6THgf8DrgWmAd+vpZDhkQJTCjAiY8fOXWzfsZOeWTPJ5VRe1Lvfxei737WvOSwOExUK7OjtpW/3roaO3dXSu7s07YVpzeyHwGeAfwZ+B7wJWGpmpZ2JCdlbMEYKBUYLsL13J0+tf5ah4b1V3xVZPGpoaIj1T69j167EIGsD5pPcDQuZ2S1mNt/MOs3sbDN7rNYyRgoFhkcLjJoxOLyXjc9vZtOWFxjeuzd+jwJjk7BRAsSDq4VCASnH8PAwmzZuZMvmzYyMNLbvlMlR8GjSL5r4K0j09Q+w5ql1DA0OsmjhvqkKFF+rlANsbIRh185eNm7YwIYNzzE4GK/RGzd3YT6pThTMKJgxWjBGc8ZIAV7d0cvgU+sYGOjnDXOPYFpXF1OnTKGtPTpFI8NDDA4O0N8/wEsv/YlNGzYwNLRnX6ENHLvLpKRcYqS7ONhQMGNgzx7WrX+W9c88w+yZ3RxxxGxmTJsGZrz+eh+vbt3Kju3b4ozlH8toxHxSRiVBPqdxsxLFXvioRbcKb9u5k97e3n2zs1agMDpSMpWR+LtB0xSQUUkd+RztOZFXLrpDKD7XBYvkFYCc8ozaSHSNiTsOUp59X6hKH7rNRQU1QFQmJXW156M5pfgeu1xirs/Moi6CAOUxFVABoqnxYhdb+119ijdHNuK6lElJnW05prTnJ7zNeOyck4fc5LcZR8N/QvnRutc3k5KOnT/vAG/Yh8o3R0aEG/brxEUXnE93fHsWlBGRaP72l1SuQdsnbPfu3fWq5hiZlNTZ2UlnZ+e4tNKHyEolTfT3Pmy/cutBJiWV9s7GPVpp+9Jgvz7c5CXXusZAFWRUUnS12a+zXOkrUJnd5V4DKM1bug4lchzYBw8L5aZJybBbtc8mN0taNpu7xo2FNoQQSSXU+oR/M6IpSJqEVgi4bDZ3ZajlOtRscZmUVOXDE1XTaHGhuZuAVmjqIEhyQZBUgUpRlEZ0BUllaJVmrkgmJU0kodUEQUZ7d1B/GY3s4WUykrwRJDkgSKqBtK5XQZIDgqQqSbPXFyQ5IEiqgrS/OwVJDsikpFoiI+0ogjDiUPffm20EmZVUJCkjjVnXashkc1eJVhQEQZILgiQHuJAkaYGk70l6XtKgpI2SboiXVzvk8dJxOJHoP9RVwAbgFOA2oiVrPpNivZqCC0lmthpYnUjaJOkE4GqCpJZmJrBjogMkdQLJp7pmNLRGDcLFNakUSYuAvwO+O8mhy4BdiS3VBawOlFQlSbpRkk2ynViS52iipu8uM7ttkrdYQRRxxe2YhnyQBpN2c/cVYOUkx2wq/hGvTf4A8AjwickKr7RypDdSlRSvPlnVCpRxBD1AtHrkR60RvxXaoqQdSVURC3oQ2ELUmzsisZ5q2fVXDyVcSAIuBRbFW+nF32cbVgMuendmttLMVG5Lu27NwIWkrBMkOSBIckCQ5IAgqY40armAIMkBQZIDgiQHBEkOCJIcECQ5IEhyQJDkgCDJAUGSA4IkBwRJDgiSHBAkOSBIckCQ5IAgyQFBkgOCJAcESQ4IkhwQJDkgSHJAkOSAIMkBQZIDgiQHBEkOCJIcECQ5IEhyQJDkgCDJAUGSA4IkB7iTJKlT0u/iFbzelHZ9moE7ScBNwEtpV6KZeFlbCABJ7wDeDrwfeEcVx4eVI5uJpLlEi+N+CBioMltYObJZKFo3bSXwHTN7vIasYeXIg0XSjcDnJjnsJKImbgbRSa+asHJkfah2ec+LgXOAoZIT/bikO83sw42pXmvgYnlPSZ8GPp9IOgq4F7gceKwxtWsd0o6kqjCzF5KvJb0e/7nRzFx2BmrBRcch67iIpFLMbDMZWHu1SIgkBwRJDgiSHBAkOSBIckCQ5IAgyQFBkgOCJAcESQ4IkhwQJDkgSHJAkOSAIMkBQZIDgiQHuJyZPVie2voU0wamNaTs/r7+upeZSUmX3HEJTGlQ4XvqX6TMGvX7Wa2HpG5g162/upUzFpwx6fH9w/1c+4tr2bRzE9/4y2+w+IjFk+Z5cvOTXH3h1QAzzWz3QVeajEbS8bOP54wjJ5bUN9THZT+8jC27t3Dvlfey5KglVZXdiOYudBzKUBT09LanueeKe6oW1CiCpBJaTRAESeNoRUEQJI3RqoIgSAJaWxAESS0vCDIuyYMgyLAkL4Igo5IGhvvdCIKMjjh87pfX8eLeF10IgoxG0uZdm/n5Ffdw1lFLXDzklElJ/3bJTeMiSLT2E2mZbO5OPPyksunlRLXCHEEmJdXCRBHWLIFB0kHQrCYyk9ckbwRJDgiSHOBKkqR3SnpM0qCkXkk/PtCyrMzWqrjpOEh6P9F6d9cDvySq+yn1fI96iFq1dlUdShmPC0mS2oCvA581s+8ldq1PqUplWfHwClauWVn3cl1IAs4EjgYKkp4EjgR+RyRtXaVMlZb3XPPCmrpXcNXaVaxcs5KlxyxlNavrW7iZtfwGXEHUGm0hWn/1zcAPgG1AzwT5llP+8tOMbUHdPn/KJ//GKj7sicCV8d+fSOTtJFor76oJyu8EuhPb0XE5R5ekV9pqPT6Zp7te5ynt5q7alSPfEP89dg0ysyFJm4B5lTJa5eU9+6q5u7TW40vy1A0vK0c+QXSyTwAeitPagQVETeAhTdqRVBVmtlvSd4AbJP2RSMxn4913pVez5uBCUsxngRFgFdBFtPbqxWbWW0MZQ8ANJJrAOh9/oHkmJFNPVXjF1bBQVgmSHBAkOSBIcsAhLUlSj6Q7Je2WtFPS9yRNLznmU5I2S9oTT4M8Ef+AVnL7TuL4D0h6JU4vSHpW0lkT1OEjZcqr7cnatMflGjzs9HOigdizgfOA54AfJPZfTtRV/ihwMvDvwF7gdqJB3OLWHR9/LjBK9FXgeuDW+PUuYE6FOnwk3p8sb25NnyPtE9lAQScRjaG9JZG2FCgAR8WvHwNuSezPxdIeqlDmD4Hekjz/C/QD100gaefBfJZDubk7h+jkJH9v6T4iSWdL6iAaTb+vuNPMCkQSlkjaJmmdpBWSpibK7E7mIfphk6F4XyWmS9oi6Y+S7pY0+WPsCQ5lSUcCryYTzGwE2BHvOxzIA1tL8q0Bngf+nOj3mj4E3JEoM1eSZyvRaPuRFerxe+BjwHuAD8b5H5FU9Q9ueRoWAmr6YawDZQ0wy8zWAmslvQzcL+m4AynMzB4FHi2+lvQI8AxwFfCFaspwJ4nqpzdeAeYkE+Np+J543zaii/7ckrxz4/1Fir/PtChOP7okz1yi5i6ZpyJmtjeeXV5UzfHgsLkzs9fM7NlJtmGi/72zJL05kf1ios/8WHzME8AlxZ2ScvHrRxN5ir9l+3KcvjuZB7gU6CjJUxFJeeDUuLyqP/QhuxF1wX8LnAW8DfgD47vgnyTqSCwnaiL/CxiMT/wCoruSeoFfxcefS9T93gtcB9zCvi743PiY24EViff4J6LfJDyW6F6N4nucXPXnSPtENlhSD9G9EH3xifwPYHpi/wKibvorRE3Wk7HU7URLOQ0Aa0lMhQMfIOosWCz4WeDsxP4HgZWJ118lmv8qNok/A86o5XOEqQoHuLsmZZEgyQFBkgOCJAcESQ4IkhwQJDkgSHJAkOSAIKkESXlJj0j6UUn6zHjS7l+bXqcwLLQ/ko4nujfi42Z2Z5x2O3A6sMSiEfTm1SdIKo+kTxONji8mGkW/i0hQ/R8TnKwuQVJ5FD1o9EuiqYhTgW+a2b+kUpcgqTKSTiSa6l4LnGnRPRJNJ3QcJuZjRHNKC4GqbxypNyGSKiDpXOBXRLOqn4+T/8JSOGEhksoQ32e3ErjVzB4A/pao8/DJNOoTJJVnBdFKadcBmNlm4DPATZIWNLsyobkrQdKFwP3ARWb2UMm+e4lug2tqsxckOSA0dw4IkhwQJDkgSHJAkOSAIMkBQZIDgiQHBEkOCJIcECQ54P8B9gMV4XoD9t4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAAFyCAYAAAD25PTWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAT5UlEQVR4nO2de5BcVZnAf98kmWaSmQlMIIGExwSGh2EEIUbkobBkcUOtVFTKisuupbKFSMm61paWhBUNumuQLfEBii7lOmvALaDKClU+ggUCloFl5WHMMKKbhERYwmsmmUlmJjNJ5ts/7u3xpunnTN++98v9flW3kr59zunT9zffuafvPf21qCpOumlKugNOZVySAVySAVySAVySAVySAVySAVySAVySAVySAUxJEpFFInK3iPSLyKiIbBaRtyfdr7iZmXQHqkVEjgI2Ao8AlwOvA6cCu5LsVyMQKxdYReQW4EJVfVfSfWk0liT1AQ8CxwMXA/8HfEdV7ypTJwfkCnZ3AANx9TOkDXhZ63VwVdXEBuwLt68A5wAfB0aBj5SpswbQhLZF9XrvliJpHHhKVS+I7PsWsExVzy9RpzCS2oCXbvvZbZx94tl17d+6zevo2dTDiuNXsOH6DQBzVXWoHm2bmTgAO4G+gn2/B64sVUFVx4Cx/GMRAeDsE8/molMvKvtiUkPHvrJxLT1/6GHNe9bwrvnvYgMbaqhdGUtT8I3A6QX7TgN21PNFhNoFrfnVGta8ew03Xri6nl2ZxJKkrwPvFJEbRaRLRK4iOC99e7oNC7XLgcYIAkOSVPU3wPuBvwF6gZuAT6vqPVNpb6pi8jRKENg6J6GqPwF+knQ/ooJWX7iauKdepiQljQJrCwQ1ApdUBflISUIQGDonJU1SgsAlVSSpIS6KS6pA0oLAJZWlcBaXFC6pBI38HFQJl1SEYoKm+qG3HrikAspF0HSuUEwHlxSh2iGu0bJcUshUzkGNkuWSmP4kIW5Rmb8sVK9ZXJyiMh1JaZpmlyOzkqwIgowOd+s2rwvWJBgQBBmNpJ5NdgRBRiPpilOv4PJTVvDsK8/Wve0/9v+x7m2aWXdXD0SkHRjkBuCImF5kH3ALUMd1d5mU9HDvw7S2zalLm4VHb3jPMMu7l0NGF0fWjbMXnEV7e/uU65f7sx6aXRcvh5BJSbWS9FjjkoqQtJRCMi8pbUKKkUlJ+e+mWCGTH2at4ZIM4JIM4JIM4JIM4JIM4JIM4JIM4JIM4JIM4JIM4JIM4JIM4JIMYFaSiNwgIioi30i6L3FjUpKILAOuBX6XdF8agTlJItIK3ANcQwZSe4JBSQQJn36qqg9VKigiORFpz28E+e7MYer2uYh8CDgXWFZlldXAF+PrUWMwE0kicgLwTeBvVXVfldXWAnMj2/ExdS9WLEXSUmA+8Ew+AyQwA3i3iFwP5FT1YLRCqcyR1rAk6WHgrQX7fgA8D3y1UNDhhBlJqrqHIBnhJCIyDPSram/xWocHZs5J9SSpfAxTxUwkFUNVL5lOfcHGIslMRlIUCxGVeUmQflEuKSTNolxShLSKckkFpFGUSypC2kS5pBKkSZRLKkNaRLmkCqRBlEsygEuqgqSjySVViWczdsrikgzgkgzgkgzgkgzgkqokyTu4LskALqkKkl4H4ZIqkLQgML5aKE7SICePR1IR0iQIPJIOIW1y8ngkhaRVEHgkpVpOnsxKsiAnTyaHO0uCIKOSrOGSDOCSDOCSDOCSDOCSDOCSDOCSDOCSDOCSDOCSDOCSDGBGkoisFpHfiMgeEXlNRNaLyOlJ96sRmJEEXEyQkPCdwGXALOAXIlKfH4xNMWbuJ6nqiuhjEfko8BpBirVfFasjIjkgF9llMnOkpUgqZG7470CZMquBwcj2UtydigOTkkSkCfgGsLFCGjXPHJkg3wa6gYvKFfLMkQkhIncA7wXeraomh69aMSNJgjC4HXg/cImqvpBwlxqGGUkEQ9xVwEpgj4gcG+4fVNXR5LoVP5YmDtcRnPwfBXZGtlUJ9qkhmIkkVa3bWd/a9MFSJGWWzEqyFE2ZlQR2RGVakhUyL8lCNGVeEqRflEsygEsKSXM0uSQDuKQIaY0ml2QAl1RAGqPJJRnAJRnAJRnAJRnAJRnAJRnAJRWQxmwpLskALilCGqMIXNIkaRUELglItyBwSSbIvKS0RxEYWsFaT1Q12GJqu95kUpKIxPZdpTjazaSkfWNjNI+NVS44BcZiaDeTkjZu3Mjs2bNjaXtkZKTubWZS0osvvkhLS0ssbY+O1v+rUpmUNDExwcTERGxt15vMT8ELqebE3+gvSGcyksoRnUI3NTVNClHVySiJY5pdDpdUBBGho6OD+fPn09bWhqqyd+9eXnvtNQYGBlxSkuRyOU455RQWLlzI7NmzaWlpYdasWQDs37+f0dFRRkZGePnll9m6dWss0+1iuKSQI488kq6uLrq6umhvby9aZu7cIFPO0UcfzezZs9myZQu7d++OvW8uCWhpaaG7u5uuri6am5snh7PCCUJ+f3t7O93d3RxxxBE88/TTjMQw7Y6S+dndzJkz6ezsZPHixcyaNQtVLXnZKL9fVWlubmbx4sWceNJJzJoZ7996ZiXlJcydO5clS5aQy+WqvqaXL5PL5ViyZAnt4TAY19S8akkisjCWHtSIiHxSRLaLyD4ReVJE3jGlhiJDV0dHx5QOsIgwb9482traDmmz3tQSSc+JyFWx9KJKRGQVcBtwM3AusAl4UETm19qWArnm5klB5STJ+vXMXLqUWfPmMXPpUmT9+nx/Jqfruebm2O5N1SLpn4Hvicj9ItIRU38q8U/AXar6A1XtAz4BjABX19JIXsec1laOOuqo8mXXr2fWqlU09fYie/fS1NvLrFWrkAcemCzT0dHBnDlzDmm7nlQtSVW/A5wFzAP6ROSKGPpTEhFpJkjl+VCkTxPh4/NL1MmJSHt+oyC9Z3Nzc8Wr4TO+/OXi+7/0pXwfaGlpobm5ufo3UyM1TUvC9GWXisj1wI9F5PfAgYIy59axf1GOBmYArxbsfxU4o0Sd1cAXp/Oisn17TfvjoOa5o4icBHwA2AU8QIGklLGW4ByWp41IHtbxsbGK93+0sxPpfXMGUe3sBILz0ujoKOPj43XobnFqkiQi1wBfIxhizlTV12PpVXHeAA4CCwr2LwBeKVahVHrP/Al+eHiYgYEBTj755JIvevCmm2ha9eZsbQe/8IXJ/w8MDDA8PHxI2/Wklin4BuCrwPWq+oEGC0JVx4GngeWRPjWFj5+otT0BxsbH2bVr1+TClKKv+773sf+++5jo7kZbW5no7mb/ffehK1dO1hsYGGBsfDy2r3LWEkkzgLMSznt6G/CfIvIU8D/Ap4E5wA9qbkkEVBkaGqK/v7/sZyVduZIDK1e+eX8oaM+ePX9uMwZqmd1dlnRiWlW9F/gM8CXgt8DbgBWqWjiZqKYtAAYHB+nr62NsbKzqWxD5cmNjY/T19TEYXmSN6xaGuctCqnqHqp6kqjlVPU9Vn5xOewcOHGDHjh288MIL7N+/P/8aRQ94fr+IMD4+zrZt2/jTjh0cOHhwOl2oiF8FJ1g80tvby759+8reqsgPh4ODg2zdupUtW7bEfgUcXNIku3fvZvPmzYyMjHDccccxZ86cojf9hoeH2blzp9/0S4r8Oaavr4958+ZxzDHH0NraCsDevXt5/fXX6e/vb3i/XFIJ+vv72bVrV9GFKI3GJRWQv6kHpdfQRcs0AnOzu7ip5uD7aqEG0NTURFNTPH+fcbSbSUknnHCCL9hPOxddeCFtJT4LTZehoaG6t5lJSblcjiNyuViuWOdyucqFaiSTkvzrmAbILyAR6n//J45lXZmfgqcxnWchmZdkAZdE+qPJJYWkWZRLipBWUS7JAC6pgDRGk0sygEsygEsygEsygEsygEsygEsqII3pPl2SAVxShDRGEbgkE7ikkLRGEbgkE7gk0h1F4JJSLwhckgkyLclCFEGGJVkRBBmWZIlMSrIURZBRSdYwIUlEOkXk+yLygoiMishWEbk5TK922GNlwf4ZBH9Q1wJbgG7gLoKUNZ9JsF8NwYQkVd0AbIjs2iYipwPX4ZJSzVxgoFwBEckB0W91tZUqm2ZMnJMKEZEu4B+A71UouhoYjGyJJrCaKolKEpFbREQrbGcU1FlEMPTdr6p3VXiJtQQRl9+Oj+WNxEzSw93XgJ4KZbbl/xPmJn8EeBz4eKXGS2WOtEaiksLsk1VloAwj6BGC7JEfCzMZZ4KkI6kqQkGPAjsIZnPHRHL+FM2/ejhhQhJwGdAVboUnf5tjWA2YmN2pao+qSrEt6b41AhOSso5LMoBLMoBLMoBLMoBLMoBLMoBLMoBLMoBLMkAmJVm7lpRJSRCIsiIrs5LyWJCVeUl50izKJUVIa1S5pCKkTZRLKkGaosolVSANolxSFSQtyiUZwCVVSZLR5JIM4JIM4JIM4JIM4JIM4JKqJMm0Ai6pCpLO++CSDOCSKpB0FIFLKksaBIFLKklaBIFLKkqaBIFLehNpEwQu6RDSKAhc0iRpFQQuCUi3IHBJqRcEGZdkQRBkVJJiRxAYlCQiORH5bZjB621J96cRmJME3Aq8nHQnGomV3EIAiMjlwHuAK4HLqyjvmSMbiYgsIEiO+2FgpMpqnjmyUUiQN60H+K6qPlVDVc8cOV1E5BbgcxWKvYVgiGsjOOhV45kj60O16T0vBc4HxgoO9FMico+qfiSe7qUDE+k9ReRTwOcjuxYCDwKrgCfj6V16SDqSqkJV/xR9LCJ7w/9uVVWTk4FaMDFxyDomIqkQVd1O8l8bahgeSQZwSQZwSQZwSQZwSQZwSQZwSQZwSQZwSQZwSQZwSQZwSQbIpCQhXanSKmHyKng9KSYqbQsnMy+pGIXikpbmkqogaWkuaQpEpTVCWCYlbXr1d7SOzKl7uwoM7xmue7uZlLT87uVwREyN76t/k6Ka9GmxcYhIOzB452N3cm7nORXLD48P84+/+DTbdm/j9r/6Fmcec2bFOs9sf5brLr4OYK6qDk2702Q0kk6bdxrnHFte0p6xPbz33ivYMbSDX1z1IMsWLquq7b0xDHeZ/DBbibyg5954jp9/6GdVC4oLl1RA2gSBSzqENAqCjJ6TijFdQXFOvzySSLcg8EialqBGfXjJtKSpCvJrdw1iKoKS+tifyXPSyPiwGUGQ0Uj63C9v4KX9L5kQBBmNpO2D21P1OagSmZT0b8tvTe0koRiZlHTG0W+pumwaFqtkUpI1XFIVJL38K5Ozu6mS1IIUlzQNGhVdmZUUjYI0TA7KYeqcJCJ/LSJPisioiOwSkfX1aFdJd8pPM5EkIlcS5Lu7EfglQd+76/06eVFTja51m9fVqyuTmJAkIjOBbwKfVdXvR57qS6hLRfnKxrX0bOqpe7smJAHnAouACRF5FjgW+C2BtN5SlUql99z0p0117+C6zevo2dTDiuNXsIEN9W1cVVO/AR8iGIl2EORfXQr8CHgD6ChTbw2HnnIauXXW7f0nfPBvqeLNngFcFf7/45G6OYJcedeWaT8HtEe2RWE7iwr2l9pqLR+t016v45T0cFdt5sjjwv9PnoNUdUxEtgEnlqqopdN77qlmdWmt5Qvq1A0rmSOfJjjYpwO/DvfNAjoJhsDDmqQjqSpUdUhEvgvcLCIvEoj5bPj0/cn1rDGYkBTyWeAAsA5oIci9eqmq7qqhjTHgZiJDYJ3LT7VOWTL1rQqrmLoslFVckgFckgFckgEOa0ki0iEi94jIkIjsFpHvi0hrQZlPish2EdkX3gZ5OvwBrej23Uj5D4rIK+H+CRF5XkTeUaYPHy3SXm3frE36ulzMl51+TnAh9jzgIuB/gR9Fnl9FMFX+GLAE+HdgP/BDgou4+a09LH8BcJDgo8CNwJ3h40Fgfok+fDR8PtregpreR9IHMkZBbyG4hvb2yL4VwASwMHz8JHBH5PmmUNqvS7R5L7CroM5/A8PADWUk7Z7Oezmch7vzCQ5O9PeWHiKQdJ6INBNcTX8o/6SqThBIWCYib4hIr4isFZHZkTbbo3UIfthkLHyuFK0iskNEXhSRB0Sk8tfYIxzOko4FXovuUNUDwED43NHADODVgnqbgBeAvyD4vaYPA3dH2mwqqPMqwdX2Y0v04w/A1cBK4O/C+o+LSNU/uGXpshBQ0w9jTZVNwJGquhnYLCI7gYdF5JSpNKaqTwBP5B+LyOPA74FrgZuqacOcJKq/vfEKMD+6M7wN3xE+9wbBSX9BQd0F4fN58r/P1BXuX1RQZwHBcBetUxJV3R/eXe6qpjwYHO5U9XVVfb7CNk7w13ukiCyNVL+U4D0/GZZ5Glief1JEmsLHT0Tq5H/Ldme4fyhaB7gMaC6oUxIRmQG8NWyv6jd92G4EU/BngHcAFwJ/5NAp+CcIJhJrCIbI/wJGwwPfSbAqaRfwWFj+AoLp937gBuAO/jwFXxCW+SGwNvIaXyD4TcKTCdZq5F9jSdXvI+kDGbOkDoK1EHvCA/kfQGvk+U6CaforBEPWs6HUfoJUTiPAZiK3woEPEkwWNBT8PHBe5PlHgZ7I468T3P/KD4k/Bc6p5X34rQoDmDsnZRGXZACXZACXZACXZACXZACXZACXZACXZACXVICIzBCRx0XkxwX754Y37f614X3yy0JvRkROI1gbcY2q3hPu+yFwNrBMgyvojeuPSyqOiHyK4Or4mQRX0e8nEFT/rwlW6otLKo4EXzT6JcGtiLcCt6vqvyTSF5dUGhE5g+BW92bgXA3WSDQcnziU52qCe0qLgaoXjtQbj6QSiMgFwGMEd1U/H+7+S03ggHkkFSFcZ9cD3KmqjwB/TzB5+EQS/XFJxVlLkBTlBgBV3Q58BrhVRDob3Rkf7goQkYuBh4FLVPXXBc89SLAMrqHDnksygA93BnBJBnBJBnBJBnBJBnBJBnBJBnBJBnBJBnBJBnBJBvh/qg5g9nNrOR8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAAFyCAYAAAD25PTWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAWe0lEQVR4nO2dfZAcZZ3HP98J2cnrBjYxiQlHQhLeQhQhRuRF8cjhhTo5fCkrHHeWL1eKlJxnXWld4ESD3hmkSnxD0aM8c0a8UqqsUKUYLBC1eDlOEGJeQC+JCXK8BLLJbrK7mU0yv/ujeza9szO7M7Pd0/2kn09VVzI9/Tz9TH/29/Qz3U//RmaGJ9sU0m6AZ2y8JAfwkhzAS3IAL8kBvCQH8JIcwEtyAC/JAbwkB3BKkqT5kr4vaZ+kAUlbJL0x7XYlzUlpN6BRJJ0CPAI8BFwJvAKcAexPs13tQK5cYJV0K3CJmb0l7ba0G5ckbQfuB04FLgP+D/immd01SpkiUKxa3QV0J9XOkOnACxbXwTUzJxbgcLh8ATgf+AgwALx/lDJrAUtpmR/XZ3cpkgaBJ8zs4si6rwErzOyiOmWqI2k68Pzt993OeaedF2v7NmzZwPrN61l16io23bAJYIaZ9cZRtzMDB+BFYHvVumeA99QrYGYloFR5LQmA8047j0vPuLThHWuM97/wyDrW/349a9++lrfMfgub2NRw3Y3g0hD8EeCsqnVnAnuS3Gkjgtb+ei1r37qWmy65MZE2uCTpy8CbJd0kaYmkawnOS99IaodZEAQOSTKz3wDvAv4G2ArcDHzCzO5OYn9ZEQRunZMws58AP0m7He0UBA5FUjsZLYraLQi8pKZIQxB4SSOoF0VpCQIvqSEaFZTUZQEvKUKtKEpbEHhJo5IFQeAlDVEdRVkRBF4SkG1B4CVlXhA4dsUhTrI6SKhF7iOpQlYFgZcEpPtFtRFyLynrgiDnklwQBDmW1IqgtGaD5HJ0t2HLhmBOQsYjqEIuI2n9ZncEQU4j6aozruLKxat46qWnmirXSHf3h31/aK1Ro+DMvLs4kNQJ9LAGmJTQTg4DtwIxzrvLpaQHtz7ItOlTmyrb6FHqO9jHymUrIaeTI2PjvDmvp7Ozs6kyjUrqnRKLl2HkcuDQLGn3NV7SGDQraKz5eq3gJY1CFgSBl1SXrAiCnA4cRiPt808tfCRFaFVQklEEPpKAbEZPlNxHUtYFQY4jKS45SXd1kNNIciF6ouRSkmt4SQ7gJTmAlzQO2jFoAC/JCbwkB/CSWqRdXR14SU7grCRJaySZpK+0fd9t3p+Tl4UkrQCuA37X1v22c2cRnIskSdOAu4EP08bUnmkJAgclESR8+qmZPTDWhpKKkjorC0G+u6ZJUxA41t1Juga4AFjRYJEbgc+2vL9WC8aMM5Ek6c+ArwJ/a2aHGyy2DpgRWU5taF9kRxC4FUnLgdnAbysZIIEJwFsl3QAUzexYtEC9zJGjkSU5FVyS9CDwuqp13wWeBb5YLagVsigIHJJkZgcJkhEOIakP2GdmW2uXapysCgKHzklJkmVB4FAk1cLM3jbeOuIWlMSteR9JDuAlOUCuJWX9XFQh15LixmeOzDFeUkz49J4Zx6f3TIi4Bg3tmLLs9JfZNGnnfPLcRtJ4aPeE/1xGUqtdXVpPY/hIapA0H5fJZSQ1QxaeZfKRNApZEAReUl2yIgi8pJpkSRB4SU7gJVWRtSgCL8kJvCQH8JIcwEtyAC8pQhYHDeAlOYGXFDM+B2vGSWqKmL8KHjKe85HPHJlhfNqaNtJKFPlkG20k64Ig591dlnJ/j0ZuI8kVQZBTSVm9slCPXEpqlrQfkfGSHMBLcgAvaQzS7urAS3ICL8kBvCQHcEaSpBsl/UbSQUl7JW2UdFba7WoHzkgCLiNISPhm4ApgIvBzSc39YKyDOHPtzsxWRV9L+gCwlyDF2q9rlZFUBIqRVU1ljszCyA7ciqRqZoT/do+yzY1AT2R5PulGJYGTkiQVgK8Aj4yRRq2lzJGQracBnenuqvgGsAy4dLSNWskc2SpJXrR1TpKkO4B3AG81s0S6r2ZVJn1V3RlJCsLg68C7gLeZ2R8T2U+T2/s8DsP5BnAtcDVwUNLccH2PmQ2Mt/JWOsJ23ZdyaeBwPcHJ/5fAi5Fl9XgrzrIgcCiSzCz2s35WvgeNhUuRFCuuCIKcShqvoHbPkcilJNfwkhzAS3IAL8kBvKQmSWNipZfkAF6SA3hJDuAlNYFP75lx0nwSw0tqgLQflfGSxiALD5s5c6ui3WTpWVofSTVIu3urxkdShKwm3PCRFJK16IniJTF+QUnf5c29pCxHUIVcS3JBEORckivkVlJcUdSOWUe5leQSXtI4aNfcvVx+mTUzMBv3l9da5c3iH47kUpKkcT+rVK90Es9A5VLS4VKJjlJp7A1brDtucinp0UceZsqUKYAwjndbivxbWRftvDS0TXWXpqE1/f39sbc3l5Ke27OHyZMnQeEkzIJDLmCop6qcV6SR75uBHUNmmASaMGy7wwPjflRqBLmUZMeOUi6XgTJmUCaQUIhIIBRQOfhE3lf5WDBAkLCCQKIcbhfUGy/5HIJHTu61Oi5hoazh74ZeoFAIBEnDoigp8hlJlW6KkQfXMBT2b4VCARUKQZSYUS6XsSDmwj/vQs064iaXkqAwFAEVhg2cCwVmzZzF7DlzmTZ9OmUzDh48xN69e+nu7g4DcULtsgmQT0mR7k6RVZOKRRYvWsS8efOYMnUqkyZPZuLEiZjBkaNHGOgfoK+/nxdfeIGdO3dSKpUIe7xh9cRNPiWFRA/uKSefzJIli1m85AxmdHYCw7sxA2bMCDLlzJo1iylTp7Bjxw56DhwIz11lkjrF51JS9ZWbKZMnsWzZuSxZvISOYgdlC848w8LCgvOVGcyY0cmyc5dRLE7i6d8+SX//IVQuYxo52IiDfI7uKpgxcUKBBQtO4/TTT2dix0TMjEKty0YKLvkUCsLM6OjoYNGi01mwYAETJ5wUDM0TEAQ5lhRIKDNjxgyWLl1KR0ex4Wt6CvvJYkeRs5cupTPsBlOfdydpXkJtaApJH5O0W9JhSY9LelNLFYV/9Z2dnczsmhlGSNW+Ri0PKoiZXTMDSRLxZ5oIaCaStkm6NplmNIak1cDtwC3ABcBm4H5Js5utywgi4ZSuU4IIov7ITBs3MnH5cjpmzmTi8uUUNm4c6v5UEF1dM+koTkrsC1Mzkv4F+LakeyR1JdOcMfkn4C4z+66ZbQc+CvQDH2qqllDG1GlT6Tql1kc5frS1cSMTV6+msHUrOnSIwtatwet77x3a5pSuLqZOn45JiQzBG5ZkZt8EXg/MBLZLuir+5tRHUgdBKs8HIm0qh68vqlOmKKmzslCV3rOjYxKTJ09mtNt/Ez7/+drrP/e5ShuYPHkyHR0dzX2gJmhqCB6mL7tc0g3AjyU9Axyt2uaCGNsXZRbB1/yXq9a/DJxdp8yNwGfHs1Pt3t3U+iRo+nuSpAXAu4H9wL1UScoY6wjOYRWmE8nDOlg6zMDAABpliGALF6KtIzOI2sKFQHBeGhgYYLA0GFOTR9KUJEkfBr5E0MWca2avJNKq2rwKHAPmVK2fA7xUq0Dd9J5h79bX10d3dzeLFi2qKnlc2rGbb6awemS2tmOf+czQ//d3d9N36GBwjymBwUMzQ/BNwBeBG8zs3W0WhJkNAk8CKyNtKoSvH2u2PgGlwRL79+/HwkkptQ6wvfOdHPnRjygvW4ZNm0Z52bLg9V9fHZQrG93d+xgsHU7si1IzkTQBeH1SeU8b5HbgPyU9AfwP8AlgKvDdpmuSwKC3t5d9+/ZxStdMJkwYfpSHbqFffTVHrr56RBV2zOje301vT09wxSEhS82M7q5IWRBm9kPgk8DngKeBNwCrzKx6MNFIXUCBnp4etm/fzuBgKYiMBrorC+/ylUolntm+nZ6eA8H6ZhvRIM5dFjKzO8xsgZkVzexCM3u85cokjhwrs+e559i1axdHBo8gQdls5Py5oRt/hiQGBwfZtWsXz+3Zw9Fjx47fqU2AXF4Frz6W/QOH2bZtO6XDJRafUftWBQIRfFnt6elhx46d7Nyxg/6BAaCAFQq1K4+BXEqqEJ2bsP/AAX63dSv9/f28dt5rmTp12tBNv7LB0SNHGBgYoK+/jxdeeJGdO3cyGN70S+QyQ4R8Sop0ZRZZVSqV2LZ9K9u3baVr5kxmz5nL1GnTMODQoUPs3fsK3d37RtyFjcpOYgieT0mUg9FYOMKD43Prghdl9u17lf0HekZMRBma9lW5E6vhEyyTIJeSNDS37qQRE+/F8Xvq5XIZKx+/shfMYA0FVe7E1qgjbpwb3cVCpE+qPptYOPNuaIZJBKkSROXjd2LNhs2VSIJ8RtKEkygUClAoBBMeqZ5mfHykVnOasSYMTTOWgi5P4XaFQvx/97mUdNqCBX7Cfta55JJL6ezsbPk8MlrX1tvb22Kt9cmlpGKxSLFYTCRNTbFYrPNO6+RSkoWXfZIYkfnHMWNC4chtvCOyWuWTeBwzn0PwmGhXRhUvyQFyKymuTqkd0ZRbSS6Ra0mu/BpZriWBG6JyLwmy//NxXlJIliPKS4pw/AJq8/j7SW0ma1HlJdWhlajy8+5SIguivKQGSLv785IaJE1RXlITpCXKS3IAL8kBvCQH8JKaJI3zkpfkAF6SA3hJDuAlOUAuJWX95+GqyaUkcOdXyCDHkoDEn9CLi1xLqtCKqHZ2eU5IkrRQ0nck/VHSgKSdkm4J06vFQitR5X/kajhnE/xBXQfsAJYBdxGkrPlknDsa9oBzAyT9vCw4IsnMNgGbIqt2SToLuJ6YJUH2RDkhqQ4zgO7RNpBUBKJPdU2vt+14SbLrc+KcVI2kJcA/AN8eY9MbgZ7I0nACqyQe1WyVVCVJulWSjbGcXVVmPkHXd4+Z3TXGLtYRRFxlOTWRD5IwaXd3XwLWj7HNrsp/wtzkDwGPAh8Zq/K6mSMbpNlzU1KkKinMPtlQBsowgh4iyB75wTCTcS5IO5IaIhT0S2APwWjuNZWoMLOa+VdPJJyQBFwBLAmX6pN/FnqkRHFidGdm681MtZa029YOnJCUd7ykMcjCVXIvyQG8JAfwkhog7S4vl5JcGxLmUhJk4+GwRsmtJGhNVBqyci0JWuv62i0q95Ig+6K8pHHQru7PSwoZz4gvaVleUoz4PA6O4PM45BQvKUJWr0R4SQ7gJTmAl+QAXpIDeElVZHHw4CU5gJdUg6xFk5dUhyyJ8pJGISuivKQxGE8a6rjwkhrEp/dsM+N5is+nUnMEn7bGEdoZVbmVFNfNuXaIyq2kOElalJcUEz6PQ87xkmIkqWjKtaS0H2lplFxLcgUvyQFyLynuLu+ESwCVFbJ+bnJOkqSipKfDDF5viKveLItyThJwG/BCEhVnVZQruYUAkHQl8HbgPcCVDWzfdObIrKRPi+JMJEmaQ5Ac931Af4PFWsocmbV84U5IUpA3bT3wLTN7oomi48ocmRVZqXZ3km4F/nmMzc4h6OKmExz0hhlv5siheki3C0z7nNRoes/LgYuAUtWBfkLS3Wb2/mSad5w0RTmR3lPSx4FPR1bNA+4HVgOPJ9O6kVS6vnbLSjuSGsLMnou+lnQo/O9OM2s4jXRs7cH/VoWnCiciqRoz203KX2faGU0+khzAS3IAL2kctOuLrpfkAF6SA3hJDpBLSVm7FTEWuZQE8U249/nu2oALUZV7SZCNRy5Hw0uK0KqopLs8L6mKLEaVl1SHLCV395JGISuivKQxyIIoL6kB0j5HeUkO4OSd2fGy+eXfMa1/atPlGunK+g72Nd+gMcilpJXfXwmTEqr8cPxVyiwLczTbg6ROoOfOX93JBQvPH3P7vsE+/vHnn2DXgV18/S+/xtLXnDtmmad2P8X1l10PMMPMesfdaHIaSWfOPJPz544u6WDpIO/44VXs6d3Dz6+9nxXzVqTW3fmBQw0qgra9uo2fXXMfK+atSLU9XlIVowlKayjuJUXIWgRV8JJCsioIvCQg24LAS2pakM8c2UYM6G0xgnzmyDZxsHSQq8bRxflHXxKmb7BvSNB919zHG1s8B7VLVC6vOKz5xRqeP/I890UiqNVHWYSf45AIu3t2DxM0Xnx6zwS4beVtNQWNJyJ8es+YOWfWOYnU6zNHtoks3rjxkhzAS6pB1qLJS3IAL6kOWYompyRJ+itJj0sakLRf0sYk95cVUc5ccZD0HoJ8dzcBvyBo+7Kk99vslYgNWzbE3gYnJEk6Cfgq8Ckz+07kre0pNakmX3hkHes3r4+9XickARcA84GypKeAucDTBNK21itUL73n5uc2x97ADVs2sH7zeladuopNbIq3cjPL/AJcQ9Dz7CHIv7oc+AHwKtA1Srm1HE8A2e5lYWyfP+WDf2sDH/Zs4Nrw/x+JlC0S5Mq7bpT6i0BnZJkf1jO/an29pdnto2U64zpOaXd3jWaOfG34/6FzkJmVJO0CTqtX0Oqn9zzYyOzSZrevKhMbrmSOfJLgYJ8FPByumwgsJOgCT2jSjqSGMLNeSd8CbpH0JwIxnwrfvie9lrUHJySFfAo4CmwAJhPkXr3czPY3UUcJuIVIFxjz9q2WGZVcPVXhKk5dFsorXpIDeEkO4CU5wAktSVKXpLsl9Uo6IOk7kqZVbfMxSbslHQ5vgzwZ/oBWdPlWZPv3SnopXF+W9KykN43Shg/UqK+5J2vTvi6X8GWnnxFciL0QuBT4X+AHkfdXEwyVPwgsBf4dOAJ8j+AibmXpDLe/GDhG8FXgJuDO8HUPMLtOGz4Qvh+tb05TnyPtA5mgoHMIrqG9MbJuFVAG5oWvHwfuiLxfCKU9XKfOHwL7q8r8N9AHrBlF0oHxfJYTubu7iODgRH9v6QECSRdK6iC4mv5A5U0zKxNIWCHpVUlbJa2TNCVSZ2e0DMEPm5TC9+oxTdIeSX+SdK+ksR9jj3AiS5oL7I2uMLOjQHf43ixgAvByVbnNwB+BPyf4vab3Ad+P1FmoKvMywdX2uXXa8XvgQ8DVwN+F5R+V1PAPbrl0WQho6oexWmUzcLKZbQG2SHoReFDS4lYqM7PHgMcqryU9CjwDXAfc3Egdzkmi8dsbLwGzoyvD2/Bd4XuvEpz051SVnRO+X6Hy+0xLwvXzq8rMIejuomXqYmZHwrvLSxrZHhzs7szsFTN7doxlkOCv92RJyyPFLyf4zI+H2zwJrKy8KakQvn4sUqbyW7Yvhut7o2WAK4COqjJ1kTQBeF1YX8Mf+oRdCIbgvwXeBFwC/IHhQ/CPEgwk1hJ0kf8FDIQHfiHBrKT9wK/C7S8mGH4fAdYAd3B8CD4n3OZ7wLrIPj5D8JuEiwjmalT2sbThz5H2gUxYUhfBXIiD4YH8D2Ba5P2FBMP0lwi6rKdCqfsIUjn1A1uI3AoH3kswWLBQ8LPAhZH3fwmsj7z+MsH9r0qX+FPg/GY+h79V4QDOnZPyiJfkAF6SA3hJDuAlOYCX5ABekgN4SQ7gJTmAl1SFpAmSHpX046r1M8Kbdv/W9jb5y0IjkXQmwdyID5vZ3eG67wHnASssuILevvZ4SbWR9HGCq+PnElxFv4dAUPyPCY7VFi+pNgoeNPoFwa2I1wFfN7N/TaUtXlJ9JJ1NcKt7C3CBBXMk2o4fOIzOhwjuKZ0ONDxxJG58JNVB0sXArwjuqn46XP0XlsIB85FUg3Ce3XrgTjN7CPh7gsHDR9Noj5dUm3UEiVDWAJjZbuCTwG2SFra7Mb67q0LSZcCDwNvM7OGq9+4nmAbX1m7PS3IA3905gJfkAF6SA3hJDuAlOYCX5ABekgN4SQ7gJTmAl+QAXpID/D/8a0YmGuNHrAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAAFyCAYAAAD25PTWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAcKUlEQVR4nO2de5QcVZ3HP7+ayQwTJhMyEt5CEJRAECGYhIeg4uoC68qqx+NrPT72yOOI6LqgkCghLAF8o+ADWSWAoC7nuLCrx+Dia4+A0fCIQIJKwlN5Z0gmyWQmmf7tH7e6p7q6qrp7pqqrinu/59RMd9W9t2/Xp3+/+6t7b90SVcWp2PLyroBTczlIJZCDVAI5SCWQg1QCOUglkINUAjlIJZCDVAI5SCVQqSCJyL4i8n0ReUFERkTkfhF5bd71ylrdeVegVYnILOAO4FfAKcBzwCuBoTzr1QlJWTpYReRy4HhVPSHvunRaZYK0FrgN2A94PfBX4Juqek1Cnl6gN7R7ENiYVT19zQD+pmmdXFUtxQZs97dLgaOA04ER4IMJeS4CNKdt37S+e5ksaQxYrarHBfZ9HVigqsfG5Alb0gzgySv+8wKO2GsueFJNWE0PeGa/CIiJq1Q8wN/nmdeKZ4776W5Y+wNWrLmOk/c7mZVnrwSYqaqb0/jupQkcgKeAtaF964B3xmVQ1VFgtPpefBhH7DWXE+f4QaE3AUOqYDwvAElQEZCuCSieh/rvVTwu/d2XWPGn67joLUs5YY8TWcnKFL92uULwO4BDQvteBTw25ZK1ErM/7GUavc6ld32Ri37771x0woUsPv6CKVclSmWC9FXgGBFZLCIHi8j7MO3SN6ZUamXixKuqAVYJQtPQfww8rXDp777IRXdcwkWv+yyLj/t0BNR0VBpIqvoH4O3Ae4EHgM8Bn1TVG9svLMJywvu02v6DVE++ag3Epau+zNI7LmXZ8UtYfOyn265COypTm4Sq/gT4STqF+VDEM9bkBxGqilCh9vtVNe2QVky7hLJ81VdYetfnWXbcBSw+5lzfAoOBXboqjSVlpiqsitZbU6UScF/qW5Oy/PdXGEDHfobFiz5ljlGFFNO2TVEOUoRqbZN5V4O1/PdXsPR3X2DZMZ9myaJPxWVOvT72QgqezBasafnqr7N01ZdYdsx5LFn4SdCABYW3lGUppMZAoA6OVghe5C9ffSUXrvoyFy/6N5Ys+ER9Xq1Q7/LSl6WQElSpt7Dlq6/kwt9/hYsXfoolrz3HQEmE4SwpG0VZE7D8nm9y4R+u4OKF/8qS136cegAa7fIykINUVfAEV9QAWv01Ll7wSZYcfTb110wJMFyb1AFpheX3fZuld3+dZUd/giVHnTXRA1EXbGh0u5SB7IYU0cuw/L6rWXrPlSw7+hyWzD+rtXz1B1OrXlV2Q4K6E758zXdYeu9VLJv/cZYcdWYoXUQfXofkIPmqATrqbJYceWbtmil5vC0ieMhADhJhQGckuDONeR3c7dxd6lq+5poJQK85PTpRCJpkZDFxshrS8jXXsPS+byQDKoCshbR8zXdKAQhKNp6Ulm54+H+49on/KgUgsNSSrv1LeQCBpZb0tv3fyCn7vY57XghPPgIzdcv/7fpTu0QEvK7QNK/AexEUk+bPGx9Ovb6lmXeXhkRkANjE+cAuGX3IduByIMV5d1ZC+uVvrqW/vy8qhf8vbEldERMmoy1py5YR3nTkqWDp5MjUdOTL5jLQv2vjAQlA8iZei9fVOGHSq06W9CGJB14Xm6ZvTb2+VgYOwASQxDRebdZrnrIXUlgFgBEnuyEVGExQVrZJdQqDkqn9brMY+LPbkkoiOyFN0Vo6rXLVNnfl04Y5SJOQdjjgcJCCCvY0FEj2QipRu1SemlosuyEFrSlsWQ3v83OB7mK2BG6v+DXMUUXoXAVbISVFb+Fj1XUdcpSdkEomB6kEshdSlMsr2EVsVaWFJCLni4iKyBXpF550WiJAZhxglBKSiCwAzgD+OKWCpmQ5nbO60kESkX7gRuCjpLm0Z0FdHZQQEmbBp5+q6u3NEopIr4gMVDfMenf18qTQgKBkPQ4i8h5gPrCgxSwXAEsbC2r+2yzKhSyUyJJE5OXA14D3q+r2FrNdBswMbPs1/6DinZIyWdLRwB7APYFfeRdwooicDfSq6ngwQ9zKkZNWTtZVJki/AF4d2nct8BDw+TCgphKvyV3kFMaqSgNJVYcxixHWJCJbgRdU9YHoXGkqvzaqGD+VoijKcqpzwHNUaSwpSqr6hikVEOPyihTZQckhpaKCtDtJKn4NM9DUo7zOnjYrIaWnzrhFBylOBXKDxalJhxXn8ooWNIDFkMokBymgVqyo0/PAwXJITaG0M0sowzbMTkgBOFVQRWyLqrITUkgNgOqmHwePufuTCizJtf/OQSqB7IUU19C3GgB00LLshdRMOc//Dqo4NclDBer6SVI5atkpTQZaB9yenZCCrqwKpsBWVdyadVKxQUQxLnAdpCg1DRrcOg4FU/7WZC+kAoXYzVSemnZKUfCCbVMOAYbdkJKsqSFokFzGksB2SGEV1AXaOe8uaBFTAdMhyyrmT8epTg5SlGoW0q6lZGNZ9kJK0VVpxtdS9kIqkRyksBrmO8RZiRv064wK0oHaTHZDglBvQtzMIHcTWQ5q1aXFZXe94IWUBh7L02k5SCWQxZCSLKL6gKtWnrHk5jg4YT2kKCsoXlhuOaRWlD80OyG1cj1UoAvd0kASkQtE5A8iMiwiz4rILSJySEqlMymLib0QTlelgQS8HrMg4THAm4FpwM9FJOJZpOmqnWHzLHrES9PjoKonB9+LyIeAZzFLrP1fVB4R6QV6A7tmBA5CSR6EXCZLCmum/39jQpoLgE2B7cmWSo6ynBwfMFJKSCLiAVcAdzRZRq39lSPrPynmfWeDitK4u5C+ARwOvC4pUdOVI6NcXjMLySHqKx0kEbkKeCtwoqq25r6SCyx821QaSGLM4Erg7cAbVPWRKZQWLjw2ZUuRXcazWksDCePi3gecBgyLyF7+/k2qOpL+x7mL2cnoLEzj/2vgqcD27hzr1BGVxpJUNeWftkCrzynP+S7AMllSB1UcVwfWQ0qGobEdsYE0HQBqOaQoFcuKwFJIjRbSYk9CyxMn05WVkKJVf8JrIBuCBn8xqMj+PTdh31pZCymvWysnI2shJaloAK2G1DqMfKFZDSlKdeDa6mlwcxwyUxBKomXl6AJL03eXrup/m0Vrg8Ky3pJiVXN1+d+n5CClIfdY7Yw02eGHpDXEM5K9kJLU7MT73UKd6AEH2yFFwajb18Z9ShnKbkiQ+6hrK7I0BJ9Mu+Ieg1BAFefayUFqpnZvb8mg/XKQIlUcKwKrIU0ShHtcXBEU5d7ymdtQleWQJn+yO3UhC5aG4IqisXdS+PuDd1uI/yeQx/Sch+7IyOgODSshCdL8QYthVyeB/w2vA91EGbhCKyGNjo0yOtoTnyA4EFjb59XDqQHx6vryRkfHUq+vlZDuuONOpvftUr8zakqxNP7XajNeN/duAtK2kfTvwrES0uNPPkHfLgFIDRMhpREQ4TQT+2r3e4jHyPbtqdfXSkgVhEpwarEy4cIQquvaBU9+oqoJValUXOCQjrwuswFhKCJeLXLTOpfm+akNBEUbIjlRdYFDapJukG5jANXVIMVAUekyrxXE8xAxWFQVrVQQ3wJrhgc+rIqf30FKRep1oV3mqxtrmQBjgMCswUFmz96DXWf0g8KWLVt49tlneXFoI6iaMB7fQYoAVctMv3/AXkhed801VXzPNa2nhwNfcRB7770PfdOn09fXR3e3OUU7d+5kZGSEbVu38vTTT7F+/Xp2jI36oASPqnU5S0pFiodiLKeiigIzZ+7GKw4+mANfcRAzBwYi8w34+wd3351d+qazYf3DbHrxRdNOiSAZLQdhJaQKyrgPp6LQ19fHofMO56CDDqKnp4eKKiJC2CaqbdPMgQHmzZvHLr293HvvPYyMjACK57vMtGVlB2tFYVzN/67ubvY/YA5z5sxhWk+PH4371z+hDf+YAj09Pcw58EBevv8BdHd3m/Yso5VVrIRUPekVhRkzBph76KH09PY2y9ZQRk9vL4cedhgzBgYMRBGywNQyJBHZJ4PPb1si8jEReVREtovIKhFZ2G4Z1UBBVRmYOZPBwUG8STyRzPM8BgcHmTHDtFUmOkwfUzs1e1BE3pd6DdqQiLwb+AqwDJgPrAFuE5E92imnCmhaTw+zZs1q2iM+7dZbmLFoAbvtNZsZixYw7dZbgnVi1qxZTPPbsizUDqQlwNUicrOIDGZSm+b6FHCNql6rqmuBM4FtwEfaLUiB/v5+Zs2alZhu2q230P/+99L94APIli10P/gA/e9/L9P++9ZamlmDg0yfblYZbQw3pq6WIanqN4EjgJcBa0XkH1OvTYJEpAezlOftgTpV/PfHxuTpFZGB6oa/vGf19z5tWg99fdMTP3eXS5dH719+Se319L4+enoShj6mqLZCcH/5spNE5GzgxyKyDtgZSjM/xfoFtTvmsv6Z0P5ngLkxeS4Alk7lQ7see7St/Vmo7eskETkAeAcwBNxKCFLBdBmmDatqBvBk1SHt2DHGyMi2xALGD5hD94ONK4iOHzCn9nrbyAhjY+kP9lXVFiQR+SjwZYyLmaeqz2VSq2g9D4wDe4b27wk8HZUhaXlPwfTHDQ0NJX7o9sVL6H//exv3L/ls7fXQxo1s27bVfGae0Z2IrAQ+D5ytqu/oMCBUdQy4G3hToE6e//6udsoyIxDCjrExhoaGEi9Cd5z2T2y56YfsnHc42t/PznmHs+WmH7LjbadV68XQ0BA7xsbwMprq1Y4ldQFHpLLu6eT1FeA6EVkN/B74JLArcG07hXi1YSJh86ZNbNy4kVmzZsVeK+1422k1KEFVKhWGhoYYHt4M+PBzju7enDMgVPVHwLnAxcB9wJHAyaoaDiYSVR1i8ASGhzezbu1axkZHm2VrKGNsdJR1a9eyedMmg0azmY1Xum4hVb1KVQ9Q1V5VXaSqq9otwxPoEvN/fOdOnnj8MR595BF2jI2Z9SR99yehzf98A2hsjEceeYQnHn+M8fHxmgvNQlb2gnsIXX4/m6CMbh9h3doHGR0d5cCD4ocqqgN8mzZtYsP69WzYsJ7tIyN4Ap7fa54FJyshmWkoZjjCE6GiMLx5E2sfvJ+RkW3stffeTJ++K319fUzr7kYJDPpt28pTT/2NDevX+8GCAVcb9MtAdkKqjCMVc3mnInj+8Pn42Bh/eWgdf163jt0GB5k9eza79vcDJlx/7rlnGdo4VAPcFWU9bppxOpLKODK+szY7CCq1yY4iXeAJm4c2MrzpRTN+5Pdua6VCtx8VNoKp+G8q4Y+bsqyEhO4E3ekPd09M6TKThsZRMdNMVCdm/wie7yb9AUGtn9Il1WDDWVJKqoybDTDhwLhvRf60kgA0k6QxCJa4dw5SOvL8qShGtYFxzBWJbyESmNmqFZKnGWvtmOe52UKpaP/9Xu4m7Bddxx9/HAMzZsQniLzjvH4SvwZmvgZnwG7evCX1+loJqbenl95mE0+a3ESm1YAjdBNZb6+7PykVNd6Oafoe4lLX2qfg3uqtl7X2S8DzXOCQlupvxwwDinF11fd1bZHWW1NdnvRUug7WdDX5X302M+yiZTmkKAVOvobD9PD+zshiSJM80Tk8pNFiSEkq1tMyHaRm0ij312L6lOQgxao41mRlCN4YCASGF2LvNE+6lspWzpI0/fGftGU3pChAdft8i8v5sdt2Q4pTM+vyu4M6dUFrL6TJurmwVXXAXdoLKU1l7A4dpDjVLCTcW955WQqp3kWJaiYTSNKSpZAmJBEzfiKVI0TrIYVVB6qtoCA7iFZDat3Fueukwqlo7ZO1kIoGIknWQiqTHKSaIpbqhIjgobbmTUQR+a+I8pKRNEznipnLEFZDl5Dru8tJxWurLIeUDKTR4iLSdACq5ZDiVCxrshhSGyByHr21GFJ5VApIIjJHRL4rIo+IyIiIrBeRZf7yahmpOC6vLLOF5mJ+UGcADwOHA9dglqw5t/3iEkLp0IR7UW2+ar5WQLqS00xBpYCkqiuBlYFdG0TkEOAsJgWprvApZe+ESgEpRjOBjUkJRKQXCN4tVn97X1yvQZLlZPTQkCSVok0KS0QOBj4OXN0k6QXApsDW5gJWYYgt9kykrFwhicjlIqJNtrmhPPtiXN/NqnpNk4+4DGNx1W2/lioWaWGV5mkyUt7u7svAiiZpNlRf+GuT/wq4Ezi9WeFJK0eWoS2qKldI/uqTLa1A6VvQrzCrR37YX8m4I2opwqumzcAV5m1JLckH9GvgMUw0N7v2PAnVyPVX21NwwY12sgWCiAwDilJAAt4MHOxv4ca//TMT23EauHMihyguTqWI7lR1hapK1NaBT8/+I5qoFJCyUxSA/KGEZTmkcshiSEkW08Z9SR0I5S2GVB45SC1KVP1eB7eOQ4cUMVTRjttyK6J0WLHrNOR/X1JVdkMqSf+d3ZCi1PIESLdKl1NA9kJK0dVlPUHSXkhJilvnrnnGtGsCOEilUFmGKtJV0NVVAmOH7T61uUPDGc6SgqoUczEouyElQWkILPJb68FuSFGKAjfpZQPSkb2QCuraomQvpJaVf9eRgxSlplbmesE7r7h2piAdsHZCClpKFVCB12K1E1KcJgPKzXHIWAW2nqDshpSkAoXo9kKKDRZahNPBoMJeSCWSg9SS2pxNlLIcJAg93496l9fC0jVZy05IgRNfBdQAqkCyE5KvpmDaifAyDOethhRWK9aUx5iSg1QCWQspzmqK2DZZC6mpCtRl5CBNSW4N1sw0ZZfWYSuzc95dUC09dDFfFbNWnVKMRRQteLAbUlhR0CqV3IfRSwdJRHpF5D5/Ba8jO/fJroO1HX0B+NuUS2ml8S9IGF6qwEFETgHeArwTOKWF9MkrR7arnNYbKo0liciemMVxPwBsazFb+ytHFsR6gioFJDHrpq0Avq2qq9vIGr1yZAsgihTh5eruRORy4DNNkh2KcXEzMCe9ZSWuHFlVxYfhFWPZtCjl3Sa1urznScCxwGjoRK8WkRtV9YNTrklFCwuqFMt7isg5wGcDu/YBbgPeDayadAUqU3FpgQUMM1beltSSVPXx4HsR2eK/XK+qbS4j3coHVhK6iCLgZBz1lSJwyERRVjQly8pOpbCksFT1UTrlawogey2pRLITUpJbCx/TSu7zwu2E1KKKckHrIGmlkF1BQdkNSSPu+It974YqnBJkL6SCu7ig7IUUpSq4gl3UOkiTUKfngztIbcndn9Q5lag9AlshlUwOUnh1/SlamWbQ72s3pIJ0+zST3ZCCKjAweyG19GykSiE6We2F1EwFuqAt5cjsVHXfCw/RP9oXc1Qm5jd4AiJmKpjXZeYx+MdUAu9FTMDgdbFly0jq9bUS0kk/+zDsklHh29Mv0kpIVx+3lKP2PTTmqB9Ci8e28e18/K5L2LD5cb5xwsXM2/2QGEvyUP//PY//kbMuPyfV+loJ6ZUzD2D+yw6LPuhPzRreOcKpt53BY8N/5fa3Xs/CvecH3J2gXlcNEJ7nQ/PYsqXVaeqtywUOERresZVTbzudB4f+wspTv8fCPY7ItT52Q4oYOh/esZVTf34mDw49zMpT/iN3QGA7pJDqAJ18TSEAgYNU0/COrZz6v2cZQH//HRbODgAKTTnWDt9I5iARAWiP18SklJjXwd2ugzV11QF6y9UsnP1qc8AzkVzkPU01ycSFrEhmk/athjS8Yyun3v6xRkBRqgFw98x2TJGAWrWExJVTnLtLRSM7R6NdHDQC8C9WzesAAAm4uOoFbkZWZmWPw3mrv8STO5+pBxQEENkemdcmskuA4QKHdPTo8F+jAcW5scj2KCJocIFDevrSgnPjg4Tgzc1Rrs4fmohX+qCsdHdzdzvIvJAQkMDrpq4uqj1ylpSyogB5gQG/qhW14+oycnn2QkpQgxWFXV0SCAcpA0VZETRYUUNUF+nqsjmdVrZJQAjIxMkXCcFqiPyauTpnSekoKtRuuIidOOG1Xu8Mw+wklQqSiPyDiKwSkRERGRKRW6ZcaJQV1T2oPiJwiHN1GQEsjbsTkXdi1rtbDPwSU/fD0ym82UVsbUfTom64/4ap1yekUkASkW7ga8B5qvrdwKG1Uyq41VW5otqoCF16x2WsWLNiSlWKUikgAfOBfYGKiNwL7AXch4H2QFymuOU9//j0Q+adF2hrqIbeXm1S5MT0LePaghewihcIyT1uWPsDVqy5jpP3O5mVrEzxq2MWlCj6BrwHc5vdY5j1V48GbgKeBwYT8l3k58tjm5Pa98/55F/ewpedC7zPf316IG8vZq28MxLK7wUGAtu+fjn7hvbHbe2mD+YZSOs85e3uWl05cm//da0NUtVREdkA7B+XUeOX9xxW1c3NKtdu+lCe1FSWlSPvxpzsQ4Df+vumAXMwLvAlrbwtqSWp6mYR+TawTESewIA5zz98c34164xKAcnXecBO4AagD7P26kmqOtRGGaPAMgIuMOX0k82TKNEC3MnmlKxSdQvZKgepBHKQSiAHqQR6SUMSkUERuVFENovIiyLyXRHpD6X5mIg8KiLb/WGQu/0HaAW3bwfSv0tEnvb3V0TkIRFZmFCHD0WU196dtXn3y2Xc7fQzTEfsIuB1wF+AmwLH340JlT8MHAZ8B9gBXI/pxK1uA37644BxzKXAYuBb/vtNwB4xdfiQfzxY3p5tfY+8T2SGgA7F9KG9NrDvZKAC7OO/XwVcFTju+dB+G1Pmj4ChUJ7fAVuB8xMgvTiV7/JSdnfHYk5O8HlLt2MgLRKRHkxv+u3Vg6pawUBYICLPi8gDInKZiEwPlDkQzIN5sMmofyxO/SLymIg8ISK3isi8dr7ISxnSXsCzwR2quhPY6B/bHegCngnlWwM8ArwR87ymDwDfD5TphfI8g+lt3yumHn8CPgKcBvyzn/9OEdmv1S9Spm4hgHYejDVZrQF2U9X7gftF5CngFyJy0GQKU9W7gLuq70XkTmAdcAbwuVbKKB0kWh/eeBrYI7jTH4Yf9I89j2n09wzl3dM/XlX1+UwH+/v3DeXZE+Pugnlipao7/NHlg1tJDyV0d6r6nKo+1GQbw/x6dxORowPZT8J851V+mruBN1UPiojnv78rkKf6LNun/P2bg3mANwM9oTyxEpEu4NV+eS1/6ZfshgnB7wEWAscDf6Y+BD8TE0hchHGRPwBG/BM/BzMraQj4jZ/+OEz4vQM4H7iKiRB8Tz/N9cBlgc+4EPNMwldg5mpUP+Owlr9H3icyY0iDmLkQw/6J/B7QHzg+BxOmP41xWff6UF/ALOW0DbifwFA48C5MsKA+4IeARYHjvwZWBN5/FTP+VXWJPwWOaud7uKGKEqh0bZKNcpBKIAepBHKQSiAHqQRykEogB6kEcpBKIAepBHKQQhKRLhG5U0R+HNo/0x+0W97xOrluoUaJyKswcyM+qqo3+vuuB14DLFDTg965+jhI0RKRczC94/Mwveg3YwCt6XhdHKRoibnR6JeYoYhXA1eq6iW51MVBipeIzMUMdd8PzFczR6LjcoFDsj6CGVM6EGh54kjacpYUIxE5DvgNZlT1s/7uv9McTpizpAj58+xWAN9S1V8B/4IJHs7Moz4OUrQuwyx/cj6Aqj4KnAt8QUTmdLoyzt2FJCKvB34BvEFVfxs6dhtmGlxH3Z6DVAI5d1cCOUglkINUAjlIJZCDVAI5SCWQg1QCOUglkINUAjlIJZCDVAL9P6A6M/7R8vpFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "kx = [0.4, 0.4, 0.1, 0.3, 0.25]\n",
    "omega = [0.1896, 0.3175, 0.4811, 0.8838, 0.2506]\n",
    "filename = [\n",
    "    \"media/holey-wvg-bands-{}-{}.mp4\".format(k, om) for (k, om) in zip(kx, omega)\n",
    "]\n",
    "for (k, om, fn) in zip(kx, omega, filename):\n",
    "    run_sim(*[k, om, fn])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video src=\"media/holey-wvg-bands-0.4-0.1896.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<video src=\"media/holey-wvg-bands-0.4-0.3175.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<video src=\"media/holey-wvg-bands-0.1-0.4811.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<video src=\"media/holey-wvg-bands-0.3-0.8838.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<video src=\"media/holey-wvg-bands-0.25-0.2506.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import IPython\n",
    "\n",
    "for i in range(len(filename)):\n",
    "    IPython.display.display(IPython.display.Video(filename[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* From the top, the first two pictures show the first two guided bands underneath the light cone at $k_x=0.4$. Note that the second guided band is propagating to the left, which is due to its negative slope (note, however, that there is a corresponding right-propagating mode at $k_x=−0.4$). Note that they are strongly (exponentially) localized to the waveguide, as they should be.\n",
    "\n",
    "* The next mode is the first leaky mode at $k_x=0.1$. As $k_x$ goes to zero, in fact, this mode actually becomes lossless, a peculiarity of symmetry related to an effect demonstrated in Phys. Rev. B. 63, 125107, 2001. However, at this non-zero $k_x$, the radiation loss is clearly visible.\n",
    "\n",
    "* The next mode is one of the many higher-order leaky modes visible in the band diagram; we arbitrarily chose the backwards-propagating mode at $k_x=0.3$, $\\omega=0.8838−0.0018i$ to plot. As can be seen from the field pattern, this mode has a very short wavelength in the material. This is short enough that it is worth checking how big the error introduced by the finite resolution is. By doubling the resolution to 40 pixels/unit, we found that this mode has shifted to $\\omegaω=0.8996−0.0021i$, or about a 2% error at the lower resolution.\n",
    "\n",
    "* Finally, we show one of the modes right along the edge of the light cone, at $k_x=0.25$, $\\omega=0.2506$. This mode is clearly not localized to the waveguide, and is just propagating through the air parallel to the waveguide — i.e. it is really part of the continuum of extended modes and its discreteness is an artifact of the finite cell and imperfect boundary conditions. For light propagating completely parallel to the boundary, PML is not very effective, so the imaginary part of ω is only -0.0008 for this field."
   ]
  }
 ],
 "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.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}