02-MPB_dispersion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MPB ModeSolver\n",
    "\n",
    "Now let's obtain the plot on the right (Figure from [3]) :\n",
    "\n",
    "![alt text](http://www.simpetus.com/files/SOI_waveguide_bands.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import meep and mpb (from meep)\n",
    "import meep as mp\n",
    "from meep import mpb\n",
    "\n",
    "# arrays\n",
    "import numpy as np\n",
    "\n",
    "# plotting\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by defining the geometry from last time :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters for the waveguide\n",
    "w = 0.50  # Si width (um)\n",
    "h = 0.22  # Si height (um)\n",
    "\n",
    "# Materials\n",
    "Si = mp.Medium(index=3.45)\n",
    "SiO2 = mp.Medium(index=1.45)\n",
    "\n",
    "sc_y = 2  # supercell width (um)\n",
    "sc_z = 2  # supercell height (um)\n",
    "resolution = 32  # pixels/um\n",
    "geometry_lattice = mp.Lattice(size=mp.Vector3(0, sc_y, sc_z))\n",
    "\n",
    "# define the 2d blocks for the strip and substrate\n",
    "geometry = [mp.Block(size=mp.Vector3(mp.inf, w, h),\n",
    "                     center=mp.Vector3(), material=Si),\n",
    "            mp.Block(size=mp.Vector3(mp.inf, mp.inf, 0.5*(sc_z-h)),\n",
    "                     center=mp.Vector3(z=0.25*(sc_z+h)), material=SiO2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the simulation\n",
    "\n",
    "Now, instead of only initializing the geometry, we will load simulation parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The k (i.e. beta, i.e. propagation constant) points to look at, in\n",
    "# units of 2*pi/um.  We'll look at num_k points from k_min to k_max.\n",
    "num_k = 20\n",
    "k_min = 0.1\n",
    "k_max = 2.0\n",
    "k_points = mp.interpolate(num_k, [mp.Vector3(k_min), mp.Vector3(k_max)])\n",
    "\n",
    "# Increase this to see more modes.  (The guided ones are the ones below the\n",
    "# light line, i.e. those with frequencies < kmag / 1.45, where kmag\n",
    "# is the corresponding column in the output if you grep for \"freqs:\".)\n",
    "num_bands = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ModeSolver object\n",
    "ms = mpb.ModeSolver(\n",
    "    geometry_lattice=geometry_lattice,\n",
    "    geometry=geometry,\n",
    "    # Add new things pertaining to simulation\n",
    "    k_points=k_points,\n",
    "    resolution=resolution,\n",
    "    num_bands=num_bands,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now run the simulation. There are many options for what to solve for and what to report, see https://mpb.readthedocs.io/en/latest/Python_User_Interface/#run-functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing eigensolver data\n",
      "Computing 4 bands with 1e-07 tolerance\n",
      "Working in 3 dimensions.\n",
      "Grid size is 1 x 64 x 64.\n",
      "Solving for 4 bands at a time.\n",
      "Creating Maxwell data...\n",
      "Mesh size is 3.\n",
      "Lattice vectors:\n",
      "     (1, 0, 0)\n",
      "     (0, 2, 0)\n",
      "     (0, 0, 2)\n",
      "Cell volume = 4\n",
      "Reciprocal lattice vectors (/ 2 pi):\n",
      "     (1, -0, 0)\n",
      "     (-0, 0.5, -0)\n",
      "     (0, -0, 0.5)\n",
      "Geometric objects:\n",
      "     block, center = (0,0,0)\n",
      "          size (1e+20,0.5,0.22)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     block, center = (0,0,0.555)\n",
      "          size (1e+20,1e+20,0.89)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "Geometric object tree has depth 2 and 8 object nodes (vs. 2 actual objects)\n",
      "Initializing epsilon function...\n",
      "Allocating fields...\n",
      "Solving for band polarization: .\n",
      "Initializing fields to random numbers...\n",
      "22 k-points\n",
      "  Vector3<0.1, 0.0, 0.0>\n",
      "  Vector3<0.19047619047619047, 0.0, 0.0>\n",
      "  Vector3<0.28095238095238095, 0.0, 0.0>\n",
      "  Vector3<0.37142857142857144, 0.0, 0.0>\n",
      "  Vector3<0.46190476190476193, 0.0, 0.0>\n",
      "  Vector3<0.5523809523809524, 0.0, 0.0>\n",
      "  Vector3<0.6428571428571428, 0.0, 0.0>\n",
      "  Vector3<0.7333333333333333, 0.0, 0.0>\n",
      "  Vector3<0.8238095238095238, 0.0, 0.0>\n",
      "  Vector3<0.9142857142857143, 0.0, 0.0>\n",
      "  Vector3<1.0047619047619047, 0.0, 0.0>\n",
      "  Vector3<1.0952380952380953, 0.0, 0.0>\n",
      "  Vector3<1.1857142857142857, 0.0, 0.0>\n",
      "  Vector3<1.2761904761904763, 0.0, 0.0>\n",
      "  Vector3<1.3666666666666667, 0.0, 0.0>\n",
      "  Vector3<1.4571428571428573, 0.0, 0.0>\n",
      "  Vector3<1.5476190476190477, 0.0, 0.0>\n",
      "  Vector3<1.638095238095238, 0.0, 0.0>\n",
      "  Vector3<1.7285714285714286, 0.0, 0.0>\n",
      "  Vector3<1.819047619047619, 0.0, 0.0>\n",
      "  Vector3<1.9095238095238096, 0.0, 0.0>\n",
      "  Vector3<2.0, 0.0, 0.0>\n",
      "elapsed time for initialization: 0.015517473220825195\n",
      "solve_kpoint (0.1,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 30 iterations.\n",
      "freqs:, 1, 0.1, 0, 0, 0.1, 0.0794202, 0.0849623, 0.297663, 0.37481\n",
      "elapsed time for k point: 0.13807034492492676\n",
      "solve_kpoint (0.190476,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 14 iterations.\n",
      "freqs:, 2, 0.190476, 0, 0, 0.190476, 0.150218, 0.160326, 0.324018, 0.383469\n",
      "elapsed time for k point: 0.05554652214050293\n",
      "solve_kpoint (0.280952,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 11 iterations.\n",
      "freqs:, 3, 0.280952, 0, 0, 0.280952, 0.219035, 0.232794, 0.36303, 0.398284\n",
      "elapsed time for k point: 0.046803951263427734\n",
      "solve_kpoint (0.371429,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 11 iterations.\n",
      "freqs:, 4, 0.371429, 0, 0, 0.371429, 0.284828, 0.301118, 0.410931, 0.417364\n",
      "elapsed time for k point: 0.04655337333679199\n",
      "solve_kpoint (0.461905,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 13 iterations.\n",
      "freqs:, 5, 0.461905, 0, 0, 0.461905, 0.34589, 0.36477, 0.440364, 0.463544\n",
      "elapsed time for k point: 0.052527427673339844\n",
      "solve_kpoint (0.552381,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 15 iterations.\n",
      "freqs:, 6, 0.552381, 0, 0, 0.552381, 0.398211, 0.424058, 0.470557, 0.516771\n",
      "elapsed time for k point: 0.05847883224487305\n",
      "solve_kpoint (0.642857,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 16 iterations.\n",
      "freqs:, 7, 0.642857, 0, 0, 0.642857, 0.436446, 0.479309, 0.513457, 0.567861\n",
      "elapsed time for k point: 0.06135964393615723\n",
      "solve_kpoint (0.733333,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 17 iterations.\n",
      "freqs:, 8, 0.733333, 0, 0, 0.733333, 0.463503, 0.529944, 0.566628, 0.615651\n",
      "elapsed time for k point: 0.06601095199584961\n",
      "solve_kpoint (0.82381,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 18 iterations.\n",
      "freqs:, 9, 0.82381, 0, 0, 0.82381, 0.485926, 0.574409, 0.623957, 0.66078\n",
      "elapsed time for k point: 0.06747317314147949\n",
      "solve_kpoint (0.914286,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 21 iterations.\n",
      "freqs:, 10, 0.914286, 0, 0, 0.914286, 0.506776, 0.611532, 0.682582, 0.698461\n",
      "elapsed time for k point: 0.07629990577697754\n",
      "solve_kpoint (1.00476,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 21 iterations.\n",
      "freqs:, 11, 1.00476, 0, 0, 1.00476, 0.527159, 0.642054, 0.725485, 0.740367\n",
      "elapsed time for k point: 0.07693672180175781\n",
      "solve_kpoint (1.09524,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 24 iterations.\n",
      "freqs:, 12, 1.09524, 0, 0, 1.09524, 0.547509, 0.667911, 0.747291, 0.775545\n",
      "elapsed time for k point: 0.08504891395568848\n",
      "solve_kpoint (1.18571,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 13 iterations.\n",
      "freqs:, 13, 1.18571, 0, 0, 1.18571, 0.568008, 0.690839, 0.766489, 0.796961\n",
      "elapsed time for k point: 0.05281376838684082\n",
      "solve_kpoint (1.27619,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 11 iterations.\n",
      "freqs:, 14, 1.27619, 0, 0, 1.27619, 0.588733, 0.711999, 0.784403, 0.816743\n",
      "elapsed time for k point: 0.046945810317993164\n",
      "solve_kpoint (1.36667,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 10 iterations.\n",
      "freqs:, 15, 1.36667, 0, 0, 1.36667, 0.609713, 0.732111, 0.80173, 0.835754\n",
      "elapsed time for k point: 0.043984174728393555\n",
      "solve_kpoint (1.45714,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 10 iterations.\n",
      "freqs:, 16, 1.45714, 0, 0, 1.45714, 0.630952, 0.751623, 0.818853, 0.854324\n",
      "elapsed time for k point: 0.04348158836364746\n",
      "solve_kpoint (1.54762,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 10 iterations.\n",
      "freqs:, 17, 1.54762, 0, 0, 1.54762, 0.652443, 0.770818, 0.83599, 0.872643\n",
      "elapsed time for k point: 0.0436398983001709\n",
      "solve_kpoint (1.6381,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 9 iterations.\n",
      "freqs:, 18, 1.6381, 0, 0, 1.6381, 0.674172, 0.78988, 0.853267, 0.890837\n",
      "elapsed time for k point: 0.04073929786682129\n",
      "solve_kpoint (1.72857,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 9 iterations.\n",
      "freqs:, 19, 1.72857, 0, 0, 1.72857, 0.696125, 0.808928, 0.870757, 0.908996\n",
      "elapsed time for k point: 0.04059171676635742\n",
      "solve_kpoint (1.81905,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 9 iterations.\n",
      "freqs:, 20, 1.81905, 0, 0, 1.81905, 0.718285, 0.828044, 0.888503, 0.927185\n",
      "elapsed time for k point: 0.0409388542175293\n",
      "solve_kpoint (1.90952,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 9 iterations.\n",
      "freqs:, 21, 1.90952, 0, 0, 1.90952, 0.740637, 0.847281, 0.906527, 0.945451\n",
      "elapsed time for k point: 0.04106616973876953\n",
      "solve_kpoint (2,0,0):\n",
      "Solving for bands 1 to 4...\n",
      "Finished solving for bands 1 to 4 after 9 iterations.\n",
      "freqs:, 22, 2, 0, 0, 2, 0.763167, 0.866673, 0.924839, 0.963831\n",
      "elapsed time for k point: 0.04073214530944824\n",
      "Band 1 range: 0.07942016624643601 at Vector3<0.1, 0.0, 0.0> to 0.7631671678430058 at Vector3<2.0, 0.0, 0.0>\n",
      "Band 2 range: 0.08496230355143654 at Vector3<0.1, 0.0, 0.0> to 0.8666725405654844 at Vector3<2.0, 0.0, 0.0>\n",
      "Band 3 range: 0.29766264003036447 at Vector3<0.1, 0.0, 0.0> to 0.9248390743934864 at Vector3<2.0, 0.0, 0.0>\n",
      "Band 4 range: 0.37480991768534555 at Vector3<0.1, 0.0, 0.0> to 0.9638311239052704 at Vector3<2.0, 0.0, 0.0>\n",
      "total elapsed time for run: 1.2828874588012695\n",
      "done\n"
     ]
    }
   ],
   "source": [
    "ms.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It solved, great! We can see from the text output lines like :\n",
    "\n",
    "`solve_kpoint (0.1,0,0):\n",
    "Solving for bands 1 to 4...\n",
    "Finished solving for bands 1 to 4 after 28 iterations.\n",
    "freqs:, 1, 0.1, 0, 0, 0.1, 0.0807582, 0.0855765, 0.3147, 0.394675\n",
    "elapsed time for k point: 0.14001059532165527`\n",
    "\n",
    "This contains our band information. Unfortunately, MPB only keeps in memory one at a time. So we only have access to the last one :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vector3<2.0, 0.0, 0.0>\n",
      "[0.7631671678430061, 0.8666725405654846, 0.9248390743934869, 0.9638311239052703]\n"
     ]
    }
   ],
   "source": [
    "print(ms.current_k)\n",
    "print(ms.freqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's do it again, this time capturing the output :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture --no-stderr output\n",
    "ms.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Initializing eigensolver data',\n",
       " 'Computing 4 bands with 1e-07 tolerance',\n",
       " 'Working in 3 dimensions.',\n",
       " 'Grid size is 1 x 64 x 64.',\n",
       " 'Solving for 4 bands at a time.',\n",
       " 'Creating Maxwell data...',\n",
       " 'Mesh size is 3.',\n",
       " 'Lattice vectors:',\n",
       " '     (1, 0, 0)',\n",
       " '     (0, 2, 0)',\n",
       " '     (0, 0, 2)',\n",
       " 'Cell volume = 4',\n",
       " 'Reciprocal lattice vectors (/ 2 pi):',\n",
       " '     (1, -0, 0)',\n",
       " '     (-0, 0.5, -0)',\n",
       " '     (0, -0, 0.5)',\n",
       " 'Geometric objects:',\n",
       " '     block, center = (0,0,0)',\n",
       " '          size (1e+20,0.5,0.22)',\n",
       " '          axes (1,0,0), (0,1,0), (0,0,1)',\n",
       " '     block, center = (0,0,0.555)',\n",
       " '          size (1e+20,1e+20,0.89)',\n",
       " '          axes (1,0,0), (0,1,0), (0,0,1)',\n",
       " 'Geometric object tree has depth 2 and 8 object nodes (vs. 2 actual objects)',\n",
       " 'Initializing epsilon function...',\n",
       " 'Solving for band polarization: .',\n",
       " 'Initializing fields to random numbers...',\n",
       " '22 k-points',\n",
       " '  Vector3<0.1, 0.0, 0.0>',\n",
       " '  Vector3<0.19047619047619047, 0.0, 0.0>',\n",
       " '  Vector3<0.28095238095238095, 0.0, 0.0>',\n",
       " '  Vector3<0.37142857142857144, 0.0, 0.0>',\n",
       " '  Vector3<0.46190476190476193, 0.0, 0.0>',\n",
       " '  Vector3<0.5523809523809524, 0.0, 0.0>',\n",
       " '  Vector3<0.6428571428571428, 0.0, 0.0>',\n",
       " '  Vector3<0.7333333333333333, 0.0, 0.0>',\n",
       " '  Vector3<0.8238095238095238, 0.0, 0.0>',\n",
       " '  Vector3<0.9142857142857143, 0.0, 0.0>',\n",
       " '  Vector3<1.0047619047619047, 0.0, 0.0>',\n",
       " '  Vector3<1.0952380952380953, 0.0, 0.0>',\n",
       " '  Vector3<1.1857142857142857, 0.0, 0.0>',\n",
       " '  Vector3<1.2761904761904763, 0.0, 0.0>',\n",
       " '  Vector3<1.3666666666666667, 0.0, 0.0>',\n",
       " '  Vector3<1.4571428571428573, 0.0, 0.0>',\n",
       " '  Vector3<1.5476190476190477, 0.0, 0.0>',\n",
       " '  Vector3<1.638095238095238, 0.0, 0.0>',\n",
       " '  Vector3<1.7285714285714286, 0.0, 0.0>',\n",
       " '  Vector3<1.819047619047619, 0.0, 0.0>',\n",
       " '  Vector3<1.9095238095238096, 0.0, 0.0>',\n",
       " '  Vector3<2.0, 0.0, 0.0>',\n",
       " 'elapsed time for initialization: 0.0075299739837646484',\n",
       " 'solve_kpoint (0.1,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 27 iterations.',\n",
       " 'freqs:, 1, 0.1, 0, 0, 0.1, 0.0794202, 0.0849623, 0.297663, 0.37481',\n",
       " 'elapsed time for k point: 0.08509659767150879',\n",
       " 'solve_kpoint (0.190476,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 14 iterations.',\n",
       " 'freqs:, 2, 0.190476, 0, 0, 0.190476, 0.150218, 0.160326, 0.324018, 0.383469',\n",
       " 'elapsed time for k point: 0.047899723052978516',\n",
       " 'solve_kpoint (0.280952,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 11 iterations.',\n",
       " 'freqs:, 3, 0.280952, 0, 0, 0.280952, 0.219035, 0.232794, 0.36303, 0.398284',\n",
       " 'elapsed time for k point: 0.03527998924255371',\n",
       " 'solve_kpoint (0.371429,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 11 iterations.',\n",
       " 'freqs:, 4, 0.371429, 0, 0, 0.371429, 0.284828, 0.301118, 0.410931, 0.417364',\n",
       " 'elapsed time for k point: 0.0352778434753418',\n",
       " 'solve_kpoint (0.461905,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 13 iterations.',\n",
       " 'freqs:, 5, 0.461905, 0, 0, 0.461905, 0.34589, 0.36477, 0.440364, 0.463544',\n",
       " 'elapsed time for k point: 0.04119586944580078',\n",
       " 'solve_kpoint (0.552381,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 15 iterations.',\n",
       " 'freqs:, 6, 0.552381, 0, 0, 0.552381, 0.398211, 0.424058, 0.470557, 0.516771',\n",
       " 'elapsed time for k point: 0.04709887504577637',\n",
       " 'solve_kpoint (0.642857,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 16 iterations.',\n",
       " 'freqs:, 7, 0.642857, 0, 0, 0.642857, 0.436446, 0.479309, 0.513457, 0.567861',\n",
       " 'elapsed time for k point: 0.050104618072509766',\n",
       " 'solve_kpoint (0.733333,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 17 iterations.',\n",
       " 'freqs:, 8, 0.733333, 0, 0, 0.733333, 0.463503, 0.529944, 0.566628, 0.615651',\n",
       " 'elapsed time for k point: 0.0529627799987793',\n",
       " 'solve_kpoint (0.82381,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 18 iterations.',\n",
       " 'freqs:, 9, 0.82381, 0, 0, 0.82381, 0.485926, 0.574409, 0.623957, 0.66078',\n",
       " 'elapsed time for k point: 0.05662083625793457',\n",
       " 'solve_kpoint (0.914286,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 21 iterations.',\n",
       " 'freqs:, 10, 0.914286, 0, 0, 0.914286, 0.506776, 0.611532, 0.682582, 0.698461',\n",
       " 'elapsed time for k point: 0.06470108032226562',\n",
       " 'solve_kpoint (1.00476,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 21 iterations.',\n",
       " 'freqs:, 11, 1.00476, 0, 0, 1.00476, 0.527159, 0.642054, 0.725485, 0.740367',\n",
       " 'elapsed time for k point: 0.06516051292419434',\n",
       " 'solve_kpoint (1.09524,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 24 iterations.',\n",
       " 'freqs:, 12, 1.09524, 0, 0, 1.09524, 0.547509, 0.667911, 0.747291, 0.775545',\n",
       " 'elapsed time for k point: 0.07437825202941895',\n",
       " 'solve_kpoint (1.18571,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 13 iterations.',\n",
       " 'freqs:, 13, 1.18571, 0, 0, 1.18571, 0.568008, 0.690839, 0.766489, 0.796961',\n",
       " 'elapsed time for k point: 0.04303407669067383',\n",
       " 'solve_kpoint (1.27619,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 11 iterations.',\n",
       " 'freqs:, 14, 1.27619, 0, 0, 1.27619, 0.588733, 0.711999, 0.784403, 0.816743',\n",
       " 'elapsed time for k point: 0.0406804084777832',\n",
       " 'solve_kpoint (1.36667,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 10 iterations.',\n",
       " 'freqs:, 15, 1.36667, 0, 0, 1.36667, 0.609713, 0.732111, 0.80173, 0.835754',\n",
       " 'elapsed time for k point: 0.03243279457092285',\n",
       " 'solve_kpoint (1.45714,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 10 iterations.',\n",
       " 'freqs:, 16, 1.45714, 0, 0, 1.45714, 0.630952, 0.751623, 0.818853, 0.854324',\n",
       " 'elapsed time for k point: 0.03248262405395508',\n",
       " 'solve_kpoint (1.54762,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 10 iterations.',\n",
       " 'freqs:, 17, 1.54762, 0, 0, 1.54762, 0.652443, 0.770818, 0.83599, 0.872643',\n",
       " 'elapsed time for k point: 0.03250932693481445',\n",
       " 'solve_kpoint (1.6381,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 9 iterations.',\n",
       " 'freqs:, 18, 1.6381, 0, 0, 1.6381, 0.674172, 0.78988, 0.853267, 0.890837',\n",
       " 'elapsed time for k point: 0.029952049255371094',\n",
       " 'solve_kpoint (1.72857,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 9 iterations.',\n",
       " 'freqs:, 19, 1.72857, 0, 0, 1.72857, 0.696125, 0.808928, 0.870757, 0.908996',\n",
       " 'elapsed time for k point: 0.029516935348510742',\n",
       " 'solve_kpoint (1.81905,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 9 iterations.',\n",
       " 'freqs:, 20, 1.81905, 0, 0, 1.81905, 0.718285, 0.828044, 0.888503, 0.927185',\n",
       " 'elapsed time for k point: 0.029439926147460938',\n",
       " 'solve_kpoint (1.90952,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 9 iterations.',\n",
       " 'freqs:, 21, 1.90952, 0, 0, 1.90952, 0.740637, 0.847281, 0.906527, 0.945451',\n",
       " 'elapsed time for k point: 0.0294950008392334',\n",
       " 'solve_kpoint (2,0,0):',\n",
       " 'Solving for bands 1 to 4...',\n",
       " 'Finished solving for bands 1 to 4 after 9 iterations.',\n",
       " 'freqs:, 22, 2, 0, 0, 2, 0.763167, 0.866673, 0.924839, 0.963831',\n",
       " 'elapsed time for k point: 0.029565095901489258',\n",
       " 'Band 1 range: 0.07942016365368101 at Vector3<0.1, 0.0, 0.0> to 0.7631671678430061 at Vector3<2.0, 0.0, 0.0>',\n",
       " 'Band 2 range: 0.08496230254981353 at Vector3<0.1, 0.0, 0.0> to 0.8666725405654854 at Vector3<2.0, 0.0, 0.0>',\n",
       " 'Band 3 range: 0.29766263809524895 at Vector3<0.1, 0.0, 0.0> to 0.924839074393488 at Vector3<2.0, 0.0, 0.0>',\n",
       " 'Band 4 range: 0.37480990708126916 at Vector3<0.1, 0.0, 0.0> to 0.9638311239052701 at Vector3<2.0, 0.0, 0.0>',\n",
       " 'total elapsed time for run: 0.9928672313690186',\n",
       " 'done',\n",
       " '']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output_text = output.stdout\n",
    "output_text.split('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyzing the results\n",
    "\n",
    "Let's parse every line that begins with 'freqs :' in the output.\n",
    "\n",
    "Notice the lines all have the same syntax `'freqs:, 22, 2, 0, 0, 2, 0.846477, 0.874081, 1.04852, 1.16543',`\n",
    "\n",
    "i.e. for all the lines that start by 'freqs:' contains all the information we need in the form `freqs:, k index, kx, ky, kz, kmag/2pi, band 1, band 2, band 3, band 4, band 5, band 6, band 7, band 8 ...`. Python's `split` function and treating strings as arrays (slicing and indexing) is handy :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "freqs:, 1, 0.1, 0, 0, 0.1, 0.0794202, 0.0849623, 0.297663, 0.37481\n",
      "['freqs:', '1', '0.1', '0', '0', '0.1', '0.0794202', '0.0849623', '0.297663', '0.37481']\n",
      "freqs:, 2, 0.190476, 0, 0, 0.190476, 0.150218, 0.160326, 0.324018, 0.383469\n",
      "['freqs:', '2', '0.190476', '0', '0', '0.190476', '0.150218', '0.160326', '0.324018', '0.383469']\n",
      "freqs:, 3, 0.280952, 0, 0, 0.280952, 0.219035, 0.232794, 0.36303, 0.398284\n",
      "['freqs:', '3', '0.280952', '0', '0', '0.280952', '0.219035', '0.232794', '0.36303', '0.398284']\n",
      "freqs:, 4, 0.371429, 0, 0, 0.371429, 0.284828, 0.301118, 0.410931, 0.417364\n",
      "['freqs:', '4', '0.371429', '0', '0', '0.371429', '0.284828', '0.301118', '0.410931', '0.417364']\n",
      "freqs:, 5, 0.461905, 0, 0, 0.461905, 0.34589, 0.36477, 0.440364, 0.463544\n",
      "['freqs:', '5', '0.461905', '0', '0', '0.461905', '0.34589', '0.36477', '0.440364', '0.463544']\n",
      "freqs:, 6, 0.552381, 0, 0, 0.552381, 0.398211, 0.424058, 0.470557, 0.516771\n",
      "['freqs:', '6', '0.552381', '0', '0', '0.552381', '0.398211', '0.424058', '0.470557', '0.516771']\n",
      "freqs:, 7, 0.642857, 0, 0, 0.642857, 0.436446, 0.479309, 0.513457, 0.567861\n",
      "['freqs:', '7', '0.642857', '0', '0', '0.642857', '0.436446', '0.479309', '0.513457', '0.567861']\n",
      "freqs:, 8, 0.733333, 0, 0, 0.733333, 0.463503, 0.529944, 0.566628, 0.615651\n",
      "['freqs:', '8', '0.733333', '0', '0', '0.733333', '0.463503', '0.529944', '0.566628', '0.615651']\n",
      "freqs:, 9, 0.82381, 0, 0, 0.82381, 0.485926, 0.574409, 0.623957, 0.66078\n",
      "['freqs:', '9', '0.82381', '0', '0', '0.82381', '0.485926', '0.574409', '0.623957', '0.66078']\n",
      "freqs:, 10, 0.914286, 0, 0, 0.914286, 0.506776, 0.611532, 0.682582, 0.698461\n",
      "['freqs:', '10', '0.914286', '0', '0', '0.914286', '0.506776', '0.611532', '0.682582', '0.698461']\n",
      "freqs:, 11, 1.00476, 0, 0, 1.00476, 0.527159, 0.642054, 0.725485, 0.740367\n",
      "['freqs:', '11', '1.00476', '0', '0', '1.00476', '0.527159', '0.642054', '0.725485', '0.740367']\n",
      "freqs:, 12, 1.09524, 0, 0, 1.09524, 0.547509, 0.667911, 0.747291, 0.775545\n",
      "['freqs:', '12', '1.09524', '0', '0', '1.09524', '0.547509', '0.667911', '0.747291', '0.775545']\n",
      "freqs:, 13, 1.18571, 0, 0, 1.18571, 0.568008, 0.690839, 0.766489, 0.796961\n",
      "['freqs:', '13', '1.18571', '0', '0', '1.18571', '0.568008', '0.690839', '0.766489', '0.796961']\n",
      "freqs:, 14, 1.27619, 0, 0, 1.27619, 0.588733, 0.711999, 0.784403, 0.816743\n",
      "['freqs:', '14', '1.27619', '0', '0', '1.27619', '0.588733', '0.711999', '0.784403', '0.816743']\n",
      "freqs:, 15, 1.36667, 0, 0, 1.36667, 0.609713, 0.732111, 0.80173, 0.835754\n",
      "['freqs:', '15', '1.36667', '0', '0', '1.36667', '0.609713', '0.732111', '0.80173', '0.835754']\n",
      "freqs:, 16, 1.45714, 0, 0, 1.45714, 0.630952, 0.751623, 0.818853, 0.854324\n",
      "['freqs:', '16', '1.45714', '0', '0', '1.45714', '0.630952', '0.751623', '0.818853', '0.854324']\n",
      "freqs:, 17, 1.54762, 0, 0, 1.54762, 0.652443, 0.770818, 0.83599, 0.872643\n",
      "['freqs:', '17', '1.54762', '0', '0', '1.54762', '0.652443', '0.770818', '0.83599', '0.872643']\n",
      "freqs:, 18, 1.6381, 0, 0, 1.6381, 0.674172, 0.78988, 0.853267, 0.890837\n",
      "['freqs:', '18', '1.6381', '0', '0', '1.6381', '0.674172', '0.78988', '0.853267', '0.890837']\n",
      "freqs:, 19, 1.72857, 0, 0, 1.72857, 0.696125, 0.808928, 0.870757, 0.908996\n",
      "['freqs:', '19', '1.72857', '0', '0', '1.72857', '0.696125', '0.808928', '0.870757', '0.908996']\n",
      "freqs:, 20, 1.81905, 0, 0, 1.81905, 0.718285, 0.828044, 0.888503, 0.927185\n",
      "['freqs:', '20', '1.81905', '0', '0', '1.81905', '0.718285', '0.828044', '0.888503', '0.927185']\n",
      "freqs:, 21, 1.90952, 0, 0, 1.90952, 0.740637, 0.847281, 0.906527, 0.945451\n",
      "['freqs:', '21', '1.90952', '0', '0', '1.90952', '0.740637', '0.847281', '0.906527', '0.945451']\n",
      "freqs:, 22, 2, 0, 0, 2, 0.763167, 0.866673, 0.924839, 0.963831\n",
      "['freqs:', '22', '2', '0', '0', '2', '0.763167', '0.866673', '0.924839', '0.963831']\n"
     ]
    }
   ],
   "source": [
    "freqs = []\n",
    "ks = []\n",
    "line_header = \"freqs:\"\n",
    "for ln in output_text.split('\\n'):\n",
    "    if ln.startswith(line_header):\n",
    "        print(ln)\n",
    "        line = ln.split(', ')\n",
    "        print(line)\n",
    "        ks.append([line[2], line[3], line[4]])\n",
    "        freqs.append([line[-4], line[-3], line[-2], line[-1]])\n",
    "        \n",
    "freqs = np.array(freqs, dtype=np.float32)\n",
    "ks = np.array(ks, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0794202 0.0849623 0.297663  0.37481  ]\n",
      " [0.150218  0.160326  0.324018  0.383469 ]\n",
      " [0.219035  0.232794  0.36303   0.398284 ]\n",
      " [0.284828  0.301118  0.410931  0.417364 ]\n",
      " [0.34589   0.36477   0.440364  0.463544 ]\n",
      " [0.398211  0.424058  0.470557  0.516771 ]\n",
      " [0.436446  0.479309  0.513457  0.567861 ]\n",
      " [0.463503  0.529944  0.566628  0.615651 ]\n",
      " [0.485926  0.574409  0.623957  0.66078  ]\n",
      " [0.506776  0.611532  0.682582  0.698461 ]\n",
      " [0.527159  0.642054  0.725485  0.740367 ]\n",
      " [0.547509  0.667911  0.747291  0.775545 ]\n",
      " [0.568008  0.690839  0.766489  0.796961 ]\n",
      " [0.588733  0.711999  0.784403  0.816743 ]\n",
      " [0.609713  0.732111  0.80173   0.835754 ]\n",
      " [0.630952  0.751623  0.818853  0.854324 ]\n",
      " [0.652443  0.770818  0.83599   0.872643 ]\n",
      " [0.674172  0.78988   0.853267  0.890837 ]\n",
      " [0.696125  0.808928  0.870757  0.908996 ]\n",
      " [0.718285  0.828044  0.888503  0.927185 ]\n",
      " [0.740637  0.847281  0.906527  0.945451 ]\n",
      " [0.763167  0.866673  0.924839  0.963831 ]]\n",
      "[[0.1      0.       0.      ]\n",
      " [0.190476 0.       0.      ]\n",
      " [0.280952 0.       0.      ]\n",
      " [0.371429 0.       0.      ]\n",
      " [0.461905 0.       0.      ]\n",
      " [0.552381 0.       0.      ]\n",
      " [0.642857 0.       0.      ]\n",
      " [0.733333 0.       0.      ]\n",
      " [0.82381  0.       0.      ]\n",
      " [0.914286 0.       0.      ]\n",
      " [1.00476  0.       0.      ]\n",
      " [1.09524  0.       0.      ]\n",
      " [1.18571  0.       0.      ]\n",
      " [1.27619  0.       0.      ]\n",
      " [1.36667  0.       0.      ]\n",
      " [1.45714  0.       0.      ]\n",
      " [1.54762  0.       0.      ]\n",
      " [1.6381   0.       0.      ]\n",
      " [1.72857  0.       0.      ]\n",
      " [1.81905  0.       0.      ]\n",
      " [1.90952  0.       0.      ]\n",
      " [2.       0.       0.      ]]\n"
     ]
    }
   ],
   "source": [
    "print(freqs)\n",
    "print(ks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot. Let's add the lightline too :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2aee299cc210>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEOCAYAAAB1g0unAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3gU5dqH70nvvQOhpdED6dTQCVKkBQt68IjlUxFUPEeOYgMUK+dYEVABC9KkKqEFpAUChAApBEIP6Qnp2Wx2d74/JgxEWoAsmzL3dXHpvDu78w7s7m/f93me3yOIooiCgoKCgkJdMDL0BBQUFBQUGg+KaCgoKCgo1BlFNBQUFBQU6owiGgoKCgoKdUYRDQUFBQWFOqOIhoKCgoJCndGbaAiC8IMgCLmCICTd4nFBEIQvBEFIFwThuCAIPfQ1FwUFBQWF+kGfK40lwLDbPB4F+Nb8eRb4Vo9zUVBQUFCoB/QmGqIo7gYKb3PKaGCZKHEAcBAEwVNf81FQUFBQuH8MGdNoAVy67jijZkxBQUFBoYFiYsBrCzcZu6mniSAIzyJtYWFtbR0UEBCgz3kpKNQZrVZLXl4eOp0OAHNzc5ydnQ08q4ZJfn4+Fy5ckI9NTEzw8fHB2tragLNqPhw5ciRfFEXX+30dQ4pGBtDquuOWQObNThRFcSGwECA4OFg8fPiw/menoHAHtFotS5YsISMjAwBra2uee+45bG1tDTyzhoUoisyZM4e3335bHmvbti0xMTH4+fkZcGbNC0EQLtz5rDtjyO2pDcCTNVlU4UCxKIpZBpyPgsJdsW3bNlkwBEFg/PjximD8Da1Wy//93//VEozu3buzf/9+RTAaKXpbaQiCsByIBFwEQcgA3gFMAURRXAD8CQwH0oEK4Cl9zUVBob5JTk7m4MGD8vHAgQNp06aN4SbUAKmsrOTRRx9l/fr18tjgwYNZs2aNIq6NGL2JhiiKj97hcRF4UV/XV1DQF/n5+WzYsEE+9vf3p2fPngacUcOjsLCQkSNHsn//fnns8ccf54cffsDMzMyAM1O4XwwZ06g3qqurycjIQKVSGXoqesXCwoKWLVtiampq6Kk0W9RqNStXrkStVgPg6OjIww8/jCDcLK+jeXLx4kWGDRtGamqqPPb6668zb948jIwUE4rGTpMQjYyMDGxtbWnTpk2T/fCKokhBQQEZGRm0bdvW0NNploiiyKZNm8jLywOk7J/o6GgsLCwMPLOGw4kTJxg2bBiZmddyWubPn8/06dMNOKvmjajRUHniRL29XpMQDZVK1aQFA6RAq7Ozs/yFpfDgOXz4MCeu+/ANHz4cDw8PA86oYbFr1y5Gjx5NSUkJAGZmZixbtoyJEycaeGbNC1EUUaenUx4XR3ncASri49GVl9fb6zcJ0QCatGBcpTncY0Pl8uXLbNmyRT4ODAyke/fuBpxRw2LVqlVMmjRJ3razs7Nj3bp19O/f38Azax5UZ2VRHndAEooDcWjz8gEwbe2N3YgRWEeEQ1RUvVyryYhGQ2Du3Ln8+uuvGBsbY2RkxHfffceiRYt49dVX6dixI8XFxUydOpV9+/YB0KtXL7788kvs7e1JTEzk//7v/ygpKcHY2Jg333xT+YXWQKioqGDVqlVotVoAPDw8GD58uIFn1XD44osvmD59OlJuC3h6erJ582a6detm4Jk1XbTFxZQfPEh5XBwVcQdQnz8PgLGzM9bh4Vj3jMA6PBzTFvVvsqGIRj0RFxfHpk2bSEhIwNzcnPz8fNRqNYsXL5bPefrpp+ncuTPLli0D4J133mHKlCmsWrUKKysrli1bhq+vL5mZmQQFBTF06FAcHBwMdUsKSEv9tWvXUlxcDEgV3xMmTFCSEZD+bmbOnMlHH30kj/n7+7NlyxZat25twJk1PXQqFZUJCfJqQpWcDKKIkZUVViEhODwyEeuInpj7+ep9R0IRjXoiKysLFxcXzM3NAXBxcQEgMjKSTz/9FAcHB44cOcKKFSvk57z99tv4+Phw5syZWoVOXl5euLm5kZeXp4iGgdm9ezfp6eny8ZgxY3BycjLgjBoG1dXVPP300/z000/yWEREBBs3blRsVOoBUaNBlZwsxyUqjx5FVKvBxATLwG64vPgi1j0jsOzSBeEB/4BpcqKhT5W9uvy+GUOGDOH999/Hz8+PQYMGMXHiRPr16yc/npKSQmBgIMbGxvKYsbExgYGBJCcn0759e3k8Pj4etVpda0zhwXPmzBl27dolH/fq1Qt/f3/DTaiBUFpayvjx49m6das8NmrUKJYvX46VlZUBZ9Z4uRa8PkD5gZrgdVkZAOYBATg+9hjWPSOwCgrCyMBeXU1ONAyFjY0NR44cYc+ePezcuZOJEycyb948+XFRFG8qaH8fz8rK4oknnmDp0qVKTrsBKS4u5vfff5eP27Rpw4ABAww4o4ZBTk4ODz30EEeOHJHHnnnmGb755htMTJSvk7uhOjPzWvD64IFrwWtvb+yGD8c6IhyrsDBMGtjKVvlXrkeMjY2JjIwkMjKSLl26sHTpUvmxTp06cfToUXQ6nSwGOp2OY8eO0aFDBwBKSkp46KGHmDNnDuHh4Qa5BwXJL2n16tVUVFQA0g+CcePGNXsRT09PZ+jQoZw9e1Yee/fdd3n77beVzL46oLlyhYqDB2tWE3FUX7gIgLGLC9ZhYZJIhEdg1rJhd4hocqJxuy0kfZKWloaRkRG+vr4AJCYm0rp1a5KSpG63Pj4+dO/evZbb55w5c+jRowc+Pj6o1WrGjBnDk08+yYQJEwxyDwoSW7duvcGI0MbGxsCzMiyHDx9m+PDhcp2QkZERCxYs4JlnnjHwzBouuvJyKo4ckbecqmoq5I2srbEKDcXp8cexCg/H3Fc/wWuNTkNyQTIHsw5yMOvgnZ9QR5qcaBiKsrIypk6dSlFRkdwnYOHChYwfP14+5/vvv2fq1Kn4+PggiiIRERF8//33AKxcuZLdu3dTUFDAkiVLAFiyZAmBgYGGuJ1mS1JSEvHx8fLxoEGDmn0mUExMDOPHj6e8pkDMwsKCFStWMGrUKAPPrGEhqtVUHjsmi0Tl8eOg0SCYmmLZoweu06dhFRYmBa/1sJUniiJnis5wMPsgB7IOcDj7MGXVUlzE37H+YnGCoX6Z3ys366eRmpoqb/E0dZrTvT5o8vPzWbRokVygFhAQQHR0dLPeelm6dClTpkxBo9EA4OTkxMaNGxWDRkDUalGlnqTiQE3ldUICYmUlGBlh0amTVC8REY5ljx4Y6clqJqssiwNZBziYLa0m8iuluEgr21aEeYYR5hlGqEcoThZOCIJwRBTF4Pu9prLSUFDgRiNCJycnRo8e3WwFQxRF5s2bx3/+8x95zNvbmy1bttBcO2eKooj63HnKD0gFdeXx8ehq6nfMfNrjMG4c1uFhWIWGYmxnp5c5FKmKiM+Ol7acsg9yoUTqq+Rk4USYZxjhnuGEeYbRwkZ/cRFFNBSaPTczIpwwYUKzNSLUarVMnz6dr776Sh7r2rUrmzdvxsvLy4Aze/BUZ2VRfuCgtJo4cBBNTg4AJl6e2A4cKGc4mbq56eX6FdUVJOQmyHGJk4UnERGxNrUm2D2Yif4TCfMMw9dB/0V9V1FEQ6HZ83cjwoceeqjZGhGqVCqeeOIJVq9eLY9FRkaybt067O3tDTizB4OU4RQvrybUNT3NjR0dsQoPwzo8AuvwMEy9vfXyJa3WqjmWd4yDWQeJz47nRN4JNKIGUyNTAt0CeSHwBcI9w+nk0glTI8O4EiiiodCsuXz5MjExMfJxjx49mm3yQVFREaNHj2b37t3yWHR0NMuWLZOdDpoaN2Q4nTxZ257j0UewDg/H3M8PQQ8p1xqdhtSCVDkmcTT3KFXaKowEIzo7d2Zy58mEeoQS6BaIpYllvV//XlBEQ6HZctWIUKfTAZLRXlQ9OYE2NjIyMoiKipJTxAGmTZvG559/3qTqU26b4dS9O64vT8UqPBzLzp31Ys+hE3WkF6VLK4mseA7nXMtw8nX0ZYLfBMI8wwhyD8LWrGG2xFVEQ6FZotPp+P3332UjQgsLCyZMmNAsq5pTUlIYOnSoXJsC8PHHHzNjxoxGnwhwpwwn56eewio8DKsePTCyrP9f8qIocq7kHIeyDhGfLYlEoaoQAG9bb4a1HUaYZxgh7iE4WzYOz67m9wnRI2vXrmXs2LGkpqYSEBBAZmYmL7/8cq39YYWGwe7duzlz5ox8PGbMGBwdHQ04I8Owd+9eRo0axZUrVwApCeDHH39k0qRJBp7ZvXHHDKexY6XgdUgIxnqI0YiiSEZphpThlH2Qw9mHyauUEizcrdzp5dWLUM9QwjzC8LTxrPfrPwgU0ahHli9fTu/evfntt99499138fLyuqlgaDSaZvmLtqGQnp7OX3/9JR/37t27lstwc2Ht2rU89thjqFQqQLJLWbNmDUOGDDHwzO6O6uxsaRVhoAynzLJM4rPjOZQtrSayy7MBcLZwJtQzlFAP6U8r21aNfuUGimjUG2VlZezbt4+dO3cyatQo3n33Xc6fP8+IESNISkpiyZIl/PHHH6hUKsrLy4mNjTX0lJslNzMibI7d5b799lteeuklOZ7j5ubGn3/+SVBQkIFndme0RUWUX5/hdLUB0fUZThHhmLbSz5d0bkXuNZHIiiejTNrWczB3IMQjhKc7P02oRyht7ds2CZH4O01ONN577z29vfY777xzy8fWrVvHsGHD8PPzw8nJiYSEhBv6LsTFxXH8+HGlH4OB0Gq1rFq1isrKSqB5GhGKosjbb7/NnDlz5DEfHx9iYmIarBW/rrKSioQEKmp6S6hSUmpnOE2ciHWE/jKc8iryOJR9iEM5hzicfZjzJecBsDWzJdg9mMc7PE6IRwi+jr4YCQ30vVSaU28v1eREw1AsX76c6dOnA/DII4+wfPlyXnzxxVrnDB48WBEMA7JlyxYuX74MSIZ7EyZMaFZGhBqNhueee44ffvhBHgsJCWHTpk246Wnr5l4QNRpUSUmUHzhA+f44qQFRdTWYmmLVrRsuL72IdURPLLvoJ8MpvzKfw9mH5dXEVZGwMbWhh3sPxvmOI8QzhADHAIyNjG//YoZCUwUXD8CZHZAeCzkn7vycOqKIRj1QUFBAbGwsSUlJCIKAVqtFEAReeOGFWudZG7h5SnMmKSmJQ4cOyceDBg3C29vbgDN6sJSXlzNx4kT++OMPeSwqKoqVK1caXDhv24CoQwccn3hCiksEBWGkhyZP+ZX5HM45LAvFueJzAFibWhPkHiSJhEcI/k7+mBg10K9MUYTCs5C+QxKKc3uguhyMTME7HAa+A++9Vi+XaqB/A/fO7baQ9MXq1at58skn+e677+Sxfv361UphVDAceXl5bNiwQT7u2LFjs+pXkpeXx4gRI2q5906ePJmFCxcarNd5dXY25fvjpAZEB+JubEDUM0JqQKSHjLaCygIO5xzmULa03XSmWMqisza1podbD8b4jCHEI4QAp4CGKxIAVaVwbrckFOnboUiqXsexLQQ+Bj4DoU1vML9a76GIRoNh+fLlvPHGG7XGxo0bxwcffGCgGSlc5aoRYXV1NQDOzs6MGjWqSQYob8a5c+cYOnQop0+flsfefPNNZs+e/UD/DrSlpVTEx0tCsX8/6nPSr3ljZ2fZDVZfDYiubjddFYqzxVITKSsTK3q492CUzyhC3EPo4NyhYYuETgfZxyWBOBMLlw6CTgOm1tCuH/ScKgmFUzu9TkOxRm9kNKd7vV9EUeT333+Xq5xNTEyYMmUK7u7uBp7Zg+Ho0aNERUWRU5OCKggCX3755Q2xNn1wrfI6jvJ9+6lMSgKtFsHSEquQYKx79sQ6oifmfvVvtJdbkVtLJK7GJKxNrenu1p0QjxCC3YPp4NzBYP5NdaY0RxKIM7FwdieUSzUfeHQBn0HQfiC0CgMTszu+lGKNrqBwB+Lj42vZYowYMaLZCMb27dsZO3YspaWlAJibm/PLL78wbtw4vVxPFEWqTp2mPG4/5fv3U3HosFx5bdmlCy7PPYt1RASW3bohmN35C+5uyCnPkTObDucclu3Crw9cB3sEN/ztJqgJYMfVxCZ2XgtgW7lA+wHSSqJdf7A13Pu4gf8NKijcGxkZGWzdulU+DgoKolu3bgac0YPj119/ZfLkyfKWnIODAxs2bKBPnz71ep3qnBzK90kiUX7gANp8KS5h1q6dVHndM0LqLWFbvx5KWWVZUuC6Jnh9sVTqtW1rakuQexAT/CZIItGQs5uuIoqQf1oKXp+JhfN7obqidgDbZyC4d4EGkhquiIZCk6O8vPwGI8Jhw4YZeFYPhs8++4wZM2bIxy1btiQmJoZOnTrd92vrKiqoOHy4Rij2UXU6HQBjFxesIyJq/oRj6ll/9hiiKJJRliGvIo7kHOFymZQ2bWdmR5B7EBP9JxLiEYKfo1/DFwmAyitw9q9r207Fl6RxZx/oPknacmrTG8wbZjq4IhoKTYqrRoQlJSWAZEQYHR3d5G1bdDodM2bMYP78+fJYp06d2Lx5M61atbqn1xR1OqpOnqRs3z4pLnHkCGJ1NYK5OVbBwdg/PAbr3r2korp6ikuIosj5kvPyKuJIzhFyKqSYjKO5I0HuQTzR8QmC3YMbdjHd9Wg1kJlQs+UUC5cPg6gDczto2xf6vCptPTm2MfRM60TT/iQpNDv++usvzp49Kx+PHTsWBwcHA85I/1RVVTF58mR+++03eaxPnz6sX7/+rk0Y5S2nffsoj4tDWyg5spr7+0v1Er16SvUS9dTVUBRFzhSdkbebjuQckftcO1s4E+wRTLC79Ke9Q/vGk/V25cK1Laezu6GqGAQj8OoOfWZIW04tgsG48X0FN74ZKyjcgvT09FoNhPr06YOvr68BZ6R/SkpKGDNmTC0vs7Fjx/LLL7/UqV2tTqWi4tAhyvfuu3HLqXcvbHr1wioiot7M/nSijtNXTtdaSVypkhx23azcCPUIlYWijV2bxiMSVaVSQd3VLafCGgdlu5bQabS0kmjbD6wavyOEIhr1yNy5c/n1118xNjbGyMiI7777jkWLFvHqq6/SsWNHQ0+vSVNUVFTLiLBdu3ZERkYabkIPgKysLKKiojh27Jg89sILL/DFF19gbHzzvX1RFFGfP0/5nr2U7dlDRXw8YlUVgpmZXractDotJ6+clGMSCTkJlKilrUMvay/6tOwjryRa2rZsPCKh00JWYo1I7LyuZsIK2vSB0Gel1YSzDzSWe6ojehUNQRCGAf8DjIHFoijO+9vj9sDPgHfNXD4VRfFHfc5JX8TFxbFp0yYSEhIwNzcnPz8ftVrN4sWLDT21Jo9Go6llRGhra8vYsWObtBFhWloaQ4cO5UJND2uQfrTMnDnzhi9eXUUF5fHxlO/eQ9mePVRfkgKvZm3a4BAdjU3fPliFhNTLllO1rprUglR5JXE096jcmc7b1ptBrQcR7B5MkHsQXjZe9329B0pxhiQQZ2Lh7C6olLbu8OwmFda1HwitQsGkabbGvYreREMQBGPga2AwkAEcEgRhgyiKKded9iKQIoriSEEQXIE0QRB+EUVRra956YusrCxcXFzkXsouLi4AREZG8umnnxIcfN81NQq3YMuWLWRmZgLXjAibss/XgQMHGDFiBAUFBQAYGxuzePFiJk+eDNSsJs6epWzPHsp376Hi8GFEtRrB0hLrsDCcJv8Dm759MbvHAPn1VGurSSpIklcSR3OPUqmRxLutfVui2kbJIuFu3chqZNTlcH7ftS2n/DRp3MYD/IbV1ExEgrWLIWf5wNHnSiMUSBdF8SyAIAi/AaOB60VDBGwF6aeRDVAIaO73wrt27arVZOd29OjRg5EjR9Ya27hxIwkJCfJxv3797rjVMWTIEN5//338/PwYNGgQEydOpF+/fnc9d4W74/jx41zvEDB48OB7zhZqDGzatIno6Gh5VWVlZcXq1asZ2rcvpbGxlO3eTfnuPVTXiKhZ+/Y4Pvoo1n37YBUcjJH5/f0KVmvVnMg/weHswxzKOcSx3GOotFITJ19HX0a3H02whyQSLpaN7Mv0qk3HVZG4eAB01WBiAa17QY8npdiEW4cmt+V0N+hTNFoAl647zgDC/nbOV8AGIBOwBSaKoqj7+wsJgvAs8CzQYJ1JbWxsOHLkCHv27GHnzp1MnDiRefPm3fmJCvdMbm4umzZtko87depEWNjf32JNh8WLF/Pcc8/J9ScuTk4sf+UV/H5fy6l//VtKh7Wywjo8HOdnn8G6d5/79nKq0lZxIu+EXHF9LO8YVdoqBAT8nfwZ7zeeYPdgerj3wNGiEbbLLcm8bstpJ1RIqzfcu0D4/0ki4R0BpvWTLdYU0Kdo3EyK/250NRRIBAYA7YFtgiDsEUWxpNaTRHEhsBAk7yk9zLVeMDY2JjIyksjISLp06cLSpUsNPaUmS1VV1Q1GhCNHjmw8gdS7QBRFZs+eXcvBuZWVNd/ZO9Di519Qe3vj+Nhj2PSPxLJHD4zuw6ajSlvF8bzjtVYSap0aAYEApwAm+E0gxCOEIPcg7M3rv8e23lFXwIX911YTeanSuLUb+AyWRKJdpEFtOho6+hSNDOD6fYKWSCuK63kKmCdKronpgiCcAwKAeO6Dq1/c98rIkSNv2LK6E2lpaRgZGckpnomJibRu3bqW95FC/SCKIhs2bJD39E1NTYmOjpbjSU2JapWK5x55lB/Xr5PHOpqb8+OggbQb/hC2A/pj5uNzz2JZpa3iWO4xDuUc4lD2IU7knaglEo8EPCKvJBqlSOh0kJN03ZZTHGjVYGwOrSMg8FEpgO3eqVlvOd0N+hSNQ4CvIAhtgcvAI8BjfzvnIjAQ2CMIgjvgD5ylEVJWVsbUqVMpKirCxMQEHx8fFi5cyPjx4w09tSbHwYMHSUm5FhobMWJEg+o8d79oy8op37uH3JgtPL90KTuKrsiPRXbtypq1a3Fqd2/219W6apLzk4nPjic+K56juUdR69QYCUYEOAXwaMCjBHtIImFnZldft/RgKc2RtpqupsOW50rjbh2lVNj2/cG7J5jVf0On5oDeREMURY0gCC8BW5BSbn8QRTFZEITnax5fAMwGlgiCcAJpO+vfoijm62tO+iQoKIj9+/ffML5r164HP5kmzKVLl9i2bZt8HBwcTNeuXQ04o/pBV1FB8YaNlG7fTsXBg1xRqXgxO4ujNR3sACZNmsT333+P2V1sP2l1WtKupBGfFc/B7IMk5CRQoakAwN/Rn2j/aMI8wxq3SFRXSiuIqyKRU7O6t3KRBKL9AMkZ1q7+PLGaM3qt0xBF8U/gz7+NLbju/zOBIfqcg0LT4e9GhF5eXgwdOtTAs7o/dJWVXPltBQWLF6MtKMC0tTdlDz3EUyt+I+06wfjXv/7Fhx9+eMfak6u2HPHZ8XKP66vFdG3s2jCy/UhCPUIJ8QhpnIFrkJxhc1Ov2XRc2A8aFRibSc6wg96VhKIBOcM2JZSKcIVGgU6nY82aNXJ/CEtLSyZMmNBojQh1KhVFK1aQv2gx2vx8rHtG4PLSS5w2Nmb08OFy3YkgCMyfP59p06bd8rUyyzLZn7mf+CxJKApUUqynhU0LBnoPJNQzlFCPUNysGvEWXnnBdVtOsVCaJY27BkDwPyWRaN0TzJpufU5DoXF+4hSaHbt27eJcTYtQgDFjxjRKI0JdVRVFK1dRsHAhmrw8rMLDcf3ff7EKCmLXrl2MHj1adug1MzPjp59+Ijo6utZrVFRXcDjnMPsz97Pv8j65M52rpSvhXuGEeYQR4hFCS9uWD/r26g+NGjLiJYFI3wFZxwARLB2lrab2A6StJ/tGfI8PAK1OJP5cIVuSs+vtNZuMaIii2CTTLa+nsbXmrS9OnTrFnj175OO+ffs2OiNCnVpN0erVFHy3EE1ODlbBwXh99inWoaEArFq1ikmTJqFWS2YIdnZ2rFu3jv79+6MTdaQVprE/cz/7M/eTkJuARqfBwtiCII8gov2j6enVk3b27RrvZ0AUofDsNZE4vwfUZSAYS9Yc/d+UhMIrEBpDzwwDotboiDtbQExSFluTcygoV2NuUn/bdE1CNCwsLCgoKMDZ2bnxfmjugCiKFBQU1Mm5tClRVFTE2rVr5eN27do1qkp7Ua2m6Pffyf9uIZqsLCyDgvD6aB5WYWHye/WLL75g+vTp8o8CT09Plq9bTolzCW/seYO4zDgKVZLPkZ+jH5M6TKKnV096uPfA3LgRpxmriuHc7mt9JopqfLQc20DXiTXOsH3BopEG6B8gqmotu0/lEZOUzfbUHEpUGqzNjBnQwZ2ozh5E+rtiPbd+rtUkRKNly5ZkZGSQl5dn6KnoFQsLC1q2bD7LcY1Gw8qVK1GpJJsKOzu7RmNEKFZXU7R2LfkLFqDJzMKye3e85s7BKiJCFgudTsfMmTP5+OOP5ee5tnal21vdeCn1JUBqPBThFUGvFr2I8IzA1crVIPdTL+h0kH0M0rdDeqzkDCtqwcxWEoeeUyU/J6d7SydubpRVadh5MpeYpGx2puVSodZib2nKkE4eRHX2oJePCxam9b8qaxKiYWpqStu2bQ09DYV6JiYmhqwsKeDZWIwIxepqijdsIP+bb6m+fBmLbl3xfO99rHv3qrUKVqvV/OOpf/Dbr9caJ1n6WOL+ijuerTwZ5zWOnl49CXAKaBzd6W5FWV7NltN26b8VNRn1noHQe/o1Z1hjU8POs5FQXFHN9tQcNidls/t0HmqNDhcbc8Z0b8Gwzh6Et3PG1Fi/75cmIRoKTY9jx45x5MgR+XjIkCENepUliiKl27aR++lnVF+8iEWXLni88zbWffrUEouCygJiTsYw85mZXDpyzZrNt7cv8xfPJ7JdJNamDVsYb4u2GjIO1awmdkg9J6CmZmIA+AyS/mvTiFdMD5i80iq2peSwOSmLuDMFaHQiXvYWPB7mTVRnT4JaO2Js9OC25RXRUGhw5OTk1DIi7Ny5M6E1AeOGSNW5c+TM/YDyvXsx9/en5bffYBMZKYvFpZJLxF6KJfZiLPGn4zk//zyq8yr5+VOemcK333zbaNOHKbpUs5LYAWf/gqqSawHsAW9JQuHRTamZuAsyrlSwJTmHLUnZHLpQiChCG6kXCrsAACAASURBVGcrpvRpR1RnD7q2tDdY/LaRvksVmioqlYqVK1ei0UgO+S4uLg3WiFBXWUn+gu8o/OEHBHNz3N98E8dHHwFjY1ILU4m9GMuOiztIL5JaqLZQtSD/k3xUGdcE49133+Xtt99ukPd3SzRqqQI7fRuc3gZ5J6Vxu5bQaYwkEu36gUUj9KoyIGfyyohJyiYmKZsTl4sBCPCwZdpAX4Z19sDf3bZBvE8U0VBoMFw1IiwslDKFrhoR3o1txoNAFEVKt28n58MP0WRmYT96NG6vzyDXQs1HCZ+y4+IOssuzMRKM6OHWg3+F/AvHXEf+Gf1P8vOkPX0jIyMWLFjAM888Y+C7qSMlmZJAnN4qda1Tl0kV2K17QfcnwHcwuPgppn93gSiKJGeWsCVZEorTuZIDQHdvB2ZGBTC0kwdtXBreVqUiGgoNhgMHDpCamiofjxw5ElfXhrX3rb5wgey5cynfvQdzPz9a/Pwxpj26sTR5Gd8d/w6tTkuvFr14MfBF+rXsh6OFIzExMYwbN46KCsnzycLCghUrVjBq1CgD381t0Gqk2MTprZJY5JyQxu1aQpcJ4DtEyngytzHsPBsZOp3I0UtX2Hwim5jkbDKuVGIkQFhbZyaFt2ZIJ3c87S0NPc3booiGQoPg4sWLtYwIQ0JC6NKliwFnVBtdZSUFixZRsGgxgpkZ7v+ZieNjj5FQcIw5G6NJL0pnQKsBvBH6Bp4214zxli5dypQpU+TtNicnJzZu3EjPnj0NdSu3pixPik2c3irFJ1TFYGQiNSEa9J4kFM28a929oNboOHC2gC3J2WxNySGvtApTY4HePi5MHeDDoA7uONs0nnobRTQUDE5ZWRmrVq2Si9tatGjBkCENw8dSFEXKYmPJmfsB1ZmZ2I0aiduMGZTZmfL2wfdYl74OL2svvhzwJZGtIms9b968efznP/+Rx7y9vdmyZQsBAQEGuJObIIqSPUfaZji9BTKPSuM27tBhpCQS7SKV2MQ9UKHWsPtUHluSc9hRU2xnaWpMpL8rwzp70D/ADTuLxplmrIiGgkG5akRYVuPo2pCMCNUXL0pbUX/txtzXB+9lS7EMCWZd+jo+j/2ccnU5T3d+mme7PouV6bXeDFqtlmnTpvH111/LY127dmXz5s14eXkZ4lauUa2SLDrS/oS0GCjNBMEIWgRLmU6+QxR32HukuKKaHSdziKmpoVBV63CwkorthnbyoI+vfortHjSG/2QqNGt27tzJ+fPn5eOxY8dib2/YX7Y6lYqChYsoWLwYwcQEtzf+jdPjj5Nedp7ZMZM5mnuUHm49mBU+Cx9Hn1rPValUTJo0iTVr1shjkZGRrFu3znD3VZ4Pp7ZIQnFmJ1SXg6k1+AwA/1mSUFi7GGZujZzcEhVbUqTU2ANnpRoKDzsLooNbMayTB6FtnTDRc7Hdg0YRDQWDkZaWxt69e+Xjfv364ePjc5tn6Be5QO/jT6jOyMBuxAjcXn+daicb5h//kp+Sf8LGzIbZvWYzuv3oG9Ifi4qKGD16NLt375bHoqOjWbZs2YNtRSuKkH+qZjWxGS7FAyLYekG3R8B/OLTpDabNy8esvrhQUE5MUjZbkrNJuFgEQFsXa6b0acewzh50bWGP0QMstnvQKKKhYBCuXLnCunXX+l63b9+evn37Gmw+FYcPk/vJp1QeOyZtRS1ZgnV4GDsv7uTD9R+SVZ7FON9xTO8xHQeLGy3ZMzIyiIqKqtUTfvr06Xz22WcPxitLq5FqJ07FSGJRWNM12bMbRL4B/lHg0VUJYt8DV1NjtyZnsyU5h7QcqadL5xZ2vDbYj6GdPfB1s2kQNRQPAkU0FB44Go2GVatWNQgjQtWpU+R9Pp+yXbswcXfHc+4c7EePJluVx8zYl9l5aSc+Dj4si1pGd7fuN32N5ORkhg0bRkZGhjz2ySef8Nprr+n3i0SjhnN/QfI6SPsDKq9ItRNt+0HES+A3DOxb6O/6TRitTuTQeakPxdbkHC4XSamxwW2cmDWiI0M6utPKqXn2GFdEQ+GBs3nz5huMCK2sHuwHsDori7wvv6J43TqMrK1xfe1VnCZNQmtuwpKUn/j22LcAvBr0KpM6TsLU6OaZLnv37mXkyJEUFUnbFCYmJvz4449MmjRJPxPXVEnGfynr4eSfUFUM5nbSSiJghNSYyNxWP9du4qiqtew9nc+W5Gx2nMylsFyNmYkRfXxcmDbQl4Ed3BpVaqy+UERD4YGSmJhIQkKCfDx06NAHakSoLS6mYNEiCn/6GXQ6nP7xD5yffQYTR0eKVEW8svUFDuccpn+r/swMnVmr5uLvrF27lkcffZSqqioAbGxsWLNmTf2nC1dXSuZ/KeulGIW6VEqD7TACOo6W0mJNlC+ze6G4spqdJ3PZmpLNrrQ8KtRabM1NGNDBjaGdPOjr54qNufI1eT3K34bCAyM7O5s//vhDPu7cuTMhISEP5Nq6qiqu/PwL+QsXoispwX7UKFxfnoppC2n75mzxWV7a8RI55Tl80PsDRrYfedvX+/bbb3nppZfQ6XQAuLm5sXnzZnr06FE/E1ZXSN5OKeulzCd1mdTqtNPD0PFhqRrbpGHZqzQWckpUbEvJYUtytuwa62or2YsP6eRBRDtnzOqx011TQxENhQeCSqVi1apVD9yIUNRqKd6wkbwvvkCTlYV13z64vfoqFtcV2MVlxvHartcwNTbl+6HfE+gWeOvXE0VmzZrF3LnX2qD5+PiwZcsW2rW7z+ZBVWVSNXbKeum/1RWSpXiX8ZJQtOmt9J24R9Jzy9iaIsUnEi9JW4ltnK14undbhnTyoHsrhyad8VSfKKKhoHdEUWT9+vUP1IhQFEXKd+8m97PPqTp1CosuXfD68EOsw8NqnbcybSUfHPyAtvZt+Xrg13jZ3Lr4rrq6mueff54ffvhBHgsJCWHTpk24ubnd20R1WskA8NhySN0EmkqwdoPAx6StJ++eYKx8TO8WyeOpiK0p2WxLzuFsfjkA3VraM2OIH0M6Na+Mp/pEeTcq6J24uDhOnjwpH48aNUqvRoSVycnkzvuIikOHMG3tTYv/zsd26NBaXxAanYZPD3/KL6m/0LdlXz7u+/Ftmx+Vl5cTHR3Nn3/+KY9FRUWxcuVKbGzuwbQvLw0Sf4XjK6WqbAsHSSg6jwPvcDBq/JXDD5oqjZb9ZwrYmpzD9lTJ48nESCCivTNP9WrDoI4N3wywMaCIhoJeuXDhAtu3b5ePQ0ND6dy5s16upVOpyP/qKwp++BFjR0fcZ72FY3Q0gmntLZ0ydRmv736dvZf38kTHJ3gt6DWMb/MlnZeXx4gRI4iPj5fHJk+ezMKFCzE1vYvtoopCSFojiUVmgtSoyHcwDPtQyn5Sgtl3TXFlNbvSctmaksOuk7mUq7VYmxkTGeDGkI7uRPq7YW+pbOnVJ4poKOiNsrIyVq9eLRsRtmzZUm9GhOXx8WTNmkX1hYs4TJiA2+szMLazu+G8jNIMpsZO5Xzxed6OeJsJfhNu+7rnzp1j6NChnD59Wh578803mT17dt22NrTVUnwi8VcpoK2rlrydhn4gWYzb3OO2VjMms6iSHak5bE3JkQPZLjbmjApswZBO7vRs74y5ibJS0xeKaCjoBZ1Ox+rVq2UjQisrK8aPH4+xcf1+mLWlpeR++hlFK1Zg2qoV3kt+xDo8/KbnHs09yrTYaWhEDd8O/pZwz5ufJ59/9ChRUVHk5OQAIAgCX375JS+++OLtJ3XVPfbYcjixCioKwNoVQp+FwEfBo+FYvjcGRFEkNauUbSk5bEvNJulyCQDtXCXrjiGd3AlsqQSyHxSKaCjohdjYWC5cuCAf68OIsDR2J9nvvYcmLw+np57C9eWpGFnefM9645mNvLP/HbxsJBvztvZtb/va27ZtY+zYsbLomZub88svvzBu3LhbP6miEBJ/gcTlkJssVWf7R0G3x8BnoJL5dBdUa3UcOlfI1hQpPpFxpRJBgB7ejrwRFcDgju60d1UaQBkCRTQU6p20tDT27dsnH0dGRtK+fft6e31NYSE5cz+g5I8/MPfzo+VXX2J5i4ZNOlHHV0e/YtGJRYR6hPJ55OfYm99evH755RcmT54spwc7ODiwYcMG+vTpc/Mn5J6EA9/A8RWgUUk24w99Bp3GgpXTfd1rc6KsSsNfaXlsS8lmZ1oexZXVmJsY0cdXalY0IMAdV1sl7mNoFNFQqFcKCwtZu3atfOzj41NvRoSiKFKyaRM5cz9AW16Oy8tTcZkyBeEWqbuVmkre3Psm2y5sY5zvON4MexPTO/za/+yzz5gxY4Z83LJlS2JiYujUqVPtE3U6qbvdgW8kWw8TC+g6EcKeB/eO932vzYWrhXbbauITaq0ORytTBnd0Z3BHd/r4umBlpnxNNSSUfw2FeqO6uppVq1bJthr29vaMGTOmXnLhq7OyyHr3Xcr/2o1lt254zpmNua/vLc/Prcjl5diXSSlIYUbwDJ7s+ORt56HT6ZgxYwbz58+Xxzp27EhMTAytWrW6dqK6HI79BgcXSPbjNh4wYBYEPQXWzvd9n00dURQ5mV3K9pptp2MZxQC0drbiHz1bM7ijB0GtHTFW4hMNFkU0FOqNzZs3k52dDYCxsXG9GBGKOh1FK1aQ++lniDqd1Jv78ccRbhNQTy1I5aXYlyhVl/LFgC9qtWG9GVVVVUyePJnffvtNHuvTpw/r16/H0dFRGijOgPhFcGQJqIrAqzuMXSwV4Cl2HrdFrdERf66Q7anSiuJykRSfCGzlwOtD/Rnc0V0ptGtEKKKhUC8cPXqUo0ePysdDhw6lRYv7s+WuOnuOrLdnUXn4CNY9I/B4/33M7mBuuOX8Ft7a+xYOFg78FPUT/k7+tz2/pKSEMWPGEBsbK4+NHTuWn3/+GUtLS7h0SNqCSlkPiFLv7PAXoFWY0pviNhRXVLPrVC7bUnL4Ky2P0ioNFqZG9PZx5eWBPvQPcMPNVmkC1RhRREPhvsnOzq5VKd2lSxeCg4Pv+fVErZbCH38k74svESws8Jw7F/uxt9/m0ok6vj32LQuOLSDQNZD5/efjYnn7FqZZWVlERUVx7NgxeeyFF17gi/mfYXzqD4j7Bi4fBnN7iHhBSpl18L7n+2rqXCyoYHuqtO0Uf66wpn7CjOFdPBnU0Z3ePi5Ymin1E40dvYqGIAjDgP8BxsBiURTn3eScSOC/gCmQL4piP33OSaF+UalUrFy5Us40cnV1ZcSIEfe81aC+cIHMN2ZSefQotoMH4T5rFqZ38HWqqK7gzb1vsv3idh72eZhZ4bMwM779llFaWhpDhw6tlRY8d/b7zBzkivBVDyi5DE7tIOoTqbZC6VFxAzqdSGJGETtSc9iekit3tPNzt+HZvu0Y1FGpn2iK6E00BEEwBr4GBgMZwCFBEDaIophy3TkOwDfAMFEULwqCoJTHNiJEUWTdunVcuXIFADMzs3s2IhRFkaIVK8j56GMEExO8PvkYuzqIz+Wyy7wc+zLpRem8Hvw6T3R84o7PiYuLY8SIEbKBorGxMYvffobJNish5iK07gUPfQ6+Q8AA3QQbMuVVGvaczmdHag4703LJL1NjbCQQWtPRblAHN1o739rDS6Hxo8+VRiiQLoriWQBBEH4DRgMp153zGPC7KIoXAURRzNXjfBTqmf3795OWliYfjxo1CheX228J3YzqnByy3nyL8r17se7ZE88P5mLq4XHH5x3JOcIrO19Bo9PwzcBv6NWi1x2fs3HjRiZOnEhlZSUAVpbmrP5nO6J0v4JFV5g0H9oPVOIV13HVtmN7ai5xZwtQa3TYWpgQ6e/GoA5uRPq5YW+lFC42F/QpGi2AS9cdZwBhfzvHDzAVBGEXYAv8TxTFZXqck0I9cf78eXbs2CEfh4WF3VjLcAdEUaTkjz/Jfv99RLUa97dn4fjoo3Xa2lpzag1zDs6hpU1LvhjwxR0rvAEWL17Mc889JzdOcrEx5Y9HTAj1NYMB30vFeMrKAp1O5MTlYnak5rAtNZfULMm2o42zFU+Et2ZgBzdC2jhhaqz8XTVH9CkaN/vkize5fhAwELAE4gRBOCCK4qlaLyQIzwLPAnh7K4FIQ1NaWlrLiLBVq1YMHjz4rl5Dc+UK2e+/T+nmGCy7dcPro3mYtWlz5+fpNHxy6BN+PfkrPb168nHfj+9Y4S2KIrNnz+add96Rx9o6CGx51gvf8bOgx5PN3uKjUq1lb7q07bTjZC55pVUYCRDU2pGZUQEM7OBOe1drJS1WQa+ikQFcVxVFSyDzJufki6JYDpQLgrAb6AbUEg1RFBcCCwGCg4P/LjwKDxCdTseaNWsoL5ea2tyLEWHprl1kzZqFtqgY11dewfnpfyKY3PmtWFxVzIy/ZnAg6wBPdHyCV4NexcTo9s/TaDS8+OKLLFy4UB7r4WXKn9/8B/eo18Gs+e6/ZxVXEnsylx2puexLz6dKo8PG3IR+fq4M7OBGf383HK2VGhSF2uhTNA4BvoIgtAUuA48gxTCuZz3wlSAIJoAZ0vbVfBQaLDt27JAzjgRBYNy4cdjdxIL8ZmjLysn96COKVq3C3NcX74ULsejQoU7PPVt0lqmxU8ksz+T9nu8zxnfMHZ9TUVHBo9Hj2PBHjDw2uEdb1vy5E1v31nW6blNCpxM5llEkC0VKzbZTKydLHg31ZlAHd0LbOin9sZsYJSUllJSU1Nvr6U00RFHUCILwErAFKeX2B1EUkwVBeL7m8QWiKKYKghADHAd0SGm5Sfqak8L9cfLkSfbv3y8fR0ZG1rkvdsWhQ2TO/A/Vly/j/MwUXKZOxaiOWVa7M3bz793/xszYjB+G/kB3t+53fE5h1gVGDu7H/uRrKbWTosfy/U/L9dpmtqFRVqVh7+k8dqTmytlORgIEt3bijagABga44aNUYzc5NBoNJ0+eJDExkbNnz95Tgsqt0GudhiiKfwJ//m1swd+OPwE+0ec8FO6fwsJC1q1bJx/7+vre2vX1OnRVVeT9938ULlmCacuWtP75J6yCgup0TVEUWZq8lM+PfI6/kz9f9P8CTxvPO1xQx4U/PmXY029xMq9aHv7Xv/7Fhx9+iFEzCHRfKqyQYxMHzhZQrRWxszChX022Uz8/Vxysmo9wNie0Wi0xMTEkJSWhUqnk8by8vHq7hlIRrnBHqqurWblypWxE6ODgUCcjwsrkZDL//W/U6WdweGQi7q+/jpF13WIIVdoq3o97nw1nNjCk9RBm95qNlekdfKzy0jj+zT+J+iyOzFIp9CUIAvPnz2fatGl1um5jRKPVkXCxiB0nc4hNzeV0rtQDpL2rNU/1asuAADeCWjsq2U7NAGNjY7KysmoJBkDbtnfOLqwrimgo3JE///xT7l531YjQ8hbNjkAyGSxcuozczz/HxMGBVgu/w+Yu7NEzyzKZ8dcMTuSf4IXAF3i+6/O3FyiNGvbOZ9eyDxj9ayklVZJgmJmZ8dNPPxEdHV3nazcWrpSr+etUHrEnc/nrlNR7wsRIIKydE4+EejMwwI02Ls03yN/U0Wg0nDol5Qt17Fjbij8wMJDLly/j6OhIt27d6NatGw4ODvzjH/+ol2sroqFwWxISEkhMTJSPo6Ki8PLyuuX5msJCMt94g/Lde7AZNBCvOXMwdnCo8/X2ZOxh5t6ZaHVa/hv5Xwa2Hnj7J1yKhw0vs3LXcZ5YV4VaIwmGnZ0d69ato3///nW+dkNGFEXSckql2MTJXBIuXkEngouNGYM7ujMgwI3evi7YWTTv1OGmjCiKZGVlkZiYyIkTJ1CpVDg7O9OhQ4daP6o6d+6Mq6sr3t7eeolVKaKhcEuysrJqGRF269aNHj163PL88gMHyHz9X2iLi3Gf9RaOjz1W5zetVqfl68SvWXRiEX6Ofnwe+Tmt7W6T4aQqgR3vw6HFfHHMgunrVXLdiKenJzExMXTt2rVuN9pAqVRriTubLwtFZrG05dC5hR0vDfBlQIAbXVvYK95OTZyysjJOnDhBYmIiubm1TTMKCgq4fPkyLa9zf7awsKB1a/1lByqioXBTVCoVq1atQqvVAuDm5sZDDz10UxEQNRryvv6aggXfYdamDa0WLcQiIKDO18qvzOeN3W9wMPsgY33HMjN0JhYmt7HNTtsMf7yGrvgyM1MC+HhdvPxQQEAAMTExev3Q6JPLRVLtRGxqDvvPFFCl0WFlZkxvHxemDfIl0t8NdzvFUrypo9VqOXXqFImJiZw+fVr+QXQ99vb2BAYGYm9/++LW+kYRDYUbuJURoanpjVsf1ZmZXJ7xOpUJCdiPHYvHW29idBeNlw5nH+Zfu/9FqbqU2b1m87DPw7c+uTQHYv4NyWupdg7g6aRe/PT7ZvnhiIgINm7ciLNz4+mgdzWIHXtSWk1cdYpt7WzFo6HeDOzgRmhbJ8xNFEvx5sT3339PVlbWDeMmJiZ07NiRwMBA2rRpY5BUaUU0FG5g3759tYwIH3744Zt+EZds20bWW7OguhqvTz7GfuTIOl9DJ+pYkryELxK+oKVtSxYMXoCfo9/NTxZFOPozbH0TqispDZvB+Pm72brtmmCMGjWK5cuX33enwAdBYbmaXWm57EzLY/d1QeyQNk68ObwDAzq40c5FsexoLoiieMO/tY+PTy3R8Pb2plu3bnTq1Alzc/MHPcVaKKKhUIvz58/X6mIXHh5Oh79Vbeuqqsj96GOu/PorFp060eLzzzC7i+2g4qpi3tr7FrsydjGk9RDe6/keNmY2Nz+54AxsnAbn90DrXuSEzeKhJ6dy5MgR+ZRnn32Wr7/+GpM6WJEYAlEUSc4sYefJXGLTckm8VIQogouNOUNqgti9lCB2s0Kr1ZKenk5iYiIajYbHH3+81uOBgYEcP36crl27EhgYiJOTk4FmeiMN81OmYBD+bkTo7e3NoEGDap1TdfYsl195laq0NJwmT8bt1VcQ7qLCOrkgmdd2vUZORQ5vhL7BYwG3CJZrq2H/l/DXR2BsDiP+S7pdL4ZGRXH27Fn5tPfee49Zs2Y1uF/lUiV2PjtPSpXYuaVVCAJ0benA9IF+9A9wpbOXEsRubuTk5MjZT1f92wCKiopwuC7L0MnJiWnTpjW49zUooqFQg1arZfXq1fIb2draupYRoSiKFP++luw5czCysKDlgm+xjYys8+uLosjKtJV8dOgjnC2dWTpsKV1db5HdlJ8Oa/4JWcekntxRn3AoLYOHeveWK1uNjIxYsGABzzzzzH3dd31yLr9cjk0cPCdVYtuam9DXz5X+AW5E+rviYmPYrQWFB09FRQVJSUkkJibeNE4BcObMGYL+5pTQEAUDFNFQqGHHjh1cvHgRuGZEaGsrtTjVlpWR/e57lGzahFVYGF4ff4ype92bLFZUV/Be3Hv8ee5PerfozYe9P8TB4ia1G6IIib/Cn6+DiRlEL4OOo4mJiWH8+PGyoFlaWrJixQpG3kUMRR9UabTEnyuUheJ8QQUAvm42/LNXW/orldjNmjNnzpCQkEBaWpqchXg9tra28vZTfXpD6RtFNBRITU0lLi5OPu7fv79sO1B5IonLr71GdUYGrtNexvnZZxHuwgb9TNEZXt31KudLzjO1+1SmdJmCkXCTL1FVCWx6BZJWQ+veMHYh2Ldg6dKlTJkyRe5B7uTkxKZNm4iIiLi/m75HsotV7ErLJfZkLnvT86lQazE3MSKivTP/7N2W/v5utHJq+MF4Bf2TlJRESkpKrTFjY2MCAgIIDAykXbt2jdILTRGNZk5BQQHr16+Xj/38/OjduzeiKHLll1/J+egjTFxcaP3TsjobDYK0HbXp7CZmH5iNpYklCwcvJMzz740ba8g4DKv/CcUZ0P8t6PMqomDER/PmMXPmTPk0b29vtmzZQsBd1IDcL1qdSOKlIimIffKanbiXvQVjurdgQIAbPdu7YGmmpMQ2VyorK8nKyrrB8TkwMFB2U2jRogXdunWjc+fOt7XgaQwootGMuZkR4cMPP4yuvIKsWW9RujkGm8hIvOZ9eFdWIBdLLjL34Fz2Z+6nh1sPPun3CW5WN9nO0ulg//8gdg7YesFTm8E7DK1WyyvTp/Hll1/Kp3bt2pXNmzff1sKkviiqULO7Joj916k8CsvVGBsJBHk78u9hAfQPcMXf3bbB7jkr6B+dTsfZs2dJTEzk5MmTGBkZ8dprr9VKh/X29qZ///4EBATg5lb37dyGzh1FQxCE/wEdkFq1HgN+FUUx8fbPUmjoiKLIH3/8IdsSGBsbEx0djXDpEuenTUd94QKur76K85SnEeq4hFZr1fyQ9AOLji/C1NiUN0Lf4BH/RzA2usmv8NJsWPscnN0FHR+Gkf8DSwdUKhVPPPEEq1evlk/t378/a9eu1VvlqyiKnMkrl+zEU3M5cvEKWp2Io5Upkf5u9A9wo5+vK/ZWSkpscyc/P5/ExESOHz9OaWmpPK7VaklJSaF792u9XgRBoO9dGHU2Fuqy0kgFNgGmQEfgZ0EQFoii+JVeZ6agVxISEjh27Jh8PHz4cCwPHuT8u+9hZGOD948/Yh0WWufXO5h1kDkH5nC+5DzD2gzj9ZDXb766ADi1FdY9D+oKGPmF1KNbECgqKmL06NHs3r1bPjU6Opply5bVe0FTtVbHoXOFbE/NZcfJHC7UBLEDPGx5vl87BgS4E9jKAWMlJbbZo1KpSE5OJjExkYyMjJue4+npiYVF87B3uaNo/K1p0p+CIHyF1MpVEY1GSlZWFps3X6um7talCx5r15G1ahVWISF4ffYppnVcTudX5vPp4U/54+wftLJtxXeDvqNni543P1lTBdvfgwNfg3tnGP8DuPoDkJGRQVRUFElJ1xo3Tps2jc8//7zegoVXytXsOpXL9tRcdqflUVqlwczEiJ7tnZnSuy0DOrjTwqFx7zcr1C/79u1j165dciLG9VhbW9OlSxcCAwNxd3c3wOwMQ51jGjVtWn0AW6D+Gs4qPFAqKytZsxLUSQAAIABJREFUuXLlNSNCR0c6/vwLxSkpOD/7LK4vT0WoQ2W1Vqdl9anV/C/hf6i0Kp7v9jxPd3761kaD+aelYHf2cQh9Dga/D6bSuSkpKQwdOrTWr7hPPvmE11577b7iBqIokp5bxo6TuexIzeHIhat24uYM7+LJwA5u9PJxwdpcCe0p3Bw7O7tagmFkZISfnx+BgYH4+PjIdUzNibv5tPwJDAbGAh/qZzoK+uSqEWFRUREAZsbGhKz5HVGlouW332Bbx94TqQWpzD4wmxP5JwjzCOPN8Ddpa3+LzmCiCIm/1NReWMCjv4F/lPzw3r17GTVqlGyOaGJiwpIlS26wVagrGq2O+POFbEuR4hMXC6Vtp05edrzU34eBHdzpotiJK1xHVVUVycnJnD59mgkTJtRa2QYEBGBubi43NOrSpQvWdew+2VSpSyB8JfCOKIqpwPeCIPwIHEWKcyg0Ivbu3St3+wII3rULZ1dXWvz3v5i1bHHH55epy/g68Wt+PfkrDuYOzOszj+Fth996NaAqhk2vSrUXbfpItRd217Kf1q5dy2OPPSa3prSxsWHNmjUMGTLkru6rSqNlf3oBm5Oy2JaSw5WKasxNjOjl48Jz/doxIMANT3tl20nhGqIocv78eRITE0lNTaW6Wuonf+7cOdq3by+fZ2pqyosvvigXuirUbaXxM7BCkL4ZjgA2gE6vs1Kod86dO8fOnTvlY//UVLr07YfbG//G6A7eUaIosu3CNj6K/4i8yjyi/aOZ2n0q9ua3yWbKPgErJkHRJRjwFvR+Fa7Lovr222956aWX0Omkt5KbmxubN2++bZOn66lQa/grLY+Y5GxiU3MprdJga27CwA5uDOvsQV8/V6zMlG0nhdpcuXKFxMREjh07RnFx8Q2PJyYm1hINQBGMv1GXQPgGYIMgCF2BQMAIaatKoZFQUlJSy4jQNb+AoZMm4ThixB2fe6nkEnPj57Lv8j46OHXgv/3/SxfXLrd/0onVsP4lsHSQay+uIoois2bNYu7cufKYj48PW7ZsuaE46u8UV1YTezKHmKRs/jqVh6pah5O1GcO7eDKsiwc92zsrfScUbkCtVpOSkkJiYiIXLly46Tlubm4EBgbSpcsd3tsKdY9piKJ4HDj+/+3deVxVZf7A8c/DIoIIbqgIiiAoAXKv+9ZmaqKGW4plOuqUZo2/1mlqaqppm/apabIps9JMy11xCdN2M3M9FxBFcUVABRFkkfU+vz8uHrmCekUuF+V5v16+4pzzHO733g73e5bn+T52jEWxg/OFCAsLLff2G5eWMmHGdJpf4Y9jz+k9LExayLdHvsXN2Y2nez3NPaH34OJ0mUOmvAw2vQi/fwgd+sP4edD0Qq+S0tJSZs6cyeeff66v69WrF+vWrcPHx6faX5mVX8zGJEui2HIwi9JySRsvNyb0bE9UhC+9OjbHRdV2Ui5jx44dbNy4scp6d3d3IiIi6NatG23btlWDNW2krt9vcN+tWUNqaioAQkrGT5pEy0uU4Sgzl/H9se9ZuHchu0/txsPFg/Gdx3N/xP20aXKFLoUFWbB0qmXei94PwtDXwPnCYLiCggJiYmKs5hwfNmwYS5curfJgMfdcKbGmdNaa0tl+JBuzhA4tPPjzgECGRrTF6N9MPchWqpWXl1fldlLXrl3ZtGmTPtlRcHAwRqORzp0719s5WOoz9YndoKSU7FywgG2HD+vr7hg8mKBqEkZucS7L9i/jm+RvOFFwAj9PP57q+RRjQsbQtJEN93PTdsLiP0FhFoz5BAz3WG3OzMzkrrvuYtu2C3N5T506lTlz5uhTyEop2XoomyU7UlmfkEFxmZmQ1p7MuiOEqPC23OSrynYo1SspKWHv3r1omsaRI0d46KGHrMp2NG3alF69euHl5UVkZKR6RnGNVNK4AZWeOkXyK6+ywdsLKr6Uu3TpwoABA6zapZxJYeG+haw9uJai8iJ6t+3Ns72f5Vb/W6sv/VGdXQtg3ZPg2Qb+vAHaGa02Hz58mKFDh3LgwAF93XPPPccrr7yCEIKTZ4tYtvM4S3ekcuR0IU0buzC+pz/39OpAhJ99yoYo1z8pJampqWiaxp49eygpKdG3aZpWpQfesGHDLv4VSg2ppHEDkVJyNjaW42+8yaZ+fSmrSBjNmzdn9OjRCCEwSzO/Hv+VBXsX8EfGH7g5u3FX0F3cG3ovXVp0sf3Fykog7hnY8RkE3Q53fw5NrOcR3717N8OGDePkyZOApRbPhx9+yIwHZ7Jp7ykWbz/Gj8mZlJslfQJb8MigEIZF+KqKscol5ebmYjKZMJlMZGdnV9kuhNCf3yn2oZLGDaL05ElOvPAieT//zK7hw8n19gIsg+ViYmIocyrjq6SvWLRvEal5qbT2aM2j3R/l7pC7ad64+dW92NkMWDoFUv+AAY/CHS+As/WhtGnTJsaMGUN+fj4Abm5uvPfx5+S26U6/N34gM68Yn6ZuzLg1iJie7Qls1bAHTCmXl56ezvfff2811W9lLVu2xGg0EhkZiZeXVx1H17CopHGdk1KSu3IVJ19/HVlaStasv3Do9Gl9e5+BfZh/bD4rf1hJQWkBBh8Dj3R7hEEBg3B1qkHV1mNbYcmfoDjf0jsqfEyVJgsXLmTq1Kl6+YUmTb3pNeN13tznjfP+Qwzs4sOEXh24vYuPmtVOsYmTk1OVhOHm5kZ4eDhGoxF/f3/1zKuOqKRxHSs9cYKMF16g4Jdfce/ZA/H44/xSqRBhYetCHj/wOC5OLgztOJRJN00iolVEzV5MStg+13JLqlkHmLwK2oRd1ETy7rvv8tRTT+nrXLxa4TXuJcxtQnmqZ3vG9fCnjVfDqAaqXL28vDxMJhPdunWz6lXXtm1b2rRpw8mTJwkKCsJoNBIaGqp3pFDqjkoa1yEpJbnLl3PyjTeR5eW0+cc/cBs9ik8+naMXIjzT6AyJzRN5+KaHGRcyDh+P6sdB2KS0CNY9YakhFTLUUg7E3XpSJrPZzF//+lfee+89fZ1rqwAm//N/PBDVm75BLdSZoFKtsrIykpOT0TSNgwcPIqXExcWFvn37WrWLjo7G09PTbvOqKLZRSeM6U5qeTsbzL1Dw22949O6N72uvcsijgG/mvYM4a/lSLncu5/YRt/N2xNs1uwVVWU6qpRxIhga3PQO3PQ0XlSovLi5mypQpLF68WF/Xtks31q9ZTbeQ9tf2+soNSUpJeno6mqaRmJio1x87T9M0+vTpY3Wi4ed35fpoiv2ppHGdkFKSs2Qpp956CyklPs8/y85+rXgh8R8U7i8kIufCbaeJ4ybWzjzax7bCNxOhvLRKddrzcnJyGBh1F9ofv+nr7hgWzboVSxrMpDSK7fLz8/XeT5mZmdW2CQwMxGAw1HFkiq1U0qjnzIWFFG7fTvb8+RRs+R3X3t3ZPDmSBWfmcerXU4QRZpUw+vfvXzsJI34prH4YvNvDxMXQKqRKk0079jJuzEhyj6fo6x5++GE++OCDBjnPgHJ5aWlpfPbZZ3oNtMqaNWuG0WjEYDDQ7Crmo1fqnkoa9Yw0mynau5eC37ZQ8NtvFO7aBaWlSA93fr8vkg8D9lCWHs8AvwE83f5pTLEmCrH0Sw8ICGDQoEHXGICEn9+En163lDOP+RI8Wlg1Sc0u5NkvNvD1yzMpP3tKX//aa6/x97//XT27UPTEUPlY8PX1xcPDg4KCAsBSdjwsLAyj0UhAQIA6bq4Tdk0aQogo4D+AMzBXSvnGJdr1ArYCE6SUy+wZU31UevIUBb/9Zvm3ZQvlFRMSuYWG0mRiDEubHWBho524uR/jnuB7mdBlAu092zNv3jx9IJOnpyfjxo27tqlRS4sgdhYkLAXjfXDX++ByoWx6TmEJH/6QwidL40hf+k/M5/IAcHZ2Zu7cuUydOrXmr63cEAoKCoiPj0fTNG699VbCw8P1bU5OThgMBtLS0jAajYSFhdHoCmX5lfrHbklDCOEMzMYy299xYLsQIlZKmVRNuzeBDfaKpb4xnztH4Y6deqIoriix4dyqFZ633kKTAQNo0q8fW4r38uKWF8kpzuEhwyNMumkSHq4eAMTFxenTowohGDduHJ6enjUPqiDL8vwi9Q8Y9CLc/DhUnPkVlZbz5e9H+PCHFE4m/kb2mrcxl1geXHp4eLBs2TJVpqEBKy8v58CBA2iaxoEDB/Q5Ukwmk1XSABg8eLC6orjO2fNKozeQIqU8BCCE+AYYBSRd1O7/gOVALzvG4lDm4mKK9+6lcKclURTu2IksKUE0aoRHzx54jx5Nk5sH4Na5M0IICkoLeHX72yw/sJyQ5iF8PPhjqxIfiYmJ/PHHH/ryoEGDCAgIqHmAmcmwcDzkn4Tx8yF8NGC5xRBrSuetuGTScs7R7uQW9qx8Q/9S8PHxYd26dfTqdcP+r1Mu48SJE2iaRkJCQrWlO44ePUpxcTFubm76OpUwrn/2TBp+QGql5eNAn8oNhBB+wBjgDi6TNIQQM4AZAB06dKj1QGuTNJspOXKUc/EmiuITOBcfT1FyMlRMJ+kWEkLziRNpMmAAHj174ORuPQ3pzpM7eW7zc6TnpzMtYhqzjLNo5HzhEj4rK4s1a9boy6GhofTv37/mAR/6yVKh1sUNpq4H/x4A5BeX8eyKBGJN6YT5NiU04xc+n/eWvltQUBAbNmwgODi45q+tXHeKi4vRNA1N0zhx4kS1bdq3b4/RaCQ8PNwqYSg3BnsmjepOKS7uNvE+8LSUsvxyZyBSyjnAHICePXtW7XrhQGVZWZyLT7iQJBITMZ89C4CThweNu3al5dQpNI6MxN1gwLVSyebKSspL+HD3h8zbMw8/Tz/mRc2jexvrqU9LSkpYsmSJXtGzefPmjBo1quZnbzvnWSrUtups6SHVzJKQk9LP8pdFuzh6uoAnBnciccm/+fTTT/Xdunfvzvr162nT5gpzbCg3HCklmzZt0kvEnHe+7LjRaKRly5aX2Fu5EdgzaRwHKo/s8gfSL2rTE/im4kuvFTBcCFEmpVxlx7hqREqJ+exZig8c4JwpnnMJlkRRlp5haeDsjFvnzngNG4Z7ZFfcIyNpFBSEsKHr6b7sffz917+TkpPC+M7j+WvPv+rPLiq//tq1a/W+7ecLEdZoLITZbJlhb8sHEDwYxn0Bjb2QUvL1tlT+uWYPzT1c+WJSJP9+dhaxsbH6rnfeeSfLli1TcxI0AKdOncLFxYUWLS70nmvcuDGhoaEkJibi4uJCaGgoRqORwMDAa+uEoVw37Jk0tgMhQohAIA24B5hYuYGUMvD8z0KIecBaRyUMKSXlOTmUpqVTmp5m+W9amuVfuuVnc0XFVgBXPz88jEYaT/4T7pFdaRwWVuVW05WUmcuYt2ces7XZNHNrxuxBs7nV/9Zq2+7YsYOEhAR9efjw4bRt2/bq32hJAayYAfvWQq/pEPUGOLtY3Y66JaQVzw9uz7SJ49myZYu+66RJk/jss89Uj5cb2Llz50hISEDTNDIyMujevTvR0dFWbfr06UNAQAARERFqAGcDZLekIaUsE0LMwtIryhn4XEq5Rwgxs2L7x/Z67csp2ruXkqPHqiSE0rQ0zBc9zHNq0gRXPz9LgujVC1c/Pxp1DMA9MhKXa7wEP3r2KM9tfg5TpomhHYfyjz7/oFnj6gc1paWlERcXpy9369aNbt26Xf2Lns2AryfAiQSIehP6zgRgb8ZZ/rJwF0dOF/DXOzszItCV4VGD2Ldvn77r3/72N15//XV1NnkDMpvNHDx4EE3TSE5O1uuXAezZs4eoqCirwoD+/v74+/s7IlSlHrDrOA0p5Xpg/UXrqk0WUsqp9ozlvPSnn6F4/34AnDw9cfX3x7V9ezz69sXVr50lMVQkCicvr1rv7SGlZEnyEt7d+S4uTi68detbDAu8dHfVwsJCli5dqvdYatu2bc26t2bEw6IJUJQL93wNXaKQUrJ4eyovxu7B292VRdP74pGfxoABt5GRYbntJoTgvffe49FHH63R+1Xqr8zMTDRNIz4+Xp/3pDJnZ2eCg4MpKipS1WQVXYMbEe778ksINzdc/fxwruPJWk4WnOTFLS/yW/pv9G/Xn5f7v0ybJpd+mCylZOXKleTm5gKW+QNiYmKu/g84OQ6W/dlSmfb+DdC2KwXFZTy3MoFVmuV21HsTjCTu+J2ho0ZxtuJBfqNGjViwYAExMTE1fs9K/VNaWsr8+fNJS0urdrufnx8Gg4GIiAjcr/KWq3Lja3BJw91ovHKjWlZuLmdVyir+vfPflJpL+UeffxDTJeaKVzG//PILKSkX6jqNGTOG5s2vcpa9HZ9beki1jbT0kGra1nI7atEujmRZbkc9fHswy5YtZfLkyXrPLC8vL1avXs3tt99+tW9XqedcXV1xcbH+0/f09CQyMhKDwUDrS/TwUxRogEmjrm1J38I7O97hwJkDdG/dnZcHvEyA15UH4qWkpPDTTz/pywMGDKBLl6uYw1tK+OVt+PE1yxwY479AunqweNsxXozdg5e7Kwsf6Eu/Ti354IMPeOyxx/R6Qb6+vsTFxREZGXm1b1epR7KysjCZTDRt2pTevXtbbTMYDKSmphIaGorBYCA4OFg9r1JsopKGnRzMOcg7O95hc9pm/Dz9eOe2d7gz4E6bnpHk5uayYsUKfbljx47ccccdtr+42QxxT8O2OWC4F0b+l4IywT+WmFi5O42bgy23o1o2ceXpp5/mrbcuDNoLDQ0lLi7u2kaYKw5TVFTEnj170DRNLzPTrFkzevXqZXXsRURE0KVLFzw8PC71qxSlWipp1LLT507zkfYRyw8sx8PFgyd7PMnEmyZajeq+nPLycpYuXcq5c+cAy22Du+++2/azwLISWDUTEpdDv1kw5BWSTxXw8MKdHM4q4IkhnfnLwGDKy0qZMmUKX331lb5rv379WLNmjRqcdZ2RUnL48GE0TWPv3r1VBt7l5ORw9OhROnbsqK9zdXVVD7eVGlFJo5YUlxezIGkBcxPmUlxWzIQuE5hpmEnzxlf3DGLDhg36A8qrLkRYnA9LJsPBH2DIyzDgUb5NyODxJRpNG7vy1QN96N+pFXl5eYwbN47vvvtO33XkyJF8/fXX6szzOpKdnY2maZhMJr3zQmVOTk507twZo9FI+/ZqBkWldqikcY3M0sy3h7/lP7v+Q0ZBBre3v50nejxBoHfglXe+SGJiItu3b9eXBw8ebPttooLTsGg8pO+GUbOh2yQ+23yYV9cl0a19Mz6Z3BOfpm6cPHmSESNGsHPnTn3XGTNmMHv27CoPR5X6S0rJwoULyc7OrrKtTZs2GI1GunbtSpMmTRwQnXIjU98S12D3qd28vf1tErISuKnFTbw64FV6+/a+8o7VyMzMtCrXcdNNN9GvXz/bds5JhQVjIDcVJiykvPMwXl2zhy9+O0JUeFvev8dIY1dnUlJSGDp0KIcOHdJ3femll3j++edV9dF6TEpJQUGB1RWnEILIyEi9s4SHhwcRERF069atZpUCFMVGKmnUQOrZVN7b9R4bj26ktXtrXh3wKtGdonESNet9cr4QYWlFJdwWLVowcuRI277IT+2Dr8Zabk1NWkGRX18eW7iLuD0n+POAQJ4bcRPOToLt27czYsQIvXaVk5MTH3/8MdOnT69RzIr9nTlzRr/95O3tzbRp06y2G41G0tPTMRqNdO7cWU2xq9QJlTSuQk5RDp8mfMqifYtwdXLlYePDTAmbUqW44NWQUrJmzRqysrKAqyxEmLrdckvKuRFMW0d20y488OlWdqfm8PxdYdx/s+UWWVxcHOPGjdOn2XR3d2fx4sVVagopjldSUkJSUhKapnH06FF9fW5uLtnZ2VbFA729vbn33nsdEabSgKmkcQXl5nJ+z/idVSmr+OHYD5SZyxgdPJpZ3WbR2uPaB0Ft376dxMREfXnEiBG2lRw/sBGW/Ak828DklRwxt2bqR7+RkVvERxO7M6yrLwBffvkl999/v96jpkWLFqxdu9b2W1+K3UkpOXbsGJqmsWfPHv2KszJ3d3eysrKskoaiOIJKGpdwJPcIq1JWsebgGk6dO0Uzt2bEdIlhXMg4gpvXzsRDx48fZ8OGC7Pcdu/eHaMtI9bjl8Cqh6D1TTBpBbuyXXlg/haklCya3oceAS2QUvLGG2/w7LPP6rsFBAQQFxdHaGhorcSvXLutW7eybds2zlTMC1+ZEILg4GD99pPqqKDUB+oorCS/JJ8NRzawKmUVWqaGk3DiZr+beSb4GW7zv83msRa2uLgQoa+vr22FCLf+D+KegY63wD0LiUs5x6PfbKWtd2PmTetNYKsmlJeX8+ijjzJ79mx9t8jISL799lvatWtXa+9BuXZZWVlVEkarVq0wGo1ERkaqeUuUeqfBJw2zNLP9xHZWp6xm49GNFJUXEegdyOM9Hic6KBofD5/af02zmRUrVuh96xs3bsz48eMvfyYpJfzwCvz6LtwUDWPnMm9bBi+tTcLg34y5U3rSytONoqIiJk2axPLly/VdBw4cyMqVK/H29q7196JcmZSS1NRUMjMz6dGjh9U2o9HIzp07ady4MRERERgMBvz8/FRvNqXearBJIy0/jdUpq4k9GEtafhqerp5Ed4pmdPBourbqatc/2l9++YWDBw/qy1csRGguh7WPwa4vocdUzMPe5V9x+5m7+TBDwtrwwT3dcG/kTE5ODqNGjeKXX37Rd50wYQLz589XczU7QG5uLvHx8WiaRnZ2Ns7OzoSFhVlVjvXz82PChAkEBwer20/KdaFBHaVSStYeWsuqlFVsO7ENgaCPbx/+r9v/MajDIBq72H8WspSUFH7++Wd9+eabb6Zz586X3qG81DLT3p4VcOtTFN38DE8sNrE+4QRT+gXwQnQ4zk6C48ePM2zYMKuH6o899hjvvvuuKkRXh0pLS9m3bx+aplmNhwFLiZjExER69eqlrxNCqGdMynWlQSUNIQTf7PuGM8VnmGWcxchOI/H19K2z17+4EGFgYCADBw689A5lJbD8ftgbC0Ne5ozxIaZ/to0dR8/w3PCbeOCWQIQQ+uxq5wvUAbz99ts8+eST6jZHHZBSkpaWhqZpJCYmUlxcXKWNm5sb4eHhqpyHct1rUEkD4IM7PqBF4xZ1/mVaVlZmVYiwadOmly9EWFYMS6dB8joY+i+OdZ7G1P9t4XjOOWZP7M6ISEuy27x5M9HR0eTk5ACWcR7z5s3jvvvuq5P3pcC6deusyrJUFhQUhNFoJDQ0VBUIVG4IDS5ptHR3TAXX7777Ti9E6OTkxLhx4y5dF6i0yDIG48AGGPY2CX4TmPrRb5RLycIH+tCro6Wv/ooVK5g4caJ+Zuvp6cmKFSsYMmRInbynhkhKWeWEIygoyCppNG/eHKPRiMFgUJ0PlBtOg0sajpCQkFClEGGHDh2qb1x6DhZPgpRNcNd7bGs5mj9/uhVvd1e+vL83nXws9Yc++ugjZs2apU+c1KZNG9avX0/37t3t/n4aGikl6enpaJrGkSNHmDlzplXJjs6dO+Pt7U1gYCBGo5EOHTqo24LKDUslDTvLzMxkzZo1+vJNN91E3759q29cUgjf3AuHfoaR/+WnJlHM/PwP2jVzZ+EDffD1dkdKyfPPP89rr72m7xYcHMyGDRsICgqy99tpUPLz8/XeT+drdoGlM0PlWRRdXFx45JFHVIcDpUFQScOOiouLrQoRtmzZklGjRlV/FlpSAIsmwJHNMPojvnUeyCNf7iCkdVO+vL83rTzdKC0tZebMmXz++ef6br169WLdunX4+NT+eJKGqLy8nOTkZEwmEwcOHNCv5Crbv39/lal3VcJQGgqVNOykukKE48ePr368RHEeLIyB1K0w9lOWlfbjb4t3YWzfjC+m9cbb3ZWCggJiYmJYv369vtvw4cNZsmSJmjOhFpw6dYqdO3eSkJCgd1aozNXVlbCwMIxGo5oKV2nQVNKwk23btrFnzx59OTo6uvpChEVnYeE4OL4D7v6M+We782KsiQHBLZkzuSdN3FzIzMzkrrvuYtu2bfpuU6dOZc6cOapHTi05fPiw1ed7XkBAAAaDgbCwMDVAUlFQScMuUlNTraZS7dGjB5GRkVUbnsuxzIWRYYLxXzD7ZDhvb9jDkLA2/PfebjR2debQoUNERUVx4MABfbfnnnuOV155RT1srYHy8nKOHDlCUFCQ1efXtWtXvvvuO8xmM97e3hgMBgwGg6oqqygXUUmjlhUUFLBs2TK9EGG7du2Iioqq2rAw2zLb3sk9yPHzefNIMB//nMwoYzveGW/A1dmJXbt2MXz4cE6ePAlYBid++OGHPPzww3X5lm4IJ06cQNM0EhISKCwsZMaMGfj6XhjY6eHhwZAhQ2jdujWBgYEqISvKJaikUYtsLkRYcBoWjILMZMwxX/HCPj++2nqQiX068OqoCJycBBs3bmTs2LHk5+cDlhHFixYtYuzYsXX9tq5bhYWFJCQkoGkaJ06csNqmaZpV0gAu3atNURSdShq16Oeff7aqNzR27FiaNWtm3Sg/E74cBdkHKZ+wiKd2+7Bi9zEevDWIZ4aFIoRg4cKFTJ06VZ84qVmzZsTGxnLLLbfU5du5LpWXl3Pw4EE0TSM5OVm/4qvMy8sLLy8vB0SnKNc/lTRqyYEDB6yqy95yyy2EhIRYN8o7CV+OhDNHKZnwNf+31YsNe9L4652d+ctAy8RO77zzDk899ZS+i7+/P3FxcYSHh9fJ+7ie7d+/n9jYWH1a28pcXFwIDQ3FaDQSGBiousgqSg2ppFELcnJyWLlypb4cFBTE7bffbt3obAbMj4az6RRN+Ibpv7jz64GTvBgdxrQBgZjNZp588knef/99fZfw8HDi4uLw9/evo3dyfWvevHmVhOHv74/RaCQ8PNy2edcVRbkslTSuUXWFCMeOHWt9JpufabnCyMsgf/w3TN3kwq5jWbw1LpKYnu0pLi5mypQpLF68WN/lllsiyaJ5AAAU/klEQVRuYfXq1ZefZ6MBMpvN+u2nqKgoq5ntfHx88PPzIzc3F4PBgNFopFWrVg6MVlFuPCppXKMNGzaQnp4OWEYFjx8/3nqwXWE2LBgNOankjlvMxDjB/pM5/PdeS6Xa3NxcxowZw48//qjvcvfdd/PVV1+pM+NKMjMz0TSN+Ph4vXOAr68vN998s1W7mJgYPD091e0nRbETlTSuQXx8PDt27NCX77zzTuv5EopyLeMwsg6QPWoB49dJjp/JZ86fejKwS2vS09MZNmwY8fHx+i5/+ctf+M9//mNVEK+hKioqIjExEU3T9ArBlZlMJgYMGGDVPVY94FYU+1JJo4ZOnTrF2rVr9eXw8HB69+59oUFxvqU0yIkETt/1OWM3NCIrr4j5f+5N36CW7Nu3j6ioKI4eParv8q9//YtnnnmmQY8RMJvNHD58GE3T2Lt3L+Xl5VXaNGnShMjISIxGY4P+rBTFEeyaNIQQUcB/AGdgrpTyjYu23wc8XbGYDzwkpTTZM6bacHEhwlatWhEdHX3hC6z0nKVa7fFtnB7+CWM2eXGmoIQFD/She4fm/P7779x1111kZ2cD4OzszGeffcaUKVMc9ZbqjSNHjvDVV19VWe/k5ESXLl0wGo0EBwer20+K4iB2SxpCCGdgNjAEOA5sF0LESimTKjU7DNwmpTwjhBgGzAH62Cum2iClJDY2ltOnTwOWQnZWhQjLimHxZDj8K6fv/IBRP7Qi95wlYRjbN2PNmjVMmDBBf3DepEkTli1bVv2o8RtccXExzs7OVoMfO3bsiJeXlz5A0tfXF4PBQNeuXfHw8HBUqIqiVLDnlUZvIEVKeQhACPENMArQk4aUckul9luBet+39I8//iAp6ULei46OpnXr1paF8jJY9mdI2Uj2HW8x8hd/8opKWfhAHyL9mzF37lwefPBBfcCZj48P69ato1evXo54Kw4hpeTw4cOYTCaSkpKIjo62qsvl5OREnz59yMvLw2g0Vl/kUVEUh7Fn0vADUistH+fyVxH3A9/aMZ5rduzYMTZu3Kgv9+zZk65du1oWzOWwaibsW0v2LS9z12/BFJSUsWh6X8LbefHSSy/xz3/+U983KCiIDRs2EBwcXMfvwjGys7P13k+5ubn6epPJVKWYY//+/es6PEVRbGTPpFHdE8qqM9oAQoiBWJLGzZfYPgOYAVx6mlQ7u7gQoZ+fH0OHDrVsNJthzaOQsJQz/f7OiG0RnCstZ+EDfQht04QHH3yQTz/9VP9d3bt3Z/369Tf8WXRxcTFJSUlomsaxY8eqbVNYWEhZWVnV+lyKotRL9vxLPQ5U6n+KP5B+cSMhRCQwFxgmpTxd3S+SUs7B8ryDnj17Vpt47MlsNrN8+XLy8vIAcHd3v1CIUEqIexp2L+BMz8eI2tmT0nIzix7oS8dmLtx9993Exsbqv+vOO+9k2bJlVoPSbjR5eXl8//33JCUl6Z0FKnN3d6dr165069aNtm3bOiBCRVFqyp5JYzsQIoQIBNKAe4CJlRsIIToAK4DJUsr9dozlmvz0008cPnxYXx47dize3t6WhLHxBdg2hzOGB7lTuxmzlHw9vS+tXEsYPDiK33//Xd9v8uTJzJ07l0aNGjnibdSZRo0aVUkYQghCQkIwGo107txZjUNRlOuU3ZKGlLJMCDEL2ICly+3nUso9QoiZFds/Bl4AWgIfVXRXLZNS9rRXTDWxf/9+fv31V335tttuu/Ac4uc3YcsH5IRPYUjiEBDw9Yy+uBVlc/PAKPbt26fv9/TTT/P666/fUOMKSkpK2Lt3L61bt7YqM+7m5kZYWBgmkwkfHx+MRiORkZF4eno6MFpFUWqDkLLO7/Zck549e8rKo7DtKScnh08++YSioiIAOnXqxMSJEy1jBDa/D5teJLfLBAan3A1Oznw9vQ+FJw4TFRVFRkYGYDnDfv/993nkkUfqJGZ7k1Jy7NgxNE0jKSmJkpISDAYDo0ePtmqXlZVFSUkJvr6+N1SiVJTrlRBiZ22clKunj5dQVlbGkiVL9ITh5eV1oRDhH5/Aphc5GzySQQfG4uTszNcz+pK6ZwejRo3Sxxg0atSIBQsWEBMT48i3Uityc3PRNA2TycSZM2estiUlJTF8+HCr226qUKCi3JhU0riEuLg4/WrhfCFCDw8P2Dkfvv0bZzsOZdDBe3F2deHr6X3Z8eN6Jk+eTElJCWBJMqtXr65aIv06Ulpayt69e9E0zeqZTmWtWrXCaDRyvV2xKopSMyppVMNkMrFz5059eejQoZY5LeKXwppHyfO/nUFHp+DSyI2vp/dl9cK5PP744/oXp6+vL3FxcVXGH1xPiouLef/99/Urrcrc3NyIiIjAaDTi5+enbj8pSgOiksZFTp48aVWIMCIiwjJie/8GWPkg+b59GXT8ARo1dmfh/b3539sv89Zbb+ntQ0NDiYuLIyAgwBHh15iU0urL383NjXbt2llNX9upUyeMRiNdunTB1dXVEWEqiuJgKmlUUlRUxJIlS/S5ufVChMd+hyV/orBlGEMyHsS1cRO+nNqd5594yKq4Xv/+/YmNjaVly5aOegtXpbS0lOTkZDRNw8/Pj4EDB1ptNxgM5OTkYDQaMRgMquy4oigqaZx3vhDh+cqzrq6uxMTE0Oj0Plg0gaImfgzLehQXD28+vSech6dM4LvvvtP3HzlyJN988w3u7u6Oegs2kVKSlpaGpmkkJiZSXFwMWCY5uu2226yqx0ZERNC1a1d1+0lRFJ1KGhW2bt3K3r179eWRI0fi43QWvrqbUpcmjMp9krLGLflwVCD3jRnOrl279LYzZsxg9uzZ9boURl5eHvHx8WiaRlZWVpXtZ8+eJSMjAz8/P32dKj+uKMrF6u+3XB26uBBhr169iOjQEj6/k7LyUmKKnyXPrS1vDvZh3PDBVvf5X3rpJZ5//vl6eTYupSQpKQmTyURKSkq1PZyaN2+OwWDAYDDQrFkzB0SpKMr1pMEnjfz8fJYuXap/ofr5+TH0ll7wZTTmgtNMLX+eDNcOPNvTlbuHDyIzMxOwnIV/8sknPPDAA44M/4p++umnKlcWrq6uhIeHYzQa6dChQ71MeIqi1E8NOmmcL0SYn58PVBQiHDUC58X3Ik8f5GH+zn7nEB4KOsu9oydRWFiot1u8eDHR0dGODN9Kfn4++fn5VgUAhRAYjUY2bdoEWCY4MhgMhIWF3fD1rxRFsY8GnTR+/PFHjhw5oi/fPXok3nEPI9N28jenv7JDRDK+aQrTJ83S56pu0aIFa9eupV+/fg6K+oLy8nL279+PpmkcOHAAX19fpk+fbtUmMjKS0tJSDAYDzZs3d1CkiqLcKBps0khOTmbz5s368m233kon0xtw8Adec36Y72UvBhZt5pnnX9LbBAQEEBcXR2hoqCNC1mVkZKBpGgkJCfq0sQDp6emcOnXqwkyCQNOmTa/rUemKotQvDTJpnDlzhlWrVunLnTp14ra81bBnJf91/hPLzbcRdnQZ//7iwsRJBoOB9evX065dO0eETEFBAQkJCWiaxsmTJ6tt06FDh2rnr1AURaktDS5pXFyI0Nvbm7GtDiG2fsEC5zF8WjKUtrv+x6L1FyZOGjhwICtXrrTMoeEAsbGxmEwmfdbAyry9vfXeTy1atHBAdIqiNCQNLmmsX7+eEydOABWFCEPBY+u7xDoP4Y38Ebj//Dob/9iit58wYQLz58/Hzc3NUSHj6upqlTBcXFwICwvDYDAQGBioej8pilJnGlTS0DSN3bt368tRES3x2/p3fnLux9NnRlG67gWS9l+YOOnxxx/nnXfeqZNBboWFhSQkJJCXl8fgwYOtthmNRrZt20b79u0xGo2EhYXRuHFju8ekKIpysQaTNMrKyvjxxx/15a4dWtAj/h/scDIwMz2anNXPcDLjwhTm77zzDk8++aRdYzKbzaSkpKBpGsnJyZjNZpycnOjXrx9NmjTR27Vt25ZHHnlE9X5SFMXhGkzScHFxYdq0aSxdupTSglxGpL5Osghi0qEo0lc+R97ZXMByK+iLL77gvvvus1ssp06dQtM04uPjKSgosNpmNptJSEigb9+++johhEoYiqLUCw0maQA0a9aMaVHdyf/yXk7QgtEJt3Mk9mVKSixF+zw9PVmxYgVDhgyp9dc+d+4ciYmJaJpGenp6tW38/PwwGo1ERETU+usriqLUhgaVNCgvQyz/M8VlZgZv7cnBDf/Wy4e0adOG9evX0717d7u89Nq1a0lKSqqy3tPTU+/95OPjY5fXVhRFqS0NKmmU4cTjZY+y7rvlHN4yX18fEhJCXFwcQUFBtfI6hYWFlqlhK4mMjNSThrOzM6GhoRgMBjp16qSqySqKct1oUElDmstJ+nUDh7es19f17t2btWvXXvNZflFREYmJiZhMJk6fPs0TTzxhVSo9JCSE4OBgOnfuTERERL2fd0NRFKU6DSZpFBUVMXbsWH749lt93fDhw1myZIlVT6WrYTabOXz4MJqmsW/fPn3GP4D9+/cTFhamLzs5Odn14bqiKEpdaDBJw83NzepqYtq0aXzyySc1muv69OnTeu+ns2fPVtnu5OTE6dOnryleRVGU+qjBJA0hBHPnzuXUqVP06NGDV1555apHUu/fv5/NmzeTmppa7fa2bdtiNBrp2rVrlWcaiqIoN4IGkzTAMgZjzZo1NZ6W9ezZs1UShoeHB5GRkRgMBqu5LBRFUW5EDSppADYljOzsbJKTk+nbt6/V1Uh4eDhxcXFIKQkJCcFoNBISEoKzs7M9Q1YURak3GlzSuJTi4mJ9Pu2jR48C0L59e/z9/fU27u7uxMTE4OfnV+OH54qiKNezBp00pJQcPXoUTdNISkqqMheFpmlWSQOgc+fOdRmioihKvdIgk8aZM2cwmUyYTCZycnKqbBdCEBISohKEoijKRRpc0ti1axdr1qypdpuPj4/e+6lp06Z1HJmiKEr91+CSRseOHa2WGzduTEREBEajkXbt2qkJjRRFUS6jwSWNFi1aEBAQgKurK0ajkS5dutS4C66iKEpD0yC/LSdPnqy6ySqKotRAgyyvqhKGoihKzdg1aQghooQQyUKIFCHEM9VsF0KIDyq2xwsh7DOZhaIoilIr7JY0hBDOwGxgGBAG3CuECLuo2TAgpOLfDOB/9opHURRFuXb2vNLoDaRIKQ9JKUuAb4BRF7UZBXwpLbYCzYQQvnaMSVEURbkG9nwQ7gdUru53HOhjQxs/IKNyIyHEDCxXIgDFQojE2g3VLloBWY4OwgYqztp1PcR5PcQIKs7a1qU2fok9k0Z1Ax5kDdogpZwDzAEQQuyQUva89vDsS8VZu1Scted6iBFUnLVNCLGjNn6PPW9PHQfaV1r2B9Jr0EZRFEWpJ+yZNLYDIUKIQCFEI+AeIPaiNrHAnyp6UfUFcqWUGRf/IkVRFKV+sNvtKSllmRBiFrABcAY+l1LuEULMrNj+MbAeGA6kAIXANBt+9Rw7hVzbVJy1S8VZe66HGEHFWdtqJU4hZZVHCIqiKIpSrQY5IlxRFEWpGZU0FEVRFJvVq6RxLWVHrrRvHcZ4X0Vs8UKILUIIQ6VtR4QQCUIIrba6v11DnLcLIXIrYtGEEC/Yum8dx/lUpRgThRDlQogWFdvq8vP8XAhx6lJjhOrJsXmlGOvLsXmlOOvLsXmlOOvLsdleCPGjEGKvEGKPEOLRatrU3vEppawX/7A8LD8IBAGNABMQdlGb4cC3WMZ39AX+sHXfOoyxP9C84udh52OsWD4CtKonn+XtwNqa7FuXcV7UPhr4oa4/z4rXuhXoDiReYrtDj00bY3T4sWljnA4/Nm2Jsx4dm75A94qfmwL77fndWZ+uNK6l7Igt+9ZJjFLKLVLKMxWLW7GMPalr1/J51NVnWZPXuhf42k6xXJaU8hcg+zJNHH1sXjHGenJs2vJZXkpdHptXG6cjj80MKeWuip/zgL1YKmtUVmvHZ31KGpcqKWJLG1v2rasYK7sfS3Y/TwLfCSF2CktpFHuxNc5+QgiTEOJbIUT4Ve5bG2x+LSGEBxAFLK+0uq4+T1s4+ti8Wo46Nm3l6GPTZvXp2BRCdAS6AX9ctKnWjs/6NAnTtZQdsakcSS2w+XWEEAOx/GHeXGn1AClluhCiNbBRCLGv4mzGEXHuAgKklPlCiOHAKizVhuvqs+QqXysa+E1KWfnMr64+T1s4+ti0mYOPTVvUh2PzatSLY1MI4YklcT0mpTx78eZqdqnR8VmfrjSupexIXZUjsel1hBCRwFxglJTy9Pn1Usr0iv+eAlZiuTS0hyvGKaU8K6XMr/h5PeAqhGhly751GWcl93DR5X8dfp62cPSxaZN6cGxeUT05Nq+Gw49NIYQrloSxUEq5opomtXd81sWDGhsf5rgAh4BALjyQCb+ozQisH+Zss3XfOoyxA5YR7v0vWt8EaFrp5y1AlAM/y7ZcGNzZGzhW8bnWyWd5Nf/fAG8s95abOOLzrPSaHbn0w1uHHps2xujwY9PGOB1+bNoSZ305Nis+my+B9y/TptaOz3pze0peQ9mRS+3roBhfAFoCHwkhAMqkpQJmG2BlxToXYJGUMq62Y7yKOMcBDwkhyoBzwD3SchTVyWd5FXECjAG+k1IWVNq9zj5PACHE11h69bQSQhwHXgRcK8Xp0GPTxhgdfmzaGKfDj00b44R6cGwCA4DJQIIQQqtY9yyWk4RaPz5VGRFFURTFZvXpmYaiKIpSz6mkoSiKothMJQ1FURTFZippKIqiKDZTSUNRFEWxmUoaiqIois1U0lCUayCEGCyEWODoOBSlrqikoSjXxgDsdnQQilJXVNJQlGtjAHYLIdyEEPOEEP8SFUOBFeVGVG/KiCjKdcoAnMJShmGulPIrB8ejKHalyogoSg1VVBbNAo4CD0opf3dwSIpid+r2lKLUXBiwHSgDygGEED5CiC+EEP4Vc0y7OjRCRallKmkoSs0ZsJS9vgf4QgjRRkqZiaWU97vAI1LKUkcGqCi1TSUNRak5A5a5FvYDTwNLhBBNgSAsZcfzHRqdotiBeqahKLVECOECzAFeAmKA7VLKnxwalKLUMpU0FEVRFJup21OKoiiKzVTSUBRFUWymkoaiKIpiM5U0FEVRFJuppKEoiqLYTCUNRVEUxWYqaSiKoig2U0lDURRFsZlKGoqiKIrN/h8nmMpT0ueIDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(ks[:,0], freqs)\n",
    "ax.plot(ks[:,0], ks[:,0]/1.45, 'k-', linewidth=3, label='SiO2')\n",
    "ax.plot(ks[:,0], ks[:,0]/1, '-', color='gray', linewidth=3, label='Air')\n",
    "ax.plot(ks[:,0], ks[:,0]/3.45, '--', color='gray', linewidth=3, label='Si')\n",
    "ax.set_xlabel('$k_x$')\n",
    "ax.set_ylabel('$\\omega$')\n",
    "ax.set_xlim([0,2])\n",
    "ax.set_ylim([0,1])\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Addendum : saving/loading files from Python\n",
    "\n",
    "Google is your friend!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test\n"
     ]
    }
   ],
   "source": [
    "# Saving text files\n",
    "string = 'test' # could be straight up output.stdout\n",
    "text_file = open(\"output.txt\", \"w\")\n",
    "text_file.write(string)\n",
    "text_file.close()\n",
    "\n",
    "# Loading text files\n",
    "text_file = open(\"output.txt\", \"r\")\n",
    "print(text_file.readline())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0794202 0.0849623 0.297663  0.37481  ]\n",
      " [0.150218  0.160326  0.324018  0.383469 ]\n",
      " [0.219035  0.232794  0.36303   0.398284 ]\n",
      " [0.284828  0.301118  0.410931  0.417364 ]\n",
      " [0.34589   0.36477   0.440364  0.463544 ]\n",
      " [0.398211  0.424058  0.470557  0.516771 ]\n",
      " [0.436446  0.479309  0.513457  0.567861 ]\n",
      " [0.463503  0.529944  0.566628  0.615651 ]\n",
      " [0.485926  0.574409  0.623957  0.66078  ]\n",
      " [0.506776  0.611532  0.682582  0.698461 ]\n",
      " [0.527159  0.642054  0.725485  0.740367 ]\n",
      " [0.547509  0.667911  0.747291  0.775545 ]\n",
      " [0.568008  0.690839  0.766489  0.796961 ]\n",
      " [0.588733  0.711999  0.784403  0.816743 ]\n",
      " [0.609713  0.732111  0.80173   0.835754 ]\n",
      " [0.630952  0.751623  0.818853  0.854324 ]\n",
      " [0.652443  0.770818  0.83599   0.872643 ]\n",
      " [0.674172  0.78988   0.853267  0.890837 ]\n",
      " [0.696125  0.808928  0.870757  0.908996 ]\n",
      " [0.718285  0.828044  0.888503  0.927185 ]\n",
      " [0.740637  0.847281  0.906527  0.945451 ]\n",
      " [0.763167  0.866673  0.924839  0.963831 ]]\n"
     ]
    }
   ],
   "source": [
    "# Saving numpy arrays\n",
    "array = freqs\n",
    "np.save('array.npy', array)\n",
    "\n",
    "# Loading numpy arrays\n",
    "array_from_file = np.load('array.npy')\n",
    "print(array_from_file)\n",
    "\n",
    "# Saving general objects : pickle, dill, etc.\n",
    "\n",
    "# MEEP/MPB likes the h5 format; there is also a package for that (see next notebook)\n",
    "# http://docs.h5py.org/en/stable/quick.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Waveguide parameters\n",
    "\n",
    "There is an inverse-solver that can return the wavevector at a specified frequency (see https://mpb.readthedocs.io/en/latest/Python_User_Interface/#the-inverse-problem-k-as-a-function-of-frequency). It needs a few parameters, since it is an iterative solver :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "f_mode = 1/1.55   # frequency corresponding to 1.55 um                                                                                                             \n",
    "band_min = 1\n",
    "band_max = 1 # Look at the lowest mode\n",
    "kdir = mp.Vector3(1)\n",
    "tol = 1e-6\n",
    "kmag_guess = f_mode*3.45\n",
    "kmag_min = f_mode*0.1\n",
    "kmag_max = f_mode*4.0\n",
    "\n",
    "k_fmode = ms.find_k(mp.NO_PARITY, f_mode, band_min, band_max, kdir, tol, kmag_guess,\n",
    "          kmag_min, kmag_max, mpb.output_poynting_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dispersion relation yields the effective index :\n",
    "\n",
    "$$ n_{eff} = \\frac{c}{v_{phase}} = \\frac{c}{\\omega/k} = \\frac{ck}{\\omega} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.3514721959617"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_fmode[0] / f_mode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is valid at this wavelength. For higher-order derivatives (group velocity, dispersion parameter, etc.), MPB will NOT give physically accurate results by default as commercial solvers because it does not include material dispersion. This can be accounted for, however, see https://github.com/abulnaga1/MPB_Simulations#material_dispersion."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mp [~/home/ELE559/anaconda3/envs/mp/]",
   "language": "python",
   "name": "conda_mp"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}