single_reflection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "644aabf8",
   "metadata": {},
   "source": [
    "## Single reflection event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "56c6ff8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from cycler import cycler\n",
    "import geo.caustic_frame as caustic_frame\n",
    "import monte_carlo_simulation as mcs\n",
    "from bandstructure.delafossite_bandstructure import delafossite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fbcf51f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = delafossite(40)\n",
    "bar_frame = caustic_frame.Frame('geo/bar.dxf')\n",
    "field = 1.0\n",
    "phi = 0\n",
    "n_inject = 1\n",
    "np.random.seed(42)\n",
    "simulation = mcs.Simulation(bar_frame, k, phi, field, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf3f8d18",
   "metadata": {},
   "source": [
    "# Horizontal Wall test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6b23e3d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial: [[22, 1], [0.25, -4.95]]\n",
      "projected: [(23, 1), -0.6390626455220065, -5.469684944405711]\n",
      "step_coords: [(0.25, -4.95), (-0.6390626455220065, -5.469684944405711)]\n",
      "intersection: (((0.5, -5.0), (-0.5, -5.0), 1), 0.1644613814393469, -4.9999999999495355, 0.0990800447124366)\n",
      "n_f_int: (22, 0.9037878613037109)\n",
      "final: (36, 0.5942767775058475) 0.6881446709423142 -5.316722294628131\n"
     ]
    }
   ],
   "source": [
    "edge = bar_frame.edges[4] \n",
    "initial_position = [0.25, -4.95]\n",
    "#initial_position = [0.8, 3]\n",
    "#initial_n_f = [22, 0.7]\n",
    "initial_n_f = [22, 1]\n",
    "print(f'initial: {[initial_n_f, initial_position]}')\n",
    "\n",
    "projected_n_f, projected_x, projected_y = simulation._update_position(initial_n_f, initial_position[0], initial_position[1])\n",
    "print(f'projected: {[projected_n_f, projected_x, projected_y]}')\n",
    "\n",
    "step_coords = ([(initial_position[0], initial_position[1]), (projected_x, projected_y)])\n",
    "print(f'step_coords: {step_coords}')\n",
    "intersection = simulation._get_sorted_intersections(step_coords)[0]\n",
    "print(f'intersection: {intersection}')\n",
    "intersection_x = intersection[1]\n",
    "intersection_y = intersection[2]\n",
    "n_f_int = simulation._get_n_f_intersection(initial_n_f, [(initial_position[0], initial_position[1]), (intersection_x, intersection_y), (projected_x, projected_y)])\n",
    "print(f'n_f_int: {n_f_int}')\n",
    "\n",
    "final_n_f = simulation._specular(n_f_int, edge)\n",
    "post_n_f, final_x, final_y = simulation._update_position(final_n_f, intersection_x, intersection_y)\n",
    "print(f'final: {final_n_f} {final_x} {final_y}')\n",
    "\n",
    "initial_r = (simulation._bandstructure.r[0][initial_n_f[0]], simulation._bandstructure.r[1][initial_n_f[0]]) + (1-initial_n_f[1])*simulation._bandstructure.dr[:, initial_n_f[0]]\n",
    "int_r = (simulation._bandstructure.r[0][n_f_int[0]], simulation._bandstructure.r[1][n_f_int[0]]) + (1-n_f_int[1])*simulation._bandstructure.dr[:, n_f_int[0]]\n",
    "ref_r = (simulation._bandstructure.r[0][final_n_f[0]], simulation._bandstructure.r[1][final_n_f[0]]) + (1-final_n_f[1])*simulation._bandstructure.dr[:, final_n_f[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "52223132",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[22, 1]\n",
      "0.09621213869628906\n",
      "0.9037878613037109\n",
      "0.6041\n"
     ]
    }
   ],
   "source": [
    "# Trying to figure out what get_n_f_intersection should be\n",
    "\n",
    "initial_n_f\n",
    "initial_position[0], initial_position[1]\n",
    "intersection_x, intersection_y\n",
    "projected_x, projected_y\n",
    "\n",
    "\n",
    "root =  np.sqrt(((intersection_x - initial_position[0])**2 + (intersection_y - initial_position[1])**2) /\n",
    "                                 ((projected_x - initial_position[0])**2 + (projected_y - initial_position[1])**2))\n",
    "f = initial_n_f[1] * (1-root)\n",
    "\n",
    "print(initial_n_f)\n",
    "print(root)\n",
    "print(f)\n",
    "print(0.7-0.137*0.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6585015e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAErCAYAAADaNFkkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABHNklEQVR4nO3deXhU9fXH8ffJigFkC5tACDuCKGJAFHAB3BfUti7gbqVVa7VqWxV/alXaqlWwbi0qrlTbWvddFhdUQFAE2QNCBFkSNllDlvP7Y24wYgKZZDIzST6v55knM/feufcMZCZnvsv5mrsjIiIiUpslxDoAERERkeqmhEdERERqPSU8IiIiUusp4REREZFaTwmPiIiI1HpJsQ5ARKQmS09P98zMzFiHISLArFmz8ty9eVn7lPCIiFRBZmYmM2fOjHUYIgKY2Yry9qlLS0RERGo9JTwiIiJS6ynhERERkVpPCY+IiIjUekp4REREpNZTwiMiIiK1XtwlPGZ2vZm5maWXs/8iM1sS3C4qtf0wM5trZtlm9nczs+hFLSIiIvEsrhIeM2sHHA/klLO/KXAbcDjQD7jNzJoEux8FLge6BLcTqz1gERERqRHiKuEBxgB/ALyc/ScA77v7BnffCLwPnGhmrYH93X2auzvwDHBGNAIWERGR+Bc3lZbNbBiwyt2/2ktvVBvg21KPVwbb2gT399xe1nVGAiMB6tevf1j37t2rGHn8WrRoEQDdunWLcSQi+7a3kvAilTFrxUYemZLN7G83ReR8B7dtxJXHdqZvZtOInE+iK6oJj5lNBFqVsWsUcDOh7qxq5e7jgHEAWVlZXptLwh9zzDEAfPDBBzGNQ6Qi9lYSXqSi3J1Pl67nocnZfLZsPU3SkjmhZyuSEqs2rLOo2Hlv3lp+8Y/P6NehKb85tjODuqSj4aI1R1QTHncfWtZ2M+sFdABKWnfaAl+YWT93X1Pq0FXAMaUetwU+CLa33WP7qogFLiIicc3dmbxwHQ9ODrXotGiYyi2nHMjwwzNIS4nMn7pbTy3i+Rk5jPtoGReOn8EhbRtx1bGdGXpgSxISlPjEu7jo0nL3uUCLksdmthzIcve8PQ59F/hzqYHKxwM3ufsGM/vezPoD04ELgQerP3IREYmlomLn7a9X8/CUpSxY/T1tm+zHXWccxM8Pa0u95MSIXmu/lEQuHdiBEf0zeOmLVTz6wVJGPjuLbi0bcuWxnTj14ANIVOITt+Ii4dkbM8sCfu3uvwwSmzuBz4Pdd7j7huD+lcBTwH7A28FNRERqoYKiYl75chWPfriUZbnb6NS8Pvf94hBO730AyYnVOx8nNSmR8/pl8IvD2vLGnNU8PCWba16YzZj3F3PFMZ0489C2pCTF25wgicuEx90zS92fCfyy1OPxwPgynjMTOCga8YmISGzsLCjiv7NW8o8PlrJq0w4ObL0/Dw/vw4kHtYp660pSYgJnHNqG0w85gPfmr+GhKdn88X9zeWDiEn51dCfO6dsu4q1MUnlxmfCIiIiUtn1XIf+aHho/s25LPn0yGnPnGT05tluLmA8cTkgwTjyoNSf0bMWHi3N5eEo2t702jwcnZ3P5oA6M6N+eBqn6cxtr+h8QEZG49f3OAp7+ZDnjP/mGjdsLOLJTM8ae25sjOjaLeaKzJzPjmG4tOKZbC6YvW89DU7L5y9sLeeSDpVwyIJNLjuxAo7TkWIdZZynhERGRuOPuvPbVd9z5xgLytuYzpHsLrhrcmT4ZTfb95DhweMdmHN6xGbO/3cRDk7MZO3EJz362glGnHMiZh7aJu2StLlDCIyIicWV53jb+79Wv+XhJHoe0bcT4i7M4uG3jWIdVKb3bNebxi7L4etVm/u/Vr7nuP1/x4qyV3HXGQXRs3iDW4dUpGkYuIlKKmTU2sxfNbKGZLTCzI2IdU12RX1jEg5OWcPzYj/gyZxN3DOvJS1cOqLHJTmkHtWnE/359JHedcRBzV23mxLEfM+b9xewsKIp1aHWGWnhERH7sAeAdd/+5maUAabEOqC6Yvmw9N788l6W52zilV2tuPa0HLfevF+uwIiohwTi/f3uO79mSu95YwAOTlvD6V99x15kHcWSn9FiHV+uphUdEJGBmjYCjgCcA3H2Xu2+KaVC13IZtu/j9f7/inHHTyC8s5slL+vLwiD61LtkprUXDevz9vEN55tJ+FLkz/LHpXPfv2azfmh/r0Go1tfCIiPygA5ALPGlmhwCzgGvcfVvpg0ovQpyRkRH1IGsDd+fFWSv581sL2LKzkCuO6cRvB3dhv5S6U7fmqK7Neffao3hocjb//Ggpkxau46aTunN2VjstVVEN1MIjIvKDJKAP8Ki7HwpsA27c8yB3H+fuWe6e1by5FngPV/a6rZw7bhq/f3EOHZs34M3fDuKPJ3avU8lOiXrJidxwQjfevmYQ3Vo15MaX5nLOuM9YvHZLrEOrdZTwiIj8YCWw0t2nB49fJJQASSVNmDCBzMxMEhISyGjfnhE3/o2THviIBau/5y9n9eK/vzqCbq0axjrMmOvcoiH/Htmfe35+MEvWbeXkBz7mnncWsmOXBjVHihIeEZGAu68BvjWzbsGmIcD8GIZUo02YMIGRI0eyYsUK3J1vc3J4/r5b6Pz9bCbfcAzn9ctQ100pZsbZWe2YdN3RDOvdhkc+WMrxYz/kg0XrYh1araCER0Tkx64GJpjZHKA38OfYhlNzjRo1iu3bt/9omxfms+D1caQ3SI1RVPGvWYNU7jv7EJ6/vD/JiQlc/OTnXDXhC1Zt2hHr0Go0DVoWESnF3WcDWbGOo6bbWVDEipycMvfllLNdfuyITs14+5pB/PPDZTw8JZtJC9dyxdGd+dXRHbUoaSWohUdERCLG3Xlv3hqOH/MRiQ3Lri2jmW0Vl5qUyG+HdGHS9UczpHtLxkxczJD7PuTtuatx91iHV6Mo4RERkYjIXreFC8fPYOSzs0hNSuDmW+8gLe3HdRvT0tIYPXp0jCKsudo2SePhEX14/vL+NKyXxBUTvmDE49NZtEazuSpKXVoiIlIl3+8s4IGJS3j60+Xsl5LIbaf14Pz+7UlOPJpurRoyatQocnJyyMjIYPTo0YwYMSLWIddYR3RqxhtXD+RfM3K4773FnPz3j7mgf3t+N7SrVmLfByU8IiJSKcXFzn9nfcs97yxiw/ZdnNs3gxuO70qzUgOSR4wYoQQnwpISE7jwiExOO/gA7nt/Ec98tpxXZ6/i9yd055y+7UjUzLcyKeEREZGwzVqxkT+9Po85KzdzWPsmPH16Pw5q0yjWYdUpTeqncNcZvRjerz23vz6Pm1+ey4TpK7j99J70zWwa6/DijsbwiIjIXpUuHtguI4OTrv4zP3v0U9Z+v5MHzu3Ni78+QslODPU4YH/+PbI/D553KBu27eIX//iM3z7/Jas3axp7aUp4RESkXHsWD1z57be8+4876F88n8nXH8Ow3m0wUxdKrJkZpx1yAJOuP5qrB3fmnXlrGPy3D3l4SjY7C1StGeIs4TGz683MzewncxnNrLeZfWZm88xsjpmdU2rfU2b2jZnNDm69oxq4iEgtdfPNN5dZPHDafx6mfqpGRcSbtJQkrj++G5OuO5qjuqZz77uLOH7MR7w5R9PY4ybhMbN2wPFAeRWptgMXuntP4ERgrJk1LrX/9+7eO7jNrtZgRUTqgFkrNpKT822Z+1Q8ML61a5rGPy/I4rnLDme/5ESu+tcXnPnIp3y+fEOsQ4uZuEl4gDHAH4AyU1B3X+zuS4L73wHrAC1TLCISYcvztnHFc7P42aOfktKo7I9ZFQ+sGQZ2SeetawZxz88OZvXmHfziH58x8pmZLM3dGuvQoi4uEh4zGwascvevKnh8PyAFWFpq8+igq2uMmZW7SIuZjTSzmWY2Mzc3t2qBi4jUIhu27eL21+Yx9P4P+XBxLtcO7cKjD/xNxQNruMQE4+y+7ZhywzFcf1xXPsnO4/gxH/F/r3xN3tb8WIcXNVHrgDWziUCrMnaNAm4m1J1VkfO0Bp4FLnL34mDzTcAaQknQOOCPwB1lPd/dxwXHkJWVVbc7NEVECK17Nf6Tb3h0ylK27SrknL4Z/G5oF1rsXw/oSmpSgooH1gJpKUlcPaQL5/bL4IFJi/nXjBxe/nIVvz66I5cN7Mh+KbV7fS6L9SAmM+sFTCI0RgegLfAd0M/d1+xx7P7AB8Cf3f3Fcs53DHCDu5+6r2tnZWX5zJkzKx17vDvmmGMA+OCDD2Iah0hFmNksd69xi3bW5M+R4mLn5S9Xcd97i/hu806GdG/BjSd1p0vLhrEOTaIge91W7n5nIe/PX0ur/etx3fFd+VmftjW6cOHePkdiPsTe3ecCLUoem9lyIMvd80ofZ2YpwMvAM3smO2bW2t1XW2hu5BnA19Udt4hITTZ1SR5/fmsB81d/T682jbjv7N4c0alZrMOSKOrcogGPXZjFjG82MPqtBfzhxTmMn/oNN518IEd3rX1DZGOe8OyNmWUBv3b3XwJnA0cBzczs4uCQi4MZWRPMrDlgwGzg19GPVkQk/i1c8z1/eWshHy7OpU3j/Xjg3N6cdvABJNTgb/VSNf06NOWVK4/kzbmrueedRVw0fgaDuqRz40nd6XlA7SkoGXcJj7tnlro/E/hlcP854LlynjM4KsGJiNRQazbv5P73F/HirJU0SE3i5pO7c+ERmdRLrt3jNqRizIxTDz6A43q05LlpOTw4eQmnPjiVMw9tww3Hd+OAxvvFOsQqi7uER0REImdbfiH//GgZj320jKJi59IBHfjN4M40TkuJdWgSh1KTErlsYAd+3qctj3yQzZOfLufNOau5fFBHfn1MJxrU4GKTNTdyEREpV1Gx879ZK/nbe4tYtyWfUw5uzR9P6E5Gs7R9P1nqvEZpydx08oGc37899767iIemZPPC599y3XFdOTurLUmJcVHVJixKeEREapmpS/K46835LFyzhUMzGvPo+X04rL1Wz5bwtWuaxt/PO5RLBmQy+s0F3PzyXJ7+dDk3n1LzBjYr4RERqSWy123hz28tZPLCdbRtsh8Pnncopx7cWot7SpUdmtGE//76CN75eg1/eXshF42fwVFdmzPq5APp1qpmlDFQwiMiUsPlbc1n7MTFPD/jW9KSE7nppO5cdKQGJEtkmRkn9WrN4ANb8OxnK/j7pCWc9MBHnNM3g+uO60rzhuUuchAXlPCIiNRQOwuKePKT5TwyJZvtBUWMODyDa4Z0oVmD+P7DIzVbalIivxzUkZ/1acvfJy/h2c9W8NrsVVx5bGcuG9ghbhNtJTwiIjWMu/P6nNXc/fZCVm3awdADW3DjSQfSuUWDWIcmdUiT+incdlpPLujfnr++vZB7313EhGkr+P2J3Rh2SJu4q+2khEdEpAaZtWIDd76xgNnfbqJH6/259+cHc2Tn9FiHJXVYx+YNGHdhFtOWrWf0mwv43b+/4slPlnPLKT3o1yF+Bssr4RERqQE2by/g5pfn8ubc1bTcP5V7f34wZ9XwdY+kdunfsRmvXjWAV2av4t53F3H2Pz/jhJ4t+ctZB9O0fuzrPinhERHZg5klAjOBVRVZiLi6bdq+i/OfmM7iNVv53dCuXH5UB9JS9PEt8SchwTirT1tOOqg1T0xdxt8nZzP8sWn86/L+MU96al7lIBGR6ncNsCDWQUAo2RnxeCjZ+ecFh3HN0C5KdiTu7ZeSyG8Gd+GJi7L4Jm8bwx+bxoZtu2IakxIeEZFSzKwtcArweKxjKUl2lqzbyj8vPIxju7eIdUgiYRnUpTlPXNR3d9Kzfmt+zGJRwiMi8mNjgT8AxbEMYuO2XQx/LJTsjLvgMI7tpmRHaqaBXdJLJT3TyYtR0qOER0QkYGanAuvcfdY+jhtpZjPNbGZubm7E49iwbRfDH59Odu5WHrswi2OU7EgNN7BLOuMv7suKDaGWnlgkPUp4RER+MAA43cyWAy8Ag83suT0Pcvdx7p7l7lnNm0d2PaEN23Yx/LFpLMvdyuMXZtW49YpEyjOgczrjL+pLzobtMUl6lPCIiATc/SZ3b+vumcC5wGR3Pz9a1y9Jdr7J28bjF2VxlJIdqWWO7Bxq6cnZsJ3zxk0jd0v0kh4lPCIicWD91vzdyc4TF/VlUBclO1I7HdkpnScv7sfKjTsY/lj0kh4lPCIiZXD3D6JVgydvaz7DH5u+O9kZ2EWVk6V2O6JTM568pC8rN+7gvMemsW7Lzmq/phIeEZEYygtadlZs2Mb4i5XsSN3Rv2MznrqkL99t2sF546ax7vvqTXriLuExs+vNzM2szHe9mRWZ2ezg9lqp7R3MbLqZZZvZv80s9nWsRUT2oiTZydmwnfEX9WWA1sSSOubwjs146pJ+rN68k3Mfq96kJ64SHjNrBxwP5OzlsB3u3ju4nV5q+93AGHfvDGwELqvGUEVEqiR3Sz7njQuSnYv7agFQqbP6dWjKU5f0Y001Jz1xlfAAYwgV/PJwnmRmBgwGXgw2PQ2cEdHIREQi6KoJX7By4w6evLgfR3ZSsiN1W78OTXn60n6s3byTy5/daxmsSoubhMfMhhFaqO+rfRxaLyj4Nc3Mzgi2NQM2uXth8Hgl0KaaQhURqZLcLfnMWL6Bq47txBGdmsU6HJG40DezKb8d0oWvvt3E6s07In7+qCY8ZjbRzL4u4zYMuBm4tQKnae/uWcBwYKyZdQozhmqtkCoisi+fLs0D0NRzkT2UvCc+yV4f8XNHNeFx96HuftCeN2AZ0AH4Kqhw2hb4wsxalXGOVcHPZcAHwKHAeqCxmZUsIdwWWFVODNVWIVVEpCI+yc5j/3pJHNSmUaxDEYkr3Vs1pFn9FD7Jzov4ueOiS8vd57p7C3fPDCqcrgT6uPua0seZWRMzSw3upxMqAz/f3R2YAvw8OPQi4NWovQARkQpyd6YuyePITukkJliswxGJKwkJxpGd05manUfoT3sEzx3Rs1UDM8sys8eDhwcCM83sK0IJzl/dfX6w74/AdWaWTWhMzxPRj1ZEZO+Wr9/Od5t3MkD1dkTKNLBzM3K35LNk3daInjdp34dEX9DKU3J/JvDL4P6nQK9ynrMM6BeN+EREKmtq0FQ/UNPQRcpUUo9q6pI8urZsGLHzxn0Lj4hIbfJpdh5tGu9HZrO0WIciEpfaNkkjs1na7sH9kaKER0QkSoqKnU+XrmdA52aEyoeJSFkGdE5n2rINFBQVR+ycSnhERKJk3neb2byjQEtIiOzDwM7pbM0vZM7KTRE7pxIeEZEoKRm/o8rKInt3RKdmmMHUJZGrx6OER0QkSj7JzqN7q4Y0b5ga61BE4lrjtBR6tWkU0Xo8SnhERKJgZ0ERny/fqNlZIhU0oHM6X+RsZFt+4b4PrgAlPCIiUTBz+UZ2FRar/o5IBQ3snE5hsTPjmw0ROZ8SHhGRKJianUdyotEvs2msQxGpEQ5r34TUpITdY9+qSgmPiEgUfJKdx6EZTaifGpf1XkXiTr3kRPpmNo3YOB4lPCIi1Wzjtl18/d1mjd8RCdOAzuksXLOFdVt2VvlcSnhERKrZtGXrcUf1d0TCVPIl4bOlVZ+eroRHRKSa5W3NByCjqZaTEAlHRrAES97WXVU+lxIeEZFqlpAQWkai2D3GkYjULEXFofdMYgRWYlHCIyJSzZKChKfkw1tEKmZ3wpNY9XRFCY+ISDVLMCU8IpXxQwtP1Zt4lPCIiFSzpEQlPCKVURR0A5e0klZF2AmPmdU3s8QqX1lEJM6YWTszm2Jm881snpldE4nzlrTwFCrhEQlLUVHoPZMQjYTHzBLMbLiZvWlm64CFwOrgA+FeM+tc5ShEROJDIXC9u/cA+gNXmVmPqp40KSH0UatByyLhiXYLzxSgE3AT0Mrd27l7C2AgMA2428zOr3IkIiIx5u6r3f2L4P4WYAHQpqrnLRlvWVikhEckHEXFxUBkWngqUuN8qLsX7LnR3TcA/wP+Z2bJVY5ERCSOmFkmcCgwvYx9I4GRABkZGfs8V6JaeEQqpSiU70Snhack2TGzLDN72cy+MLM5ZjbXzOaUPqaqzOx6M3Mz+0k5UjM71sxml7rtNLMzgn1Pmdk3pfb1jkQ8IlI3mVkDQl/ornX37/fc7+7j3D3L3bOaN2++z/OVtPBo0LJIeEreMwkRmKUVzip2E4DfA3OB4ipfeQ9m1g44Hsgpa7+7TwF6B8c2BbKB90od8nt3fzHScYlI3RK0WP8PmODuL0XinBq0LFI5u6elR6lLq0Suu79W5SuWbwzwB+DVChz7c+Btd99ejfGISB1jZgY8ASxw9/sjdV4NWhapnFhNS7/NzB43s/PM7KySW5UjAMxsGLDK3b+q4FPOBZ7fY9vooKttjJml7uVaI81sppnNzM3NrWzIIlI7DQAuAAaX6iI/uaonTdCgZZFKifag5RKXAN2BZH7o0nKgQk2+ZjYRaFXGrlHAzYS6sypyntZAL+DdUptvAtYAKcA44I/AHWU9393HBceQlZWlTx8R2c3dpwIRWLXnx9TCI1I5kRy0HE7C09fdu1X2Qu4+tKztZtYL6AB8FWpNpi3whZn1c/c1ZTzlbODl0gOl3X11cDffzJ4EbqhsnCIikbZ7WrrG8IiEpbCkhSfKS0t8GokCXHty97nu3sLdM909E1gJ9Ckn2QE4jz26s4JWn5L+9zOAryMdp4hIZe2elq6ERyQsxSUtPBFYLj2cFp7+wGwz+wbIJ9Ts6+5+cJWjKIeZZQG/dvdfBo8zgXbAh3scOsHMmgcxzQZ+XV0xiYiEK1GztEQqJZItPOEkPCdW+WoVELTylNyfCfyy1OPllFH11N0HRyM2EZHKKJlSWzIAU0QqpmTcW1Snpbv7ij23mVmHKkcgIlLLlTTHq4VHJDwFRZGbll7hhMfMbt1jUyJwIaEBxyIiUo6SD2tNSxcJT8l7JjkxnCHHZQunS2tbqfvJwFHAf6ocgYhILVfyYV1QpC4tkXCUjOGJ6qBld7+v9GMzuweYUeUIRERquZKER11aIuEp6dJKiUALT1XO0ARYW+UIRERquZJvp2rhEQlPyXsmqi08ZjaXUGVlCE3/zgTWl2yvzunpIiI1WXJCSZeWWnhEwlFYkvAkRHcMz6lVvpqISB20e5aWWnhEwlKwe9ByFFp4zMw85CfT0ksfU+VIREQqwMzqAzvdvSjWsVSUpqWLVM4Pg5ajM4ZnipldbWYZpTeaWYqZDTazp4GLqhyJiEgZzCzBzIab2Ztmtg5YCKw2s/lmdq+ZdY51jPvyQ5eWWnhEwhHtOjwnApcCzweFBjcB9QjV4XkPGOvuX1Y5EhGRsk0BJgI3AV+7ezGAmTUFjgXuNrOX3f25GMa4VwkJRmKCqQ6PSJiiWofH3XcCjwCPmFkykA7scPdNVb66iMi+DXX3gj03uvsG4H/A/4LPpriWlGBq4REJU0FRMQkW5aUlAIIPndVVvqpERe/evWMdgkiVlSQ7wWLCo4D2hD67di9gXFZCFG+SExM0S0skTAXFxREZvwNhJjxSs4wdOzbWIYhE0gTg98BcoMY1lSQl2u4BmCJSMYVFTnIEWndACY+I1By57v5arIOorKQEtfCIhKuwKAYtPGaWCvyMUMHB3c9z9zsiEolIDTBhAowaBTk5kJEBo0fDiBGxjqrOuM3MHgcmAfklG939pdiFVHHJiaY6PCJhKij2iNTggfBaeF4FNgOzKPVhI1JXTJgAI0fC9u2hxytWhB6Dkp4ouQToTmjx4pLMwYEakvAkaNCySJgKCosjMkMLwkt42rr7iRG5qkgNNGrUD8lOie3bQ9uV8ERFX3fvFusgKisp0ShQ4UGRsBQWe0TW0YLwFg/91Mx6ReSqIjVQTk542yXiPjWzHrEOorKSEowijeERCUthsUdkHS0Ir4VnIHCxmX1DqEtr95TQiEQiEucyMkLdWGVtl6joD8zWZ5CIVEY4Cc9J1RYFYGa3A5cDucGmm939rTKOOxF4gFCl58fd/a/B9g7AC0AzQuOMLnD3XdUZs9Qto0f/eAwPQFpaaLtERY3uUneHCH1RFakzEgzcI9MyWuG3X7B4aGPgtODWeG8LilbSGHfvHdzKSnYSgYcJJV89gPNKNXHfHTy/M7ARuCzCsUktNmHuBDLHZpLwpwQyx2YyYe6EnxwzYgSMGwft24NZ6Oe4cRq/Ey3uvmLPG+F1y8dUsTuG1lkWCYcBkRr6Fs609GsItcCUzIh4zszGufuDkQmlQvoB2e6+LIjpBWCYmS0ABgPDg+OeBm4HHo1ibFJDTZg7gZGvj2R7QajpZsXmFYx8PTT9akSvH2czI0YowYkVM7t1j02JwIVAhxiEEzYnlCiLSMUlmOFEuYWHUIvJ4e5+q7vfSqg//fKIRPGD35jZHDMbb2ZNytjfBvi21OOVwbZmwCZ3L9xju8g+jZo0aneyU2J7wXZGTRoVo4ikHNtK3XYBfYH/RPoiZnaimS0ys2wzuzFS53UHU8YjEh6DSBUoD2cMjwFFpR4XBdsqfgKziUCrMnaNItQacyehL0J3AvcRWqU9osxsJDASIEOjTQXI2Vz2NKvytktsuPt9pR+b2T3AjEheo1S3+XGEvjh9bmavufv8qp7b3dWhJRKmSHYDh5PwPAlMN7OXg8dnAE+EczF3H1qR48zsMeCNMnatAtqVetw22LYeaGxmSUErT8n2smIYB4wDyMrK0hxRIaNRBis2/3Q4WkYjJcRxrgmwNsLnLLPbHKh6wkNoAKaIVFysBi3fT6jFZUNwu8Tdx0YkCsDMWpd6eCbwdRmHfQ50MbMOZpYCnAu85qF/jSnAz4PjLiJUGVpkn0YPGU1actqPtqUlpzF6iKZfxRMzmxt0ec8xs7nAcqBnyfYIXaa8bvM9YxlpZjPNbGZubu6eu8tU7K4uLZEwmcVg0DKAu88iNOW7OtxjZr0JfRFaDvwKwMwOIDT9/GR3LzSz3wDvEhqwON7d5wXP/yPwgpndBXxJmK1PUneVDEweNWkUOZtzyGiUwegho38yYFli7tRYB1CiMi3F7mGOARARjMgNWt5nwmNmU919oJltgR9dtaTo1/6RCMTdLyhn+3fAyaUevwX8ZMp60ATdLxKxSN0zotcIJThxyszMQ8otg2GRazopr9u8yjRoWSR8CQmh905EzrWvA9x9YPCzobvvX+rWMFLJjojIXkwxs6vN7EeDqswsxcwGm9nThLqxI6HMbvNInNjdNS1dJGwWsS6tCo/hMbO7K7JNRCTCTiQ0K/R5M/vOzOab2TJgCXAeMNbdn4rEhYJJDyXd5guA/5TqNq/auVGXlki4Ql8SotSlVcpxhMbJlHZSGdtERCLG3XcCjwCPmFkykA7scPdN1XS9MrvNq37eUBE1Eam40CytyJyrImN4rgCuBDruMROiIfBJZMIQEdk3dy8AVsc6jsooVpeWSNgMozhCGU9FWnj+BbwN/AUoXXV0i7tviEgUIiK1nJaWEAmfWaQ6tCqQ8Lj7ZmAzob5yEZGYM7NUd8+PdRzh0CwtkfAlmEVvlpaZTQ1+bjGz74PblpLHkQlDRCQsj5jZcbEOIhxaWkKkciLVpVXZaekNNS1dRGLF3S8DMszs72aWHut4KqLInUStLSESlsQEozhC89LDmZb+CzNrGNy/xcxeMrNDIxKFiEgYzOwEoAPQGXjMzM6IbUT7ll9QTGpShT9yRQRITUogvzAyy6WH8+77P3ffYmYDgaGElm74R0SiEBHZCzNLMLObS21qTWhpmZPd/Uzg2BiFViHuTn5hEalJibEORaRGSU1KpLDYKSyqetITTsJTFPw8BRjn7m8CKVWOQERkH9y9mFJrabn7UyUrmgf+L/pRVVxhsVPsUC9ZLTwi4Sh5z+yKcsKzysz+CZwDvGVmqWE+X0SkKuaY2W1m9pPPHXeP6wkUOwtC3xfVwiMSnpJu4J0F0U14ziZUbv2EoMJpU+D3VY5ARKRimhJa2+o7M3vVzO40s1/EOqiKKBmDoBYekfDUSw59ScgvLNrHkftW4aUl3H27mS0FTggGDH7s7u9VOQIRkQpw97MhVIMH6An0AvoB/41lXBWhFh6RyklNjlwLT4UTHjO7BrgceCnY9JyZjXP3B6schYhIBQUFB78IbjVCSQtPqlp4RMJS8iUhqi08wGXA4e6+DXavlP4ZoIRHRGQv8oNvp2rhEQlPSTdwfpTH8Bg/zNQiuK8qWiIi+7Az+HaqFh6R8JR8SSjpFq6KcFp4ngSmm9nLhBKdYYRq8YiIyF6UfDutpxYekbDsbuGJQPHBcAYt329mHwADCS1eeom7f1nlCEREajm18IhUTkxaeMysHnAMMAgoBpLMbIG776xyFCIitdgPY3iU8IiEo+Q9E4kWnnDefc8Qmgr6d+AhoAfwbJUjAMzsdjNbZWazg9vJZRzTzsymmNl8M5sXzBqr8PNFRGKlZIZJSU0REamYH+rwRLFLCzjI3XuUejzFzOZXOYIfjHH3v+1lfyFwvbt/ESxiOsvM3nf3+RV8vohITKiFR6Ryfqi0XPUurXDefV+YWf+SB2Z2ODCzyhFUkLuvdvcvgvtbgAVAm2hdX0SkstTCI1I5qRFs4Qkn4TkM+NTMlpvZckI1ePqa2Vwzm1PlSOA3ZjbHzMabWZO9HWhmmcChwPRwn29mI81sppnNzM3NjUDYIiJ7t7vwoFp4RMLywxie6E5LP7EqFzKziUCrMnaNAh4F7iQ0++tO4D7g0nLO0wD4H3BtqQUDK/x8dx8HjAPIysrySr4cEZEK09ISIpUTycVDw5mWvqIqF3L3oRU5zsweA94oZ18yoWRngruXLHGBu6+tyPNFRGLBLFSjtdj1HUskHEXFofdMQgTKHMdF+6qZtS718Ezg6zKOMUKFDhe4+/3hPl9EJFaapKUAsHH7rhhHIlKzbNpRAPzwHqqKuEh4gHtKjQU6FvgdgJkdYGZvBccMAC4ABpcx/bzM54uIxIMmackAbNxWEONIRGqWTcGXhMbBe6gqwik8eDXwnLtvrPJV9+DuF5Sz/Tvg5OD+VMpZu6u854tIOSZMgFGjICcHMjJg9GgYMSLWUcWUmd0LnAbsApYSqia/KRLnbhx8O92kFh6RsGzcHpsWnpbA52b2HzM70Uo6pUWkZpkwAUaOhBUrwD30c+TI0Pa67X1C9cYOBhYDN0XqxE3rhz6sNyjhEQnLhm2h90zJe6gqKpzwuPstQBdC42guBpaY2Z/NrFOVoxCR6Bk1CrZv//G27dtD2+swd3/P3QuDh9OAtpE69+4ure3q0hIJRyS7tMIaw+PuDqwJboVAE+BFM7unypGISHTk5IS3vW66FHi7vJ3h1vPa3aW1TS08IuGISZeWmV1jZrOAe4BPgF7ufgWhgoQ/q3IkIhIdGRnhba9FzGyimX1dxm1YqWNGEfpCV24fn7uPc/csd89q3rz5Pq+bkpRA/ZREtfCIhGnj9l2kJCaQllL1GlbhFB5sCpy1Zz0edy82s1OrHImIRMfo0aExO6W7tdLSQttruX3VAzOzi4FTgSFBi3bENKmfomnpImHauG0XTeonE4lhw+GM4bmtvOKD7r6gypGISHSMGAHjxkH79mAW+jlunGZpmZ0I/AE43d237+v4cDVJU8IjEq6N2wsi0p0F4bXwiEhtMWJEnU9wyvAQkAq8H3ybnObuv47UyRunJatLSyRMm7bvisiAZVDCIyICgLt3rs7zN0lLIWdDxBuORGq1jdsL6NKiQUTOFS+VlkVEarUmacm7a4qISMVs3LZr9yzHqlLCIyISBU3qp7BlZyGFRVVf9VmkLnB3Nu0ooGn9yHRpKeEREYmCkoGXJYshisjefb+zkKJij9igZSU8IiJRUDLwUutpiVTMD1WWlfCIiNQYJd9SN2jFdJEKKRnz1iRCs7SU8IiIREHJ4oerN++IcSQiNcPqzTuB0Pi3SFDCIyISBV1bNqRd0/14ZMpSDVwW2YeiYufhKdkc0KgePVrvH5FzKuEREYmClKQERp18IIvWbuH5GVqoVWRv/jvzW+Z99z03nnwg9ZKrvo4WKOGJWxPWriXzs89I+OADMj/7jAlr18Y6JBGpohN6tuKIjs247/3FGrwsUo7vdxZw77uLyGrfhNMObh2x8yrhiUMT1q5l5KJFrMjPx4EV+fmMXLRISY9IDWdm3HpaD77fUcDYiUtiHY5IXHpocjYbtu/ittN6RmTR0BJKeOLQqGXL2F784z7+7cXFjFq2LEYRiUikHNh6f4YfnsGz01awZO2WWIcjEleW5W7lyU++4ezD2tGrbaOInjsuEh4zu93MVpnZ7OB2cjnHLTezucExM0ttb2pm75vZkuBnk+hFH3k5+flhbReRmuW647pRPyWRO96Yj7vHOhyRuDH6zQWkJiVywwndIn7uuEh4AmPcvXdwe2svxx0bHJNVatuNwCR37wJMCh7XWBmpqWFtF5GapWn9FK4d2pWPl+QxacG6WIcjEhc+WLSOSQvXcfXgzjRvGPm/d/GU8FTFMODp4P7TwBmxC6XqRnfsSFrCj/9r0hISGN2xY4wiEpFIu+CI9nRqXp+73pxPfmFRrMMRiamComLufGM+mc3SuHhAZrVcI54Snt+Y2RwzG7+XLikH3jOzWWY2stT2lu6+Ori/BmhZrZFWsxEtWzKuWzfap6ZiQPvUVMZ168aIljX6ZYlIKcmJCdx6Wk+Wr9/OU58sj3U4IjH17GcrWJq7jVtO6UFqUmSmoe8pqVrOWgYzmwi0KmPXKOBR4E5CCc2dwH3ApWUcO9DdV5lZC+B9M1vo7h+VPsDd3czK7RQPEqWRABkZGZV6LdEwomVLJTgitdzRXZszpHsLHpyczVl92lZLM75IvNuwbRdjJy5mUJd0hhzYotquE7UWHncf6u4HlXF71d3XunuRuxcDjwH9yjnHquDnOuDlUsetNbPWAMHPcjvF3X2cu2e5e1bz5s0j+RJFRMI26pQDyS8s4m/vLop1KCIxcf/7i9i2q4hbT+0R0Wnoe4qLLq2SZCVwJvB1GcfUN7OGJfeB40sd9xpwUXD/IuDV6otWRCRyOjZvwMVHZvKfWd8yd+XmWIcjElULVn/Pv6bncEH/9nRp2bBarxUXCQ9wTzDdfA5wLPA7ADM7wMxKZmy1BKaa2VfADOBNd38n2PdX4DgzWwIMDR6LiNQIVw/pQtO0FP70+jxNU5c6w9254/X57L9fMtcO7VLt14vaGJ69cfcLytn+HXBycH8ZcEg5x60HhlRbgCIi1Wj/esn8/oRu3PjSXF6fs5rTDzkg1iGJVLt3563hs2XruXNYTxqnRWZF9L2JlxYeEZE67RdZ7eh5wP785a0FbNymdbakdtu8vYC73lxAt5YNOa9fdCYQKeEREYkDiQnGn07vSe6WfI4f+xGTF2rtPKmdPlycy/FjP2T15p3cfnpPkhKjk4oo4RERiRNZmU155aoBNE1L4dKnZvLHF+ewZWdBrMMSiYht+YXc/PJcLho/g4b1knn5yiM5olOzqF0/LsbwiIhIyEFtGvHa1QMY8/4Sxn20lKnZefztF4dE9Q+DSKTN+GYDN/z3K77duJ3LB3Xg+uO7US+5egoMlkctPCIicSY1KZEbT+rOf399BMmJxnmPTeNPr89jZ4GWoJCaZWdBEX9+awHnjPsMgH+PPIJRp/SIerIDauEREYlbh7VvylvXDOLutxfy5CfL+XBxLvef3Zve7RrHOjSRfZq7cjPX/Wc2S9ZtZcThGdx88oHUT41d2qEWHhGRUszsejNzM0uPdSwAaSlJ/GnYQTx32eHs3FXEWY98wt/eXcSuwuJYhyZSpoKiYsa8v5gzHvmE73cW8PSl/Rh9Zq+YJjughEdEZDcza0eointOrGPZ08Au6bzzu6M4q09bHpqSzRkPf8LCNd/HOiyRH1mydgtnPfIpD0xawumHHMB71x7N0V3jYxknJTwiIj8YA/yB0ELGcWf/esn87ReH8NiFWazbspPTHpzKIx9kU1Qcl+FKHVJU7Iz7aCmnPDiVVZt28I/z+zDmnN40SkuOdWi7aQyPiAhgZsOAVe7+1b4WMDSzkcBIgIyM6BRNK+24Hi05rP3R3PLKXO55ZxET56/lvrN70yG9ftRjEclZv50b/vsVM5Zv4PgeLfnzWb1Ib5Aa67B+QgmPiNQZZjYRaFXGrlHAzYS6s/bJ3ccB4wCysrJi0rzStH4KDw/vw2tffcetr87jpAc+4vcndOfCI9qTHKVCblK3FRQV8/yMHP769kISE4z7zz6EMw9tU60rnleFEh4RqTPcfWhZ282sF9ABKGndaQt8YWb93H1NFEMMi5kxrHcb+ndsxh//N4c735jPM58t5/rju3Fqr9YkJMTnHx6p2dydt+au4b73FrEsbxuDuqRz988O5oDG+8U6tL1SwiMidZ67zwValDw2s+VAlrvnxSyoMLTcvx5PXtyXKYvWcc87i/jt81/yzw+X8scTuzOoS3rcfuOWmueT7Dzufmchc1ZupmvLBjx2YRZDD2xRI37HlPCIiNQCZsbg7i05umsLXvtqFfe9t5gLx8/giI7N+ONJ3VW7R6pk7srN3PPuQj5ekscBjepx788P5qw+bUmsQa2ISnhERPbg7pmxjqGyEhOMMw9ty8m9WvOv6Tk8NDk0hf2kg1pxwwnd6NS8QaxDlBrkm7xt/O29Rbw5ZzVN0pK55ZQDOb9/+5hUSq4qJTwiIrVQalIilwzowC+y2vHYR8t4/ONlvDd/LWdnteWaIV1p1aherEOUOLbu+508MGkJ//78W5ITE7h6cGcuP6oj+9eLn2nm4VLCIyJSizVITeJ3x3XlgiPa89DkbCZMX8FLX6zi4gGZXHl057iqkyKxt3lHAeM+Wsr4qcspKCpm+OEZ/GZwZ1o0rPkJshIeEZE6IL1BKref3pPLBnbg/vcXM+6jZTw/PYcrjunMJQMya2QXhUTOzoIinvlsOY98sJRN2ws4/ZADuP74rrRvVntqOynhERGpQ9o1TWPMOb0ZeVRH7nlnIXe/s5CnPv2Ga4d25ReHtSVJNXzqlKJi53+zVjJm4mJWb97JUV2b84cTunFQm0axDi3ilPCIiNRBB7benycv6cf0Zeu5+52F3PTSXMZP/YZbTu0RN2sfSfWauiSPO9+Yz6K1WzikXWPuO/sQjuwUF2vmVou4SOXN7HYzW2Vms4PbyWUc063U/tlm9r2ZXVvR54uIyE8d3rEZ/7viSP55wWHsKirmovEzuOTJGWSv2xrr0KSaLMvdyi+f/pzzn5jO9oJCHhnRh1euPLJWJzsQXy08Y9z9b+XtdPdFQG8AM0sEVgEvV/T5IiJSNjPjhJ6tOKZbc575dAV/n7SEE8Z+xAX923PNkC40qZ8S6xAlAjZvL+Dvk5fw9KfLqZecyI0ndefiI+vO+K14SnjCMQRY6u4rYh2IiEhtkZqUyOVHdeSsPm24//3FPPPZcl7+chXXDu3C+f21RldNVVhUzL9m5DDm/cVs2lHAuX3bcd1x3WjeMP4W+KxO8fTb+xszm2Nm482syT6OPRd4vjLPN7ORZjbTzGbm5uZWOWgRkdqmWYNURp/Zi7euGUSvNo340+vzOWHsR0xZuA73mKyVKpX04eJcTnrgY259dR7dWjXkjasH8pezDq5zyQ6AReuXdx+rFE8D8gAH7gRau/ul5ZwnBfgO6Onua4NtLSv6/NKysrJ85syZlXg1IhJpZjbL3bNiHUe4avvniLszeeE67npzAd/kbeOors255ZQD6dqyYaxDk73IXreV0W/OZ8qiXNo3S+Pmkw/k+B4ta8SaV1Wxt8+RqHVplbdK8Z7M7DHgjb0cchLwRUmyE5x79/0KPF9ERCrIzBhyYEsGdWnOs9NW8MDExZz0wMeMODyDa4d2pekr/4VRoyAnBzIyYPRoGDEi1mHXWZu272LsxCU8N20F+yUncvPJ3bnoyExSk+rGOJ29iYsxPGbW2t1XBw/PBL7ey+HnsUd3VpjPFxGRMKUkJXDZwA6ceWgbxk5czITpOeQ/8yx3vfkgyfk7QgetWAEjR4buK+mJqoKiYiZMW8GYiUvYsrOAc/tlcN1xXUlvUPe6rsoTFwkPcI+Z9SbUJbUc+BWAmR0APO7uJweP6wPHlezf1/NFRCSymtZP4Y5hB3F+//Y07XHZD8lOie3bQy0+SniiZsqiddz1xnyW5m5jQOdm/N+pPejeav9YhxV34iLhcfcLytn+HXByqcfbgGYVfb6IiFSPri0b4hvXlr0zJye6wdRR327Yzh1vzOf9+WvpkF6fxy/MYsiBLWr9OJ3KiouER0REah7LyAh1Y+3h++atSSkoqjP1XaItv7CIxz5axkNTsjGMP57YncsGdiAlKZ4mXscf/euIiEjljB4NaWk/2pSfUo9b+p0Xmsa+aF2MAqu9Pl6Sy0ljP+Zv7y3m2G4tmHj90VxxTCclOxWgfyEREamcESNg3Dho3x7MoH17Usc/ztn330iiGZc8+Tm/enYmqzbt2Pe5ZK9Wb97BVRO+4IInZlDszlOX9OXR8w+jTeP9Yh1ajRG1OjzxqLbXzxCpSVSHp3bJLyzi8Y+/4cHJSzCMq4d05pcDO6olIkwFRcWMn/oND0xaQlGxc9WxnRl5VEd1F5YjLurwiIhI3ZGalMhVx3ZmWO8DuOP1+dzzziJenLWSO4cdxIDOtXuRykj5bOl6bn31a5as28qQ7i247bSeZDRL2/cTpUxKtUVEpNq0bZLGuAuzePLivhQWOSMen85v/vUFazbvjHVocWvdlp1c+8KXnPfYNHYUFPH4hVk8cXFfJTtVpBYeERGpdsd2b8ERnZrx6AdLefTDpUxZuI7fHdeVi5d/StL/3aJKzYQW+XzmsxWMeX8x+YXFXD24M1ce05n9UtR9FQlKeEREJCrqJSfyu+O6clafNtz22jzm3PMIhe8+RFJBfuiAOlypedaKDdzyyjwWrP6eQV3S+dPpPenYvEGsw6pVlPCIiEhUtW9Wnycv7svOUT+jXkmyU6KOVWrO3ZLP3e8s5MVZK2m1fz0eGdGHkw5qpeKB1UAJj4iIRJ2Zsd+aVWXu85wcavuf+8KiYp6dtoL731/MzoIifnV0R347uAv1U/VnubroX1ZERGKjnErN6xq3YOnSPI7sVDtnc01ftp7bXpvHwjVbGNQlndtO60nnFuq+qm6apSUiEjCzq81soZnNM7N7Yh1PrVdGpeaievvx6PGXMvyx6Vz1ry9Yvbn2FC1c+/1OrnnhS84ZN40tOwv5x/l9eObSfkp2okQtPCIigJkdCwwDDnH3fDNrEeuYar2ScTqjRu2epZU4ejQ3nn0ujT9cyqMfLGXygnVcPaQzlw3sQGpSzZyttKuwmKc+/YYHJi6hoMg1+ypGVGlZFVJF4kKsKy2b2X+Ace4+MZzn6XOk+pReDbxjen1uO70nR3dtHuuwwjJ1SR63vfY1S3O3Mbh7C249tQeZ6fVjHVattbfPEXVpiYiEdAUGmdl0M/vQzPrGOqC6rl3TNB67MIsnL+mLAxeNn8HIZ2by7YbtsQ5tn1Zt2sGVE2Zx/hPTKShynrgoi/EX91WyE0Pq0hKROsPMJgKtytg1itDnYVOgP9AX+I+ZdfQymsHNbCQwEiAjI6P6AhYAju3WgiM7NePxj7/hocnZDL3/Q648pjO/Ojr+1pQqWUPsocnZOM71x3Xlcq19FRfUpaWmaJG4EAddWu8Ad7v7lODxUqC/u+fu7Xn6HImu7zbtYPRbC3hzzmoymqZx66k9GNqjZazDAmDKonX86bV5LF+/nRN7tuKWUw+kbRMtBxFNWjxURGTfXgGOBaaYWVcgBciLaUTyEwc03o+Hh/dheL88bnttHr98ZiZdWzZgUJfmDOySzuEdmpKWEp0/bTt2FTFj+QY+XpzL1Ow8Fq7ZQsfm9Xnm0n4cVcPGGtUFSnhERELGA+PN7GtgF3BRWd1ZEh8GdE7n7WsG8cKMHN6dt5Znp63gianfkJKYQJ/2jRnUpTmDuqTT84BGJCZEpoxhcbEzf/X3fLQkl6lL8pi5fCO7iopJSUwgK7MJt5/Wg+GHtyclScNj45G6tNQULRIXYt2lVVn6HIkPOwuK+Hz5Bj5eksfHS/JYsPp7ABqnJTOgczqDOqczsEt62F1MqzbtYOqSXD5eksenS9ezYdsuALq3asigLukM7NKcfplNNcU8TtSYLi0zuxq4CigC3nT3P5RxzInAA0Ai8Li7/zXY3gF4AWgGzAIucPdd0YpdRERip15yYtCqE+pKyt2Sz6dL8/hocR5Ts3N5c85qADqk1w8lKp3TOaJTMxrWS/7RebbmFzJt6Xo+XpLLx9l5LMvdBkCLhqkc0y3UajSgczotGtaL7guUKoubhKciRb/MLBF4GDgOWAl8bmavuft84G5gjLu/YGb/AC4DHo3eKxARkXjRvGEqw3q3YVjvNrg72eu28vGSPKZm5/HirJU889kKEhOMQ9s1Zv/9knF3tuYX8mXOJgqLnf2SEzm8Y1OG98tgUJfmdG3ZQAt61nBxk/AAVwB/dfd8AHdfV8Yx/YBsd18GYGYvAMPMbAEwGBgeHPc0cDtKeERE6jwzo0vLhnRp2ZBLB3ZgV2ExX+RsZOqSPD7OzmPywtCfm4PbNmLkUR0Z2CWdw9o3qbGVnaVs8ZTwlBT9Gg3sBG5w98/3OKYN8G2pxyuBwwl1Y21y98JS29uUdZHS9TOA/GCAYm2WTu2faaLXWDt0i3UAUjekJCXQv2Mz+ndsxg0ndGPWig0UO/TNbBrr0KQaRTXhiVTRr6pw93HAuCCemTVxkGQ49Bprh7ryGmMdg9RNh7VXolMXRDXhcfeh5e0zsyuAl4IEZ4aZFRP6Vlu66NcqoF2px22DbeuBxmaWFLTylGwXERERiau1tF4hVPSLvRT9+hzoYmYdzCwFOBd4LUiSpgA/D467CHg1GkGLiIhI/IunhGc80DEYU/MCQdEvMzvAzN4CCFpvfgO8CywA/uPu84Ln/xG4zsyyCY3peaIC1xwX6RcRh/Qaawe9RhGRKqjThQdFRKpKhQdF4sfeCg/GUwuPiIiISLVQwiMiIiK1Xp1KeMzsF2Y2z8yKzazcKb5mdqKZLTKzbDO7MZoxVpWZNTWz981sSfCzSTnHFZnZ7OD2WrTjrIx9/b+YWaqZ/TvYP93MMmMQZpVU4DVebGa5pf7vfhmLOCvLzMab2bry6l9ZyN+D1z/HzPpEO0YRqZ3qVMIDfA2cBXxU3gGllq84CegBnGdmPaITXkTcCExy9y7ApOBxWXa4e+/gdnr0wqucCv6/XAZsdPfOwBhCy43UGGH87v271P/d41ENsuqeAk7cy/6TgC7BbSSqli4iEVKnEh53X+Dui/Zx2O7lK4LFR18gtMZXTTGM0NIaBD/PiF0oEVWR/5fSr/1FYIjVrMVvavrv3j65+0fAhr0cMgx4xkOmEaqv1To60YlIbRZPS0vEi/KWr6gpWrr76uD+GqBlOcfVCyrbFhJaw+yVaARXBRX5f9l9jLsXmtlmQiUKasqSDBX93fuZmR0FLAZ+5+7flnFMTVXWv0EbYHXZh8ferFmz8sxsRQUOjaflQeIllniJAxRLWeIlDqh4LO3L21HrEp69LV/h7rWiGOE+lujYLahjVF7dgfbuvsrMOgKTzWyuuy+NdKwSca8Dz7t7vpn9ilCL1uAYx1SnuXvzihwXT8uDxEss8RIHKJZ4jgMiE0utS3j2tnxFBZW3fEXc2McSHWvNrLW7rw66AspadR53XxX8XGZmHwCHAvGc8FTk/6XkmJVmlgQ0IrTsSE2xz9fo7qVfz+PAPVGIK5ri/v0nIjVTnRrDU0FlLl8R45jC8RqhpTWgnCU2zKyJmaUG99OBAcD8qEVYORX5fyn92n8OTI704rPVbJ+vcY/xLKcTqjhem7wGXBjM1uoPbC7VRSsiUmm1roVnb8zsTOBBoDnwppnNdvcTzOwA4HF3PzkY+1GyfEUiML7U8hU1wV8JrTR/GbACOBsgmIb/a3f/JXAg8M9ggdYEQmN44jrhKe//xczuAGa6+2uElhN5NlheZAOhhKHGqOBr/K2ZnU5o7NUG4OKYBVwJZvY8cAyQbmYrgduAZAB3/wfwFnAykA1sBy6JTaTVIp6WzoiXWOIlDlAsZYmXOCACsWhpCREREan11KUlIiIitZ4SHhEREan1lPCIiESRmV1tZguDZW5iPsvOzK43Mw8mMMTi+vcG/x5zzOxlM2scgxhivpyQmbUzsylmNj/43bgmFnGUiifRzL40szdiHEdjM3sx+B1ZYGZHVPZcSnhERKLEzI4lVE36EHfvCfwtxvG0A44HcmIYxvvAQe5+MKFimjdF8+JxtJxQIXC9u/cA+gNXxXhZo2uIj1mgDwDvuHt34BCqEJMSHhGR6LmC0KzIfAB3L7NOVhSNAf4AxGz2iru/5+6FwcNphGovRVNcLOni7qvd/Yvg/hZCf9jbRDsOADNrC5xCqNZXzJhZI+AoQjNwcfdd7r6psudTwiMiEj1dgUFmNt3MPjSzvrEKxMyGAavc/atYxVCGS4G3o3zN8pYziRkzyyRUDHZ6jEIYSygRLo7R9Ut0AHKBJ4PutcfNrH5lT1an6vCIiFS3fSz9kgQ0JdRl0ZdQzayO1VUgcx+x3EyoO6vaVWTJHzMbRahbZ0I0YopXZtYA+B9wrbt/H4Prnwqsc/dZZnZMtK+/hySgD3C1u083sweAG4H/q+zJROKSmX3q7kcGgxiHu/sje+6LwDX2A94BBrt7URXPlQJMDM5VuK/jpXbax9IvVwAvBQnOjKD4Zzqhb7FRi8XMehH69vyVmUGoG+kLM+vn7muiFUepeC4GTgWGxKA6etwsZ2JmyYSSnQnu/lIsYiBUef90MzsZqAfsb2bPufv5MYhlJbDS3Utaul4klPBUirq0JG6VSmgaA1eWs6+qLiX0B6hKyQ6E+peBScA5VY5KaqtXgGMBzKwrkEIMVqN297nu3sLdM909k9Aflj7Vkezsi5mdSKj75HR33x7t6xMnywlZKPN8Aljg7vdH+/ol3P0md28b/F6cS2iJnlgkOwS/j9+aWbdg0xCqsAySEh4Ji5n1DaaP1jOz+sH0yYP2OCYzmEI4IZhG+KKZpQX7rjOzr4PbtcG2+mb2ppl9FWw/J9i+NTjlX4FOZjbbzO7dY19558wMrv1YEON7QWvOnkYQrDcWPOfrUue9wcxu3+M1PWVmi4PXNtTMPjGzJWbWL3jaK8E5RcoyHugY/J69AFxUw9Z7qw4PAQ2B94P3+D+iefGgNbZkSZcFwH9itJzQAOACYHDw7zA7aGWp664GJpjZHKA38OfKnkhLS0jYzOwuQk2d+xFqbvzLHvszgW+Age7+iZmNJ5SVTwGeIjR+wQgNyDsf6Aic6O6XB89v5O6bzWyruzcIzveGux9U6hol+w4r55wbCa3HlOXus83sP8Br7v5cqXOkADnu3qpU3LuvY2Y3AA3c/fZgXzahgYTzCH0r/Aq4jNAinpe4+xnBFNc17t68Kv/GIiISWWrhkcq4AzgOyALKK5z2rbt/Etx/DhgY3F52923uvhV4CRgEzAWOM7O7zWyQu28OI5byzgnwjbvPDu7PAjL3eG46sCmMa30TdAUUE0p6JgXfzueWnDvoGttlZg3DOK+IiFQzJTxSGc2ABoSaoeuVc8yeTYflNiW6+2JCI/HnAneZ2a2RCBLIL3W/iJ8O0t/BT+O3UveT93K+4lKPi/c4dyqwM6xIRUSkWinhkcr4J6FpgROAu8s5JqNUCfDhwFTgY+AMM0sLaimcCXxsZgcA24PupnsJJT+lbSGUXJWlzHNW5EW4+0Yg0cxKJz3tzay5mSUQKniVWJFzlTCzZkCeuxeE8zwREaleSngkLGZ2IVDg7v8iNJi4r5kNLuPQRYRKoy8AmgCPBlVEnwJmEBpr87i7fwn0IjRFdzZwG3BX6RO5+3rgk2BQ8r177CvvnBX1HqFusRLrgWeAmcDXwIVm1imM8x0LvBnG8SIiEgUatCwRV9Yg43hlZn2A37n7BZGI28xeAm4MuulERCROqIVH6rSghWhKMLuqSoJZX68o2RERiT9q4REREZFaTy08IiJSp5jZfhZavDUiLbtm9pGZaammOKeER0REapUKJDJaUqYOUsIjIiI1npn918z+aWbTgJv2cbiWlKmDlPCIiEht0AtY6+793f2u8g4KJhd0dPflFTxvZ+A+oHtwG06olMUNwM3BMV8DfSsZt0SJEh4REanRguKhTQkte7MvWlKmjlLCIyIiNV1PYHqw8jlBtfQnzaytmY03s9LLxGhJmTpKCY+IiNR0vYA5JQ/cPRfIIdQV9dvSS71oSZm6SwmPiIjUdD9KeMysAdARKHT3rWUcryVl6iAVHhQRkVojqIczDvgTcDbwubt/sMcxWlKmDlKhJBERqTWCcTyXBg/vLeeYL8xMS8rUMWrhERERkVpPY3hERESk1lPCIyIiIrWeEh4RERGp9ZTwiIiISK2nhEdERERqPSU8IiIiUusp4REREZFa7/8Bv9rlxN6/4nsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 648x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['figure.figsize'] = [9, 5]\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "x, y = [edge.xs, edge.ys]\n",
    "plt.plot(x, y, color='k')\n",
    "plt.xlabel(r'x position ($\\mathrm{\\mu m}$)')\n",
    "plt.ylabel(r'y position ($\\mathrm{\\mu m}$)')\n",
    "plt.plot([initial_position[0]], [initial_position[1]], 'bo')\n",
    "plt.plot(intersection_x, intersection_y, 'go')\n",
    "plt.plot(projected_x, projected_y, 'co')\n",
    "plt.plot(final_x, final_y, 'ro')\n",
    "ax1.set_xlim(-1, 1)\n",
    "ax1.set_ylim(-6, -4)\n",
    "ax1.set_aspect('equal')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(simulation._bandstructure.r[0], simulation._bandstructure.r[1])\n",
    "plt.xlabel(r'$r_x\\ (\\mathrm{\\mu m})$')\n",
    "plt.ylabel(r'$r_y\\ (\\mathrm{\\mu m})$')\n",
    "plt.plot(simulation._bandstructure.r[0][initial_n_f[0]], simulation._bandstructure.r[1][initial_n_f[0]], 'ko')\n",
    "plt.plot(simulation._bandstructure.r[0][initial_n_f[0]+1], simulation._bandstructure.r[1][initial_n_f[0]+1], 'ko')\n",
    "plt.plot(simulation._bandstructure.r[0][final_n_f[0]], simulation._bandstructure.r[1][final_n_f[0]], 'ro')\n",
    "plt.plot(simulation._bandstructure.r[0][final_n_f[0]+1], simulation._bandstructure.r[1][final_n_f[0]+1], 'ro')\n",
    "ax2.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "eff2ba64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAFFCAYAAAAgpt3YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA68UlEQVR4nO3dd3yV5fnH8c+VBBJG2GGFDWEFkBG2CAgyZbgHLmpLba2jttY6QEatW1HrLP4cFavWwRaQjYooKAIJe++hgAgykty/P87BpjGBBJLznCfn+3698uKMO+d8OSe5cp37foY55xAREREJd1FeBxARERHJCzUtIiIi4gtqWkRERMQX1LSIiIiIL6hpEREREV9Q0yIiIiK+ENKmxczKmdn7ZrbazFaZWcds95uZPWtm681suZm1DmU+EQl/qiMikSsmxM/3DDDdOXe5mRUHSma7vy+QFPxqD7wY/FdE5BTVEZEIFbKZFjMrC1wAvArgnDvhnDuYbdgg4E0X8AVQzsyqhSqjiIQ31RGRyBbK5aG6wD7gNTP7xszGmVmpbGMSgW1Zrm8P3iYiAqojIhEtlMtDMUBr4Dbn3GIzewb4KzD8bB7MzIYBwwBKlSrVpnHjxgUWVETO3tKlS/c75xIK6eFVR0QiQG51JJRNy3Zgu3NucfD6+wSKTVY7gJpZrtcI3vYLzrlXgFcAUlJS3JIlSwo2rYicFTPbUogPrzoiEgFyqyMhWx5yzu0GtplZo+BNPYC0bMMmATcEt/7vABxyzu0KVUYRCW+qIyKRLdR7D90GjA9u8b8RGGpmtwA4514CpgH9gPXAUWBoiPOJSPhTHRGJUCFtWpxzy4CUbDe/lOV+B9waykwi4i+qIyKRS0fEFREREV9Q0yIiIiK+oKZFREREfEFNi4iIiPiCmhYRERHxBTUtIiIi4gtqWkRERMQX1LSIiIiIL6hpEZECk5npvI4gIkWYmhYRKRDr9/7IlS8v8jrGWTtw5ASBg+mKSLhS0yIi5+REeibPzV5Hv2cWsm7vj17HOWvbD/7EkHGL2bz/iNdRRCQXalpE5Kx9s/UAA577lCc/WUuv5CrMuqur15HOWmK5EqzYfojeYxfw8vwNpGdkeh1JRLJR0yIi+XbkeDqjJqdy6Yuf88Oxk4y7IYV/XNuahPhYr6OdtQqlivPJXV25oGECD3+8mkHPf8bKHYe8jiUiWahpEZF8mbdmL72eXsBrn23m+g61mfnHC+jZtIrXsQpE1bJxvHJ9G14c0pq9h48z6PnPePjjVRw7meF1NBEBYrwOICL+8P2RE4yZksZH3+ygfkIp3r+lIyl1Kngdq8CZGX2bV6NT/Ur8fdoqXp6/kekrd/Pwpc3pVL+S1/FEIppmWkTktJxzTFy2g55PzWfK8p3c3iOJaXd0KZINS1ZlSxbj0ctb8Pav2wNw7T8Xc8/7yzl09KTHyUQil2ZaRCRX2w8c5YEJK5m3Zh8ta5bj0cta0KhqvNexQqpTg0rMuPMCxs5axz8XbmTOmr2MHphMn2ZVMTOv44lEFDUtIvILGZmONxdt5vEZawB4cEBTbuhYh+ioyPwjHVcsmr/2bczFLapxzwfL+d34r+nVtAqjBzWjatk4r+OJRAwtD4nI/1i75zCXv/Q5oyan0bZOBWb+8QKGdq4bsQ1LVs0SyzLx1s7c27cx89fu46Kn5vP24q06ErBIiGimRUQAOJ6ewfNzN/DivPWUjo1h7FUtGdSyupZAsomJjuK3XevTO7kq9320gvs+WsHEZTt4+NLm1Eso7XU8kSJNMy0iwtIt39P/2U95dvY6Lm5RnVl3dWVwq0Q1LKdRp1Ipxv+6PY9d1oJVu36gzzMLeX7uek7qoHQihUYzLSIR7Mfj6Tw2fTX/+mIL1cuW4PWhbenWqLLXsXzDzLiybU26NU5g1KQ0Hp+xhsnf7uSxy1vQokY5r+OJFDmaaRGJUHNW7+Gip+bzry+2cFOnOsz84wVqWM5S5fg4nh/Smleub8OBoycY/PxnPDQ1jaMn0r2OJlKkaKZFJMLs//E4oyanMfnbnTSqEs8LQ1rTqlZ5r2MVCb2Sq9KhfkUe+Xg1/1y4iempu/n7Jc3pkpTgdTSRIkEzLSIRwjnHB0u30/Op+cxYuZu7LmrI5NvOV8NSwMrEFePvlzTn3WEdKBYVxfWvfsmf3vuWA0dOeB1NxPc00yISAbZ9f5T7PlrBwnX7Saldnkcua06DypF1kLhQa1+vItPu6MI/5qznpfkbmL92Lw8OSObiFtW0gbPIWVLTIlKEZWQ6XvtsE0/OXEuUwZhByQxpX5soHXMlJOKKRfPn3o3o17waf/1wObf9+xsmfLODMYObUb1cCa/jifhOyJsWM9sMHAYygHTnXEq2+8sCbwG1gvmecM69FuqcIn63atcP/PWD5Xy7/RA9GlcuMn8o/VhDmlYvw0e/7/xzA9nr6QXc06eRGkiRfPJqpqW7c25/LvfdCqQ55waYWQKwxszGO+e0ICySB8dOZvDcnHW8PH8jZUsU47lrWhXFJQnf1ZDoKOPXXer9fFC64RNTmbBsJ49qqU4kz8JxQ1wHxFugwpYGvge036BIHize+B39nlnI83M3MLhVIrPu6sqA8yLuqLZhXUNqVijJm79qx5NXnMeGfT/S75lPeWbWOk6k66B0ImfixUyLA2aamQNeds69ku3+fwCTgJ1APHCVc06/zSKn8cOxkzzy8WreXryVmhVK8K+b2xXl3Wx9X0PMjMva1KBrowRGT07j6VlrmbpiJ49c1oLW2ptLJFdezLSc75xrDfQFbjWzC7Ld3xtYBlQHWgL/MLMy2R/EzIaZ2RIzW7Jv375CjiwSvmak7uaip+bzzpdb+U2Xusy484Ki3LBAAdUQ8L6OVCody7PXtOLVG1M4fCydy178nJGTUjlyPGwmhkTCSsibFufcjuC/e4GPgHbZhgwFPnQB64FNQOMcHucV51yKcy4lIaFIF2iRHO09fIzfj1/Kb/+1lPIlizPh1s7c378pJYsX7Z0CC6qGBB8jLOpIjyZV+OSurtzQoTZvLNpMr6cXMHfNXs/yiISrkDYtZlbKzOJPXQZ6ASuzDdsK9AiOqQI0AjaGMqdIOHPO8e5XW+n55HxmrdrL3b0bMfm28yPiXDdFuYaUjo1h1KBmvH9LR0oUj2boa19x5zvf8N2Px72OJhI2Qv2RrArwUXCjwBjgbefcdDO7BcA59xIwBnjdzFYABtxzmr0ERCLK5v1HuPfDFSza+B3t6lbgkUubUy+htNexQqnI15A2tSsw9fbzeWHuBl6Yt575a/cxYkBTBrfUWbdFzDnndYZzlpKS4pYsWeJ1DJFCk56RybhPN/H0J2spHh3Fvf2acHXbmmF5jA8zW5r92Cl+EI51ZO2ew9zzwXK+2XqQrg0TeOiSZtQoX9LrWCKFLrc6Eo67PItIFit3HGLQ85/xyMer6dowgVl/6sq17WuFZcMiBathlXjev6UTIwc05avN39Pr6QW89tkmMjL9/2FT5GwU7S32RHzspxMZjJ29lnELN1GhVHFeuq41fZpV8zqWhFh0lHFT57r0bFqFByasZNTkNCYu28mjl7WgUVUdlE4ii2ZaRMLQ5+v30+eZBbw8fyNXtKnBrD92VcMS4WqUL8lrN7XlmatbsvX7o1z83EKemrmG4+kZXkcTCRnNtIiEkUNHT/LQtDTeW7KdOhVL8vZv2tOpfiWvY0mYMDMGtUykS1ICY6ak8eyc9UxdsYtHL2tBSp0KXscTKXSaaREJA845pq3YRY+n5vPB1zu4pWt9pt95gRoWyVGFUsV5+qqWvD60LcdOZnL5S4sYPmElh4+d9DqaSKHSTIuIx3YfOsbwiSv5JG0PzRLL8PrQtjRLLOt1LPGBbo0qM/OPF/DkzLW89vkmZq3aw98GN6NHkypeRxMpFJppEQmR8eOhTh2Iigr8+9ZbjvGLt3DRU/NZuG4f9/VrzITfd1bDIvlSKjaGEQOa8uHvOlEmrhg3v7GEP7z9NfsO66B0UvRopkUkBMaPh2HD4OjRwPUtW+CmmzMp1+t7eg0sy8OXNqd2xVLehhRfa1WrPJNvO5+X52/guTnrWbhuPw/0b0LiiW3Mnz+fbt260bFjR69jipwTNS0iIXD//f9tWE7JOBGNLWnO+EnROtKpFIjiMVHc1iOJvs2rce+Hy7n9mXfZ9+4DkJlO8eLFmT17thoX8TUtD4mEwNatOd/+3Z4YNSxS4BpULs27wzrSpfReMtJPkpGRwfHjJ5g1e47X0UTOiZoWkUK254djlK54Isf7atUKcRiJGFFRxh9vuIS4uFgsKhoXFc2Hu8vq7NHia1oeEikkJzMyeePzzTz9yVpKdq7KsektOHn8v58TSpaEhx7yMKAUeR07dmTO7NnMmzeP8vVb8t62Egx97St6Na3CiAFNdR4j8R01LSKFYPHG7xgxMZU1ew7TvVECI+9I4tPpUdx/f2CpqFatQMMyZIjXSaWo69ix48/bsfwqPZNxn27kudnr6fnUfG67MIlfd6lLbEy0xylF8kZNi0gB2nv4GI9MW82H3+wgsVwJXrm+DRc1rYKZUXuImhTxVvGYKH7frQGDWiYyZnIaj89YwwdLtzNqUDJdkhK8jidyRmpaRApAekYmb32xhSdnruV4eiZ/6N6AW7s3oERxfYKV8JNYrgQvXd+GeWv2MnJSKte/+iX9m1fjgYubUK1sCa/jieRKTYvIOVq65XsemJDKql0/0CWpEqMGJlMvobTXsUTOqFujyky/syL/XLCRf8xdz9w1e7mjRxJDO9eleIz205Dwo6ZFJB/GrxjP/bPvZ+uhrSTG16RF/C2kbmhBtbJxvDikNX2aVdUuzOIrccWiua1HEoNbJTJqchoPf7ya/yzdzuhByTr3lYQdtdIieTR+xXiGTR7GlkNbcDi2H97KxztG0T45jVl3daVv82pqWMS3alYoybgbU3j1xhSOp2dw7T8Xc/u/v2HPD8e8jibyM3POeZ3hnMXHx7s2bdp4HUOKuC86fMHxEr88n0vsT7F0+KKDB4nC0/z585c651K8zpFfqiP/lekcOw8eY+fBnzCDGuVLUrVMHOrJJVRyqyOaaRHJg/QMx/G4nE9Al9vtIn4VZUaN8iVoUaMs8XHF2PLdEVbsOMQPx056HU0iXJHYpqVRo0bMmzfP6xhSBGVmOt5dso1Hp68mJnMV6ez7xZja5Wrr5y8Lvy6RqY7kzDnHJ2l7GDU5jR0Hf6Jbq0T+2q8xlePjvI4mRVhudUQzLSK5WLH9EJe8+Dn3friChlXiebjHw5Qs9r9HEC1ZrCQP9dBhbaXoMjN6JVdl1l1d+UP3BkxZvoseT8zn9c82kZ6R6XU8iTBFYqZFpCAdPHqCJ2auYfzirVQsFcvTV53H4JaJmHWkWrm4n/ceqlW2Fg/1eIghzXXEOCn6ShSP5s+9G3Fp60QenJTKyMlpvLtkO38bnEyb2hW8jicRokhsiJuSkuKWLFnidQzxucxMx/tfb+eRj1dz8OgJbuhYh7t6NaRMXDGvo/mKmflyQ1zVkbxzzvHxyt2MmZLGrkPHuKJNDe7p25hKpWO9jiZFRG51RDMtIkDqzkOMmJjK0i0HaFO7PGMGtadp9TJexxIJS2ZGv+bV6NowgefmrGfcwo3MSN3N3X0ac227WkRH+XO7Jgl/alokoh366SRPf7KWNxdtpnzJ4jx+eQsua12DKBVdkTMqFRvDX/s25vI2iQyfkMrwCSt576ttjBncjJY1y3kdT4ogNS0SkZxzfPTNDv4+bTXfHznOdR1q86eLGlG2pJaCRPKrQeV43v5NeyYv38XfpqRxyQufcXXbmvyld2PKlyrudTwpQtS0SMRZvfsHRkxI5cvN39OyZjleH9qWZollvY4l4mtmxsDzqnNh48o8M2st//fZZj5euZt7+jTmqpSamr2UAhHSpsXMNgOHgQwgPaeNbMysGzAWKAbsd851DV1CKcoOHzvJ2FnreP3zzZSJi+GRS5tzpYqp76iOhLfSsTHc378pl7epyfCJK7n3wxW889U2/jaoGc1r6MOBnBsvZlq6O+f253SHmZUDXgD6OOe2mlnlkCaTIsk5x6Rvd/LQ1FXs+/E4V7etxV96N9K0tb+pjoS5RlXjeXdYByYs28FDU1cz8PlPGdK+Fnf3aqxlWDlr4bY8dC3woXNuK4Bzbq/HecTn1u05zIiJqSza+B3NE8vyyg0p2kCw6FMdCRNmxiWtatCjSRWemhnY4H3ait2BjXe1wbuchVAfEdcBM81sqZkNy+H+hkB5M5sXHHNDbg9kZsPMbImZLdm375eHVpfIduR4Og9PW0XfZxaStusH/ja4GRNu7ayGpWhQHfGZMnHFGDkwmSm3daFupVL85f3lXPHyIlJ3HvI6mvhMSA8uZ2aJzrkdwenaT4DbnHMLstz/DyAF6AGUABYB/Z1za0/3uDoolJzinGPaisBBr3b/cIwrU2pwT5/GVNRBr0KmsA8upzrib5mZjg+CB3E8EDyI4x8vakjZEloykv8Ki4PLOed2BP/da2YfAe2ABVmGbAe+c84dAY6Y2QLgPOC0xUYEYMO+Hxk5KZWF6/bTtFoZnh/Smja1y3sdSwqY6oi/RUUZV6TUpFfTqjwxcw1vLNrMlOW7uK9fYy5plejbE25KaIRsecjMSplZ/KnLQC9gZbZhE4HzzSzGzEoC7YFVocoo/nT0RDqPTV9Nn7ELWLbtIKMGJjPpD53VsBRBqiNFR9mSxRgzuBmTbj2fGuVLcNd733LVy1+wevcPXkeTMBbKmZYqwEfBLjoGeNs5N93MbgFwzr3knFtlZtOB5UAmMM45l70giQCBpaAZqXsYMyWNHQd/4tLWidzbtwkJ8VoKKsJUR4qY5jXK8uHvOvHekm08On01/Z/9lJs61eHOnknE67xfko1OmCi+tHn/EUZOTmXemn00rhrP6EHNaFdXZ5oNBzphopytA0dO8NiMNbzz1VYSSsdyf/8mDDyvupaMIlBYbNMicq6OnczghbnreWn+RorHRDH84qbc2LE2MdGh3hFORApa+VLFefjS5lzVtibDJ6zkjneW8c6X2xg9KJmkKvFex5MwoKZFfGNW2h5GTk5l+4GfGNSyOvf1a0KVMnFexxKRAtayZjkm3NqZf3+5lcdnrKHvMwu5+fy63N4jiVKx+rMVyfTuS9jb9v1RRk1OZdaqvSRVLs2/f9OBjvUreh1LRApRdJRxXYfa9G1WlUenr+blBRuZuGwnwy9uSr/mVbVkFKHUtEjYOnYyg5fnb+SFeeuJjjLu69eYoZ3rUkxLQSIRo2LpWB67/DyualuL4RNWcuvbX9MlqRIjByZTP6G01/EkxNS0SFiau2YvIyelsuW7o/RvUY0H+jehWtkSXscSEY+0qV2eSX/ozPjFW3li5hr6jF3Ab7rU4w8XNqBkcf0pixR6pyWsbD9wlDFT0piRuod6CaV46+b2nJ9UyetYIhIGYqKjuLFTHfo1r8bDH6/ihXkbfl4y6p1cRUtGEUBNi4SF4+kZjFu4iefmrAPg7t6N+HWXusTGRHucTETCTUJ8LE9d2ZKrg0tGt7y1lG6NEhg5IJk6lUp5HU8KkZoW8dzCdft4cGIqG/cfoU9yVYYPaEpiOS0FicjptatbgSm3n8+bi7bw9Cdr6fX0Am7pWo/fd29AXDF94CmK1LSIZ3Yd+om/TVnF1BW7qF2xJK8PbUu3RpW9jiUiPlIsOoqbz6/LxS2q8fdpq3h2zno+/GYHIwck07NpFa/jSQFT0yIhdyI9k9c+28Qzs9eRken400UN+c0F9fTJSETOWpUycTxzdSuubluLERNX8us3l9CzSWUeHJBMzQolvY4nBURNi4TU5xv2M2JiKuv3/kjPJlV4cEBTFRQRKTAd61dk2h1deO2zTYydtY6eT83n1u4NGKYPRkWCmhYJiT0/HOOhqauY9O1OalYowas3ptCjiaZuRaTgFYuOYtgF9RlwXnX+NnUVT32ylg++3s7Igcl01xK0r6lpkUJ1MiOTNz7fzNhZ6ziRkcntPZL4fbf6+sQjIoWuWtkSPH9ta65uG9jYf+hrX9E7uQrDL25KjfKa4fUjNS1SaBZv/I4RE1NZs+ewdkcUEc90SUrg4zu7/HxYhflPzee2C5N0WAUfUtMiBW7v4WM8Mm01H36zg8RyJXj5+jb0aqoDP4mId2Jjorm1ewMGt0pkzOQ0Hp+xhg+WbmfUoGS6JCV4HU/ySE2LFJj0jEze+mILT85cy7H0DG7tXp8/dE+iRHF9khGR8JBYrgQvXd+GeWv28uCkVK5/9Uv6N6/GAxfrVCF+oKZFCsTSLd8zfEIqabt+oEtSJUYNTKaeTmYmImGqW6PKzLizIq8s2Mjzc9czd81e7uiRxNDOdSkeo5Oyhis1LXJO9v94nEc/Xs1/lm6nWtk4XhjSmr7NdNp4EQl/ccWiub1HEpe0SmTU5DQeDtay0YOS6VRf5zwLR2pa5KxkZDre/nIrj09fzdETGfy2az1uvzCJUrH6kRIRf6lZoSTjbkxh9qo9jJycyrX/XMzA86pzf/8mVCkT53U8yUJ/YSTflm07yPAJK1mx4xAd61VkzOBkGlSO9zqWiMg56dGkCp0bVOKFeRt4af4G5qzey509k7ixUx2KRWvJKByoaZE8O3DkBI/NWM07X20joXQsz17TigEtqmkpSESKjLhi0dx1UUMubZXIyMmp/G3qKv6zJLBk1L5eRa/jRTw1LXJGmZmOd5ds49Hpqzl8LJ2bO9fljp5JxMcV8zqaiEihqFOpFK/d1JaZaXsYPTmNq175gktbJXJvvyYkxMd6HS9iqWmR01qx/RAPTFzJt9sO0q5uBcYMakajqloKEpGiz8zonVyVLkmVeH7uel5ZsJFP0vbwp14Nua5DbWK0ZBRyalokR4eOnuTxmasZv3grFUvF8vRV5zG4ZaKWgkQk4pQsHsPdvRtzWesaPDgplZGT03hvyXbGDE6mTe0KXseLKGpa5H9kZjre/3o7j3y8moNHT3Bjxzr88aKGlC2hpSARiWz1Ekrz5q/a8fHK3YyZksZlLy7iijY1uKdvYyqV1pJRKKhpkZ+l7jzEiImpLN1ygDa1yzNmUHuaVi/jdSwRkbBhZvRrXo2uDRN4ds46Xl24iRmpu7m7T2OubVeL6CjNRhcmNS3CD8dO8tTMtby5aDPlSxbn8ctbcFnrGkTpl09EJEelYmO4t28TrmhTg+ETUhk+YSXvfbWNMYOb0bJmOa/jFVkh34rIzDab2QozW2ZmS04zrq2ZpZvZ5aHMF0mcc3z49XYufGI+byzazJD2tZnzp25ckVJTDYuELdUQCScNKsfz9m/a8+w1rdjzwzEueeEz7v1wOQeOnPA6WpHk1UxLd+fc/tzuNLNo4FFgZugiRZbVu39gxIRUvtz8PefVLMdrN7WleY2yXscSySvVEAkbZsbA86rTvVECz8xax2ufb+bjlbu5p09jrtKHwAIVrvtr3QZ8AOz1OkhRc/jYScZMSaP/s5+ybu9hHrm0OR/9rpMaFilqVEMk5OLjivHAxU2ZdnsXGlaJ594PV3DJi5+zYvshr6MVGV40LQ6YaWZLzWxY9jvNLBG4BHgx5MmKMOccE5ftoMeT8/m/zzZxZUpN5vypG1e3q6VPAeI3qiES1hpVjefdYR14+qrz2HHgJwY+/ykPTFjBoaMnvY7me14sD53vnNthZpWBT8xstXNuQZb7xwL3OOcyT3dMkGCxGgZQq1atwszre+v2HGbExFQWbfyO5olleeWGFG0oJn5WIDUEVEek8JgZl7SqwYWNq/D0J4EdHaat2M1f+zbmcu3ocNbMOefdk5uNBH50zj2R5bZNwKl3sxJwFBjmnJuQ2+OkpKS4JUty3R4vYh05ns6zs9fx6qebKBUbw929G3GNdsmTQmZmS51zKSF6rpEUQA0B1REpXNkPKTF6UDLJ1bUsn5vc6khIZ1rMrBQQ5Zw7HLzcCxiddYxzrm6W8a8DU85UbOR/OeeYtiJw8KPdPxzjypQa3NOnMRV18CPxOdUQ8avk6mX5z2878kHw4J0DnvuUGzrW4a5eDSmj87jlWaiXh6oAHwWnbGOAt51z083sFgDn3EshzlPkbNj3IyMnpbJw3X6aVivD80Na6TDTUpSohohvRUUZV6TUpFfTqjwxcw1vLNrMlOW7uK9fYy5ppdOk5IWny0MFRdO6cPREOv+Ys55/LtxIXLFo/tyrEUPa19IJvSTkQrk8VJBURyTU/ueEtHUqMHpwMo2r6ijkECbLQ1LwnHPMSN3DmClp7Dj4k06dLiLiE81rlOWj33Xi3SXbeHT6avo/+yk3darDnT2TiNeSUY7UtPjY5v1HGDk5lXlr9tGoSmAXu/b1KnodS0RE8igqyrimXS36JFflsRlr+L/PNjH5253c378JA8+rriWjbNS0+NCxkxm8MG8DL83fQPHoKB7o34QbO9WhmJaCRER8qXyp4jx8aXOualuT4RNWcsc7y3jny22MHpRMUpV4r+OFDTUtPjN71R5GTk5l2/c/Mahlde7r14QqZeK8jiUiIgWgZc1yTLi1M//+ciuPz1hD32cWcnOXutx+YRKlYvUnW6+AT2z7/iijJqcya9VeGlQuzdu/aU+n+pW8jiUiIgUsOsq4rkNt+jaryqPTV/Py/I1M/GYnwy9uSr/mVSN6yUhNS5g7djKDVxZs5Pm564mOMu7t25ihnetSPEZLQSIiRVnF0rE8dvl5XNW2Jg9MSOXWt7+mS1IlRg5Mpn5Caa/jeSLfTUvwgE7HnHMZhZBHspi7Zi8jJ6Wy5buj9G9RjQf6N6Fa2RJexxI5Z6ojInnXpnYFJv+hM299sYUnZ66lz9gF/KZLPf5wYQNKFo+suYcz/m/NLAq4GhgCtAWOA7Fmth+YCrzsnFtfqCkjzI6DPzF6ciozUvdQr1Ip/nVzO7okJXgdS+SsqY6InJuY6Chu6lyX/i2q8/C0VbwwbwMTl+1kxICm9GpaJWKWjPKyxjAXqA/cC1R1ztV0zlUGzge+AB41s+sKMWPEOJ6ewfNz19PjyXnMX7uPu3s34uM7u6hhkaJAdUSkACTEx/LUVS15d1gHSsfG8Nt/LWXo61+xef8Rr6OFxBmPiGtmxZxzpz2fdl7GFKaicCTLT9ftZ8TElWzcf4TeyVUYMSCZxHJaChL/yelIlqojIgXvZEYmb3y+mbGz1nEiI5Nbutbn993qE1cs2uto5+ysj4h7qoiYWQpwP1A7+H0WuNu18LLQ+N2uQz/xtymrmLpiF7UrluS1oW3p3qiy17FECpTqiEjBKxYdxa+71GPAedV5aOoqnp29jo++2c7IAcn0aFLF63iFIj9b8IwH7gZWAJmFEydynEjP5LXPNvHM7HVkZDruuqghwy6oVyQ6ZJHTUB0RKWBVysTx7DWtuLpdTUZMTOXmN5bQs0llHhyQTM0KJb2OV6Dy07Tsc85NKrQkEeTzDfsZMTGV9Xt/LLI/WCK5UB0RKSSd6ldi2u1dfv5A3POp+dzavUGR+kCcn6blQTMbB8wmsOU/AM65Dws8VRG154djPDR1FZO+3UnNCiUYd0MKPZsWzSk8kVyojogUouIxUfy2a30GtqzO36au4qlP1vLB19sZOTC5SGx6kJ+mZSjQGCjGf6d1HaBicwbZN5a6vUdSkdlYSiSfVEdEQqBa2RI8f21rrm67jwcnpjL0ta+KxE4e+Wla2jrnGhVakiLqy03fM2LiSlbvPkzXhgmMGphMnUqlvI4l4hXVEZEQ6pKUwMd3dmHcwk08N2cdPZ6cx20XJvHrLnWJjfHfB+f8HAv+czNrWmhJiph9h49z17vLuPLlRRw+ls7L17fh9aFt1bBIpFMdEQmx2Jhobu3egFl3daVrw4TAiRjHLuTTdfu9jpZv+Zlp6QAsM7NNBNaif95VsVCS+VR6RubPh1o+lp7Brd3rc2v3yDvUskguVEdEPFKjfElevj7l51PEXPfqYvo3r8YDF/vnFDH5+Uvap9BSFBFLt3zP8AmppO36IeJPaiWSC9UREY91b1SZjndW/PlkvHPX7OWOHkm+OBlvnpsW59yW7LeZWd2CjeNP3/14nEc+Xs1/lm6napk4XhjSmr7NIvv04SI5UR0RCQ9xxaK5vUcSl7RKZNTkVB4O/g0bPSiZTvUreR0vV3luWsxsRLabooEbgIgtOBmZjre/3Mrj01dz9EQGv+1aj9svTKJUrJaCRHKiOiISXmpWKMm4G9syK20PIyencu0/FzPwvOrc378JVcrEeR3vF/Lz1zXr2ZiKARcA7xVsHP9Ytu0gwyesZMWOQ3SsV5HRg5JJqhLvdSyRcKc6IhKGejatwvlJlXhh3gZemr+BOav3cmfPJG7sVIdi0eGzZHTGEybm+o2BU81/mdMJjUItlCc6O3DkBI/NWM07X20joXQs9/dvwsDzqmspSCQotxOd5TI2IuuISDjbvP8ID05KZf7afTSuGs/oQc1oV7dCSDOc9QkTT6M8sOccvt9XMjMd7y7ZxqPTV3P4WDo3d67LHT2TiI8r5nU0ET+LqDoi4gd1KpXi9aFtmZm2h9GT07jy5UVc2iqRe/s1ISE+1tNs+dmmZQWBI1dCYDfFOsB3p24vyrssrth+iAcmruTbbQdpV6cCowcn07hqGa9jifhOJNcRET8xM3onV6VLUiWen7ueVxZs5JO0PfypV0Ou61CbGI+WjPIz03JxoaUIU4eOnuSJmWt4a/EWKpYqzlNXnsclrRK1FCRy9iKujoj4WcniMdzduzGXtq7BgxNTGTk5jfeWbGfM4GTa1A7tkhHkoWkxM3MBv9hVMeuYgo3lrcxMx/tfb+eRj1dz8OgJbuxYhz9e1JCyJbQUJHI2IrGOiBQl9RNK86+b2zFtxW7GTEnjshcXcUWbGvy1b2Mqlg7dklFeZlrmmtkHwETn3NZTN5pZceB84EZgLvD6mR7IzDYDh4EMID37RjZmNgS4h8C08WHgd865b/P0PykgqTsPMWJiKku3HKB1rXKMubkdydXLhjKCSFEUUXVEpCgyM/q3qEa3Rgk8O2cdry7cxIzU3dzdpzHXtqtFdFThf+7IS9PSB/gV8O/gQaAOAnEEjq8wExjrnPsmH8/Z3TmX2wkPNgFdnXMHzKwv8ArQPh+PnS/j9+zh/o0b2Xr8ODWKx9L2YHG+WbibciWL89jlLbi8dQ2iQvAmiESAIltHRCJNqdgY7u3bhMtb12DExFSGT1jJe19tY8zgZrSsWa5QnztfuzybWTGgEvCTc+5gvp8s8Akp5TTFJuvY8sBK51zimcaeza6K4/fsYdiaNRzNzPzvc2Y4Bh4rzRvdW1K2pJaCRM7GmXZ5Lkp1RCTSOeeY9O1OHpq6in0/HqdLme+pc2ILfXv1oGPHjmf9uLnVkXxt/uucO+mc23U2hebUQwAzzWypmQ07w9ibgY/P8nnO6P6NG/+nYQFw0cayCulqWEQKUVGqIyKRzswY1DKR2X/qykUVD/HWiF8zZtSDdO1+IQs//azAny/U+yyd75xrDfQFbjWzC3IaZGbdCRSbe3J7IDMbZmZLzGzJvn378hVi16Gf2HLseI73bT2e8+0iEjbCoo6IyH/FxxWj7sktmMsAl8nJEye46e+vM2f1Hs72ILY5OaumxczOalNh59yO4L97gY+Adjk8dgtgHDDIOffdaR7rFedcinMuJSEhIU/Pf/REOk9/spbuT8wj5lhmjmNqxXp74ByRSOHXOiIiOevWrRuxxYsTHR1NbGxxytVrya9eX8IN//cla3YfLpDnONuZlhfM7KL8fIOZlTKz+FOXgV7AymxjagEfAtc759aeZbZfyMx0fPj1di58Yj7PzF5HjyZVeLJhA0pG/e9/v2RUFA/Vq1dQTysip+erOiIip9exY0dmz57NmDFjmDtnDovG/o7hFzfl220H6fvMAu7/aAXf/Xhuqxnncu6hm4HzgNF53CCuHoFPRRDYa+lt59xDZnYLgHPuJTMbB1wGnDqWwy92Z8zJ6TagW7L5e8ZMSePb7YdoUaMswy9uSts6gQPiZN17qFZsLA/Vq8eQKlXO9HQichr5PPeQL+qIiJy9A0dO8Mzsdfzriy2ULBbNbT0acGOnOsTGROf6PbnVkbNqWsysN9AFaA0cB95wzk3I9wMVkJyKzbbvj/LI9NVMXb6LKmViuadPYwa3TNQuzCKFLK9Nix/qiIgUnPV7D/PQ1FXMXbOP2hVLcm/fJvROrpLjUebPae8hM4sys/uy3FQN+D/nXD/n3CVA97P9TxS0H4+n89j01fR4aj6zV+3hjh5JzP1zNy7VMVdEPOWnOiIiBa9B5XheG9qON37VjuLRUdzy1lKu+ecXrNxxKM+PkeeZFjP73DnXKZf7yjjnfsjzsxawlJQUt/jLr3h/6TYen7GW/T8e55JWifylTyOqlS3hVSyRiHS6mZZwryOaaREJjfSMTP795Vae+mQtB386yRVtavDn3o2oHB8H5F5H8nPCxOVm9iAwxjn3P7veeFloAI4cT2fAc5+StusHWtcqx7gbUwr9qHwiclbCto6ISOjEREdxfcc6DGyZyHOz1/HGos1MXb6L33dvwM3n1839+/LxHBWArsDvzGwxsBxY7pz7z7lFP3cb9x+hwk8nee6aVlzcoprOwiwSvsK2johI6JUtUYwHLm7KkA61eXjaKh6fsYa3F2/NdXyemxbn3JXw87EVkoHmBI6P4HmxqVImjtl/6kpcsdy3RBYR74VzHRER79StVIpXbkjh8/X7GTN1Va7j8jPTAoBz7jjwdfArLFSOj1XDIuIj4VhHRMR7nRpUYspt5xNzZ873h/ow/iIiIiK5ij7Nnr5qWkRERMQX1LSIiIiIL6hpEREREV9Q0yIiIiK+oKZFREREfEFNi4iIiPiCmhYRERHxBTUtIiIi4gtqWkRERMQX1LSIiIiIL6hpEREREV9Q0yIiIiK+oKZFREREfEFNi4iIiPiCmhYRERHxBTUtIiIi4gtqWkRERMQX1LSIiIiIL6hpEREREV9Q0yIiIiK+EPKmxcw2m9kKM1tmZktyuN/M7FkzW29my82sdagzikj4Ug0RiVwxHj1vd+fc/lzu6wskBb/aAy8G/xUROUU1RCQChePy0CDgTRfwBVDOzKp5HUpEfEM1RKSI8qJpccBMM1tqZsNyuD8R2Jbl+vbgbSIioBoiErG8WB463zm3w8wqA5+Y2Wrn3IL8PkiwWA0DqFWrVkFnFJHwVSA1BFRHRPwm5DMtzrkdwX/3Ah8B7bIN2QHUzHK9RvC27I/zinMuxTmXkpCQUFhxRSTMFFQNCT6G6oiIj4S0aTGzUmYWf+oy0AtYmW3YJOCG4B4AHYBDzrldocwpIuFJNUQksoV6eagK8JGZnXrut51z083sFgDn3EvANKAfsB44CgwNcUYRCV+qISIRLKRNi3NuI3BeDre/lOWyA24NZS4R8QfVEJHIFo67PIuIiIj8gpoWERER8QU1LSIiIuILalpERETEF9S0iIiIiC+oaRERERFfUNMiIiIivqCmRURERHxBTYuIiIj4gpoWERER8QU1LSIiIuILalpERETEF9S0iIiIiC+oaRERERFfUNMiIiIivqCmRURERHxBTYuIiIj4gpoWERER8QU1LSIiIuILalpERETEF9S0iIiIiC+oaRERERFfUNMiIiIivqCmRURERHxBTYuIiIj4gpoWERER8QU1LSIiIuILalpERETEF0LetJhZtJl9Y2ZTcrivlpnNDd6/3Mz6hTqfiIQ/1RGRyOTFTMsdwKpc7nsAeM851wq4GnghZKlExE9UR0QiUEibFjOrAfQHxuUyxAFlgpfLAjtDkUtE/EN1RCRyxYT4+cYCfwHic7l/JDDTzG4DSgE9c3sgMxsGDAOoVatWgYYUkbA2FtURkYgUspkWM7sY2OucW3qaYdcArzvnagD9gH+ZWY4ZnXOvOOdSnHMpCQkJhZBYRMKN6ohIZAvl8lBnYKCZbQbeAS40s7eyjbkZeA/AObcIiAMqhTCjiIQ31RGRCBaypsU5d69zroZzrg6BjePmOOeuyzZsK9ADwMyaECg2+0KVUUTCm+qISGTz/DgtZjbazAYGr/4J+I2ZfQv8G7jJOee8SycifqA6IhIZQr0hLgDOuXnAvODlEVluTyMw/SsiclqqIyKRx/OZFhEREZG8UNMiIiIivqCmRURERHxBTYuIiIj4gpoWERER8QU1LSIiIuILalpERETEF9S0iIiIiC+oaRERERFfUNMiIiIivqCmRURERHxBTYuIiIj4gpoWERER8QU1LSIiIuILalpERETEF9S0iIiIiC+oaRERERFfUNMiIiIivqCmRURERHxBTYuIiIj4gpoWERER8QU1LSIiIuILalpERETEF9S0iIiIiC+oaRERERFfUNMiIiIivqCmRURERHwh5E2LmUWb2TdmNiWX+680szQzSzWzt0OdT0TCn+qISGSK8eA57wBWAWWy32FmScC9QGfn3AEzqxzqcCLiC6ojIhEopDMtZlYD6A+My2XIb4DnnXMHAJxze0OVTUT8QXVEJHKFenloLPAXIDOX+xsCDc3sMzP7wsz6hCyZiPjFWFRHRCJSyJoWM7sY2OucW3qaYTFAEtANuAb4p5mVy+XxhpnZEjNbsm/fvoKOKyJhSHVEJLKFcqalMzDQzDYD7wAXmtlb2cZsByY550465zYBawkUn19wzr3inEtxzqUkJCQUZm4RCR+qIyIRLGRNi3PuXudcDedcHeBqYI5z7rpswyYQ+HSEmVUiMM27MVQZRSS8qY6IRDbPj9NiZqPNbGDw6gzgOzNLA+YCdzvnvvMunYj4geqISGQw55zXGc5ZSkqKW7JkidcxRAQws6XOuRSvc+SX6ohI+Mitjng+0yIiIiKSF2paRERExBfUtIiIiIgvqGkRERERX1DTIiIiIr6gpkVERER8QU2LiIiI+IKaFhEREfEFNS0iIiLiC2paRERExBfUtIiIiIgvqGkRERERX1DTIiIiIr5QJM7ybGb7gC3n+DCVgP0FEKcgKEvuwimPsuSskXMu3usQ+WVmh4E1XucICqf3U1lypiw5K6gstZ1zCdlvjCmAB/ZcTv+x/DKzJTmdBtsLypK7cMqjLDkzsyVeZzhLa8LpNVSWX1KWnEVSFi0PiYiIiC+oaRERERFfUNPyX694HSALZcldOOVRlpyFU5b8CKfcypIzZclZxGQpEhviioiISNGnmRYRERHxhYhtWsxsjJktN7NlZjbTzKrnMKalmS0ys9Tg2Ku8yhIcd6OZrQt+3VhIWR43s9XBPB+ZWblcxv0x+LqsNLN/m1mch1nKmdn7wbGrzKxjQWfJT57g2Ggz+8bMpniVxcxqmtlcM0sLvld3eJUlOK6Pma0xs/Vm9tfCyHI6eXk9LODZYMblZtbawyxDghlWmNnnZnZeYWTJa54sY9uaWbqZXe5lFjPrFqyZqWY236ssZlbWzCab2bfBMUMLKUucmX2Z5XlG5TAm1szeDf78LjazOh5muSv4ui03s9lmVrtAntw5F5FfQJksl28HXsphTEMgKXi5OrALKOdRlgrAxuC/5YOXyxdCll5ATPDyo8CjOYxJBDYBJYLX3wNu8iJL8L43gF8HLxcvjPcoP3mC998FvA1M8SoLUA1oHbwcD6wFmnqUJRrYANQLvkffFkaWM+Q84+sB9AM+BgzoACz2MEunU7/jQN/CypKfn5Xg+zgHmAZc7uFrUw5IA2oFr1f2MMt9p37mgQTge6B4IWQxoHTwcjFgMdAh25jfE/z7AVwNvFtIr0tesnQHSgYv/66gskTsTItz7ocsV0sBv9i4xzm31jm3Lnh5J7CXwA9lyLMAvYFPnHPfO+cOAJ8AfQohy0znXHrw6hdAjVyGxgAlzCwGKAns9CKLmZUFLgBeDX7PCefcwYLOktc8wUw1gP7AuMLIkdcszrldzrmvg5cPA6sINJwhzwK0A9Y75zY6504A7wCDCjrL6eTx9RgEvOkCvgDKmVk1L7I45z4P/q7D6X8XQ5In6DbgAwK10Mss1wIfOue2BscVSp48ZnFAvJkZUJpA05JOAQv+TP4YvFos+JX9b8UgAh/iAN4HegRzhTyLc26uc+5o8GqB/fxGbNMCYGYPmdk2YAgw4gxj2xH4hLjBoyyJwLYs17dTCH+AsvkVgU+d/8M5twN4AthKYPbpkHNuphdZgLrAPuC14HLMODMrVchZTpcHYCzwFyAzBDnOlAWA4DRxKwKfiLzI4sXPb65O83qEPGce35ubOcN7XNh5zCwRuAR4MRQ5TpeFwCx4eTObZ2ZLzewGD7P8A2hC4IPbCuAO51yh/O4Hl52XEWgaP3HO5frzG/wgcQio6FGWrArs57dINy1mNssC21xk/xoE4Jy73zlXExgP/OE0j1MN+Bcw9Gx/GAsqS0E4U5bgmPsJfFoYn8P3lyfQ0dclsGxWysyu8yILgRmf1sCLzrlWwBHgrLeXKIDX5mJgr3Nu6dlmKKgsWcaUJvDp+M5ss3ohzxIOCuL1CGUWM+tOoOjf43GescA9hfUHOZ9ZYoA2BGY0ewPDzayhR1l6A8sI1MKWwD/MrExh5HDOZTjnWhKYtWhnZs0K43kKMkvwb0MK8HhBPXHEfwG1gJW53FcG+JpCWr/NaxbgGuDlLNdfBq4ppAw3AYsIrkfmcP8VwKtZrt8AvOBRlqrA5izXuwBTC/H9OVOehwl8Ot8M7AaOAm95kSU4phgwA7irsF6TPL4uHYEZWa7fC9xbmJnO5vXI/ntF4FxE1bzIEhzTgsDsbsMweG02BX+uNwM/EviEPdijLH8FRmW5/ipwhUdZpgJdslyfA7QLwfs1AvhztttmAB2Dl2MInAPIvMgSvL0ngSW1AtvmqFD/I+H8RXAD2+Dl24D3cxhTHJhNoLv2OkuFYNEoH/zaBFQohCx9CGzglnCaMe2BVALbshiBNdTbvMgSHLeQwEn6AEYCjxfS+5SnPFnGd6PwNsTNy/tkwJvA2MLIkM8sMQQ2Hq/LfzfETS7MXGfzehD45J51Q9wvPcxSC1gPdAqH1ybb+NcpvA1x8/LaNAnW5lPb1K0EmnmU5UVgZPByFWAHUKkQsiQQ3MkAKBGsexdnG3Mr/7sh7nuF9B7lJUsrAg13UoE+d2H8h/zwRWCqbyWwHJgMJAZvTwHGBS9fB5wkMPV36qulF1mC138VLGLrCSxVFcbrsp7Amuip/++pX4DqwLQs40YBq4O5/wXEepilJbAk+PpNoBD2qspPnizju1F4TcsZswDnE9g4bnmWcf08fJ/6Edj7YgNwf2G8LmfImePrAdwC3BIcY8DzwYwrgBQPs4wDDmS5f4mXr0228a9TeE1LnrIAdxNolldSSB8s8/g+VQdmBn9eVgLXFVKWFsA3wSwrgRHB20cDA4OX44D/BH8nvwTqeZhlFrAny+s2qSCeW0fEFREREV8o0hviioiISNGhpkVERER8QU2LiIiI+IKaFhEREfEFNS0iIiLiC2paRERExBfUtIiIiIgvqGmRkDOzEmY238yiC+CxipvZguDZpkUkQqiORCY1LVLg8lBEfkXgtPIZ5/pczrkTBA7nfdW5PpaIhA/VEcmJmhYpEGb2HzN72cy+IHAyvNMZAkwMfl8dM1uZ5XH+bGYjs9y32sxeN7O1ZjbezHqa2Wdmts7M2gW/bULwMUXEx1RH5EzUtEhBaQ7scc51cM79LbdBZlacwPkwNufxcRsATwKNg1/XEjgfyJ+B+4JjVgJtzzK3iIQP1RE5LTUtcs7MLI7AWahH52F4JeBgPh5+k3NuhXMuk8CZpWe7wAmzVgB1AILTwyfMLD4/uUUkfKiOSF6oaZGCkAwsds6lA5hZgpm9ZmY1zOz/zKxYlrE/ETgTaVaW5XKxbPcdz3I5M8v1TAKnpT8lFjh2tv8BEfGc6oickZoWKQjNCZyiHADn3D5gK4Hp2Nudcyez3HcAiA5+qjqldrBARQEXAPnaG8DMKgL7sz6PiPiO6oickZoWKQj/U2zMrDRQD0h3zv2Yw/iZBNaTT/kOeBNYQmBd+QYzq5+P5+8OTM1vaBEJK6ojckYWWNYTKRjB4xy8AowCrgS+cs7NyzamNfBH59z1ZlYHmOKca3YOz/kh8Ffn3NqzDi4iYUN1RHKjA+lIgQquR/8qePXxXMZ8bWZzC+qgUMAEFRqRokN1RHKjmRYRERHxBW3TIiIiIr6gpkVERER8QU2LiIiI+IKaFhEREfEFNS0iIiLiC2paRERExBfUtIiIiIgvqGkRERERX/h/znCe3xy5+8EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 648x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "\n",
    "\n",
    "for plot in [1, 2]:\n",
    "    plt.subplot(1, 2, plot)\n",
    "    plt.plot(simulation._bandstructure.r[0], simulation._bandstructure.r[1])\n",
    "    plt.plot(simulation._bandstructure.r[0], simulation._bandstructure.r[1], 'k.')\n",
    "    plt.xlabel(r'$r_x\\ (\\mathrm{\\mu m})$')\n",
    "    plt.ylabel(r'$r_y\\ (\\mathrm{\\mu m})$')\n",
    "\n",
    "    plt.axhline(int_r[1], color='k')\n",
    "    plt.plot(initial_r[0], initial_r[1], 'bo')\n",
    "    plt.plot(simulation._bandstructure.r[0][initial_n_f[0]+1], simulation._bandstructure.r[1][initial_n_f[0]+1], 'co')\n",
    "    plt.plot(int_r[0], int_r[1], 'go')\n",
    "\n",
    "    plt.plot(ref_r[0], ref_r[1], 'go')\n",
    "    plt.plot(simulation._bandstructure.r[0][final_n_f[0]], simulation._bandstructure.r[1][final_n_f[0]], \n",
    "            color='red', linestyle='dashdot', marker='o')\n",
    "    plt.plot(simulation._bandstructure.r[0][final_n_f[0]+1], simulation._bandstructure.r[1][final_n_f[0]+1], \n",
    "            color='red', linestyle='dashdot', marker='o')\n",
    "\n",
    "    if plot == 1:\n",
    "        plt.xlim(-3.3, -2)\n",
    "        plt.ylim(4.5, 6)\n",
    "\n",
    "    if plot == 2:\n",
    "        plt.xlim(2, 3.3)\n",
    "        plt.ylim(4.5, 6)\n",
    "\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd43db0b",
   "metadata": {},
   "source": [
    "# Vertical wall test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6b23e3d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial: [[22, 0.7], [-0.4, 0]]\n",
      "projected: [(23, 1), -1.0223438518654047, -0.3637794610839975]\n",
      "step_coords: [(-0.4, 0), (-1.0223438518654047, -0.3637794610839975)]\n",
      "intersection: (((-0.5, -4.9), (-0.5, 4.9), 0), -0.49999999991366717, -0.05845312999872237, 0.11583077479401131)\n",
      "n_f_int: (22, 0.5875219868730923)\n",
      "final: (16, 0.4124780131269027) -0.8667187934839317 0.1559054833217101\n"
     ]
    }
   ],
   "source": [
    "edge = bar_frame.edges[6] \n",
    "initial_position = [-0.4, 0]\n",
    "#initial_position = [0.8, 3]\n",
    "initial_n_f = [22, 0.7]\n",
    "#initial_n_f = [22, 1]\n",
    "print(f'initial: {[initial_n_f, initial_position]}')\n",
    "\n",
    "projected_n_f, projected_x, projected_y = simulation._update_position(initial_n_f, initial_position[0], initial_position[1])\n",
    "print(f'projected: {[projected_n_f, projected_x, projected_y]}')\n",
    "\n",
    "step_coords = ([(initial_position[0], initial_position[1]), (projected_x, projected_y)])\n",
    "print(f'step_coords: {step_coords}')\n",
    "intersection = simulation._get_sorted_intersections(step_coords)[0]\n",
    "print(f'intersection: {intersection}')\n",
    "intersection_x = intersection[1]\n",
    "intersection_y = intersection[2]\n",
    "n_f_int = simulation._get_n_f_intersection(initial_n_f, [(initial_position[0], initial_position[1]), (intersection_x, intersection_y), (projected_x, projected_y)])\n",
    "print(f'n_f_int: {n_f_int}')\n",
    "\n",
    "final_n_f = simulation._specular(n_f_int, edge)\n",
    "post_n_f, final_x, final_y = simulation._update_position(final_n_f, intersection_x, intersection_y)\n",
    "print(f'final: {final_n_f} {final_x} {final_y}')\n",
    "\n",
    "initial_r = (simulation._bandstructure.r[0][initial_n_f[0]], simulation._bandstructure.r[1][initial_n_f[0]]) + (1-initial_n_f[1])*simulation._bandstructure.dr[:, initial_n_f[0]]\n",
    "int_r = (simulation._bandstructure.r[0][n_f_int[0]], simulation._bandstructure.r[1][n_f_int[0]]) + (1-n_f_int[1])*simulation._bandstructure.dr[:, n_f_int[0]]\n",
    "ref_r = (simulation._bandstructure.r[0][final_n_f[0]], simulation._bandstructure.r[1][final_n_f[0]]) + (1-final_n_f[1])*simulation._bandstructure.dr[:, final_n_f[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "793a0b8d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-2.6826267075816213, 5.507046406439765) (-2.6826267075816213, 105.50704640643977)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-2.6826267075816213, -4.282143153892547)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print((int_r[0], int_r[1]), (int_r[0], int_r[1]+100))\n",
    "simulation._calc_int_of_two_lines([(int_r[0], int_r[1]), (int_r[0], int_r[1]+100)], [(-3.65356602408424, -4.618027135479627), (-5.128135499854627, -5.128135499854627)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "52223132",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[22, 0.7]\n",
      "0.16068287589558233\n",
      "0.5875219868730923\n",
      "0.6041\n"
     ]
    }
   ],
   "source": [
    "# Trying to figure out what get_n_f_intersection should be\n",
    "\n",
    "initial_n_f\n",
    "initial_position[0], initial_position[1]\n",
    "intersection_x, intersection_y\n",
    "projected_x, projected_y\n",
    "\n",
    "\n",
    "root =  np.sqrt(((intersection_x - initial_position[0])**2 + (intersection_y - initial_position[1])**2) /\n",
    "                                 ((projected_x - initial_position[0])**2 + (projected_y - initial_position[1])**2))\n",
    "f = initial_n_f[1] * (1-root)\n",
    "\n",
    "print(initial_n_f)\n",
    "print(root)\n",
    "print(f)\n",
    "print(0.7-0.137*0.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6585015e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAFFCAYAAAB42U7BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABC80lEQVR4nO3dd3xV9f3H8dcnmyUrYZMESJDlJAwFrFvcWlur4sKRWqvVLqvSaqvF1vZXR+tExFFxrypSN8pQwKAiUwgbZC+ZIePz++Pe2AQTSMjNPTfJ+/l43EfuPffknM+F3HM/9zs+X3N3RERERErFBR2AiIiIxBYlByIiIlKOkgMREREpR8mBiIiIlKPkQERERMpRciAiIiLlxExyYGadzWyCmc01szlmdkMF+5iZ/dPM8s3sKzM7MohYRURE6rOEoAMoowj4tbt/bmbNgBlm9p67zy2zz6lAdvg2AHg4/FNEREQiJGZaDtx9tbt/Hr6/DZgHdNxrt7OBpz1kKtDCzNpHOVQREZF6LWaSg7LMLBM4Api211MdgRVlHq/k+wmEiIiI1EAsdSsAYGZNgVeAG9392xocJxfIBWjSpEnfHj16RChCqU9mzJixwd3Tgo5DRCSWxFRyYGaJhBKDse7+agW7rAI6l3ncKbzte9x9FDAKICcnx/Py8iIcrdQHZrYs6BhERGJNzHQrmJkBjwPz3P2eSnZ7A7g0PGthILDV3VdHLUgREZEGIJZaDgYBlwCzzOzL8LZbgXQAd38EGA+cBuQDO4Hh0Q9TRESkfouZ5MDdJwO2n30c+Hl0IhIREWmYYqZbQURERGKDkgMREREpR8mBiIiIlKPkQERERMpRciAiIiLlKDkQERGRcmJmKqOI1H+pqamemZkZdBgiElZZCXklByISNZmZmaiUuUjsqKyEvLoVREREpBwlByIiIlKOkgMREREpR8mBiIiIlKPkQERERMpRciAiIiLlKDkQERGRcpQciIiISDkqgiQiIgDMWLaJl/JWsquwuMbHapaSwLABGfRsf1AEIpNoU3IgItKAuTufLNrIAx/m8+nijTRLTqB106QaH3fdtgKembqcE3u24efHZXFEessIRCvRouRARKQBcnc+nL+Of32Yz5crttCmWTK/P70nFw1Ip3FSzT8atu4s5MlPlvLEJ0s496FPGJTVmuuOy2Zg11aYWQRegdQmJQciIg1IcYnz39mreXDCIuat/pZOLRvx53P68KO+nUhJjI/YeZo3TuSGE7O5ckgXnp22jFETl3DhY1Ppm9GS647L4tiD05QkxDAlByIiDUBhcQmvf7GKhz9exOL1O+ia1oR//Pgwzjq8A4nxtTc2vWlyArnHdOPSozJ5MW8Fj368mOFPfkbvDgdx3XFZnNK7HXFxShJijZIDEZF6bHdhMS/NWMkjHy1i1ZZd9Gx/EA9edCRD+7QjPoofyimJ8Vx6VCYX9Ev/Lkn52djPyWrTlGuP7cZZh3UgoRaTFKkeJQciIvXQzj1FPDttOaMmLmbdtgKOSG/Bnef05riD2wTanJ+UEMf5/TpzXt9OvDVrNQ9NyOdXL87k3vcX8LMfZHFe344kJ0Sue0MOjJIDEZF6ZOuuQp7+ZCljpixh885Cju7Wmvt+cjhHdWsdU3388XHGWYd14IxD2vPB/HU88OFCbn1tFv/8YCFXH9OVC/t3jsjASDkw+pcXEakHikucpz9dyj3vLmBbQRHH9whNIeybEdtTCOPijJN6teXEnm2YnL+BBz7M585xc3n4o3xuObUnPzyyY0wlNQ1FTCUHZjYGOANY5+59Knj+WOA/wJLwplfd/Y6oBSgiEoNmrdzKra/NYtaqrQzJTuXmU3vQu0PzoMOqFjNjSHYaQ7LT+GzpJu4aP49fvzSTVz5fyZ3n9KFbWtOgQ2xQYm30x5PA0P3sM8ndDw/flBiIBMzMWpjZy2Y238zmmdlRQcfUUGwvKOJPb87h7Acns3rrbv554RE8fUX/OpcY7K1fZiteueZo/nxOH2at2sqp903i3vcWsDsClRulamKq5cDdJ5pZZtBxiEi13A+87e4/MrMkoHHQAdV37s47c9bwxzfmsnbbboYNSOe3p/SgeaPEoEOLmLg44+KBGZzcuy1/HjeP+z9YyBszv2HkOX04Ois16PDqvVhrOaiKo8xsppn918x6Bx2MSENmZs2BY4DHAdx9j7tvCTSoem7l5p1c9VQe1zzzOS0aJ/LKz47mz+ccUq8Sg7LaNEv5rkWkxJ2LRk/jly98yYbtBUGHVq/FVMtBFXwOZLj7djM7DXgdyK5oRzPLBXIB0tPToxagSAPTBVgPPGFmhwEzgBvcfUfpDnovRkZhcQljJi/hvvcXAnDraT0YPqhLrRYwiiXHdE/jnRuP4YEP83l04iI+nL+OW07twfk5nVVEqRaYuwcdQznhboVxFQ1IrGDfpUCOu2/Y1345OTmel5cXmQClXjGzGe6eE3QcdZWZ5QBTgUHuPs3M7ge+dfc/VLS/3osH5vPlm7n11VnMX7ONE3u24Y9n9aZTy4bbe5O/bhu3vjab6Us2kZPRkrt+eAjd2zYLOqw6qbJrYJ1KOc2snYXntJhZf0Lxbww2KpEGbSWw0t2nhR+/DBwZYDx12tixY8nMzCQuLo7MzExGP/E0I16bxXkPf8KWnYU8cnFfHrs0p0EnBgBZbZrxQu5A/vajQ1m0fjun3T+Ju9+ez56n/g2ZmRAXF/o5dmzQodZZMdWtYGbPAccCqWa2ErgdSARw90eAHwE/M7MiYBdwgcda04dIA+Lua8xshZkd7O5fAycAc4OOqy4aO3Ysubm57Ny5E4Bly5aR+9NcWg+9juuuHs6vTz6YpskxdckOlJlxfk5nTuzZlrvGz2PVg6MpfucBKAyPRVi2DHJzQ/eHDQsu0Doq5roVaoOaMqUy6laoOTM7HBgNJAGLgeHuvrmiffVerFxmZibLli373vb2HTvzzcrlAURUt+zu2JmUb1Z+/4mMDFi6NOrx1BWVXQOVhopIjbj7l4ASrBpavrziBGBNRR948j0pq1dVuN2XL0fDFauvTo05EBGpb0pKnFdmrCTxoLQKn9cMjyqq5N9pbfM03p2zhobQSh5JSg5ERAIyc8UWznvkE3790kz6nHMNKY0alXu+cePGjBw5MqDo6piRI6Fx+YGaxSmNGHPa1eT+ewaXjplO/rptAQVX9yg5EBGJsvXbCrjp5Zmc89AUVmzaxd9/dCifjfkjox97jIyMDMyMjIwMRo0axTANpquaYcNg1KjQGAMzyMggfvRj/PbpO7j9zF58uWILQ++bxJ3j5vLt7sKgo415GpAoDZoGJEZXQ38vFhaX8NQnS7n//YXsLipm+KAuXH98Fs1S6md1w1iycXsB//fuAp7/bDmtGidx09CD+XFfFVDSgEQRkQBNXLCeP705h0Xrd/CD7mncdmYvrTQYRa2bJvOXHx7CsAHp3P7GHH73yizGTlvO7Wf2jvllrYOgbgURkQjau5DR/Y+O4aqn8rh0zHSKS5zHL8vhyeH9lBgEpE/H5rx8zVHc95PDWfvtbs57+BN+9cKXrP12d6hokoooAWo5EBGJmIoKGf3yumtpf8Yv+P0NP+WKwZkkJ8QHHKWYGecc0ZGTerXlwQn5jJ60hPjnn+Wu//6LxILdoZ0aeBEljTmQBk1jDqKrvr8XKytk1KlzZ1ZUUsdAgrd0ww6a9sgideOa7z9Zz4soacyBiEgtq6yQ0aqVKmQUyzJTm8CmtRU/2UCTOo05EBGpoQVrtzH8ienENUut8HkVMqoDKvk/2tS6HSs27YxyMMFTciAicoDWfbubm1/5iqH3TSRv2WYuvf4WGu9ViEeFjOqICoooFSanMHLwxZzwj4/587i5bNm5J6Dgok/dCiIi1bSjoIhHJy7msYmLKSop4fKjQ/UKWjY5hRN6tmHEiBEsX76c9PR0Ro4cqUJGdUHp/9GIEaGuhPR0EkeO5LdnnEf8e1/z+JQlvJi3guuPz+aSozJISazfA0s1IFEaNA1IjK66/l4sKi7hhbwV3PveQjZsL+D0Q9tz0ykHk9G6SdChSS2bv+Zb/jJ+Ph8vWE+nlo347SkHc+ahHep8ESUNSBQROUDuzgfz1vHXt+eTv247/TJb8tilfTkiXcVzGooe7Q7iqSv6M3nhBu4aP48bnv+Sxycv4ZZTe3JUt9ZBhxdxSg5ERPZh5oot3DV+HtOWbKJrahMevaQvJ/dqi1nd/sYoB2Zwdirjrh/Ma1+s4v/e/ZoLH5vKCT3acMtpPchq0yzo8CJGyYGISAVWbNrJ3975mjdnfkPrJknceXZvLuifTmK8xnE3dHFxxnl9O3H6oe0ZM2UJD09YxCn3TeIn/Tpz44nZtGmWEnSINabkQESkjK27Cnngw4U89cky4uLguuOy+OkPumpxJPmelMR4rj02i5/kdOZfH+bzzNRlvP7FKq75QTeuHtKVRkl1d9CikgMREUIrJo6duoz7P1jIll2FnHdkJ359cnfaN28UdGgS41o3TeaPZ/XmsqMzufu/87nnvQU8O205vz3lYM49omOdHLSo5EBEGjR35/156/jL+Hks3rCDo7u1ZsTpPendoXnQoUkd0yW1CY9c0pfpSzYx8q25/PqlmTzxyRJGnNarzg1aVHIgIg3W7FVb+fNbc5m6eBNd05ow+tIcTujZRoMNpUb6d2nFa9cO4o2Z3/C3t+dz4WNTOalXW245tQdd68hqnEoORKTBWbN1N39/52te/WIlLRsnccfZvblQgw0lguLiQis/Du3TjscnL+HhjxZx8r0TuXhgBjeckE3LJklBh7hPSg5EpMHYUVDEox8vYtSkxZSUQO6Qrlx7XBbNG2mwodSOlMR4fn5cFufndObe9xfw9KdLefXzlVx/fDaXHp0Rs0t4KzkQkXqvuMR5ecYK/u/dBazfVsAZh7bnd0N70LlV4/3/skgEpDVL5q5zD+HyozO5a/w8Ro6fx7+nLuPmU3twap92MdeVpeRAROq1SQvXM/Ktecxfs40j0lvwyMV96ZuhyoYSjO5tm/Hk8P5MXLCeu8bP49qxn5OT0ZIRp/eMqYqbMZUcmNkY4Axgnbv3qeB5A+4HTgN2Ape7++fRjVJE6oJF67fz53FzmfD1ejq3asQDFx3B6Ye0j7lvaNIwHdM9jUFZqbyUt4J/vLeAcx/6hDMP68Atp/agQ4vgp8/GVHIAPAk8ADxdyfOnAtnh2wDg4fBPEQmQmcUDecAqdz8j6Hg+X76Zyx6fDga3ntaDy47OjNm+XWm44uOMC/qnc+ZhHb4bCzN18Uaeu3ogWW2CndUQU0Nz3X0isGkfu5wNPO0hU4EWZtY+OtGJyD7cAMwLOggIJQaXPj6dlk2SePvGY8g9ppsSA4lpTZIT+NXJB/PmdYNxhwsfm0r+uu2BxhRTyUEVdARWlHm8MrxNRAJiZp2A04HRQccyY1koMWjdNIkXfjqQjjHQPCtSVdltm/F87gDc4YJRU8lfty2wWOpaclBlZpZrZnlmlrd+/fqgwxGpz+4DbgJKggxixrJNXDZmOqlNk3g+d6DKHkudlNUmlCAAXDBqGgvXBpMg1LXkYBXQuczjTuFt3+Puo9w9x91z0tLSohKcSENjZqUDiGfsY59aT9Tzlm7i0senk9Ysmedzj1JiIHVaKEEYiFmoi2FBAAlCXUsO3gAutZCBwFZ3Xx10UCIN2CDgLDNbCjwPHG9mz5TdobYT9c+WhloM2hyUwnNXD6Rd87q/XK5IVpumPHf1QMyMC0dFP0GIqeTAzJ4DPgUONrOVZnalmV1jZteEdxkPLAbygceAawMKVUQAd7/F3Tu5eyZwAfChu18crfOXJgZtD0rh+VwlBlK/ZLVpyvO5A4mPCyUIX6+JXoIQU1MZ3f3C/TzvwM+jFI6IxLDpSzZx+RPTaXdQCs/lDqTtQUoMpP7plhZKEC58bCoXPjaVZ68eQI92B9X6eWOq5UBE6i53/yhaNQ6mLd4YSgyah1oMlBhIfdY1rSnP5x5FYrxx0WPTmLf621o/p5IDEalTpi7eyPAnP6N9ODFoo8RAGoAuqU14PvcokuLjuOixqbWeICg5EJE649NFGxn+xGd0aNGI53IH0qaZEgNpOEIJwkBSEuO56LGpzP2m9hIEJQciUid8umgjVzz5GZ1aNuK5q5UYSMOUWTZBGD2VOd9srZXzKDkQkZg3e9VWhj85nU4tG/Hs1QNJa5YcdEgigcloHUoQGifGM2z0NFZs2hnxcyg5EJGY9/KMlRjGc7lKDEQglCA8c9UAtuws5M2vvon48ZUciEjMm5y/gf5dWpHaVImBSKmuaU3p0a4ZU/I3RPzYSg5EJKat2bqb/HXbGZyVGnQoIjFncFYqny3dzO7C4ogeV8mBiMS00m9FR2e1DjgSkdgzKDuVPUUl5C3dHNHjKjkQkZg2JX8DrZok0TMKVeFE6pr+ma1IjDcmR7hrQcmBiMQsd2dy/gaO7taauDgLOhyRmNMkOYEj0ltGfNyBkgMRiVn567azbluBxhuI7MPgrFRmf7OVzTv2ROyYSg5EJGaVNpUOUnIgUqlBWa1xh08Xb4zYMZUciEjMmpK/kYzWjencqnHQoYjErEM7taBpckJExx0oORCRmFRUXMLUxRvVaiCyH4nxcQzs2iqi4w6UHIhITJq5civbC4o03kCkCgZlpbJs486IlVJWciAiMWlK/gbM4Kiuqm8gsj+lSXSkWg+UHIhITJqcv4HeHQ6iZZOkoEMRiXlZbZrSpllyxMYdKDkQkZizo6CIL5Zv1ngDkSoyMwZnpfLJoo2UlHiNj6fkQERizvSlmygsdo03EKmGQVmpbNqxh3lrvq3xsZQciEjMmbJwA0kJcfTLbBV0KCJ1xqAIjjtQciAiMWfKoo3kZLQkJTE+6FBE6ox2zVPIatOUKfk1L4ak5EBEYk7+um306dg86DBE6pxDOjYnf932Gh9HyYGIxJw4M7TMksiBiYvAJ7uSAxGJOfFxRnEERlyLNDTFJU5CBLKDmEoOzGyomX1tZvlmdnMFz19uZuvN7Mvw7aog4hSR2hVvRpGSA5FqKy5xIrG6eULNDxEZZhYPPAicBKwEPjOzN9x97l67vuDu10U9QBGJmvh4o8SVHIhUV2AtB2bWJPxBHmn9gXx3X+zue4DngbNr4TwiEiFm1tnMJpjZXDObY2Y3ROK48aZuBZEDUexOXASaDvabHJhZnJldZGZvmdk6YD6wOnwx+LuZZdU4ipCOwIoyj1eGt+3tPDP7ysxeNrPOETq3iByYIuDX7t4LGAj83Mx61fSgGnMgcmCKS5z4KA1InAB0A24B2rl7Z3dvAwwGpgJ3m9nFNQ+lSt4EMt39UOA94KnKdjSzXDPLM7O89evXRyk8kYbF3Ve7++fh+9uAeVSc1FeLkgORAxNKDmqeHVRlzMGJ7l6490Z33wS8ArxiZok1jgRWAWVbAjqFt5U9Z9nKDqOBv1V2MHcfBYwCyMnJ0VVGpJaZWSZwBDCtpsdSciByYIpLnPgIDEjcb3pRmhiYWY6ZvWZmn4eb9WeZ2Vdl96mhz4BsM+tiZknABcAbZXcws/ZlHp5F6FuKiATMzJoS+rJwo7t/u9dz1W7Fi48zijUgUaTaQi0HNc8OqjNbYSzwW2AWUFLjM+/F3YvM7DrgHSAeGOPuc8zsDiDP3d8AfmFmZxHq59wEXB7pOESkesIth68AY9391b2fP5BWPA1IFDkwxR795GB9+AO61rj7eGD8XttuK3P/FkJjH0QkBpiZAY8D89z9nkgdV90KIgemuMRJTIzOmINSt5vZaOADoKB0Y0XfFESkwRgEXALMMrMvw9tuDSf6B0zJgciBieaAxFLDgR5AIv/rVnBAyYFIA+XukyHyyyAoORA5MJEakFid5KCfux9c81OKiOybBiSKHJhIDUisTtvDJ5EobiIisj9qORA5MCUBDEgcCHxpZksIjTkwwMMFiUREIkazFUQOTFEAUxmH1vhsIiJVEBenVRlFDkRJiRNnUUwO3H3Z3tvMrEuNIxAR2UtivLG7MOLlVETqvcKSEhKi2XJgZrfttSkeuBRQgiAiEZUQF0dRcVHQYYjUOUXFTmIEVl6qTrfCjjL3E4FjgBdrHIGIyF4S443CYnUriFRXYbGTEM3kwN3/Ufaxmf0NmF7jCERE9pIQF0dRiboVRKqrqKSExAgUOqhJetESWFvjCERE9pKYEEeRWg5Eqi3q3QpmNotQRUQITWPMBDaWbteURhGJlMQ4Y0+xWg5EqmtPcQkJEWg5qM6YgzNqfDYRkSpIiDe1HIgcgKLiEhKjsbaCmZmHfG8qY9l9ahyJiNQKM2sC7Hb34qBjqaqEeI05EKmukhKnxIlIy0FV0osJZna9maWX3WhmSWZ2vJk9BVxW40hEJCLMLM7MLjKzt8xsHTAfWG1mc83s72aWFXSM+5MUH6fZCiLVVBhOqCMx5qAqRxgKFAPPmdk34QvMYmAhcCFwn7s/WeNIgjJ2LGRmQlxc6OfYsUFHJFJTE4BuwC1AO3fv7O5tgMHAVOBuM7s4yAD3JyHOKNSYA5FqKU2oIzFbYb/dCu6+G3gIeMjMEoFUYJe7b6nx2YM2dizk5sLOnaHHy5aFHgMMGxZcXCI1c6K7F+690d03Aa8Ar4TfyzErIV6zFUSqqyicUCdEYMxBtY7g7oXuvrpeJAYAI0b8LzEotXNnaLtIHVWaGJhZjpm9Zmafm9lXZjbLzL4qu0+sSoy375pIRaRqotpyUK8tX1697SJ1y1jgt8AsoE590ibExeEeubXpRRqC0kG8Ua2QWC+lp4e6EiraLlL3rXf3N4IO4kAkJoQSgsLiEuLj4gOORqRuKPqu5SC6RZCSgfMIFT/67vfc/Y4aRxGUkSPLjzkAaNw4tF2k7rvdzEYDHwAFpRvd/dXgQqqa0nnahcUlpCQqORCpitLCYdHuVvgPsBWYQZkLTZ1WOuhwxIhQV0J6eigx0GBEqR+GAz0ILZRW2q3gQMwnB6XztDUoUaTqSt8vkRiQWJ3koJO7D63xGWPNsGFKBqS+6ufuBwcdxIEo7TPVoESRqiud/hutIkilPjGzQ2p8RhGJlk/MrFfQQRyIxDi1HIhUV1FJMLMVBgOXm9kSQt0KhhZcEollA4Ev6+J7trQgu1IDkapzD71jjOgmB6fW+Gz7YWZDgfuBeGC0u/91r+eTgaeBvsBG4CfuvrS24xKpo+psN2D4GodmMYpUXekyRx6BtLrK3QrhhZdaAGeGby32tRhTdZlZPPAgoSSkF3BhBU2iVwKb3T0LuBe4O1LnD4IqN0ttcvdle9+oZuGzoIRbRyPyDUikoShNpj0CTW7Vmcp4A3A1/xvp/IyZjXL3f9U8DAD6A/nuvjh8vueBs4G5ZfY5G/hj+P7LwAOlq0ZGKIaoUeVmqW1mdttem+KBS4EuET7PPlv8DkTpNx+t9ypSdaXJdEkEPhGr8y3iSmCAu9/m7rcR6s+8uuYhfKcjsKLM45XhbRXu4+5FhKZWto5gDFGjys0SBTvK3PYA/YAXI3mCKrb4VVtpuq/kQKTqvhurE4Hvy9UZc2CEVmcsVRzeFpPMLBfIBUiPwYqHqtwstc3d/1H2sZn9DZge4dNUpcWv2iI5sEqkoShNDqLdcvAEMM3M/mhmfyS09OvjNQ/hO6uAzmUedwpvq3AfM0sAmhMamPg97j7K3XPcPSctLS2CYVbf2Fljybwvk7g/xZF5XyZjZ42ttEJzDOYxUn+0BNZG+JhVafGrttJrmwYkilTd/5Lp6A5IvAe4AtgUvg139/tqHMH/fAZkm1kXM0sCLgD2rgv/BnBZ+P6PgA9jfbzB2FljyX0zl2Vbl+E4y7YuI/fNXE773VgaNy6/ryo3SySVrsJYuiIjsBToXXZ1xijFkWtmeWaWt379+ir9TklJ6ZgDZQciVVVaGDGqAxJDJ/QZhMonR5y7F5nZdcA7hAY2jXH3OWZ2B5AXXkDmceDfZpZPKEG5oDZiiaQRH4xgZ2H5wQU7C3cyvmAEo0YNU+VmqU1nROEc+23xc/dRwCiAnJycKl22SndSaiBSdZEckLjf5MDMJrv7YDPbRvm2itKCKgfVPIwQdx8PjN9r221l7u8Gfhyp80XD8q0VDyJYvnU5w25UMiCRVzqDZ19TjS1yX8m/a/EjlBRcAFxU04P+r86B0gORqvpf8bAodCu4++Dwz2buflCZW7NIJgb1VXrzigcRVLZdJAImmNn1Zlbuj8zMkszseDN7iv91z9VIeNZQaYvfPOBFd59T0+OWfDddoaZHEmk4IlnnoMpjDszsewWHKtom5Y08YSSNE8sPLmic2JiRJ2hwgdSaoYRmEz1nZt+Y2VwzWwwsBC4E7nP3JyN1Mncf7+7d3b2bu0f0D1sDEkWqo7RbIYoDEoGTKthW6yWV67phhwxj1JmjyGiegWFkNM9g1JmjGHaI+hOkdrj7bnd/yN0HARnACcCR7p7h7le7+xcBh7hfpRc3DUgUqbpIvl2qMubgZ8C1QNe9Rjg3A6ZELpT6a9ghw5QMSCDcvRBYHXQc1aVeBZHqKx2jE63ZCs8C/wX+AtxcZvs2d99U8xBERMr7X50DpQciVVX6bolEt8J+kwN330qoTPGFNT6biATCzJLdvSDoOKrqf90KAQciUodYNAckmtnk8M9tZvZt+Lat9HHNQxCRKHjIzCoaNxSTtLaCSPV9160QiWPtb4dKpjI201RGkbrD3a8E0s3sn2aWGnQ8+6O1FUSq77u1FSJQBak6Uxl/bGbNwvd/b2avmtkRNY5ARGqdmZ1CaKnmLOAxMzsn2Ij2bU9RCWaQGK/kQKSqkhJCH+l7iktqfKzqTGX8g7tvM7PBwImEShk/UuMIRCTizCzOzG4ts6k9oZLkp7n7ucBxAYVWJbuLSkhOiNNURpFqSE6IB2B3YfF+9ty/6iQHpWc7HRjl7m8BSTWOQEQizt1LKLO2grs/Wbqsctgfoh9V1RUUFpOSGB90GCJ1Skpi6CO9oCi6LQerzOxR4CfAeDNLrubvi0h0fWVmt5vZ996n7h7Tg4l3F4ZaDkSk6pLiw8lBlFsOzidUP/0Ud98CtAJ+W+MIRKS2tCK0ENI3ZvYfM7vTzOrEwmUFRcXfNZGKSNWYGckJcRFpOajyks3uvtPMFgGnhAc3TXL3d2scgYjUCnc/H0I1DoDewCFAf+ClIOOqioKiku+aSEWk6lIS46ObHJjZDcDVwKvhTc+Y2Sh3/1eNoxCRWhMufvR5+FYn7C5Uy4HIgUhOiIvIgMQqJwfAlcAAd98B363I+Cmg5EBEIkotByIHJlItB9V59xn/m7FA+L7mGYlIxBUUlajlQOQAhMYcRLfl4Algmpm9RigpOJtQrYOYt6mwMOgQRKQadhcW06JRYtBhiNQ5yYlx7C6MYsuBu98DDAc2ARuA4e5+X40jiIJlBQWMXbs26DBEpIpC3QpqORCprpSE+Ii0HFSnfHIKcCyhymrHAseGt8W8EndGLF68/x1FJCaEBiRqzIFIdUW95QB4mtB0qH8CDwC9gH/XOIIoWV5QZ1arFWnwCopKSFbLgUi1RarloDpjDvq4e68yjyeY2dwaRxAl6cnJQYcgIlVUoJYDkQOSnBhHQZRbDj43s4GlD8xsAJBX4wiiIM6MkV27Bh2GiFTR7qISkjWVUaTakhPi2R3lloO+wCdmtjz8OB342sxmAe7uh9Y4mlqSkZzMsLZtgw5DRKrA3dlTVEKKpjKKVFtKhFoOqpMcDK3x2QLSKlFTokTqitICLmo5EKm+5IT46FZIdPdlNT5bJcysFfACkAksBc53980V7FcMzAo/XO7uZ9VWTCISjOISB0LdgSJSPXFm372HanScCMQSCTcDH7h7NvBB+HFFdrn74eGbEgOReqhxUjyJ8cbWXSpeJlJdW3btoUXjpBofJ1aSg7OBp8L3nwLOCS4UEQmSmdGicRJbdu4JOhSROmfLzkJaNK55V3p1iiBdb2Yta3zGirV199Xh+2uAykYPpphZnplNNbNzaikWEakCM/u7mc03s6/M7DUzaxGpY7dsnMjmHWo5EKmuzTv30DLKLQdtgc/M7EUzG2pWvQ5BM3vfzGZXcDu77H7u7kBlHSYZ7p4DXATcZ2bd9nG+3HAikbd+/frqhCoiVfMeofonhwILgFsideAWjZPYpJYDkWrbvGNPdFsO3P33QDahxZYuBxaa2V37+oDe6/dPdPc+Fdz+A6w1s/YA4Z/rKjnGqvDPxcBHwBH7ON8od89x95y0tLSqvkwRqSJ3f9fdi8IPpwKdInXsVupWEDkgm3cW0qpJlMcchL/VrwnfioCWwMtm9rcaxvEGcFn4/mXAf/bewcxamlly+H4qMAioMxUaReq5K4D/RupgLZsksnmnuhVEqqO4xPl2d2FEBiRWeSqjmd0AXEpoRcbRwG/dvdDM4oCFwE01iOOvwItmdiWwDDg/fM4c4Bp3vwroCTxqZiWEkpq/uruSA5FaZGbvA+0qeGpEuNUPMxtB6MvC2EqOkQvkAqSnp1fpvKUDEt2davZgijRYW3cV4h4as1NT1SmC1Ar44d71Dty9xMzOqEkQ7r4ROKGC7XnAVeH7nwCH1OQ8IlI97n7ivp43s8uBM4ATwi2LFR1jFDAKICcnp0oTsFs2TqSw2NleUESzFBUxE6mKTTtCXXGRGJBYnSJIt+/juXk1jkRE6hQzG0qoxfAH7r4zkscuvbht2Vmo5ECkikrH6bSM9pgDEZEyHgCaAe+Z2Zdm9kikDlyaHGzWoESRKisdpxPtbgURke+4e1ZtHbtlk9DFTYMSRaquNJmOdp0DEZGoaPFdt4JaDkSqqvT9EtU6ByIi0VL6zad0gJWI7N+mHYUkxBlNk2veKaDkQERiTvNGiZipW0GkOrbs3EPLJkkRmf6r5EBEYk58nNG8UaK6FUSqIbSuQmRm9yg5EJGY1LJxkloORKph887IVEcEJQciEqNaNE5k046CoMMQqTM27VDLgYjUc306NGf6kk0s2bAj6FBEYt6MZZvIX7edPh2aR+R4Sg5EJCZdf0IWyQnxjHxLS6iI7EtJifOnN+fS9qBkrhjcJSLHVHIgIjGpTbMUrjs+i/fnrWPigvVBhyMSs179YhVfrdzKzaf2oEkEpjGCkgMRiWHDB2WS0boxd46bS2FxSdDhiMSc7QVF3P32fA7v3IKzD+sYseMqORCRmJWcEM/vT+/FwnXbGTt12f5/QaSBeXBCPuu3FXD7mb2Ii4vc8uZKDkQkpp3Ysw1DslO5570FqpgoUsayjTt4fNISfnhkR45IbxnRYys5EJGYZmb84Yxe7NhTzL3vLQg6HJGYcdf4eSTEG78b2iPix1ZyICIxr3vbZlw8IJ2x05Yxf823QYcjErgp+Rt4Z85afn5cFm0PSon48ZUciEid8MuTunNQo0TueHMu7h50OCKBKSou4Y4359K5VSOujNDUxb0pORCROqFF4yR+dVJ3Plm0kXfnrg06HJHAPDd9OV+v3caI03qSkhhfK+dQciAidcZF/dPp3rYpI9+ax+7C4qDDEYm6LTv3cM97Cziqa2tO6d2u1s6j5EBE6oyE+DhuO6M3yzftZMyUJUGHIxJ1972/kK27CrntzF4RWZq5MkoORKROGZydykm92vLAh/ms+3Z30OGIRM3Ctdv499RlXDQgnZ7tD6rVcyk5EJE6Z8RpPSkqdv72ztdBhyISFe7OHePm0iQpnl+ddHCtn0/JgYjUOZmpTRg+OJOXZ6zkyxVbgg5HpNZ9OH8dkxZu4MYTu9OqSVKtn0/JgYjUSdcdl0Vas2SufPIz3p69JuhwRGqFu/PiZyu44fkvyWrTlEuOyojKeWMiOTCzH5vZHDMrMbOcfew31My+NrN8M7s5mjGKSGxplpLIs1cNoH2LFK55Zga/euFLtu4qDDoskYhZ9+1urnoqj5te+Yo+HQ/iyeH9SIyPzsd2TCQHwGzgh8DEynYws3jgQeBUoBdwoZn1ik54IhKLsts247VrB3HDCdn8Z+Y3nHLvRC3vLPXCuK++4eT7JjI5fwO3ndGLZ68aSKeWjaN2/phIDtx9nrvvb2RRfyDf3Re7+x7geeDs2o9ORGJZYnwcvzypO69dezRNUxK4dMx0fv/6LHYUFAUdmki1bd6xh+uf+4Lrnv2CjNZNGH/DEK4Y3CWiKy5WRUwkB1XUEVhR5vHK8DYRCZCZ/drM3MxSg4zj0E4tGHf9YK4a3IWx05Zz2j8n8dnSTUGGJFItH85fy8n3TeTt2av5zcndeeWao+iW1jSQWKKWHJjZ+2Y2u4JbrXz7N7NcM8szs7z169XMKFIbzKwzcDKwPOhYAFIS4/n9Gb14/uqBlLhz/qOfctd4VVOU2LZtdyG/e/krrngyj9ZNknj954O47vhsEqI0vqAiCdE6kbufWMNDrAI6l3ncKbytsvONAkYB5OTkaJUWkdpxL3AT8J+gAylrQNfWvH3DMYwcP49RExfz0dfruOf8w+nTsXnQoYmU8+mijfzmpZms3rqLnx3bjRtPzCY5oXbWS6iOutSt8BmQbWZdzCwJuAB4I+CYRBqscKvfKnefGXQsFWmSnMBd5x7Ck8P7sXVXIec8OIX7319IYXFJ0KGJsLuwmD+9OYcLH5tKYrzx0jVH8buhPWIiMYAothzsi5mdC/wLSAPeMrMv3f0UM+sAjHb309y9yMyuA94B4oEx7j4nwLBF6j0zex+oaHWXEcCthLoU9neMXCAXID09PaLxVcWxB7fh3Rt/wO1vzObe9xfwwfy1/OPHh5HdtlnUYxEB+HLFFn714pcsXr+Dy47K4Hen9qBxUkx8HH/HGsK66Dk5OZ6Xlxd0GBKDzGyGu1daW0MqZmaHAB8AO8ObOgHfAP3dvdKKREG/F8fPWs3vX5/N9oIifnvywVwxuAvxUR4FLg1XQVEx//ogn4c+yqfdQSn8/ceHMSgr0HG8lV4DYytVEZE6wd1nAW1KH5vZUiDH3TcEFlQVnHZIe/pltuLW12Yxcvw8Xpqxgt+e0oMTe7ap1RXupGErKXH+M3MV/3h3ASs37+JHfTtx25m9OCglMejQKqXkQEQalLRmyYy6pC9vz17D39/5mqufzqNvRkt+N7QH/bu0Cjo8qUfcnY++Xs/db89n/ppt9O5wEH/54SEMyU4LOrT9UnIgIjXm7plBx1AdZsaph7TnpF5teWnGSu57fwHnP/opx/dow01DD6ZHu9pdDlfqvxnLNnP32/OZvmQTGa0b888Lj+CMQ9pHvZjRgVJyICINVkJ8HBf2T+ecwzvy5CdLefijfE69fxLnHt6RX57Unc6toleuVuqHhWu38bd3vua9uWtJbZrMnef04YJ+naO2JkKkKDkQkQavUVI8Pzu2Gxf1T+fhjxfxxJQlvPnVNwwbkMH1x2fRumly0CFKjFu1ZRf3vreAVz9fSZOkBH5zcneuGNwl5mYhVFXdjFpEpBY0b5zIzaf24PKjM7n/gwX8e+oyXspbwdXHdOWqIV1pmqxLppS3ecceHpyQz9NTlwFw5eAuXHtsFi2bJAUcWc3oL11EZC/tmqfwlx8eylVDuvKPd7/mvvcX8u9Pl3H98VlcOCA9ZgrVSHB27ilizOQlPPrxYnbsKeK8Iztx40nd6diiUdChRYSSAxGRSnRLa8pDw/oyc8UW7n57Pn98cy6jJy/h1yd35+zDOtaZwWUSOYXFJTw/fTn3f5DPhu0FnNSrLb895WC617OiWkoORET247DOLRh71QAmLdzA3W/P55cvzGTUxCX84fSeHB1wERuJDnfnnTlr+et/57F04076Z7bi0UuOpG9G/Zz+quRARKQKzIxjuqcxOCuVcbNW87e353PR6Gmc1Kstt57Wky6pTYIOUWrJnG+2cue4uUxdvInsNk154vJ+HHtwWr0unKXkQESkGuLijLMO68DJvdryxJSlPPDhQk6+92MuPzqT647Ppnmj2K16J9WzflsB/3j3a17IW0GLRonceU4fLuzXOdCllKNFyYGIyAFISQxNfzyvb0fueXcBoycv4ZXPV/Grk7pzQQP5AKmvdhcW88SUpTw4IZ/dhcVcOagL15/QsBI/JQciIjXQplkKfz3vUC45KoM73pzL71+fzb8/Xcbvz+hZJ8rkyv+4O2/PXsNd/53Hik27OLFnW249rQdd05oGHVrUKTkQEYmA3h2a83zuQN6Zs5a7xs/jksenc2LPNtx6Ws8G+eFS18xetZU7xs1l+pJNHNy2Gc9cOYDB2Q13sKmSAxGRCDEzhvZpx3E90sLjEfI5+d6JXHZ0Jr84PpvmjRtOs3Rdse7b3fzfu1/z0oyVtGycxMhz+/CTHHULKTkQEYmw5IR4rvlBN847shP3vLeAMVOW8OrnK/nlSd25qH96g//giQW7C4t5fPISHpqQz57iEq4e0pWfH5fVoMYV7IuSAxGRWpLWLJm//PAQLhmYwZ3j5nLbf+bw70+X8YczenFMd41HCIK7M37WGu4aP49VW3Zxcngqaqamopaj5EBEpJb16nAQz149gPfmrmXk+HlcOiY0HuH2M3tr5ccoWrB2G394fTbTlmyiR7tmPHvVABWxqoSSAxGRKDAzTu7djh8cHBqP8M8PFnLiPR9z3XFZ5P6gq9ZrqEXbC4q4//0FPDFlKU1TEhh5bh8u6JdOvMpfV0rJgYhIFJWORzj78A7cOW4u/3hvAa9+sYo/ndVbXQ0R5u68NWs1fx43jzXf7uaCfp25aWgPWtXxFROjQcmBiEgA2jdvxEPD+jJxwXpuf2MOl46Zzql92vGHM3rRoZ6s7BekReu3c/t/5jA5fwO9OxzEQxcfyZHpLYMOq85QciAiEqBjuqfx9o1DeGziYh6YkM/HC9bzixOyuWJQF5ISNKuhunbtKeZfHy7ksUmLSUmM509n9ebigRnqQqgmJQciIgFLTojnuuOzOfvwjvzpzbn89b/zeXnGSu44uzdHd9OAuapwd96du5Y73pzLqi27+OGRHbnl1J6kNUsOOrQ6SWmpiEiM6NyqMaMvy+Hxy3IoKCrmosem8YvnvmDdt7uDDi2mLd+4kyue/Iyf/nsGTZMTePGnR3HP+YcrMagBtRyIiMSYE3q2ZVBWKg99tIhHPl7Eh/PX8cuTunPZURkqoFTG7sJiHvl4EQ99tIjEOOP3p/fksqMzSdS/UY3FxL+gmf3YzOaYWYmZ5exjv6VmNsvMvjSzvGjGKCISTSmJ8fzqpO68e+Mx9M1oyZ3j5nLGvyaTt3RT0KHFhAnz13HyvRO57/2FnNK7HR/+5liuGtJViUGExErLwWzgh8CjVdj3OHffUMvxiIjEhMzUJjw5vB/vzFnDHW/O5UePfMp5R3bi5lN7NMhm8xWbdnLnuLm8O3ct3dKaMPaqAQxSIaOIi4nkwN3nQahIiIiIlBda0Kk9x3RP418f5jN60mLenbOGX57UnUsbSFfD7sJiHv14MQ99lE+cGTcNPZirBnfVjI5aEhPJQTU48K6ZOfCou48KOiARkWhpnJTA74b24Ed9O/HHN+Zwx7i5vPDZCv50dm8Gdm0ddHi15v25a7lj3FyWb9rJ6Ye2Z8RpPVULopZFLTkws/eBdhU8NcLd/1PFwwx291Vm1gZ4z8zmu/vESs6XC+QCpKenH1DMIrJvZnY98HOgGHjL3W8KOKQGoVtaU56+oj/vzFnLnePmcsGoqZx1WAdGnN6TtgelBB1exCzdsIM7xs3lw/nryG7TVGshRFHUkgN3PzECx1gV/rnOzF4D+gMVJgfhVoVRADk5OV7Tc4tIeWZ2HHA2cJi7F4STdomSUFdDO37QPY2HPw7Navhg3lp+cUI2w+t4AaVde4p5cEI+oyYuJikhTrMQAlBnuhXMrAkQ5+7bwvdPBu4IOCyRhuxnwF/dvQBCSXvA8TRIjZJCsxrOO7Ijd46by1/+O58X81bwp7P6MDi7bn3Ldnfenr2GP78VWk753CM6csupPWhTj1pD6oqYSMPM7FwzWwkcBbxlZu+Et3cws/Hh3doCk81sJjCdUBPm28FELCJAd2CImU0zs4/NrF/QATVkGa2bMPqyfoy5PIeiEufix6fxs2dmsGrLrqBDq5L8ddu5dMx0fjb2c5qlhAoZ3fuTw5UYBCQmWg7c/TXgtQq2fwOcFr6/GDgsyqGJNGj7GitE6PrRChgI9ANeNLOu7l6uG0/jf6Lr+B5tObpbKqMnhdZqmPD1Oq47Lourj4nNZaG3FxTxrw8W8vjkJTRKCq2FMGxAeoOYgRHLbK/3cb2Uk5PjeXmqmSTfZ2Yz3L3SwltSOTN7G7jb3SeEHy8CBrr7+sp+R+/F6Fq1ZRd/HjeX/85eQ2brxtx+Zm+O6xEbQ0PcnTdmfsNd4+ex9tsCzs/pxE1De5DatOHVbghSZdfAmGg5EJE66XXgOGCCmXUHkgAVKIshHVs04uGL+zJpYWhZ6OFPfsaQ7FSG9mnHkKw00ls3jmo87s6i9duZtHADb321mrxlmzmkY3MeubgvR2g55Zii5EBEDtQYYIyZzQb2AJft3aUgsWFIdhpv33AMT0xZwpOfLGXEa7MBSG/VmCHZqQzJTuWobqk0b5QY8XNv3F7A5PwNTF64gcn5G1i9NbSIVEbrxow8tw8X9EvXcsoxSMmBiBwQd98DXBx0HFI1SQlx/PQH3cg9piuL1u9g8sL1TM7fwOtfrGLstOXEGRzaqQXHZKcyODuNI9JbHNDUwd2FxeQt3cyk/PVMWrCBuau/BeCglAQGZ6dyfVYaQ7JT6dwquq0WUj1KDkREGhAzI6tNU7LaNOXyQV0oLC7hi+VbmLxwPZPyN/DAhHz++WE+TZLiGdi1NUPCyUK3tCYVlrgvKXHmr9nG5Pz1TFq4gelLNlFQVEJivHFEekt+c3J3BmencUjH5mohqEOUHIiINGCJ8XH079KK/l1a8auTD2brzkI+XbyBSeFugA/mh8pXtG+ewuCsVIZ0TyMrrSlzV3/LpIXrmZK/gQ3b9wCQ3aYpFw1IZ0h2KgO6tKZJsj5i6ir9z4mIyHeaN05kaJ/2DO3THoDlG3cyKX89kxdu4J05a3hpxsrv9k1tmsSgrFQGZ6UyODuV9s213kF9oeRAREQqld66McNaZzBsQAbFJc4H89by65dm8uTwfhzRuSVx6iqol5QciIhIlcTHGSf3bses3hXVxZL6RCWoREREpBwlByIiIlKOkgMREREpR8mBiIiIlKPkQERERMpRciAiIiLlKDkQERGRcpQciIiISDlKDkRERKQcVUgUkaiZMWPGBjNbVsXdU4ENtRlPFcVKHBA7scRKHKBYKlKdODIq2qjkQESixt3TqrqvmeW5e05txlOX4oDYiSVW4gDFUltxqFtBREREylFyICIiIuUoORCRWDUq6ADCYiUOiJ1YYiUOUCwVqXEc5u6RCCSm5eTkeF5eXtBhSAwysxmx0EcoIhJL1HIgIiIi5cREcmBmfzez+Wb2lZm9ZmYtKtlvqJl9bWb5ZnZzlMMUkSgzs+vD14Y5Zva3GIjn12bmZpYaYAxVul7W4vlj4jpsZp3NbIKZzQ3/fdwQVCzheOLN7AszGxdwHC3M7OXw38g8MzvqQI4TE8kB8B7Qx90PBRYAt+y9g5nFAw8CpwK9gAvNrFdUoxSRqDGz44CzgcPcvTfwfwHH0xk4GVgeZBxU4XpZW2LsOlwE/NrdewEDgZ8H/JlwAzAvwPOXuh942917AIdxgDHFRHLg7u+6e1H44VSgUwW79Qfy3X2xu+8Bnid04RCR+ulnwF/dvQDA3dcFHM+9wE1AoAO1qni9rC0xcx1299Xu/nn4/jZCH4Idg4jFzDoBpwOjgzh/mTiaA8cAjwO4+x5333Igx4qJ5GAvVwD/rWB7R2BFmccrCegPQUSiojswxMymmdnHZtYvqEDM7GxglbvPDCqGSlR2vawtMXkdNrNM4AhgWkAh3EcocSwJ6PylugDrgSfCXRyjzazJgRwoahUSzex9oF0FT41w9/+E9xlBqKlobATOlwvkhh8WmNnsmh4zAupiac3aFAtxHBzw+Ru0fV0XCF2fWhFqMu4HvGhmXb2WpljtJ5ZbCXUpREW0r5d1mZk1BV4BbnT3bwM4/xnAOnefYWbHRvv8e0kAjgSud/dpZnY/cDPwhwM5UFS4+4n7et7MLgfOAE6o5M2/Cuhc5nGn8LbKzjeK8FzP+lTSUnFEPoYgz9/Q7eu6YGY/A14NXw+mm1kJoYRyfTRjMbNDCH0jm2lmELr2fG5m/d19TTRjKRPT5ez7ellbqnUdrm1mlkgoMRjr7q8GFMYg4CwzOw1IAQ4ys2fc/eIAYlkJrHT30haUlwklB9UWE90KZjaUUJPMWe6+s5LdPgOyzayLmSUBFwBvRCtGEYm614HjAMysO5BEAC1N7j7L3du4e6a7ZxK6AB9ZW4nB/lTxellbYuY6bKFM7XFgnrvfE0QMAO5+i7t3Cv9tXAB8GFBiQPhvcoWZlbaIngDMPZBjxURyADwANAPeM7MvzewRADPrYGbjAcIDcK4D3iE08ORFd58TVMAiUuvGAF3DXYLPA5dF+VtyrKrwehkNMXYdHgRcAhwf/nf4MvztvaG7HhhrZl8BhwN3HchBGkSFRDPLDXczKA7FEXMxiIjEmgaRHIiIiEjVxUq3goiIiMSIepkcmNmPw+U0S8ys0tHwtV0G1Mxamdl7ZrYw/LNlJfsVl+kzi8jgnv29NjNLNrMXws9PC88TjrgqxHG5ma0v8/qvqqU4xpjZusqmtFrIP8NxfmVmR9ZGHCIidUG9TA6A2cAPgYmV7RClMqA3Ax+4ezbwAZVPKdnl7oeHb2fV9KRVfG1XApvdPYtQ5be7a3reA4wD4IUyr7+2Kow9CQzdx/OnAtnhWy7wcC3FISIS8+plcuDu89z96/3sFo0yoGcDT4XvPwWcE+HjV6Yqr61sbC8DJ4SnBkU7jqhw94nApn3scjbwtIdMBVqYWfvoRCciElvqZXJQRdEoA9rW3VeH768B2layX4qZ5ZnZVDM7JwLnrcpr+26f8PSkrUDrCJy7unEAnBduyn/ZQovbBCEmy8KKiAQhahUSI60q5UWDjqPsA3d3M6tsakiGu68ys67Ah2Y2y90XRTrWGPUm8Jy7F5jZTwm1ZhwfcEwiIg1anU0O9ldetAoiUgZ0P+Vf15pZe3dfHW6irnBVOXdfFf652Mw+IrSASE2Sg6q8ttJ9VppZAtAc2FiDcx5QHO5e9pyjgb9FOIaqiqmysCKxwswaAW8Dx7t7cQ2PlQS8Hz5W0f72l+A05G6FaJQBfQO4LHz/MuB7LRpm1tLMksP3UwlV/TqgcpdlVOW1lY3tR4RKfka66MV+49irX/8sglsP/Q3g0vCshYHA1jJdQiL1Vnjg8L5cQWiNixolBhBaQpjQ4Oyf1PRYUsvcvd7dgHMJ9RkXAGuBd8LbOwDjy+x3GrCA0Lf0EbUQR2tCb4SFhLLlVuHtOcDo8P2jgVnAzPDPKyN07u+9NuAOQvXYIbRAyEtAPjAd6FpL/xf7i+MvwJzw658A9KilOJ4DVgOF4b+NK4FrgGvCzxuhmRWLwv8POUH/HeumW23dwu/9R4GpwO/3s+8nQGb4fiYwu8xzvwH+WOa5+YRmBi0gtFrkicCU8DWwf3i/w8peh3WLzZsqJIqINDBmNp/Qugi37We/JGC5u7cLP84Exrl7n/Dj3wBN3f2P4efyCXWLziHUcjiTUCJ+FjDc3c8Jt1Sscfe0WnlxEhENuVtBRKTBMbMUoBWhFrz9SQW2VOPwSzy0imUJoQThAw99A51FqGUBD3VP7DGzZtWJW6JLyYGISMPSG5jm4QGBZpZmZk+YWadwJdHEMvvuItQFWVbZeiiJez1XUOZ+SZnHJZQfAJ8M7D7QFyC1T8mBiEjDcgjwVekDd18PLAf+AfzC3QvLPLcZiA+3NpTKCCcUccAxwP4GNJZjZq2BDWXPI7FHyYGISMNSLjkws6ZAV6DI3bdXsP+7wOAyjzcCTwN5hErVX2pm3apx/uOAt6obtESXBiSKiDRQ4Rono4A/AecDn7n7R3vtcyTwS3e/ZO8BiQd4zleBm919wQEHLrWuzhZBEhGRmgmPO7gi/PDvlezzuZlNqEI9hP0Kz354XYlB7FPLgYiIiJSjMQf1hJl9Ev7Zwsyurei5CJyjkZl9HKlvEGY2MdysKSIiMUTJQT3h7keH77YArq3kuZpSGVURkQZAyUEtMrN+4aWIU8ysiZnNMbM+e+2TaWbzzWysmc0LL1vcOPzcr8xsdvh2Y3hbEzN7y8xmhrf/JLy9dJTxX4FuZvalmf19r+cqO2Zm+NyPhWN8N7zYyt6GEV4fIvw7s8sc9zdm9se9XtOTZrYg/NpONLMpZrbQzPqHf+318DFFRCSGqEm3Frn7Z2b2BvBnoBHwjLvPrmDXgwmtqTDFzMYA15rZBGA4MIBQ0ZFpZvYxoSlH37j76QBm1nyvY90M9HH3w/c+iZn1reSYm4Fs4EJ3v9rMXgTOA54p87tJhNZfWFrFl58F/JhQa8NnwEWEpkOdBdwKnENoGlS/Kh5PRESiRC0Hte8O4CRCiy1VthzxCnefEr7/DKEP0cHAa+6+Izz3+FVgCKEypCeZ2d1mNsTdt1YjlsqOCaGyp1+G788gXOq0DJVRFRFpIJQc1L7WQFOgGd8vQ1pq7ykjlU4hCU8BOpLQh+yfzWyfC6dUQ9myp8V8v1VJZVRFRBoIJQe171HgD4SWL727kn3Szeyo8P2LgMnAJOAcM2tsZk0ILUM9ycw6ADvd/RlC85KP3OtY2wglIhWp8JhVeREqoyoi0nAoOahFZnYpUOjuzxIaKNjPzI6vYNevgZ+b2TygJfCwu39OaF306cA0YLS7f0Go9Ol0M/sSuJ3QeIbvuPtGYEp4wOHf93qusmNWlcqoiog0ACqCFLBIlCONFpVRFRFpGNRyIFUWbnlQGVURkXpOLQciIiJSjloOREREpBwlByIiIlKOkgMREREpR8mBiIiIlKPkQERERMpRciAiIiLlKDkQERGRcpQciIiISDn/D07esEhC+BfUAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 648x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['figure.figsize'] = [9, 5]\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "x, y = [edge.xs, edge.ys]\n",
    "plt.plot(x, y, color='k')\n",
    "plt.xlabel(r'x position ($\\mathrm{\\mu m}$)')\n",
    "plt.ylabel(r'y position ($\\mathrm{\\mu m}$)')\n",
    "plt.plot([initial_position[0]], [initial_position[1]], 'bo')\n",
    "plt.plot(intersection_x, intersection_y, 'go')\n",
    "plt.plot(projected_x, projected_y, 'co')\n",
    "plt.plot(final_x, final_y, 'ro')\n",
    "ax1.set_xlim(-1, 1)\n",
    "ax1.set_ylim(-2, 2)\n",
    "ax1.set_aspect('equal')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(simulation._bandstructure.r[0], simulation._bandstructure.r[1])\n",
    "plt.xlabel(r'$r_x\\ (\\mathrm{\\mu m})$')\n",
    "plt.ylabel(r'$r_y\\ (\\mathrm{\\mu m})$')\n",
    "plt.plot(simulation._bandstructure.r[0][initial_n_f[0]], simulation._bandstructure.r[1][initial_n_f[0]], 'ko')\n",
    "plt.plot(simulation._bandstructure.r[0][initial_n_f[0]+1], simulation._bandstructure.r[1][initial_n_f[0]+1], 'ko')\n",
    "plt.plot(simulation._bandstructure.r[0][final_n_f[0]], simulation._bandstructure.r[1][final_n_f[0]], 'ro')\n",
    "plt.plot(simulation._bandstructure.r[0][final_n_f[0]+1], simulation._bandstructure.r[1][final_n_f[0]+1], 'ro')\n",
    "ax2.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "eff2ba64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAFFCAYAAADl1EdJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABILUlEQVR4nO3dd3hUZfrG8e+ThBB6l96r1AARpaggrF2wYnfX1WVdG3bX3tYugl3R1S0WQEWwICoodlSQ0HvvBJEOCUme3x8Z9hdjAknIzJlM7s91zcWZOe+cczNhXp6c95zzmrsjIiIiEovigg4gIiIiEi4qdERERCRmqdARERGRmKVCR0RERGKWCh0RERGJWSp0REREJGZFtNAxs+pm9o6ZLTCz+WbWM896M7OnzWyJmc0ys26RzCciIiKxJSHC+3sKmOjuZ5tZIlAxz/qTgNahx5HAC6E/RURERIosYkd0zKwacAzwTwB3z3D3rXmaDQL+4zmmAtXNrH6kMoqIiEhsieTQVXMgDXjNzGaY2StmVilPm4bA6lzP14ReExERESmySA5dJQDdgGvc/Qczewr4O3BXcTZmZkOAIQCVKlXq3q5duxILKuGzcOFCANq2bRtwEgmX6dOnb3b3OkHnKKratWt7s2bNCly/aUc6G7fvpX61JGpXLh+5YGWQ+gkpyX4kkoXOGmCNu/8Qev4OOYVObmuBxrmeNwq99jvuPhIYCZCSkuLTpk0r2bQSFn379gVgypQpgeaQ8DGzlUFnKI5mzZpxoH4kO9v56+vT+WLBJp7/y1H0aF4zgunKFvUTUpL9SMSGrtx9A7DazPaX6P2BeXmavQ9cErr66ihgm7uvj1RGEZGCxMUZwwZ3oUnNilz5xs9s3L436EgiUgiRvo/ONcAbZjYLSAYeMrMrzOyK0PoJwDJgCfAycGWE84mIFKhqUjlevLg7uzMyufKNn8nIzA46kogcREQvL3f3VCAlz8sv5lrvwFWRzCQiUhRt6lbhsbM7c/WbM3jwo3ncN6hj0JFE5AB0Z2QRkSI6tXMDLu/TnH9/v5L3ZqwJOo6IHIAKHRGRYvj7Se04snlNbhs7m3nrtgcdR0QKoEJHRKQYEuLjePaCblSvkMgVr09n2+59QUcSkXyo0BERKaY6Vcrz/EXdWL9tD0NHzyA724OOJCJ5qNARETkE3ZrU4J7TOjBlYRojJi8OOo6I5KFCR0TkEF14ZBPO7t6IpycvZvL8jUHHEZFcVOiIiBwiM+Mfp3ekY8OqXDc6lRWbdwUdSURCVOiIiJSApHLxvHBhd+LjjCten87ujMygI4kIKnREpASV9ZNxG9esyNPndWXhxh3cNnY2OfdAFZEgqdARkRKxZNNOBr/0fdj3Y2Y3mpmbWe0C1jcxs0/NbL6ZzTOzZmEPlcsxbepw4x/aMD51Hf/6bkUkdy0i+YjoFBAiEnsyMrN56culPPP5Eiokxod1X2bWGDgeWHWAZv8BHnT3z8ysMhDxCamu7NuK1NXbePCj+XRsWI0jmmmmc5Gg6IiOiBTbjFW/ctoz3zDss0Uc36Euk244Nty7HA7cAuQ7JmRm7YEEd/8MwN13uvvucIfKKy7OePLcLjQOzXS+STOdiwRGhY6IFNmu9Ezu+2AuZ77wHdv37uOVS1J49oJu1KlSPmz7NLNBwFp3n3mAZm2ArWY21sxmmNnjZhbew0wFqJpUjhcv6s7OvZrpXCRIGroSkSKZsnATd7w3h7Vb93BJz6bcfEJbqiSVK5Ftm9kkoF4+q+4Abidn2OpAEoCjga7kDG+NBv4E/LOA/Q0BhgA0adKkWJkPpG29nJnOr3lLM52LBEWFjogUypZdGTzw4Tzem7GWlnUq8c4VPUkp4XNP3H1Afq+bWSegOTDTzAAaAT+bWQ9335Cr6Rog1d2Xhd43DjiKAgoddx8JjARISUkJyyVSp3VpwMzVW3nlm+VU+HUpvn4effv2pWfPnuHYnYjkoUJHRA7I3Xl/5jru+2AeO/bu49r+rbmqX0vKJ0RuRMjdZwOH7X9uZiuAFHffnKfpT0B1M6vj7mnAccC0iAUtwN9PascXX3/D7X+5BvMsyicmMnnyZBU7IhGgc3REpEBrft3Npf/6iaGjUmlSsyIfXnM0N/yhTUSLnIMxsxQzewXA3bOAm4DJZjYbMODlIPNBzkznPStsxLMyyc7KIiMjgylTpgQdS6RM0BEdEfmdrGznP9+v4PFPFgJwz2ntuaRnM+LjLOBkOdy9Wa7lacDluZ5/BnQOINYBnXLCAJ549CHS0zOw+ASOPTbsV6iJCCp0RCSPRRt3cOu7s5ixaivHtqnDg2d0pFGNikHHKvV69uzJF59/zoOvvMP0jPqsTWwcdCSRMkGFjogAkJ6ZxXNfLOWFKUuoXD6BEecmMyi5AaGTf6UE9OzZk/FHHsWFr0zl7vFz6dqkBq0Oqxx0LJGYpnN0RITpK7dwytPf8PTkxZzauQGTbjiW07s2VJETBvFxxlPndaVCYjxXv/kze/dlBR1JJKap0BEpw3amZ3L3+Dmc/eL37MnI4l+XHsHwc5OpVTl8N/4TqFs1iWHndGHBhh08+NH8oOOIxDQNXYmUUZ8v2Mgd781hw/a9/KlXM246vi2VyqtLiJR+7Q7jL0c35+Wvl9OrZS1O6lQ/6EgiMUm9mkgZs3lnOvd9MI8PZq6jbd0qPH9hN7o2qRF0rDLp5hPa8ePyLdzy7iw6NqxG45o66VukpGnoSqSMcHfenb6GAU9+ySdzNnDDH9rwwTV9VOQEKDEhjmfO7wYO146awb4szYclUtJU6IiUAau37OaSV3/kxrdn0qpOZSYM7cO1/VuTmKAuIGhNalXkoTM7MWPVVp78bFHQcURijoauRGJYVrbz2rfLGfbpIuIMHhjUgQuPbEpclNz4T3Kc1qUB3y3dzAtTltKzRS2OaVMn6EgiMSPihU5ojpodQBaQ6e4pedZXA14HmoTyPeHur0U6p0hpN3/9dv7+7ixmrtlG/3aH8cDpHWlQvULQsaQAd5/agekrf+WGMalMGHo0h1VJCjqSSEwI6rh1P3dPzlvkhFwFzHP3LkBfYJiZJUY0nUgptndfFo9/soDTnvmGNb/u4Znzu/LKH1NU5ES5ConxPHtBN3amZ3LD6JlkZ4dlMnWRMicaB+gdqGI5dyqrDGwBMoONJFI6/LDsF05+6mue+2Ipp3dtyKQbjuW0Lrq7cWnRpm4V7j2tA98s2cwLXy4NOo5ITAjiHB0HPjUzB15y95F51j8LvA+sA6oA57q7LkUQOYDte/fxyMcLePOHVTSuWYH/XtaDo1vrPI/S6NwjGvPt0l948rNFHNm8JinNagYdSaRUC+KITh937wacBFxlZsfkWX8CkAo0AJKBZ82sat6NmNkQM5tmZtPS0tLCHFkken0ydwN/ePJLRv24ir8c3ZxPrjtGRU4pZmY8dEZHGlavwNBRqWzdnRF0JJFSLeKFjruvDf25CXgP6JGnyaXAWM+xBFgOtMtnOyPdPcXdU+rUUacuZc+mHXu58o3p/PW/06lRMZFxV/XmjlPaUzFRF1OWdlWSyvHsBV3ZtGMvt7wzC3edryNSXBEtdMyskplV2b8MHA/MydNsFdA/1KYu0BZYFsmcItHM3Rn90yoGDPuSSfM3cfMJbfngmj50blQ96GhSgjo3qs6tJ7bj03kb+e/UlUHHESm1Iv2rX13gvdCJkQnAm+4+0cyuAHD3F4EHgH+Z2WzAgFvdfXOEc4pEpRWbd3Hb2Nl8v+wXejSvySNndqJFncpBx5IwuaxPc75b+gv/+HA+3ZvWoEODakFHEil1IlrouPsyoEs+r7+Ya3kdOUd6RCQkMyubV75ZzvDPFpEYH8dDZ3TivCMa68Z/Mc7MeOKcLpz01Fdc8+YMPrimjyZeFSmiaLy8XERymbN2G4Oe+5ZHPl7AsW3qMOnGY7ngyCYqcsqImpUSeeq8rqz4ZRd3jc870i8iB6NfDUSi1J6MLEZMXsQrXy+nZqVEXryoGyd2rB90LAnAUS1qcc1xrXlq8mJ6t6zNWd0bBR1JpNRQoSMShb5bspnb3pvNyl92c94RjbntpMOpVrFc0LEkQNf2b83UZb9w1/g5JDepTkudmyVSKBq6Eoki23bv45Z3ZnLBKz9gwJt/OZJHzuqsIkeIjzOeOq8r5RPiuPrNGezdlxV0JJFSQYWOSBRwdybMXk//J7/k3Z/XcsWxLZl43TH0alk76GgSRepVS2LY4C7MX7+dhyfMDzqOSKmgoSuRgG3Ytpe7xs/hs3kb6diwKv+69Ag6NtRlxJK/49rV5bI+zfnnN8vp2bI2J3asF3QkkaimQkckINnZzls/reKRCQvYl53N7Se348+9m5MQrwOtcmC3ntiOn1Zs4ZZ3ZtKxYVUa1agYdCSRqKVCRySCvv/+e6ZMmUKrLj14e1VFflyxhV4ta/HwmZ1oWqtS0PGklEhMiOOZ87tyytPfMHRUKqOGHEU5Fcgi+VKhIxIh33//Pf3792dvejrEJdD8kkd47G9nc073RoTuFi5SaE1rVeKhMztx7VszGDFpETef8LspAUUEnYwsEjH/fvcj9uxNx7OzITuL8xrvZHBKYxU5UmwDuzTgvCMa8/yUpXyzWDPliORHhY5ImG3cvpdr35rB+I3ViUsoR1x8PEnlEzn1hAFBRyu1zOxGM3Mzy/eyNDN7zMzmmtl8M3vaYriavOe0DrSqU5nrRqeStiM96DgiUUeFjkiY7MvK5pWvl3HcE1OYOHcDt1wykMmTJ/GPBx5g8uTJ9OzZM+iIpZKZNSZnPrxVBazvBfQGOgMdgSOAYyMWMMIqJMbz7AXd2LF3HzeMSSU724OOJBJVdI6OSBj8sOwX7h4/l4Ubd9CvbR3uHdghdLJxG/oe3SfoeKXdcOAWYHwB6x1IAhIBA8oBGyMTLRht61XhntM6cPt7s3npq2X8rW/LoCOJRA0VOiIlaNOOvTwyYQFjZ6ylYfUKjLy4O39oX1fn4ZQQMxsErHX3mQV9pu7+vZl9Aawnp9B51t3zvbuemQ0BhgA0adIkPKEj5Pwejfl26Wae+HQhPZrXpHvTGkFHEokKKnRESkBmVjavT13JsE8XkZ6ZzdX9WnFVv1ZUSIwPOlqpY2aTgPzugncHcDs5w1YHen8r4HBg/8yXn5nZ0e7+dd627j4SGAmQkpJSqsd8zIyHz+zErDVbufatGUy49mhNHSKCCh2RQzZ95RbuHDeX+eu3c3Tr2tw3sAMtNOFisbl7vmdpm1knoDmw/2hOI+BnM+vh7htyNT0DmOruO0Pv+xjoCfyu0Ik1VZPK8cz53Tj7he+49d1ZvHBRNx1NlDJPJyOLFNPmnenc/PZMznrhe7buzuCFC7vxnz/3UJETJu4+290Pc/dm7t4MWAN0y1PkQM5JyseaWYKZlSPnROQyMzFUcuPq3HJiWybO3cDrP+R7vrZImaIjOiJFlJXtvPnjKh6fuIDdGVlccWxLrjmuFZXK6+sUFDNLAa5w98uBd4DjgNnknJg80d0/CDJfpF3epwXfLf2FBz6cR/cmNWjfoGrQkUQCo55ZpAhmrPqVu8bPYc7a7fRsUYsHTu9Aq8OqBB2rTAod1dm/PA24PLScBfw1oFhRIS7OGHZOF0566muufutnPrymDxUT1d1L2aShK5FC+HVXBreNncWZL3zHpu3pPH1+V978y5EqciRq1apcnhHnJbN88y7uHj836DgigVGJL3IA2dnO6GmreXTiAnbszeTyPs0ZOqANlTVMJaVAr5a1uaZfK57+fAm9W9XijK6NDv4mkRij3lqkALPXbOPO8XOYuXorPZrX5IFBHWlbT0dwpHS5tn9rpi7bwh3vzaFLo+o6WV7KHA1dieSxdXcGd46bzcDnvmHtr3sYfm4XRg85SkWOlEoJ8XE8dX4yiQlxXPPWDNIzs4KOJBJRKnREQrKznTHTVnPcsC9584dV/LFnMz6/6VjO6NpI9yKRUq1+tQo8cXYX5q7bzsMTFgQdRySiNHQlAsxdt427x89l+spf6d60Bg8MOlKX5EpMGdC+Ln/u3ZxXv11Or5a1OL5DfjefFok9KnSkTNu2Zx/DP1vEf75fQY2KiTx+dmfO6taIuDgdwZHYc+tJbflxxS/c/M4sOjasRoPqFYKOJBJ2GrqSMsndGfvzGvoP+5L/fL+Ci45qyuc39uWclMYqciRmlU+I59nzu5GZlc21b80gMys76EgiYadCR8qcBRu2c+5LU7lhzEwa1ajA+1f34f5BHTUBopQJzWpX4qEzOzFt5a88NXlx0HFEwi6iQ1dmtgLYAWQBme6ekk+bvsAIoByw2d2PjVxCiWU79u5jxKTF/Ou7FVRNSuCRMzsxWEdwpAwalNyQb5ds5tkvlnBUi1r0blU76EgiYRPEOTr93H1zfivMrDrwPHCiu68ys8Mimkxikrvz/sx1PPjRfNJ2pnPeEU245YS21KiUGHQ0kcDcO7ADP6/aynWjU/l46NHUrlw+6EgiYRFtQ1cXAGPdfRWAu28KOI+Ucos37uCCl39g6KhU6lZN4r0re/PwmZ1U5EiZVzExgWcv6Mq2Pfu4YcxMsrM96EgiYRHpQseBT81supkNyWd9G6CGmU0JtbmkoA2Z2RAzm2Zm09LS0sIWWEqnXemZPDxhPic99TXz1m/nH6d3ZNxVvUluXD3oaCJRo129qtx9anu+WpTGy18vCzqOSFhEeuiqj7uvDQ1JfWZmC9z9qzx5ugP9gQrA92Y21d0X5d2Qu48ERgKkpKToVxEBcoapJszewAMfzmPD9r0MTmnErSe2o5YOy4vk68Ijm/Dd0s08/slCjmhek25NagQdSaRERfSIjruvDf25CXgP6JGnyRrgE3ffFTqP5yugSyQzSum1NG0nl7z6I1e9+TM1KyXy7t968djZXVTkiByAmfHwmZ2pVy2Ja9+awbY9+4KOJFKiIlbomFklM6uyfxk4HpiTp9l4oI+ZJZhZReBIYH6kMkrptDsjk8cmLuDEEV+Runor9w3swPtX96Z7U/1mKlIY1SqU4+nzu7Jh215uGzsLdx0kl9gRySM6dYFvzGwm8CPwkbtPNLMrzOwKAHefD0wEZoXavOLueYshKePeeAOaNYO4OKdug0ySL17I81OWclqXBnx+Y1/+2KsZCfHRdp69SHTr1qQGN53QlgmzN/Dmj6uCjiNSYiJ2jo67LyOfYSh3fzHP88eBxyOVS0qXN96AIUNg924AY9P6BDa/1457HmvK3YMrBx1PpFQbcnQLvlv6C/d/MI/uTWvQrp7me5PST7/2Sqly2+0eKnL+X/a+eF4doSJH5FDFxRlPDu5C1QrluPrNGezOyAw6ksghU6EjpcakeRtZXcAR9VU60i5SImpXLs+Ic5NZmraTe9+fG3QckUOmQkei3uotu7n83z9x+X+mkVQjPd82TZpEOJRIDOvdqjZX92vFmGlrGJ+6Nug4IockiCkgRApl774sXvpyGc9PWUJ8nHH7ye1I6pDI367gN8NXFSvCgw8Gl1MkFg3t35qpy37h9rGz6dKoOs1qVwo6kkix6IiORKUvFm7ihBFfMXzSIga0r8vkG49lyDEtueTiOEaOhKZNwSznz5Ej4cILg04sElsS4uN46ryuJMTHcfVbP5OemRV0JJFiUaEjUWXNr7v563+ncelrPxEfZ7x+2ZE8d0E36ler8L82F14IK1ZAdnbOnypyRMKjQfUKPHFOF+as3c6jHy8MOo5IsWjoSqJCemYWr3y9nGc+XwzAzSe05fKjm1M+IT7gZCJl2x/a1+VPvZrx6rfL6dWyFgPa1w06kkiRqNCRwH29OI17xs9l2eZdnNihHned1p6G1Ssc/I0iEhG3ndyOn1Zs4aZ3ZvLx0KN/c4RVJNpp6Eoi5o3ZbzD1qKl82fdLmo1oxnM/vMZVb/zMxf/8kSx3/nXpEbx4cXcVOSJRpnxCPM9e0I19mdkMfSuVzKzsoCOJFJoKHYmIN2a/wZAPhpBeIR0MVm5byTUf/42xC0Zx4x/a8Ml1x9C37WFBx5QoZmb3mtlaM0sNPU4uoN2JZrbQzJaY2d8jnTNWNa9diX+c0ZEfV2zh6c+XBB1HpNBU6EhE3DH5Dnbv++0tjd3Sia8ximv6tyapnM7FkUIZ7u7JoceEvCvNLB54DjgJaA+cb2btIx0yVp3RtRFnd2/EM58v5rulm4OOI1IoKnQkIlZty//Wxet2rIlwEolxPYAl7r7M3TOAUcCggDPFlPsGdqB57UpcNyqVX3bmfwNPkWiiQkfCal9WNq98vYwEr5Pv+ibVdEtjKZKrzWyWmb1qZjXyWd8QWJ3r+ZrQa1JCKpVP4Nnzu7F1zz5ufHsm2dkedCSRA1KhI2Hzw7JfOPXpb/jHR/M5qtZVVEj47UnGFctV5MH+uqWx/D8zm2Rmc/J5DAJeAFoCycB6YFgJ7G+ImU0zs2lpaWmHurkyo32Dqtx1yuFMWZjGP79ZHnQckQPS5eVS4jbt2MsjExYwdsZaGlavwEsXd+f49ifz5pyWXPbGZaQnpdO0elMe7P8gF3bS3f7k/7n7gMK0M7OXgQ/zWbUWaJzreaPQawXtbyQwEiAlJUWHJorgoqOa8u2SX3h04gKOaF6T5MbVg44kki8VOlJiMrOyeX3qSoZ9uoi9mVlc1a8lV/drTYXEnBONL+x0IS9PfRmAKVOmBJhUSiMzq+/u60NPzwDm5NPsJ6C1mTUnp8A5D7ggQhHLFDPj0bM6c/LTX3PNWz/z0bVHUzWpXNCxRH5HQ1dSIqav3MLAZ7/l3g/mkdykOp9cdww3n9Duf0WOSAl4zMxmm9ksoB9wPYCZNTCzCQDunglcDXwCzAfGuPvcoALHumoVy/H0+V1Zt3Uvt42djbsOikn00REdOSSbd6bz6McLeHv6GupXS+L5C7txUsd6mFnQ0STGuPvFBby+Djg51/MJwO8uPZfw6N60Bjcd35ZHJy6gT6vanN9DFxhIdFGhI8WSle28+eMqHp+4gN0ZWfz12BZce1xrKpXXPymRsuavx7Tgu6Wbuff9uXRrUoO29aoEHUnkfzR0JUWWunorpz/3LXeNm0OHBtWYeN3R3HbS4SpyRMqouDjjycHJVEkqx9Vv/syejKygI4n8jwodKbRfd2Vw29hZnPH8t2zcvpenz+/Km385klaH6bc3kbKuTpXyjDg3mSVpO7nvA50WJdFDv4LLQWVnO6OnrebRiQvYsTeTy3o3Z+iA1lTRFRYikkuf1rW5sm9LnvtiKb1a1WZglwZBRxJRoSMHNnvNNu4cP4eZq7fSo3lNHhjUUePvIlKg6we0YeqyLdw+djZdGlWjaa1KQUeSMk5DV5Kvbbv3cee42Qx87hvW/rqH4ed2YfSQo1TkiMgBJcTH8fT5XYkzuOatGWRkZgcdSco4FTryG9nZzphpq+k3bApv/rCKP/ZsxuQbj+WMro10ybiIFErD6hV4/JwuzFqzjccmLgg6jpRxGrqS/5m7bht3j5/L9JW/0r1pDR4YdCTtG1QNOpaIlEIndKjHH3s25ZVvltOrVS2Oa1c36EhSRqnQEbbv3ceTny7iP9+voEbFRB4/uzNndWtEXJyO4IhI8d128uH8tOJXbhwzk4+HHkO9aklBR5IyKOJDV2a2InQb91Qzm3aAdkeYWaaZnR3JfGWJuzP25zUc98SX/Pv7FVx4ZFM+v7Ev56Q0VpEjIocsqVw8z1zQlfTMbIaOmkFWtqaIkMgL6ohOP3ffXNBKM4sHHgU+jVyksmXBhu3cPW4uP67YQpfG1XntT0fQqVG1oGOJSIxpWacyDwzqyI1vz+SZzxdz3YA2QUeSMiZah66uAd4Fjgg6SKzZsXcfIyYt5l/fraBqUgKPnNmJwTqCIyJhdFb3Rny7dDNPT17MUS1qcVSLWkFHkjIkiKuuHPjUzKab2ZC8K82sIXAG8ELEk8Uwd2d86lr6D/uSV79dzuCUxnx+Y1/O69FERY6IhN0DgzrSrFYlho6awZZdGUHHkTIkiEKnj7t3A04CrjKzY/KsHwHc6u4HvPmCmQ0xs2lmNi0tLS1MUWPD4o07uODlHxg6KpW6VZN478rePHxmJ2pUSgw6moiUEZXKJ/DMBV35dfc+bnp7Ju46X0ciI+KFjruvDf25CXgP6JGnSQowysxWAGcDz5vZ6flsZ6S7p7h7Sp06dcIbupTalZ7JwxPmc9JTXzNv/Xb+cXpHxl3Vm+TG1YOOJiJlUIcG1bjzlMP5fMEm/vnN8qDjSBkR0XN0zKwSEOfuO0LLxwP3527j7s1ztf8X8KG7j4tkztLO3ZkwewMPfDiPDdv3MjilEbee2I5alcsHHU1EyriLj2rKt0s28+jEBfRoXpPOjaoHHUliXKRPRq4LvBe6w24C8Ka7TzSzKwDc/cUI54k5S9N2cu/7c/l68Wba16/Kcxd2pXvTmkHHEhEBwMx47KwunPz011zy8OucXm8rJwzoT8+ePYOOJjEqooWOuy8DuuTzer4Fjrv/KdyZYsXujEye/XwJL3+9jKRy8dw3sAMXHtmEhHjN8iEi0aVaxXJc3mYfl911I7OyM3ns4YeYPHmyih0Ji2i9vFwKyd35ZO5GHvhwHmu37uHMrg257eTDqVNFw1QiEr02LPwZsjPx7GzSMzKYMmWKCh0JCxU6pdiKzbu494O5TFmYRtu6VRg95CiO1P0pRKQU6Nu3L0nly7MnPR23BDqlqMiR8FChUwrt3ZfF81OW8uKXS0mMj+POUw7nj72aUU7DVCJSSvTs2ZPJkyfz3kefMmZtVUavqsDJ2a77ekmJU6FTykyev5F7P5jL6i17GJTcgNtPPpy6VTVRnoiUPj179qRnz550+WEld7w3h/9OXckfezULOpbEGBU6pcTqLbu574O5TJq/iVaHVebNvxxJr5a1g44lInLILujRhEnzNvLQhPn0bqV+TUqWCp0ot3dfFiO/WsZzXywhPs647aR2XNq7OYkJGqYSkdhgZjx6dmdOGP4V149OxR1MI1hSQopc6IRu9LfX3bPCkEdy+WLhJu59fy4rf9nNKZ3rc+cph1O/WoWgY4kcTJyZxauPkKI4rEoSD5/ZiSte/5mErbtpVKNi0JEkRhz0sICZxZnZBWb2kZltAhYA681snpk9bmatwh+zbFm7dQ9//e80Ln3tJ+LN+O9lPXjugm4qciQqZWdn8+abb3LKKadw2GGHAXREfYQUw4kd63NWt0as3bqXnemZQceRGFGY8Y8vgJbAbUA9d2/s7ocBfYCpwKNmdlEYM5YZ6ZlZPPfFEvoPm8KXi9K4+YS2fHzd0RzdWnN5SfTq168fS5cu5eGHH2bDhg0As9RHSHHdM7A9ifFxLNm0k90ZKnbk0BVm6GqAu+/L+6K7bwHeBd41s3IlnqyM+WbxZu4eP4dlm3dxQoe63H1aBxpW1xEciX6TJk2iXLnfdwEl3UeY2b3AX4C00Eu3u/uEPG0aA/8hZ7oZB0a6+1OHum+JnKpJ5Wh5WCXmrdvOgx/N58EzOgUdSUq5gxY6+4scM0sB7gCaht5nOau9c36FkBTO+m17+MeH8/lo9nqa1qrIa5ceQb+2hwUdS6TQ9hc506ZN48EHHwQ43MxmEZ4+Yri7P3GA9ZnAje7+s5lVAaab2WfuPq+E9i8RUDWpHA2qV+CNH1Yx4PC69GunPlGKrygnI78B3AzMBrLDE6fsyMjM5rVvl/PU5MVkZTs3/KENQ45pQVK5+KCjiRTLhRdeyOOPP864ceOWAqcFkcHd1wPrQ8s7zGw+0BBQoVPKNKpRgYr1qnDLu7P45LpjqFkpMehIUkoV5RrlNHd/392Xu/vK/Y+wJYth3y3dzMlPf83DHy+gV8taTLrhWK7t31pFjpRqderUYeDAgQAZYewjrjazWWb2qpnVOFBDM2sGdAV+KOEMEgFxZgw/N5ltu/dx+9jZuHvQkaSUKsoRnXvM7BVgMpC+/0V3H1viqWLUxu17efCj+bw/cx2Na1bglUtSGNC+btCxRErEfffdx+WXXw5Q08zO3P96UfoIM5sE1Mtn1R3AC8AD5Jx78wAwDPhzAdupTM75Qde5+/YD7G8IMASgSZMmhY0pEXJ4/arceHwbHv54AWN/XstZ3RsFHUlKoaIUOpcC7YBy/P/QlQMqdA5iX1Y2//5uBSMmLSYjK5tr+7fmyr4tdQRHYsprr73GggULAKry/0NXReoj3H1AYdqZ2cvAhwWsK0dOkfPGwYosdx8JjARISUnRIYModPnRLZi8YBP3vD+XI1vU1P11pMiKUugc4e5tw5YkRv24fAt3j5/Dgg07OLZNHe4b2IFmtSsFHUukxP30008sXLgQM1vh7peW9PbNrH7oHByAM4A5+bQx4J/AfHd/sqQzSOTFxxnDzunCSU99zY1jZvLWX47SxJ9SJEU5R+c7M2sftiQxJm1HOjeMTmXwS9+zY28mL13cnX9deoSKHIlZvXr1Yt68sJ7z+5iZzQ5d0dUPuB7AzBqY2f7LzHsDFwPHmVlq6HFyOENJ+DWuWZF7TmvPD8u38M9vlgcdR0qZohzROQpINbPl5Jyj879LR8OSrJTKzMrm9akrGfbpIvZmZnFVv5Zc1a8VFRM1rZjEtqlTp5KcnAzQMe/l5SWxfXe/uIDX1wEnh5a/Ce1XYszZ3Rsxaf5GHv9kIUe3qU27elWDjiSlRFH+9z0xbClixPSVW7hr3Fzmrd/O0a1rc+/ADrSsUznoWCIRMXHiRACaNWu2iIAuL5fYZWY8dEYnThjxNdeNSmX81b0pn6DzHOXgCl3o5HeZqJk1L9k4pdMvO9N55OMFvD19DfWqJvH8hd04qWM9TNPvShnStGnT/YsZ+/sL9RFSkmpVLs+jZ3Xisn9P48nPFnHbSYcHHUlKgUIXOmZ2d56X4oFLgDLbkWVlO2/+uIrHJy5gd0YWfz22Bdce15pK5TVMJWXP/fffv3+xfqi/KPN9hJS8/ofX5fweTRj51TL6t6tLj+Y1g44kUa4o/yPvyrVcDjgGGFOycUqP1NVbuWvcHGav3UbPFrW4f1AHWtetEnQskcBUqvS/E+2zgQzKeB8h4XPnKYfz3dLN3DAmlY+HHk2VJE23KAUrytDVsNzPzewx4McSTxTlft2VwWOfLGDUT6upU7k8T52XzMAuDTRMJWXejTfeCMBNN9200d0fKat9hIRfpfIJPDk4mXNe/I77P5jH4+d0CTqSRLFDGWOpAWwsqSDRLjvbGT1tNY9OXMCOvZlc1rs5Qwe01m8SIgUrU32ERFb3pjW4sm8rnv1iCQPa1+WEDvndUFukaOfozCbnLqeQc/lmM+CX/a/H8mXms9ds487xc5i5eis9mtXk/tM76NJGkTw6deq0/8hm+1C/0Iwy0kdIMK7t35opizZx29jZdGtSgzpVygcdSaJQUY7onBq2FFFq2+59PPHpQl7/YSW1KiXy5OAunNG1oYapRPLx4Yc5MzI0a9ZsCbq8XCIgMSGO4YOTOfWZb/j7u7N45Y8p6p/ldw5a6JiZeY4CZyG2GPuXlZ3tvPPzGh75eAFbd2fwx57NuP4PbahWQcNUInm5O2aW7+Xl+8VaHyHRo3XdKtx6Yjvu/3Aeo35azfk9NDmr/FZhpoD4wsyuMbPf/Osxs0QzO87M/g38sTA7M7MVoVu4p5rZtHzWX2hms0JtvjOziJ9hNnfdNs556XtueWcWzWpV5INr+nDvwA4qckQK0K9fP5555hlWrVr1m9eL00eIFMefejWjd6taPPDhPFb+suvgb5AypTBDVycCfwbeCt38ayuQRM49Mj4FRrj7jCLss5+7by5g3XLgWHf/1cxOImdW4SOLsO0ieWPjRu5YtoxV6ek0SizPEVsTmfH1BqpXTOSxsztzdrdGmjxO5CAmTpzIq6++yvnnn8/y5csBOpjZMorfR4gUSVyc8cQ5XThh+FdcPzqVMX/tSUJ8UaZylFh20ELH3fcCzwPPm1k5oDawx923lnQYd/8u19OpQKOS3sd+b2zcyJCFC9mdnQ3A6ox01pTfy8A+h/HvfslUq6gjOCKFkZSUxJVXXsmVV17Jvn37SExMXETOLyxbg84mZUf9ahV44PSODB2VyktfLeOqfq2CjiRRokglr7vvc/f1h9CBOfCpmU03syEHaXsZ8HEx93NQdyxb9r8i53/h4o3UmpkqckSKqVy5cgD7VORIEAYlN+TUzvUZ/tki5qzdFnQciRKRPrbXx927AScBV5nZMfk1MrN+5BQ6txa0ITMbYmbTzGxaWlpakUKs37aHlXvT8123Kj3/10VEJPr94/SO1KqcyHWjU9m7LyvoOBIFilXomFmxblbg7mtDf24C3gN65LPtzsArwCB3/+UA2xrp7inunlKnTp1C7X93RibDP1tEvyemkLA3O982TcrrPgwih6q4fYTIoapeMZEnzunCkk07eWziwqDjSBQo7hGd583sD0V5g5lVMrMq+5eB44E5edo0AcYCF7v7omJm+53sbGfsz2s47okveWryYvofXpdhbVpRMe63f/2KcXE82KJFSe1WpCwrch8hUlKObl2HP/VqxqvfLufbJQVd+yJlRbEKHXe/DGhiZk+bWe1Cvq0u8I2ZzSRn/puP3H2imV1hZleE2twN1CKnk8z3EvSimrZiC2c8/y03jJnJYVXL8/YVPXnugm5c27IJI9u2pWn58hjQtHx5RrZty4V16x7qLkXKvGL2ESIl5tYT29GiTiVuensm23bvCzqOBKhYc12Z2QlAc6AV8LKZ/dvdxx3oPe6+DPjdfXHc/cVcy5cDlxcnU16rt+zmkYkL+GjWeupWLc+Tg7twenLD31wufmHduipsRMKgOH2ESEmqkBjPiHOTOfP577j7/Tk8dV7XoCNJQAp1RMfM4szs9lwv1QdedfeT3f0MoF9Y0hXDzvRMHpu4gP5Pfsnk+RsZ2r81X9zUlzN1TxyRsMnOzuahhx7K/VLU9hFSdnRuVJ1r+7dmfOo6Ppi5Lug4EpBCFTrunk2uua7c/V+hIzT73VXSwYoqK9sZ/dMq+j4+heenLOWUTvX54qa+XP+HNlRMPJRJ2kXkYOLi4v431xVEZx8hZdOVfVuS3Lg6d46bw4Zte4OOIwEoyjk6s8zsHjP73XvcfXsJZiqyXemZnPbMN9z67mya1KzAuKt6M/zcZOpXqxBkLJEypXPnztx33335rgu6j5CyKyE+juHnJpORmc3N78wkO9uDjiQRVpRCpyZwHrDOzMab2QNmdk6YchXJss272LZnH8+c35V3/9aL5MbVg44kUuZs2bKFUaNGAXSJtj5CyrbmtStxxymH8/Xizfx3aoHzU0uMKvSYjrsPhv/dH6MD0Imc++C8HZ5ohVe3ahKTbzyWpHLxQUcRKbPGjBkDgJnNAu4jivoIkQuPbMKk+Rt5+OP59G5Vm1aHVQ46kkRIkS8vd/d0d//Z3f/t7jeHI1RRHValvIockejh0dZHiJgZj53VmQrl4rlhTCr7svK/aazEHk3vKiIiZcJhVZN46IxOzFqzjWc+XxJ0HIkQFToiIlJmnNSpPmd2a8hzXyxhxqpfg44jEaBCR0REypR7B3agXtUkbhgzk90ZmUHHkTBToSMipYKZ3Wtma0PTw6Sa2ckHaBtvZjPM7MOC2kjZVTWpHMMGd2HFL7t4aML8oONImKnQEZHSZLi7J4ceEw7Qbiig/8GkQEe1qMXlfZrz+tRVfLFwU9BxJIxU6IhITDGzRsApwCtBZ5HoduPxbWlbtwq3vDOLLbsygo4jYaJCR0RKk6vNbJaZvWpmNQpoMwK4BdD1w3JASeXiGX5uMlt3Z3DHe7Nx112TY5EKHRGJGmY2yczm5PMYBLwAtASSgfXAsHzefyqwyd2nF3J/Q8xsmplNS0tLK8G/iZQW7RtU5cbj2/LxnA28N2Nt0HEkDDTbpYhEDXcfUJh2ZvYykN+Jxr2BgaETlZOAqmb2urtfVMD+RgIjAVJSUvTrfBn1l6Nb8Pn8Tdwzfi49mtekUY2KQUeSEqQjOiJSKphZ/VxPzwDm5G3j7re5eyN3b0bO3HyfF1TkiOwXH2cMG9yFbHduelsTf8YaFToiUlo8ZmazQ3Np9QOuBzCzBmZ2oCuwRA6qcc2K3DOwA1OXbeHVb5cHHUdKkIauRKRUcPeLC3h9HfC7e+q4+xRgSnhTSSw5p3sjPpu3kccmLuTo1nVoW69K0JGkBOiIjoiICDkTfz58ZieqVkjgutGppGdmBR1JSoAKHRERkZDalcvzyJmdmb9+O8M/Wxx0HCkBKnRERERyGdC+Luf3aMxLXy3lpxVbgo4jh0iFjoiISB53ntKexjUqcv3oVHbs3Rd0HDkEKnRERETyqFQ+gScHd2Hd1j088OG8oOPIIVChIyIiko+UZjX5W9+WjJm2hk/nbgg6jhSTCh0REZECDO3fhg4NqnLb2Nmk7UgPOo4UgwodERGRAiQmxDHi3GR2pGdy29hZmvizFFKhIyIicgCt61bh1hPbMWn+Jkb/tDroOFJEES90zGxF6DbuqWY2LZ/1ZmZPm9kSM5tlZt0inVFERCS3S3s1o1fLWtz/4TxW/rIr6DhSBEEd0enn7snunpLPupOA1qHHEOCFiCYTERHJIy7OeOKcLsTHGTeMmUmWJv4sNaJx6GoQ8B/PMRWonmfWYhERkYhrUL0CDwzqyPSVv/Lil0uDjiOFFESh48CnZjbdzIbks74hkHsQdE3oNRERkUANSm7AKZ3rM/yzRcxZuy3oOFIIQRQ6fdy9GzlDVFeZ2THF2YiZDTGzaWY2LS0trWQTioiI5MPMePD0jtSslMj1o1PZu08Tf0a7iBc67r429Ocm4D2gR54ma4HGuZ43Cr2Wdzsj3T3F3VPq1KkTrrgiIiK/Ub1iIo+f04XFm3by+CcLg44jBxHRQsfMKplZlf3LwPHAnDzN3gcuCV19dRSwzd3XRzKniIjIgRzbpg6X9GzKP79ZzrdLNgcdRw4g0kd06gLfmNlM4EfgI3efaGZXmNkVoTYTgGXAEuBl4MoIZxQRETmo2046nBa1K3HT2zPZtkcTf0arhEjuzN2XAV3yef3FXMsOXBXJXCIiIkVVITGe4ecmc+YL33HP+DmMOK9r0JEkH9F4ebmIiEip0KVxda49rjXjUtfx4ax1QceRfKjQEREROQRX9WtJl8bVueO9OWzYtjfoOJKHCh0REZFDkBAfx/DBXUjPzOLmd2Zq4s8oo0JHRETkELWoU5k7TmnP14s389+pK4OOI7mo0BERESkBFx3ZhGPb1OGhCfNZmrYz6DgSokJHRESkBJgZj5/dmaRy8dwwOpV9WdlBRxJU6IiIiJSYw6om8dAZnZi5ZhvPfr4k6DiCCh0REZESdXKn+pzZtSHPfrGE1NVbg45T5qnQERERKWH3DupAvapJXD86ld0ZmUHHKdNU6IiIiJSwqknlePyczizfvIuHJywIOk6ZpkJHREoFM7vXzNaaWWrocXIB7aqb2TtmtsDM5ptZz0hnFQHo1bI2l/dpzn+nruSLhZuCjlNmqdARkdJkuLsnhx4TCmjzFDDR3duRM7fe/MjFE/mtm05oS5u6lbnlnVn8uisj6DhlkgodEYkZZlYNOAb4J4C7Z7j71kBDSZmWVC5n4s+tuzO4Y9xs3TU5ACp0RKQ0udrMZpnZq2ZWI5/1zYE04DUzm2Fmr5hZpYI2ZmZDzGyamU1LS0sLW2gp2zo0qMb1f2jDhNkbGJe6Nug4ZY4KHRGJGmY2yczm5PMYBLwAtASSgfXAsHw2kQB0A15w967ALuDvBe3P3Ue6e4q7p9SpU6fE/z4i+/31mJakNK3B3ePmsnbrnqDjlCkqdEQkarj7AHfvmM9jvLtvdPcsd88GXgZ65LOJNcAad/8h9PwdcgofkUDFxxlPDk4m252bxswkO1tDWJGiQkdESgUzq5/r6RnAnLxt3H0DsNrM2oZe6g/Mi0A8kYNqUqsid5/Wnu+X/cKr3y4POk6ZoUJHREqLx8xstpnNAvoB1wOYWQMzy30F1jXAG6F2ycBDEU8qUoDBKY0ZcHhdHvtkIQs37Ag6TpmgQkdESgV3v9jdO7l7Z3cf6O7rQ6+vc/eTc7VLDZ1309ndT3f3X4NLLfJbZsYjZ3WiSvkErhudSnpmVtCRYp4KHRERkQiqXbk8j5zVmfnrtzNi0uKg48Q8FToiIiIR9of2dTk3pTEvfbmUn1ZsCTpOTFOhIyIiEoC7TmtPwxoVuGFMKjvTNfFnuKjQERERCUDl8gkMH5zM2l/3cMUTb/Hwww/z/fffBx0r5iQEHUBERKSsSmlWk5MO284Lt1yOZWdRvnwikydPpmdPzUVbUnRER0REJED1dy/HszLJzs4iIyODKVOmBB0ppqjQERERCVD/4/qRVL48WBzEJXDssccGHSmmaOhKREQkQD179uTzzyfz5L/f45tdh7HUGtAr6FAxREd0REREAtazZ0/GvPAox/c7hn98NJ8lm3TX5JIS8ULHzOLNbIaZfZjPuiZm9kVo/SwzOzm/bYiIiMQaM+PxczpTqXwCQ0elkpGZHXSkmBDEEZ2hwPwC1t0JjHH3rsB5wPMRSyUiIhKww6ok8ciZnZi7bjtPfrYo6DgxIaKFjpk1Ak4BXimgiQNVQ8vVgHWRyCUiIhItju9Qj/N7NOalr5by/dJfgo5T6kX6iM4I4BagoONx9wIXmdkaYAI5sxDny8yGmNk0M5uWlpZW0jlFREQCc+cp7WlasyI3jkll2559Qccp1SJW6JjZqcAmd59+gGbnA/9y90bAycB/zSzfjO4+MjRDcUqdOnXCkFhERCQYlconMOK8rmzckc5d4+YEHadUi+QRnd7AQDNbAYwCjjOz1/O0uQwYA+Du3wNJQO0IZhQREYkKyY2rc13/1rw/cx3jU9cGHafUilih4+63uXsjd29GzonGn7v7RXmarQL6A5jZ4eQUOhqXEhGRMulvfVvSvWkN7nxvDmt+3R10nFIp8PvomNn9ZjYw9PRG4C9mNhN4C/iTu3tw6URERIKTEB/HiHOTceCG0TPJytZ/iUUVSKHj7lPc/dTQ8t3u/n5oeZ6793b3Lu6e7O6fBpFPREQkWjSuWZF7B3bgxxVbeOmrpUHHKXUCP6IjIiIiB3ZWt4ac0qk+T366iNlrtgUdp1RRoSMiIhLlzIwHz+hI7crlGTp6BnsysoKOVGqo0BERESkFqldMZNjgLixL28WDE+YFHafUUKEjIiJSSvRuVZu/HN2c16eu4vMFG4OOUyqo0BERESlFbjqhLe3qVeGWd2axeWd60HGingodERGRUqR8QjxPndeV7XszueWdWeguLAemQkdERKSUaVuvCn8/sR2fL9jEGz+sCjpOVFOhIyIiUgr9qVczjm5dm398NI8lm3YGHSdqqdAREREpheLijCfO6UKFcvFcN3oGGZnZQUeKSip0RKRUMLN7zWytmaWGHicX0O56M5trZnPM7C0zS4p0VpFIqVs1iYfP7MyctdsZMWlR0HGikgodESlNhoemh0l29wl5V5pZQ+BaIMXdOwLx5EwiLBKzTuxYj3NTGvPCl0v5cfmWoONEHRU6IhJrEoAKZpYAVATWBZxHJOzuPq09TWpW5PrRqWzfuy/oOFFFhY6IlCZXm9ksM3vVzGrkXenua4EngFXAemCbJgeWsqBS+QRGnJvMhu17uXvcnKDjRBUVOiISNcxsUujcmryPQcALQEsgmZwiZlg+768BDAKaAw2ASmZ20QH2N8TMppnZtLS0tHD8lUQipmuTGlx7XGvGpa5jfOraoONEjYSgA4iI7OfuAwrTzsxeBj7MZ9UAYLm7p4XajQV6Aa8XsL+RwEiAlJQU3XVNSr2r+rXky0WbuHPcHFKa1aRh9QpBRwqcjuiISKlgZvVzPT0DyO/4/CrgKDOraGYG9AfmRyKfSDRIiI9j+LnJZGc7N4xOJStb9bsKHREpLR4zs9lmNgvoB1wPYGYNzGwCgLv/ALwD/AzMJqePGxlQXpFANK1ViXsGduCH5Vt4+etlQccJnIauRKRUcPeLC3h9HXByruf3APdEKpdINDqneyO+WLCJYZ8upE+r2nRsWC3oSIHRER0REZEYY2Y8dEYnalZKZOioGezJyAo6UmBU6IiIiMSgGpUSeeKcLixN28XDH5fdU9VU6IiIiMSoo1vX4bI+zfnP9yv5YsGmoOMEQoWOiIhIDLv5hLa0q1eFm9+Zxead6UHHiTgVOiIiIjEsqVw8I85LZvueffz93dm4l61LzlXoiIiIxLh29apyy4ltmTR/I2/9uDroOBGlQkdERKQM+HPv5vRpVZsHPpzHsrSdQceJGBU6IiIiZUBcnDFscBfKl4vjutGp7MvKDjpSRKjQERERKSPqVk3i4TM6MWvNNp6atDjoOBER8ULHzOLNbIaZ5TchH2Y22MzmmdlcM3sz0vlERERi2Umd6nNO90Y8P2UJP63YEnScsAviiM5QCphkz8xaA7cBvd29A3BdBHOJiIiUCfcM7ECjGhW5fnQq2/fuCzpOWEW00DGzRsApwCsFNPkL8Jy7/wrg7mXz7kYiIiJhVLl8AsPPTWbd1j3c+/7coOOEVaSP6IwAbgEKOgOqDdDGzL41s6lmdmLEkomIiJQh3ZvW4OrjWjP257V8OGtd0HHCJmKFjpmdCmxy9+kHaJYAtAb6AucDL5tZ9QK2N8TMppnZtLS0tJKOKyIiEvOuPa4VyY2rc/vY2azbuifoOGERySM6vYGBZrYCGAUcZ2av52mzBnjf3fe5+3JgETmFz++4+0h3T3H3lDp16oQzt4iISExKiI9jxLnJZGY7N46ZSXZ27N01OWKFjrvf5u6N3L0ZcB7wubtflKfZOHKO5mBmtckZyloWqYwiIiJlTbPalbjntPZ8v+wXXvkm9v7LDfw+OmZ2v5kNDD39BPjFzOYBXwA3u/svwaUTERGJfYNTGnNCh7o8/slC5q3bHnScEhVIoePuU9z91NDy3e7+fmjZ3f0Gd2/v7p3cfVQQ+URERMoSM+PhMztTo2IiQ0fNYO++rKAjlZjAj+iIiIhI8GpWSuSJc7qweNNOHvl4QdBxSowKHREREQHgmDZ1uLR3M/713QqmLIyNW9mp0BEREZH/ufXEdrSpW5mb35nFLzvTg45zyFToiIiIyP8klYtnxLld2bZ7H7eNnY176b7kXIWOiIiI/Eb7BlW5+YS2fDpvI6N/Wh10nEOiQkdERER+57I+zenVshb3fTCP5Zt3BR2n2FToiIiIyO/ExRnDBnchMSGO60ansi+roGkqo5sKHREREclX/WoVeOiMTsxcvZVnJi8OOk6xqNARERGRAp3SuT5ndWvEs18sYfrKLUHHKTIVOiJSapjZNWa2wMzmmtljBbQ50cwWmtkSM/t7pDOKxKJ7B7anYY0KXDc6lR179wUdp0hU6IhIqWBm/YBBQBd37wA8kU+beOA54CSgPXC+mbWPaFCRGFQlqRzDByez9tc93PfBvKDjFIkKHREpLf4GPOLu6QDunt9tW3sAS9x9mbtnAKPIKY5E5BClNKvJVf1a8c70NUyYvT7oOIWmQkdESos2wNFm9oOZfWlmR+TTpiGQ+6Yfa0KviUgJuLZ/a7o0rs5tY2ezftueoOMUSkLQAUrC9OnTN5vZykPcTG1gc0nkKQExn8XMivvWmP9siimasrQt7hvNbBJQL59Vd5DTX9UEjgKOAMaYWQs/hNu2mtkQYEjoabqZzSnutsIomn62uYU9VzH7iTL7eRVTsXM1uLdkg+RR7H4kr5godNy9zqFuw8ymuXtKSeQ5VMpSsGjKoyz5M7NpxX2vuw84wHb/BowNFTY/mlk2OZ10Wq5ma4HGuZ43Cr1W0P5GAiP3546WzzA35Soa5SqaaM5VUtvS0JWIlBbjgH4AZtYGSOT3v4n+BLQ2s+ZmlgicB7wfyZAiEl1U6IhIafEq0CI0vDQK+KO7u5k1MLMJAO6eCVwNfALMB8a4+9zAEotI4GJi6KqEjAw6QC7KUrBoyqMs+QtLltBVVBfl8/o64ORczycAE4qxi2j6DHNTrqJRrqKJ+VxW2qdfFxERESmIhq5EREQkZpXZQsfMHjCzWWaWamafmlmDfNokm9n3odvNzzKzc4PKEmr3RzNbHHr8MUxZHg/dYn+Wmb1nZtULaHd96HOZY2ZvmVlSgFmqm9k7obbzzaxnSWcpSp5Q23gzm2FmHwaVxcwam9kXZjYv9LMaGlSWULuonJrBzG40Mzez2gWsD/v3Ls/+CtsfZIXapJpZ2E+4jqZ+Ks/+Cvvvb4WZzQ7lL7ErekogV0S/F2Z2Tqg/yDazAq+2CuDzKmyuon9e7l4mH0DVXMvXAi/m06YN0Dq03ABYD1QPKEtNYFnozxqh5RphyHI8kBBafhR4NJ82DYHlQIXQ8zHAn4LIElr3b+Dy0HJiOH5GRckTWn8D8CbwYVBZgPpAt9ByFWAR0D6gLPHAUqBF6Gc0MxxZipG9MTknLq8EauezPiLfuzz7PGh/EFq3M8KfVdT0U0X99xdatyK/n3GQuYL4XgCHk3OPmilAygHaRfrzOmiu4n5eZfaIjrtvz/W0EvC7k5XcfZG7Lw4trwM2AYd8z57iZAFOAD5z9y3u/ivwGXBiGLJ86jlXrgBMJec+JPlJACqYWQJQEVgXRBYzqwYcA/wz9J4Md99a0lkKmyeUqRFwCvBKOHIUNou7r3f3n0PLO8i5CqnE7xJcyM8lWqdmGA7cQv7fOYjQ9y63QvYHERdN/VSeXIXtsyIqWr8X7j7f3ReGcx/FUchcxfq8ymyhA2BmD5rZauBC4O6DtO1BTgW5NKAsQdza/s/Ax3lfdPe15EyouIqco1zb3P3TILIAzcm5YdxroaGiV8ysUpizHCgPwAhy/vPMjkCOg2UBwMyaAV2BHwLKEnVTM5jZIGCtu888QLNAcheyb0oys2lmNtXMTg93pkLmCvrnfKDvggOfmtl0y7kjdiSVmu9FLkF+XgUp1ucV05eX2wFuJ+/u4939DuAOM7uNnHtv3FPAduoD/yXnvh3F+s+rpLKUhINlCbW5A8gE3sjn/TXIqaKbA1uBt83sInd/PdJZyPk33A24xt1/MLOngL8DdxU1S0nkMbNTgU3uPt3M+hYnQ0llybWdysC7wHV5fiuPeJZIO1Bu4HZyhhciroT6g6buvtbMWgCfm9lsdz+kX8SiqZ8qSq5Qm4P9++sT+rwOAz4zswXu/lUU5CpxhclVCIF8XuEQ04WOH+B28nm8Qc59N373pTWzqsBH5PwgpgaYZS3QN9fzRuSMZZZ4FjP7E3Aq0N9DA6N5DACWu3taqP1YoBdQ5EKnBLKsAda4+/4jFe+QU+gUSwnk6Q0MNLOTgSSgqpm97u6/u/9LBLJgZuXIKXLecPexRc1QglmKNDVDSSkot5l1IqdQn2k58yk1An42sx7uviFX0xL73hUmVz4K7JtCR1Zx92VmNoWcI3aHVOhEUz9VlFyF+S7k+rw2mdl75AyDHNJ/3NH6vSjCz/FA24j451UIxfu8DnYST6w+CJ1kHFq+BngnnzaJwGRyfhMOOktNck4ArhF6LAdqhiHLicA8oM4B2hwJzCXn3Bwj52Tga4LIEmr3NdA2tHwv8HiYfk6FypOrfV/CdzJyYX5OBvwHGBGODEXMkkDOianN+f+TCDuEM1cR/w4rKPhk5LB/7/LsszD9QQ2gfGi5NrCY8J/EGjX9VDH+/VUCquRa/g44MQpyBfa94MAn/Ub88ypkrmJ9XmEPHa0Pcn7LnQPMAj4AGoZeTwFeCS1fBOwDUnM9koPIEnr+Z2BJ6HFpmD6XJeSMge7/+74Yer0BMCFXu/uABaHc/93f6QaUJRmYFvr8xhGmqzwKmydX+76Er9A5aBagDznj7LNytTs5wJ/TyeRc+bWUnCOkJf65HMLfYQWhQieI712eLIXpm3oBs8np6GcDl0VDroA+r8J8F1qEPquZ5PySFvZ/f9H6vQDOIOdIeDqwEfgkSj6vg+Yq7uelOyOLiIhIzCrTV12JiIhIbFOhIyIiIjFLhY6IiIjELBU6IiIiErNU6IiIiEjMUqEjIiIiMUuFjoiIiMQsFToScWZWwcy+NLP4EthWopl9FZpFXUTKCPUjUlgqdKTEFaLj+TMw1t2zDnVf7p5BzjQd5x7qtkQkeqgfkZKiQkdKhJm9bWYvmdlU4LaDNL8Q2D+zbzMzm5NrOzeZ2b251i0ws3+Z2SIze8PMBpjZt2a22Mx6hN42LrRNESnF1I9IOKjQkZLSCdjo7ke5+z8KamRmiUALd19RyO22AoYB7UKPC8iZw+km4PZQmznAEcXMLSLRQ/2IlDgVOnLIzCyJnFmL7y9E89rA1iJsfrm7z3b3bHIml5vsORO0zQaaAYQOXWeYWZWi5BaR6KF+RMJFhY6UhA7AD+6eCWBmdczsNTNrZGavmlm5XG33AEl53m+5lsvlWZeeazk71/NsIPeJg+WBvcX9C4hI4NSPSFio0JGS0AmYtf+Ju6cBq8g5VHytu+/Lte5XID7029t+TUOdWhxwDFCkqyjMrBawOfd+RKTUUT8iYaFCR0rCbzooM6sMtAAy3X1nPu0/JWd8fL9fgP8A08gZJ7/EzFoWYf/9gI+KGlpEoor6EQkLyxmmFCkZoftQjATuAwYDP7n7lDxtugHXu/vFZtYM+NDdOx7CPscCf3f3RcUOLiJRQ/2IlCTdHElKVGh8/c+hp48X0OZnM/uipG70BYxT5yQSO9SPSEnSER0RERGJWTpHR0RERGKWCh0RERGJWSp0REREJGap0BEREZGYpUJHREREYpYKHREREYlZKnREREQkZqnQERERkZj1f5CO3pzspiyEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 648x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "\n",
    "\n",
    "for plot in [1, 2]:\n",
    "    plt.subplot(1, 2, plot)\n",
    "    plt.plot(simulation._bandstructure.r[0], simulation._bandstructure.r[1])\n",
    "    plt.plot(simulation._bandstructure.r[0], simulation._bandstructure.r[1], 'k.')\n",
    "    plt.xlabel(r'$r_x\\ (\\mathrm{\\mu m})$')\n",
    "    plt.ylabel(r'$r_y\\ (\\mathrm{\\mu m})$')\n",
    "\n",
    "    plt.axvline(int_r[0], color='k')\n",
    "    plt.plot(initial_r[0], initial_r[1], 'bo')\n",
    "    plt.plot(simulation._bandstructure.r[0][initial_n_f[0]+1], simulation._bandstructure.r[1][initial_n_f[0]+1], 'co')\n",
    "    plt.plot(int_r[0], int_r[1], 'go')\n",
    "\n",
    "    plt.plot(ref_r[0], ref_r[1], 'go')\n",
    "    plt.plot(simulation._bandstructure.r[0][final_n_f[0]], simulation._bandstructure.r[1][final_n_f[0]], \n",
    "            color='red', linestyle='dashdot', marker='o')\n",
    "    plt.plot(simulation._bandstructure.r[0][final_n_f[0]+1], simulation._bandstructure.r[1][final_n_f[0]+1], \n",
    "            color='red', linestyle='dashdot', marker='o')\n",
    "\n",
    "    if plot == 1:\n",
    "        plt.xlim(-3.3, -2)\n",
    "        plt.ylim(4.5, 6)\n",
    "\n",
    "    if plot == 2:\n",
    "        plt.xlim(-4, -1)\n",
    "        plt.ylim(-6, -4.5)\n",
    "\n",
    "\n",
    "\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.8.2 64-bit",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}