08. Bent Waveguide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bent waveguides (layout and simulation)\n",
    "\n",
    "Step by step process:\n",
    "1. Draw a bent waveguide and save it in a .gds file.\n",
    "2. Load the bent waveguide using Meep.\n",
    "3. Setup simulation environment.\n",
    "4. Simulate FDTD and visualize results.\n",
    "5. Compute loss and reflection of the bend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using MPI version 3.1, 1 processes\n",
      "Meep version: 1.12.0\n"
     ]
    }
   ],
   "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\n",
    "\n",
    "# Debug info\n",
    "print(\"Meep version:\", mp.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pya\n",
    "import numpy as np\n",
    "\n",
    "SIM_CELL = pya.LayerInfo(0, 0)\n",
    "Si = pya.LayerInfo(1, 0)\n",
    "MEEP_SOURCE1 = pya.LayerInfo(10, 0)\n",
    "MEEP_PORT1 = pya.LayerInfo(20, 0)\n",
    "MEEP_PORT2 = pya.LayerInfo(21, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ring_radius = 8 # um\n",
    "ring_width = 0.5 # um\n",
    "pml_width = 1.0 # um\n",
    "straight_wg_length = pml_width + 0.2 # um\n",
    "\n",
    "# Simulation resolution\n",
    "res = 100        # pixels/μm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1. Drawing a bent waveguide and saving into a temporary .gds file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Produced file /tmp/tmpssrr48zg.gds.\n",
      "-rw-r--r-- 1 tlima grad 2.3K Feb 25 20:50 /tmp/tmpssrr48zg.gds\r\n"
     ]
    }
   ],
   "source": [
    "from zeropdk.layout import layout_arc, layout_waveguide, layout_path, layout_box\n",
    "from tempfile import NamedTemporaryFile\n",
    "\n",
    "# Create a temporary filename\n",
    "temp_file = NamedTemporaryFile(delete=False, suffix='.gds')\n",
    "filename = temp_file.name\n",
    "\n",
    "# Clean existing file\n",
    "!rm {filename}\n",
    "\n",
    "# Instantiate a layout and a top cell\n",
    "layout = pya.Layout()\n",
    "layout.dbu = 0.001\n",
    "TOP = layout.create_cell(\"TOP\")\n",
    "\n",
    "# Unit vectors\n",
    "ex = pya.DVector(1, 0)\n",
    "ey = pya.DVector(0, 1)\n",
    "\n",
    "# Draw circular bend\n",
    "layout_arc(TOP, Si, - ring_radius*ey, ring_radius, ring_width, 0, np.pi/2)\n",
    "\n",
    "# Extend the bend to avoid discontinuities\n",
    "layout_waveguide(TOP, Si, [0*ex, - straight_wg_length*ex], ring_width)\n",
    "layout_waveguide(TOP, Si, [-1*ring_radius*ey + ring_radius*ex, \n",
    "                           -straight_wg_length * ey - ring_radius*ey + ring_radius*ex], ring_width)\n",
    "\n",
    "# Add the ports as 0-width paths\n",
    "port_size = ring_width * 4.0\n",
    "\n",
    "# Source port\n",
    "layout_path(TOP, MEEP_SOURCE1, [-port_size/2*ey - 0.2*ex, port_size/2*ey - 0.2*ex], 0)\n",
    "# Input port (immediately at the start of the bend)\n",
    "layout_path(TOP, MEEP_PORT1,   [-port_size/2*ey, port_size/2*ey], 0)\n",
    "# Output port (immediately at the end of the bend)\n",
    "layout_path(TOP, MEEP_PORT2,   [-1*ring_radius*ey + ring_radius*ex - port_size/2*ex, \n",
    "                               -1*ring_radius*ey + ring_radius*ex + port_size/2*ex], 0)\n",
    "\n",
    "# Draw simulation region\n",
    "layout_box(TOP, SIM_CELL, \n",
    "           -1.0*ring_radius*ey - straight_wg_length * (ex + ey), # Bottom left point \n",
    "            1.0*ring_radius*ex + (straight_wg_length + port_size / 2) * (ex + ey),  # Top right point\n",
    "           ex)\n",
    "\n",
    "# Write to file\n",
    "layout.write(filename)\n",
    "print(f\"Produced file {filename}.\")\n",
    "!ls -lh {filename}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2. Load gds file into meep\n",
    "\n",
    "### Visualization and simulation\n",
    "\n",
    "If you choose a normal filename (not temporary), you can download the GDSII file from the cluster (see Files in MyAdroit dashboard) to see it with your local Klayout. Otherwise, let's get simulating:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdsII_file = filename\n",
    "CELL_LAYER = 0\n",
    "SOURCE_LAYER = 10\n",
    "Si_LAYER = 1\n",
    "PORT1_LAYER = 20\n",
    "PORT2_LAYER = 21\n",
    "\n",
    "t_oxide = 1.0\n",
    "t_Si = 0.22\n",
    "t_SiO2 = 0.78\n",
    "\n",
    "oxide = mp.Medium(epsilon=2.25)\n",
    "silicon=mp.Medium(epsilon=12)\n",
    "\n",
    "lcen = 1.55\n",
    "fcen = 1/lcen\n",
    "df = 0.2*fcen\n",
    "nfreq = 25\n",
    "\n",
    "cell_zmax =  0\n",
    "cell_zmin =  0\n",
    "si_zmax = 10\n",
    "si_zmin = -10\n",
    "\n",
    "# read cell size, volumes for source region and flux monitors,\n",
    "# and coupler geometry from GDSII file\n",
    "# WARNING: Once the file is loaded, the prism contents is cached and cannot be reloaded.\n",
    "# SOLUTION: Use a different filename or restart the kernel\n",
    "\n",
    "si_layer = mp.get_GDSII_prisms(silicon, gdsII_file, Si_LAYER, si_zmin, si_zmax)\n",
    "\n",
    "cell = mp.GDSII_vol(gdsII_file, CELL_LAYER, cell_zmin, cell_zmax)\n",
    "src_vol = mp.GDSII_vol(gdsII_file, SOURCE_LAYER, si_zmin, si_zmax)\n",
    "p1 = mp.GDSII_vol(gdsII_file, PORT1_LAYER, si_zmin, si_zmax)\n",
    "p2 = mp.GDSII_vol(gdsII_file, PORT2_LAYER, si_zmin, si_zmax)\n",
    "\n",
    "\n",
    "sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),\n",
    "                              size=src_vol.size,\n",
    "                              center=src_vol.center,\n",
    "                              eig_band=1,\n",
    "                              eig_parity=mp.NO_PARITY,\n",
    "                              eig_match_freq=True)]\n",
    "\n",
    "# Display simulation object\n",
    "sim = mp.Simulation(resolution=res,\n",
    "                    default_material=oxide,\n",
    "                    eps_averaging=False,\n",
    "                    cell_size=cell.size,\n",
    "                    boundary_layers=[mp.PML(pml_width)],\n",
    "                    sources=sources,\n",
    "                    geometry=si_layer,\n",
    "                    geometry_center=cell.center)\n",
    "\n",
    "# Delete file created in previous cell\n",
    "\n",
    "import os \n",
    "temp_file.close()\n",
    "os.unlink(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3. Setup simulation environment\n",
    "\n",
    "This will load the python-defined parameters from the previous cell and instantiate a fast, C++ based, simulation environment using meep. It will also compute the eigenmode of the source, in preparation for the FDTD simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000456095 s\n",
      "Working in 2D dimensions.\n",
      "Computational cell is 11.4 x 11.4 x 0 with resolution 100\n",
      "     prism, center = (5.11221,-2.88779,0)\n",
      "          height 20, axis (0,0,1), 228 vertices:\n",
      "          (7.75,-8,-10)\n",
      "          (7.748,-7.81,-10)\n",
      "          (7.745,-7.715,-10)\n",
      "          (7.741,-7.62,-10)\n",
      "          (7.729,-7.43,-10)\n",
      "          (7.713,-7.24,-10)\n",
      "          (7.703,-7.146,-10)\n",
      "          (7.692,-7.051,-10)\n",
      "          (7.666,-6.863,-10)\n",
      "          (7.652,-6.769,-10)\n",
      "          (7.636,-6.675,-10)\n",
      "          (7.619,-6.581,-10)\n",
      "          (7.601,-6.488,-10)\n",
      "          (7.582,-6.395,-10)\n",
      "          (7.562,-6.302,-10)\n",
      "          (7.54,-6.209,-10)\n",
      "          (7.518,-6.117,-10)\n",
      "          (7.494,-6.025,-10)\n",
      "          (7.469,-5.933,-10)\n",
      "          (7.443,-5.841,-10)\n",
      "          (7.416,-5.75,-10)\n",
      "          (7.388,-5.659,-10)\n",
      "          (7.359,-5.569,-10)\n",
      "          (7.297,-5.389,-10)\n",
      "          (7.231,-5.211,-10)\n",
      "          (7.196,-5.122,-10)\n",
      "          (7.16,-5.034,-10)\n",
      "          (7.123,-4.947,-10)\n",
      "          (7.085,-4.859,-10)\n",
      "          (7.046,-4.773,-10)\n",
      "          (7.006,-4.686,-10)\n",
      "          (6.965,-4.601,-10)\n",
      "          (6.879,-4.431,-10)\n",
      "          (6.835,-4.347,-10)\n",
      "          (6.79,-4.263,-10)\n",
      "          (6.743,-4.18,-10)\n",
      "          (6.696,-4.098,-10)\n",
      "          (6.598,-3.934,-10)\n",
      "          (6.548,-3.854,-10)\n",
      "          (6.444,-3.694,-10)\n",
      "          (6.391,-3.616,-10)\n",
      "          (6.336,-3.537,-10)\n",
      "          (6.281,-3.46,-10)\n",
      "          (6.225,-3.383,-10)\n",
      "          (6.168,-3.307,-10)\n",
      "          (6.11,-3.232,-10)\n",
      "          (6.051,-3.157,-10)\n",
      "          (5.991,-3.083,-10)\n",
      "          (5.93,-3.01,-10)\n",
      "          (5.806,-2.866,-10)\n",
      "          (5.742,-2.795,-10)\n",
      "          (5.678,-2.725,-10)\n",
      "          (5.613,-2.656,-10)\n",
      "          (5.547,-2.588,-10)\n",
      "          (5.48,-2.52,-10)\n",
      "          (5.412,-2.453,-10)\n",
      "          (5.344,-2.387,-10)\n",
      "          (5.275,-2.322,-10)\n",
      "          (5.205,-2.258,-10)\n",
      "          (5.134,-2.194,-10)\n",
      "          (4.99,-2.07,-10)\n",
      "          (4.917,-2.009,-10)\n",
      "          (4.843,-1.949,-10)\n",
      "          (4.768,-1.89,-10)\n",
      "          (4.693,-1.832,-10)\n",
      "          (4.617,-1.775,-10)\n",
      "          (4.54,-1.719,-10)\n",
      "          (4.463,-1.664,-10)\n",
      "          (4.384,-1.609,-10)\n",
      "          (4.306,-1.556,-10)\n",
      "          (4.146,-1.452,-10)\n",
      "          (4.066,-1.402,-10)\n",
      "          (3.902,-1.304,-10)\n",
      "          (3.82,-1.257,-10)\n",
      "          (3.737,-1.21,-10)\n",
      "          (3.653,-1.165,-10)\n",
      "          (3.569,-1.121,-10)\n",
      "          (3.399,-1.035,-10)\n",
      "          (3.314,-0.994,-10)\n",
      "          (3.227,-0.954,-10)\n",
      "          (3.141,-0.915,-10)\n",
      "          (3.053,-0.877,-10)\n",
      "          (2.966,-0.84,-10)\n",
      "          (2.878,-0.804,-10)\n",
      "          (2.789,-0.769,-10)\n",
      "          (2.611,-0.703,-10)\n",
      "          (2.431,-0.641,-10)\n",
      "          (2.341,-0.612,-10)\n",
      "          (2.25,-0.584,-10)\n",
      "          (2.159,-0.557,-10)\n",
      "          (2.067,-0.531,-10)\n",
      "          (1.975,-0.506,-10)\n",
      "          (1.883,-0.482,-10)\n",
      "          (1.791,-0.46,-10)\n",
      "          (1.698,-0.438,-10)\n",
      "          (1.605,-0.418,-10)\n",
      "          (1.512,-0.399,-10)\n",
      "          (1.419,-0.381,-10)\n",
      "          (1.325,-0.364,-10)\n",
      "          (1.231,-0.348,-10)\n",
      "          (1.137,-0.334,-10)\n",
      "          (0.949,-0.308,-10)\n",
      "          (0.854,-0.297,-10)\n",
      "          (0.76,-0.287,-10)\n",
      "          (0.57,-0.271,-10)\n",
      "          (0.38,-0.259,-10)\n",
      "          (0.285,-0.255,-10)\n",
      "          (0.19,-0.252,-10)\n",
      "          (0,-0.25,-10)\n",
      "          (0,0.25,-10)\n",
      "          (0.202,0.248,-10)\n",
      "          (0.304,0.244,-10)\n",
      "          (0.405,0.24,-10)\n",
      "          (0.607,0.228,-10)\n",
      "          (0.708,0.22,-10)\n",
      "          (0.809,0.21,-10)\n",
      "          (0.909,0.2,-10)\n",
      "          (1.01,0.188,-10)\n",
      "          (1.11,0.175,-10)\n",
      "          (1.211,0.161,-10)\n",
      "          (1.311,0.145,-10)\n",
      "          (1.41,0.129,-10)\n",
      "          (1.51,0.111,-10)\n",
      "          (1.609,0.091,-10)\n",
      "          (1.709,0.071,-10)\n",
      "          (1.808,0.05,-10)\n",
      "          (1.906,0.027,-10)\n",
      "          (2.005,0.003,-10)\n",
      "          (2.103,-0.022,-10)\n",
      "          (2.2,-0.049,-10)\n",
      "          (2.298,-0.076,-10)\n",
      "          (2.395,-0.105,-10)\n",
      "          (2.492,-0.135,-10)\n",
      "          (2.588,-0.166,-10)\n",
      "          (2.684,-0.199,-10)\n",
      "          (2.779,-0.232,-10)\n",
      "          (2.874,-0.267,-10)\n",
      "          (2.969,-0.303,-10)\n",
      "          (3.063,-0.34,-10)\n",
      "          (3.157,-0.378,-10)\n",
      "          (3.25,-0.417,-10)\n",
      "          (3.436,-0.499,-10)\n",
      "          (3.527,-0.542,-10)\n",
      "          (3.619,-0.586,-10)\n",
      "          (3.709,-0.631,-10)\n",
      "          (3.799,-0.677,-10)\n",
      "          (3.889,-0.724,-10)\n",
      "          (3.978,-0.772,-10)\n",
      "          (4.154,-0.872,-10)\n",
      "          (4.328,-0.976,-10)\n",
      "          (4.414,-1.03,-10)\n",
      "          (4.499,-1.085,-10)\n",
      "          (4.583,-1.14,-10)\n",
      "          (4.667,-1.197,-10)\n",
      "          (4.75,-1.255,-10)\n",
      "          (4.833,-1.314,-10)\n",
      "          (4.915,-1.374,-10)\n",
      "          (4.995,-1.434,-10)\n",
      "          (5.076,-1.496,-10)\n",
      "          (5.155,-1.559,-10)\n",
      "          (5.234,-1.623,-10)\n",
      "          (5.312,-1.687,-10)\n",
      "          (5.389,-1.753,-10)\n",
      "          (5.465,-1.82,-10)\n",
      "          (5.54,-1.887,-10)\n",
      "          (5.615,-1.956,-10)\n",
      "          (5.689,-2.025,-10)\n",
      "          (5.762,-2.095,-10)\n",
      "          (5.834,-2.166,-10)\n",
      "          (5.905,-2.238,-10)\n",
      "          (5.975,-2.311,-10)\n",
      "          (6.044,-2.385,-10)\n",
      "          (6.113,-2.46,-10)\n",
      "          (6.18,-2.535,-10)\n",
      "          (6.247,-2.611,-10)\n",
      "          (6.313,-2.688,-10)\n",
      "          (6.377,-2.766,-10)\n",
      "          (6.441,-2.845,-10)\n",
      "          (6.504,-2.924,-10)\n",
      "          (6.566,-3.005,-10)\n",
      "          (6.626,-3.085,-10)\n",
      "          (6.686,-3.167,-10)\n",
      "          (6.745,-3.25,-10)\n",
      "          (6.803,-3.333,-10)\n",
      "          (6.86,-3.417,-10)\n",
      "          (6.915,-3.501,-10)\n",
      "          (6.97,-3.586,-10)\n",
      "          (7.024,-3.672,-10)\n",
      "          (7.128,-3.846,-10)\n",
      "          (7.228,-4.022,-10)\n",
      "          (7.276,-4.111,-10)\n",
      "          (7.323,-4.201,-10)\n",
      "          (7.369,-4.291,-10)\n",
      "          (7.414,-4.381,-10)\n",
      "          (7.458,-4.473,-10)\n",
      "          (7.501,-4.564,-10)\n",
      "          (7.583,-4.75,-10)\n",
      "          (7.622,-4.843,-10)\n",
      "          (7.66,-4.937,-10)\n",
      "          (7.697,-5.031,-10)\n",
      "          (7.733,-5.126,-10)\n",
      "          (7.768,-5.221,-10)\n",
      "          (7.801,-5.316,-10)\n",
      "          (7.834,-5.412,-10)\n",
      "          (7.865,-5.508,-10)\n",
      "          (7.895,-5.605,-10)\n",
      "          (7.924,-5.702,-10)\n",
      "          (7.951,-5.8,-10)\n",
      "          (7.978,-5.897,-10)\n",
      "          (8.003,-5.995,-10)\n",
      "          (8.027,-6.094,-10)\n",
      "          (8.05,-6.192,-10)\n",
      "          (8.071,-6.291,-10)\n",
      "          (8.091,-6.391,-10)\n",
      "          (8.111,-6.49,-10)\n",
      "          (8.129,-6.59,-10)\n",
      "          (8.145,-6.689,-10)\n",
      "          (8.161,-6.789,-10)\n",
      "          (8.175,-6.89,-10)\n",
      "          (8.188,-6.99,-10)\n",
      "          (8.2,-7.091,-10)\n",
      "          (8.21,-7.191,-10)\n",
      "          (8.22,-7.292,-10)\n",
      "          (8.228,-7.393,-10)\n",
      "          (8.24,-7.595,-10)\n",
      "          (8.244,-7.696,-10)\n",
      "          (8.248,-7.798,-10)\n",
      "          (8.25,-8,-10)\n",
      "          dielectric constant epsilon diagonal = (12,12,12)\n",
      "     prism, center = (-0.6,0,0)\n",
      "          height 20, axis (0,0,1), 4 vertices:\n",
      "          (-1.2,-0.25,-10)\n",
      "          (-1.2,0.25,-10)\n",
      "          (0,0.25,-10)\n",
      "          (0,-0.25,-10)\n",
      "          dielectric constant epsilon diagonal = (12,12,12)\n",
      "     prism, center = (8,-8.6,0)\n",
      "          height 20, axis (0,0,1), 4 vertices:\n",
      "          (7.75,-9.2,-10)\n",
      "          (7.75,-8,-10)\n",
      "          (8.25,-8,-10)\n",
      "          (8.25,-9.2,-10)\n",
      "          dielectric constant epsilon diagonal = (12,12,12)\n",
      "subpixel-averaging is 84.1535% done, 0.753601 s remaining\n",
      "subpixel-averaging is 84.3803% done, 0.741163 s remaining\n",
      "subpixel-averaging is 84.3803% done, 0.741198 s remaining\n",
      "time for set_epsilon = 13.0522 s\n",
      "-----------\n",
      "MPB solved for omega_1(2.2349,0,0) = 0.682374 after 23 iters\n",
      "MPB solved for omega_1(2.10209,0,0) = 0.645189 after 9 iters\n",
      "MPB solved for omega_1(2.10199,0,0) = 0.645161 after 4 iters\n",
      "MPB solved for omega_1(2.10199,0,0) = 0.645161 after 1 iters\n"
     ]
    }
   ],
   "source": [
    "sim.reset_meep()\n",
    "\n",
    "# Could add monitors at many frequencies by looping over fcen\n",
    "# Means one FDTD for many results!\n",
    "mode1 = sim.add_mode_monitor(fcen, df, nfreq, mp.ModeRegion(volume=p1))\n",
    "mode2 = sim.add_mode_monitor(fcen, df, nfreq, mp.ModeRegion(volume=p2))\n",
    "\n",
    "# Let's store the frequencies that were generated by this mode monitor\n",
    "mode1_freqs = np.array(mp.get_eigenmode_freqs(mode1))\n",
    "mode2_freqs = np.array(mp.get_eigenmode_freqs(mode2))\n",
    "\n",
    "sim.init_sim()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verify that the structure makes sense.\n",
    "\n",
    "Things to check:\n",
    "- Are the sources and ports outside the PML?\n",
    "- Are dimensions correct?\n",
    "- Is the simulation region unnecessarily large?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAFtCAYAAAD8jSoDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3wU1f3/8dfZBIJgAkFAEEQQUCI35VaKXESlglJQ+Klf+xVRQC14qaWooIAgWGiLYCtqVS4K1WqtWEQEyldRoBTCHQIhAgJCuAshIZiEJOf3x2Z2N4EkO7s7O7O7n+fjsY+Yzc7MWY5578lnzpxRWmuEEEJEF5fdDRBCCBF6Eu5CCBGFJNyFECIKSbgLIUQUknAXQogoJOEuhBBRSMJdCCGikIS7EEJEoXi7GxBOSikFXAXk2N0WIYQwIRE4ok1cdRpT4Y472A/b3QghhAhAIyDT3xdHTLgrpcYCA4GWwE/AWuB5rXWGid3kAAx+bzDdr+se+kaWOJh1kFdWvcKA6wewKGMRw9oP46F2D1l2vIqkn0xn1PJRNE1uyvTe06letXrY23C+4DyjV4xm/5n9zLhjBil1U8LeBoD52+YzZ/Mc6Q/pD49A+yPjxwwe+/wxXuzxItfUuiaoNmT9lMX4lePJL8xnyI1DuLnxzaV+vvq71Sx4eAGYrThorSPiASwDHgZaAe2AL4CDQA0T+0gC9Dtr39FW2nRkk2Yimonoyd9OtvRYFVl/eL1Ompqku87pqrPzsm1pQ3Zetu46p6tOmpqk1x9eb0sbtNZ68reTpT+09IevYPrD+B3fdGRTUG3IPJupE3+fqJmInrt57iVf887adzSggSRtJjPNvNhJD6BuyRvuYWKbsIR753c6ayaiR3wxwtLjVESCxCvSgyRUpD+8gu2PUIS7P8GudeDhHjFlmUuoWfL1dHkvUEolAAk+TyVa2iLg9vm3k3okFYDh7YdbfbhLSs1MpfeC3rSu15pl/7uMxATL3/ZFcvJz6PNBH9JOpLFi8Ao6N+wc9jYATFk1hfErxzO512TG9RhnSxukP7ykP9yOZB+h5RstySnIYW7/uTxy0yMhP0ZEToUsmfUyA1ijtU6r4KVjgbM+D0tPpt4+/3a+2v8Vna+y5xcHnPE/rgSJl/SHl/SHWziCHYjMsgzwBnAAaFTJ6xJwl2KMR0MsKsvc9v5tmono296/LWT1OLPkT3+vaPjTPxSkP7xC2R+B/o77W4rxFTM1d+B14BDQNIBtLam5+wa71qE72WKGBIlXtAVJoKQ/vELdH4H8jgcS7FrHQM29pBTzOnAPcIvWer/NTQK8pZjbmt7G/z30f7a0wQl/asqf/l7SH17SH25hK8X4iJhwx12K+RUwAMhRStUvef6s1vonOxokwe4mQeIl/eEl/eFmR7ADkVOWwf1nyaUeD5vYR8jKMmVLMb7CVZaRP/29ovFP/0BIf3hZ2R/+/o4HWorxFfVlGa21srsNBhmxu8kI0Uv6w0v6w822EXuJiJwKaScJdjcJEi/pDy/pDze7gx0k3E2RYHeTIPGS/vCS/nBzQrCDhLvfJNjdJEi8pD+8pD/cnBLsIOHuFwl2NwkSL+kPL+kPNycFO0i4V0qC3U2CxEv6w0v6w81pwQ4S7hWSYHeTIPGS/vCS/nBzYrCDhHu5JNjdJEi8pD+8pD/cTuaedGSwA5FzEVMoHvh5EVNFFyj5IxQXMckFMV7RfkGMv6Q/vOzuD+N3vPor1YO6QMkfgV7EJCP3MmTE7iYjRC/pDy/pD7fcglwAzl8477wRu8HMJ0GkP6hk5B7siN0QzMjd7hGJ1jJC9CX94SX94Zadl63bvdVOMxE9YeUEy48nI/cgyYjdTUaIXtIfXtIfbkZ/7D29F4AB1w8Iexv8JeGOBLtBgsRL+sNL+sPNtz/evOvNsB/frIhZOMwqEuxuEiRe0h9e0h9uZfsj3uX86HR+Cy1wMOsgm49u5tdf/JoNRzbQ6apO/LH3H9l8dHNI9p9+Kr3U14qknUhj5JKRNKvdjGm3TWPP6T0haYMZuQW5PLn0Sfad3sebd71JvCs+ZP8WZszePJu3Nr7FiI4juLPFnba0QfrDS/rD7VL9YeZ3PFgHsw4GtJ3S7hONMUEplQScZQxQze7WCCGEH/KAaQDU1Fpn+7tZTI7cB1w/gEUHFzGi4wiGtx8e8v2nn0rnwYUP8reBfyOlTkrI9y+EsFc4f8dfXv4yi1hkeruYDPdFGYuY3Nf6GmJKnRTaN2hv6TGEEPax+nd8yqopLMowH+wQo7NlhrUfZtvJISGE8IdxMntY+2EBbR+T4f5Qu4fsboIQQpTLd5ZSoHkVk+Futev7PMihV91ffXXsCI0aub8KIcSlhGr6aUzW3K1W5eSPNMqBgpM/lnr+2DHIzLSpUUIIxwvldQUychdCCAcI9QVjEu5CCGEzK64ElnAXQggbWbXEg4S7BS4UXbC7CUKICGDl2j0S7iGWmplKVl6W3c0QQjic1YuySbiHkLF6XSSsGCeEsE84VtuUcA8R32VJa1WrZXdzhBAOFa5llCXcQ6DsetNKKbubJIRwoHCujy/hHiQn3EhACOF84b7xiYR7ECTYhRD+sOOOVhLuAZJgF0L4w65bFUq4B0CCXQjhDzvvQSvhbpIEuxDCH3bfXFzC3QQJdiGEP+wOdpBw95sEuxDCH04IdpBw94sEuxDCH04JdojRm3UUFxdTXFzs12tTM1O544M7aF23NV8+8CU1qtSofFvt/Vr6tarkoSku1hdtVvbiJ7kYSojI4aRghxgNd5fLhctV+R8tnmA3O2JX3q+XPo7C5TIf3MYHhdbuDwYj/H0/BOQDQYjwc1qwQ4yG+/fff0/tgtqe77XWF4XijtM7eGz1YzSv2ZwZN83gwHcHKg1OYz/XFhZSFSgsLOS7tDTP8xcuXAdU5cKFC+zYsduzP+Pnl112GS6Xi6KiIurUqUNSUhJKKc/r/PlAMvZnPMp+AEj4i1hQWFhIfHx44s2JwQ4xGu6PPPwIHK/gBQ2BwcAJ2Pq3rXQp6GJq/4dcUB04/eNp2rRp4/sToBEnThynbdu2F21XpUoVlFIUFhZSr149kpOTcblcVKlShYSEBG644QaqVKlC1apVady4MU2bNiUhIYFmzZpRs2ZNqlevTnJycqkPhEspLi6+KPgl9EW0KC4uJj8/PyzhPnvzbN7a+Jbjgh1iNNwr5BPs/A0oCN+hL1zw3uTj2LFjHDt2rNTP169ff8nt4uPjqVKlComJiTRu3Jjq1avTtm1batWqRZs2bbjyyitp0qQJtWvXJjExsdy/AHxDXwJfRKpw/r/r1GCHCAx3pdRI4FmgAbATeEZrvTokO7cx2MtTXj3dqLuD+0/QwsJCfvrpJ06cOAHAqlWrSu2nRo0aJCUl0bx5c1q0aEGLFi1o2bIl1157LU2aNCEpKemSoS+BL0T5RnQc4chghwgLd6XU/cBrwEjgP8DjwFKl1A1a6x+C2rkDgx1Kh7jvf5enbG3dqL3n5uaSm5vL0aNHWb269Gdh/fr1ad68OS1btqRDhw60a9eO5s2bU7du3YsCX2tNcXGxhL0QwPD2w+1uQrkiKtyBUcAcrfXsku+fUUrdAYwAxga8V4cGeyCMD4DyPgh8A9mYfWOUgNasWcPs2e5/2uTkZFq0aEHbtm3p0KEDbdu2JSUlheTkZOLi4krts6ioSMJeCIeJmHBXSlUFOgDTyvzo30DXcrZJABJ8nrp4LmMUBbs/jJF8WWVD/8yZM6SmppKamuoJ/CuvvJKUlBS6dOlCjx49uPHGG2nQoEGpsPcd2fs7u0cIEXoRE+5AHSCOi+e5HAfql7PNWOClcvcYY8FekUuFvhH4xs+OHz/O8ePH+eabb5g2bRpJSUnccMMNdO3alZ49e9KhQwcaNmxYKuyNmr3L5ZJRvRBhFEnhbig77FSXeM4wFZjh830icBiQYPdDRYFfXFxMdnY269atY926dcyYMYPLL7+cTp060bt3b3r16sWNN95ItWrVSu1PRvVChEckhfspoIiLR+n1KGfWutY6H8g3vveMHOsAv0SCPQBlA9837M+dO8fKlStZuXIlAE2aNOFnP/sZ/fr145ZbbqFRo0alRvVGrV6CXojypZ9MD2i7iAl3rXWBUmoT0Bv4zOdHvYFFpnZ2BxLsIVJR2B84cIADBw7w8ccfc/nll9OxY0cGDBjAbbfdRps2bSTohahEamYqo5aPCmjbiAn3EjOABUqpjcB/gceAxsBfTe3lDPABEuwWuFTYA5w7d45vvvmGb775BoBOnTrRv39/+vfvT9u2bS8KeqnRi1hnrEbbNLkpO9lpevuIGiZprT8GngEmAFuBHsCdWuuDpnb0byIu2I1RrfGIlKmHvmHvOzLfsGED48ePp127dnTu3JkpU6awfft2AOLi4jwnco0TskLEEt9lxqf3nh7QPiJt5I7W+k3gzWD2USe5DtXjqlsWGq5jmVBcjCvOxdUNG3qeP3IkjqIid3hdddXVpbZRSlFQUEBxcTEul4sff/yx1HIE5U1hvBTfWS7Gtk5QNuiN8s2GDRs8Yd+9e3cGDRrEgAEDaNKkyUVz8qVsI6Jd2ftHbDm4JaD9RFy4h8KHH3xI9+u6W7Z/1bQxHDvJFVdcwXfffed5vlmzBI4cgSuvrF/qecOFCxcoLCwkLi6OzMxMsrOzPYF/6NAhDh06RGFhIadPn2bXrl2e/z5w4IBn9gpU/EFQdnqjXcoL+tWrV7N69WrGjBnDbbfdxsMPP0zfvn2pUaOGZ1sp24hoFcobA8VkuCdUSyg1RS/UCnyWAEjwOY6RRUqpSx7f97mkpCS/jqW15syZM+Tl5bF3715yc3NJT0/n8OHD7Nmzh8zMTA4ePEhOTg4XLlzwaz57uF0q6PPy8liyZAlLlizh6quv5u677+bhhx+mffv2nvq8nIQV0STUd3yLyXB3qrLBWt73vrV2pRS1a7vXpr/qqqsA6Nu3b6ntTp48yenTp9mzZw8HDhwgLS2NjIwM9u3bx6FDhxwV+GWDHuDQoUO8/vrrvP7663Tr1o2hQ4dyzz33UKtWrVLbSMiLSGXFrTwl3B0k0Nvsla2tG1+N0kXdunWpW7cu119/fantcnNz+eGHH9i1axc7d+5k69at7Ny5k4MHD5Kfn3/JwPfdv9V8j+NyuSguLmbNmjWsWbOGF154gXvvvZfhw4fTtm1bT9ukZCMijVX3aJZwt0C4SxuV3WXJN/SN/46Li6NGjRqkpKSQkpLCoEGDAPfywT/88APbtm1j8+bNbNq0iS1btnDs2DFbw944oWoc89ixY57RfJ8+fXjiiSe48847S5VsJOSF01kV7CDhHnI5+Tmcz8viSrsb4qO88C8b+kop4uPjufbaa7n22mu55557ADh79iw7duxgy5YtrFq1ig0bNnDw4MGLQt3lcllexrnUaH7ZsmUsW7aMdu3a8fjjj3P//fd7SlW+yxML4SRWBjtE2Dx3p8vJz6HPB30oLC60uyl+MULP5XIRFxfnqVkb88uLiorQWlOzZk26devGU089xSeffEJ6ejobN27kL3/5C/369aN+ffeKEL5z0sMRqGVH89u2bWPkyJGkpKQwduxYDh8+7Bm9FxcXe14vhN2sDnaQcA8ZI9jTTqRRq1otu5sTFN/A972YqKioiOLiYi677DI6dOjAU089xeLFi9m9ezdff/01Y8aMoVOnTlSpUqXci5esUPY4J06cYNq0aVx33XUMHTqU3bt3ey7+kpAXdgtHsIOEe0j4BvuKwSuoElfF7iaF1KVG92VH9r169WLq1Kmkpqayc+dO3n33Xe6++25q167t+XAwWBX0vsdxuVz89NNPzJs3jzZt2vDggw+yefNmCXlhq3AFO0i4B61ssHdu2NnuJoVFeSN7gBYtWjB8+HA+++wzMjIy+PTTT3n00Ue5+mr3Vbm+5RSrSje+IV9YWMgHH3xAhw4d+J//+R+2bNkiIS/CLpzBDhLuQYnVYC/Ld2QP7hG0UcKpU6cOAwcO5J133iE9PZ2lS5fy2GOP0ahRo0tevBRqZZct+Pjjj2nfvv0lQ94pyzSI6BPuYAcJ94BJsJdPKeUp4Rij+uLiYmrUqEGfPn14++232bVrF//85z954IEHqFWrVqmgt6JsU17IDxkyxFOTV0p5Sk1ChIodwQ4S7gGRYPef72qWvkGfmJjIoEGD+PDDD8nIyGDu3LncfvvtxMfHW1qfLxvy8+fPp02bNgwbNoxDhw55Sk1GiUmIYNgV7CDhbpoEe+DKC/p69erxyCOPsGLFCrZv386kSZNo2bIlYN1qkGVr8nPnzqVVq1a88MILZGVleUpMUo8XgbIz2EHC3RQJ9tApG/RGOSQlJYUJEyawY8cOli5dyv3330+1atUsOwnrG/I5OTlMnTqVG264gdmzZ3ueLzvbR4jK2B3sIOHuNwl26xg1emPmTVFREfHx8fTp04ePPvqI9PR0Jk2axDXXXGNZbd435I8ePcqjjz7KTTfdxLJlyzwfRFKPF/5wQrCDhLtfJNjDxwh6wFO2adKkCRMmTGD37t0sWLCA7t27e34O1oS8UoqtW7fSt29ffvnLX/Ldd99JPV5UyinBDhLulZJgt0/Zsk21atV48MEHWbVqFatXr+bBBx+katWqloS8sdYOwBdffEHr1q156aWXyMvLIy4uTko14iJOCnaQcK+QBLsz+I7mjdJIt27dWLBgAWlpaYwePZq6deuGPOR9yz8XLlzg5ZdfplWrVixcuLBUqUYIpwU7SLiXS4LdmYzSiFGyadGiBX/605/YuXMn06ZNo2HDhiEPed/9ff/99wwaNIj+/fuTkZEhs2qEI4MdJNwvKdhgv1D3Cg4nur/6ql8fGjZ0fxXBKVuyqVu3Ls8//zy7du3ij3/8Iw0aNLAk5I3ZOosXL6Zt27a8+uqrnmPICdfYk3YizZHBDhLuFwnFiD1j2d+4+nfur742boTDh91fRWgYJRsj5JOSknj22WfJyMiwJOR9b+lXUFDA6NGj6dKlC+vWrSv1V4WIDSOXjHRksIOEeylSiolcZUM+MTGRZ599lt27d1sS8r77Wr9+PT//+c8ZN24chYWFUouPIc1qN3NksIOEu4cEe3SoaCT/yiuvkJycHPKQN2bVvPLKK9x0002sXbtWavExYlbfWY4MdpBwByTYo9GlRvIvvPACe/bsYcyYMaWmUAZ7xavvrJq0tDRuvvlmXnjhBYqLi2UUH+VqVK1hdxPKFfPhLsEe3cqG/BVXXMHUqVPZuXMngwcPBrxz2oMNed9R/NSpU+nQoQNpaWkyihe2iOlwl2CPHWVDvnnz5syfP5/Vq1dz6623ljpRGgzfUfzWrVvp0KEDf/7znz3PyShehEvMhrsEe2zyDfni4mK6devGV199xccff0zz5s1DVo83RvEFBQU888wz9OnTh6NHjxIXFycjeBEWMRnu5wvOS7DHOOMKU+MOTPfddx/bt29n/PjxVKlSpVSJJVC+o/jly5fTqlUrFi9eXOo+tEJU5nzB+YC2i8lwH71itAS7ACh1B6bLLruMl19+me3bt9OvX7+QrT5pnFg9c+YM/fv3Z9SoUZ79SplGVCQnP4fRK0YHtG1Mhvv+M/sl2EUpvvX4li1bsnjxYv75z3/SuHHjkMyq8R2lz5w5k5///OccOHCAuLg4ubJVXJJROt5/Zn9A28dkuM+4Y4YEu7iIUY83SjWDBg1i165d/Pa3vwUIyQlXcI/Y161bR7t27fj000/lylZxEd9zgjPumBHQPmIy3FPqptjdBOFgvqWaGjVqMGPGDL799ltSUlJCNop3uVxkZ2fz//7f/5MyjSil7GSPQPMqJsNdCH/4lmp69OjBjh07mDx5sueOUcGM4suWaW699VZOnjzpKdOI2BTKWXwS7kJUwLdUExcXx7hx49iyZQs/+9nPQnaFa3x8PCtXruSmm25i/fr1pUpDInaEenq2hLsQfjBG6UVFRbRr145169YxefJkIPhavLHYWGZmJl26dGHevHme0pDU4WODFdfdSLgLYYLvRUjjxo1j9erVtGzZMugQ9p1XP3ToUKnDxxCrLqiUcBfCJN9RfLdu3di+fTtPPfXURT83y7cMM3PmTO68807OnTsndfgoZuWV8hLuFpi9ebbdTRBhYIRulSpV+Mtf/sLChQs993INdsqky+Vi6dKldOzYkb179xIXF0dhYWGIWi6cwOolUCTcQ2zKqim8tfEtu5shwsR3nZp77rmHbdu20bt376BPthYXFxMfH09GRgadOnVi9erVxMfHywg+SoRjbSsJ9xCasmoK41eOZ0THEXY3RYSRsU5NUVERDRo04N///jdTp04FgjvZapxozcrKokePHnzyyScykyYKhGvRQgn3EDGCfXKvyQxvP9zu5ggb+J5sHTNmDCtXrqR+/fpBlWl8T7Ted999zJo1yzOTRgI+8oRzNdqICHelVBOl1Byl1H6l1E9KqX1KqUlKqap2tw1KB/u4HuPsbo6wke/J1ltuuYUtW7bQq1evoGbT+Ib4U089xdixYwFkqmSECfcy4xER7kBL3G19HGgF/Bb4NfB7OxsFEuzi0oyTrfXr1+frr7/m+eef9/wsmJOtSimmTZvGsGHDPPuSgHc+O+4fERHhrrVeprV+RGv9b63191rrz4HpwEA72yXBLiriW6aZNm0aH3zwAQkJCUGVaYwa/ty5c+nXrx9FRUUyF97h7LoxUESEezlqAqcreoFSKkEplWQ8gJDdplyCXfjD98Ycv/rVr9iwYQNNmzYNug7vcrlYsmQJt9xyC2fPnpW58A5l5x3fIjLclVLNgKeAv1by0rHAWZ/H4VAcX4JdmGWMrtu0acOWLVvo06dPUHd7MqZKrlmzhp49e5KdnS0B7zB238rT1nBXSk1USulKHh3LbHMVsAz4RGtd2dVCU3GP8I1Ho2DbLMEuAmWEb82aNVm6dClPPPFEUDNeCgsLiY+PZ9u2bXTr1o2DBw9KwDuE3cEOEB/2I5Y2C/ioktccMP6jJNhXAv8FHqts51rrfCDfZ/uAGmmQYBfBMurwLpeLWbNm0bhxY8/J1kCmNxoBv2PHDrp06cK3337LddddR1FREXFxcVa8BVEJJwQ72BzuWutTwCl/XquUaog72DcBj2itwzpFQIJdhIrL5UJrjdaa5557jmbNmnHvvfd6Tpaanf1iBPyxY8fo2bOnBLyNnBLsECE195IR+zfAIWA0UFcpVV8pVT8cx5dgF6Hme1XroEGDWLt2LVdeeaWnlm6WcTWrEfDfffedlGjCzEnBDhES7sAvgObArbhPih71eVhKgl1YyVgQrEuXLqSmppKSkuIZiZtllHsk4MPPacEOERLuWuv3tNbqUg8rjyvBLsIhPj6ewsJCGjduzNq1a7nllls8I3GzJODDz4nBDhES7naQYBfhZKz4WKtWLb7++mv69u0b8Fz4sgEvs2isk1uQ68hgBwn3S5JgF3YwZtIopfjyyy956KGHAp4Lb9Tujx07Rr9+/WQevEWeXPqkI4MdJNwvIsEu7OQ7W+b999/nN7/5TcBz4Y3afVpaGj169PBcySpr0YTOvtP7HBnsIOFeigS7cAJjqiTAa6+9xksvvRTwvnwvdPJdi0YCPjTevOtNRwY7SLh7SLALJzEuaNJaM3HixJAE/Jo1a7jrrruA0h8gInCt67W2uwnlknBHgl04k1FrD1XAu1wuli9fztChQz37l4CPXjEf7hLswsmUUp4QDjbgjVk08+bN89zwQ8I9esV0uEuwi0gSqoA3bvjx+uuvy1rwUSxmw12CXUSSUJZojNH6008/zT/+8Q/PVbIiusRkuM/fNl+CXUScUAa84f7772f16tWei6iE88zfNj+g7WIy3OdsniPBLiJS2Rr8hAkTPM8Hsi+A/v37s3fvXpkD70BTVk1hzuY5AW0bk+E+rP0wCXYR8bTWTJo0yXOhk9mlCoxtsrKy6NevH+fOnZM58A5ilI6HtR8W0PYxGe4PtXvI7iYIERTfkfprr73GkCFDAlqLxlimICMjg/vuuw+QOfBO4HtOMNC8islwFyIa+M5Tf++99xgwYEBA68Ebc+CXLl3KqFGjAJkiaadQTfaQcBcigimlPGWUzz77zLNcsNmAN6ZIzpw5k3nz5skUSZuEchafhLsQEc6okyul+Oyzzzw3/AikBg8wdOhQ1q9fL6tIhlmop2dLuAsRBYyRdq1atVi+fLnnln1mZ9EYrx80aBAnT56UGTRhYsV1NxLuQkQJ42Kkq6++mkWLFgW0dowxgyYzM5P7778fkBOsVrPqgkoJdwuknUizuwkiRhm37PvZz37GJ598EtA+jJOyK1eu9JxgldG7Nay8Ul7CPcRSM1MZuWSk3c0QMcy42nTQoEH84Q9/AMxf5GQsRzBz5kw+/fRTqb9bwOolUCTcQyg1M5XeC3rTrHYzu5siYlxcXBxaa5577jmeeOKJgC5yMgwdOpQDBw5I/T2EwrG2lYR7iBjB3rpea2b1nWV3c4TwmDVrFn369AnoIieXy0V2djYPPPCA53upvwcnXIsWSriHgG+wL/vfZdSoWsPuJglRag78Rx99RNOmTU0HvPH6devWyQVOIRDO1Wgl3INUNtgTExLtbpIQHsYUyZo1a7Jo0SISEhJMl1aM18+cOZPPP/9cLnAKULiXGZdwD4IEu4gExsnQNm3aMG/evKD2NWTIEI4cOSL1d5PsuH+EhHuAJNhFJDHC+IEHHuD5558HCKj+npWV5bkHa6AnaGONXTcGkt4JgAS7iERGGE+bNo1bb7014Pr78uXLee211wCkPFMJO+/4JuFukgS7iGRGKeXDDz+kQYMGppcoMLZ/7rnnSEtLk/JMBey+laff4a6UamRlQyKBBLuIdMbJ0CuvvJK///3vgPnZL0opLly4wODBgykqKpLyzCXYHexgbuSeppQabFlLHE6CXUQL4wRrz549mTp1KmCufm5cELV161bGjx8PSHnGlxOCHcyF+wvAG0qpT5VSV1jVICeSYBfRJi4uDoAxY8bQu3fvgOrvAFOnTmXt2rVSninhlGAHE+GutX4TaAckAzuVUv0ta5WDSLCLaP/QrwwAAB9/SURBVGWE8fvvv0+9evVMh7NRq3/sscc868fH8gVOTgp2MHlCVWu9X2t9KzAF+FQptV0ptdn3YU0z7SHBLqKZUX9v0KABb7/9tuntjfLMzp07mTBhAhC7q0c6LdghgNkySqlrgEHAaWDRJR5RQYJdxAKjnHL33Xfz1FNPAebq777lmXXr1sVkecaJwQ5g6kaLSqlHgVeB/wNaa61PWtIqm0mwi1hilFdeffVVVqxYwe7du03d6MN47W9+8xvWr18fU7NnZm+ezVsb33JcsIO5qZDLgD8AT2qtB0qwCxEdlFIUFRVRpUoV3n33XcDc9EijPJOamsr06dOB2Jk949RgB3NlmTigrdZ6vlWNsZsEu4hVxvTIbt26MXnyZCCw8syLL75IRkZGzJRnRnQc4chgB3OzZXprrQ9b2Rg7SbCLWGdMjxw3bhxdunQxffWqy+WioKCAZ599FjB/96dINLz9cLubUK7YKY5VQIJdCDdjtP3Xv/7V9NRG48Ng8eLFLFy40FPuEfaI+XCXYBfCy5ge2a5dOyZNmuR5zqxnn32WvLw8z+3+RPjFdLhLsAtxMaM8M3bsWFq1amWqPGOcXP3+++/5/e9/D8Tu3He7xWy4S7ALUb7i4mLi4uJ46623AHOzZ4wwnzZtWkydXHWaiAt3pVSCUmqrUkorpW4MZB/pJ9Ml2IWogFGe6d69O7/97W89z5nZ/sKFC4wePRqIjZOrVkk/mR7QdhEX7sAfgSPB7GDU8lES7EJUwgjzKVOm0LhxY1PlGWOk/sUXX/Dll1+Wulm38F9qZiqjlo8KaNuICnelVF/gF8DoYPbTNLmpBLsQlTBmu1SvXp2ZM2cGtD24a/dGLV74zygdN01uGtD2EfOvrZS6EngXGAyc93ObBKVUkvEAEgGm954uwS6EH4zZLgMHDqRfv36mQlprjVKK7du3M2fOHCB2rlwNlu85wem9pwe0j4gId+UeArwH/FVrvdHEpmOBsz6PwwDVq1YPdROFiFrGydTp06dTtWpVU+UVY9sJEyZw5swZmRrph7KTPQLNK1vDXSk1seTEaEWPjsBTQBIw1eQhpgI1fR4xf6tAIcwyTq5ef/31jBkzxvOcme2PHj3KH//4R0CmRlYklLP47B65zwJSKnmkAbcCXYB8pVQhsLdk+41KqffL27nWOl9rnW08gBzL3okQUcwI87Fjx9KiRYuATq7OmjWLH374QUbv5Qj19Gxbw11rfUprvbuSRx7wNO67QN1Y8rizZBf3Ay/a1HwhYoZxcrVatWq88sorprd3uVycO3eOiRMnAuZvyh3trLjuxu6Ru1+01j9ordOMB/BdyY/2RfNiZkI4iTHivvfee7n99ts9J0z9YYzeFyxYQHp6Oi6XS8ozJay6oDIiwj3S5Bbk2t0EISxhjLiNdWfMjMBdLheFhYVMmzbNkrZFIiuvlI/IcNdaH9BaK631VrvbUlZOfg5PLn3S7mYIYQljxN21a1eGDBniec4fxkh9/vz5bN68OeZH71YvgRKR4e5UOfk59PmgD/tO77O7KUJYxijFTJgwwfTUSOOD4A9/+IMlbYsU4VjbSsI9RIxgTzuRxpt3vWl3c4SwjHFy9dprr+V3v/sdYH70/o9//CNmR+/hWrRQwj0EfIN9xeAVtK7X2u4mCWEpI8x/97vfkZycLKN3P4VzNVoJ9yCVDfbODTvb3SQhLGeM3q+44grPyo+BjN43bdoUM6P3cC8zLuEeBAl2EcuMMH/66adp2LBhQKP3V1991ZK2OY0d94+QcA+QBLuIdcbo/fLLL+eZZ54BzI/eP/nkk6if927XjYEk3AMgwS6EmxHmjz32WECj98LCQv70pz9Z1Tzb2XnHNwl3kyTYhfAyRu9JSUkBj97//ve/c+jQIVwuV1QtS2D3rTwl3E2QYBfiYkaYP/744wGN3vPy8njjjTeA6Fkx0u5gBwl3v0mwC3Fpxug9MTEx4NH73LlzOXXqVFSsGOmEYAcJd79IsAtRsWBH7ydPnuSjjz4CInv07pRgBwn3SkmwC1E539H7008/Dfg/ejdG6u+88w7FxcXExcVZ1k4rOSnYQcK9QhLsQvjPCPOhQ4dSt25dv0fgRrjv2LGDJUuWAJE3endasIOEe7kk2IUwxxi916lTh4cffhjwf/RuLEZmnFj1d514J3BisIOE+yVJsAsRGN957wkJCaZH78uXL2f79u0opSJi9J52Is2RwQ4S7heRYBcicEYoN2/enPvuuw/wf/RuvO7dd98FIuNWfCOXjHRksIOEeykS7EKEzq9//WvA//q575IEZ86ciYhpkc1qN3NksIOEu4cEuxChYVxp2rVrV3r06AH4X0NXSnH8+HEWLlwIOP/E6qy+sxwZ7CDhDkiwCxFqRig/+uijAW0/b948AMdPi6xRtYbdTShXzIe7BLsQoWeE8r333kuTJk38Lq8Yr/vPf/7Dpk2bAOeP3p0qpsNdgl0I6xQVFZGQkMCwYcMA8ydW33vvPauaFhNiNtwl2IWwlhHSgwcPplq1aqZPrP7rX/8iNzc36laLDJeYDPfzBecl2IWwmDEt8pprrmHAgAGAuYuaDh8+zJdffglExrRIq5wvOB/QdjEZ7qNXjJZgFyKMHnnkEcB8/dw4sRpJV6yGUk5+DqNXjA5o25gM9/1n9kuwCxEGxkj9tttuIyUlxe/tjJH6ypUr2b9/P0qpmBu9G6Xj/Wf2B7R9TIb7jDtmSLALESZFRUXEx8fzwAMPAOZOrObl5bFo0SIgtmbN+J4TnHHHjID2EZPhnlLX/xGEECI4Rpjfe++9xMfHm15v5tNPPwWcP+c9VMpO9gg0r2Iy3IUQ4WOUVFq2bEmvXr08z1XGCPc1a9awbds2IPpH76GcxSfhLoSwnBHKv/rVr0xtZ4z6P//8cyC6Z82Eenq2hLsQwnJGSN95553UqlXL9BWrRrhHa2nGiutuJNyFEJYz5rzXq1ePu+66y/NcZYxw37hxY9SWZqy6oFLCXQgRVgMHDjT1emPUv3jxYiC6SjNWXikv4W6B2Ztn290EIRzHGKn37t2bq6++2nRpxpgSGS2lGauXQJFwD7Epq6bw1sa37G6GEI5jlGYSExO58847Af/mvPuWZrZv3w5EfmkmHGtbSbiH0JRVUxi/cjwjOo6wuylCONo999wD+B/SxofAV199BUR2aSZcixZKuIeIEeyTe01mePvhdjdHCEcyQrp79+40btzY7+2ipTQTztVoJdxDwDfYx/UYZ3dzhHC04uJiqlevTp8+fQBzpZlNmzZx6NChUs9FinAvMy7hHiQJdiEC07dvX8D/0oxSinPnzvHtt9+a2s4J7Lh/hIR7ECTYhTDPmDXTvXt36tSpY3p744KmSFkG2K4bA0m4B0iCXYjAGLNmrrjiCrp37+55rjJGGSY1NZW8vDy/V5e0k513fHP+v44DSbALERwjqI2rVc04ePAgW7ZsAZxdmrH7Vp4S7iZJsAsRPGPU3bNnT6pWrer3ydFImRJpd7BDhIW7UuoupdR6pdRPSqlTSqmF4Ty+BLsQoWGUYZo3b86NN95Y6rmKGGFuhLsTSzNOCHaIoHBXSg0CFgDzgHbAzcCH4Tq+BLsQoVVUVAS4b8EH5sJ9w4YNHD582HG333NKsEOEhLtSKh74M/Cs1vqvWuvvtNYZWut/huP4EuxChJ4R5r/4xS8Ac1Mic3Nz2bRpE+Cc0oyTgh0iJNyB9kBDoFgptUUpdVQptVQp1aqijZRSCUqpJOMBJJo9sAS7ENYwwr19+/bUr1/f9HbGfHcncFqwQ+SE+7UlXycCU4B+wBngW6VU7Qq2Gwuc9XkcNnNQCXYhrGOUVJKSkujUqZPnucoYI/W1a9cC9tfdnRjsYHO4K6UmKqV0JY+OPu18RWv9qdZ6E/AIoIF7KzjEVKCmz6ORv22TYBfCekYpxijN+MMI9127dnHkyJFSz4VbbkGuI4MdIN7m488CPqrkNQfwllN2GU9qrfOVUt8D5a4+pLXOB/KN7/29ok2CXYjwMH4nu3btCpgL6ZycHLZu3cpVV12F1tqWK1afXPokB7IOOC7YweaRu9b6lNZ6dyWPPGAT7pC+3thWKVUFaAIcDGWbJNiFCB+jpJKSkkKTJk1Mb2fU3e0aue87vc+RwQ4RUnPXWmcDfwUmKaV+oZS6HjDuiPFJqI4jwS5E+BUXF3PZZZfRubM7IM3U3detWwfYV3d/8643HRnsECHhXuJZ3CWcBcAG4BrgVq31mVDsXIJdCHsYQW2sM2Nmm927d3P69Gnb5ru3rtc67Mf0V8SEu9b6gtZ6tNb6Sq11kta6t9Z6Zyj2LcEuhH2MkXqHDh0AcyWWEydOsGvXLtPbxYKICXerSLALYS8j3Fu1akWDBg1Mb2fcV1XCvbSYDncJdiHs5zvf3cw6M8ZrjCtVI2V993CJ2XCXYBfCOYz57oFczLRjxw7A/ouZnMbuee62mL9tPnN2zZFgF8IhjDC/6aabAP9KLMZr9uzZw4kTJ6hXr551DbTR/G3zA9ouJj/q5myWYBfCSYxwb9euHdWqVTNVP8/KymLv3r1A9NXdp6yawpzNcwLaNibDfVj7YRLsQjiIEe6NGzfmmmuu8Xs7oxSzbds2ILrC3SgdD2s/LKDtYzLcH2r3kN1NEEKUUVxcTFxcHK1auRd7NXOCdOPGjVY1yxa+5wQDzauYDHchhPMYo26j7m7mpGpGRgbgHslH+ug9VJM9JNyFEI7gO98dzJ1U3bt3L1lZWX5v51ShnMUn4S6EcAQj3FNSUkwvJ3D8+HEOHDgARG64h3p6toS7EMJRGjduzNVXX216u++//96C1oSHFdfdSLgLIRzBGK1Xr16dZs2aeZ6rjDFjJlLXmLHqgkoJdwuknUizuwlCRCTjStWWLVsC5mbM7Nu3z/Q2drPySnkJ9xBLzUxl5JKRdjdDiIjWurV7KV2zV6pC5MyYsXoJFAn3EErNTKX3gt40q93M7qYIEZF8L2YCc+G+b98+srOzI2KNmXCsbeX8f4UIYQR763qtmdV3lt3NESIiGeF+/fXXU6VKFVPbZmVlcfr0acDZdfdwLVoo4R4CvsG+7H+XUaNqDbubJERES05OpmbNmqa2ycvLc/x0yHCuRivhHqSywZ6YkGh3k4SIWMbIvU6dOp7SjJm13Y8fP25d44IU7mXGJdyDIMEuROgZo+6GDRv6vY0R7sba7k4budtx/wgJ9wBJsAthDWM6ZIsWLQBzUxuNJQicNB3SrhsDSbgHQIJdCOs1atQIMDdjxrifqlNmzNh5xzdn/AtEEAl2Iazlu8YMmCux5ObmltqHney+laeEuwkS7EKEz+WXX+73a40PgMOHD3Pq1Cni4+29g6jdwQ4S7n6TYBciPIxRd7NmzUxPh8zOziY3N9fWsowTgh0k3P0iwS5E+FWrVs10SBcWFnquUrVjxoxTgh0k3CslwS5EeBkj9+TkZM/9VP2d615YWMjevXttCXcnBTtIuFdIgl0IexjBXLduXdPb5ufnh7o5lXJasIOEe7kk2IWwX1xcnOltjCUIwsWJwQ4S7pckwS6EvYyRuzEd0swSBAcPHrSuYWXM3jzbkcEOEu4XkWAXwn5GuNeuXdv0tkZZJhxz3d/a+JYjgx0k3EuRYBfCWcws+2t8IBQWFlrVnIuM6DjCkcEOEu4eEuxCOIcx6jYWDzMz8yU9PZ38/PywjNyHtx9u+TECJeGOBLsQTmXmjkyGvLw8z+JjsSzmw12CXQjnMnMRk29ZJpylGaeK6XCXYBfC2apWrWp6G62149Zzt0PMhrsEuxDOZdTLa9asafpeqmfOnOHMmTOA827aEU4xGe7pJ9Ml2IWIAFdddRVXXHGFqW1OnDjBsWPHgOgI9/ST6QFtF5PhPmr5KAl2ISJAUVFRQAFt95K/oZKamcqo5aMC2jYmw71pclMJdiEiQHx8fEDL9164cMGC1oSXUTpumtw0oO1jMtyn954uwS6Egxk19ypVqpg6qWpst3fvXiByyzK+5wSn954e0D5iMtyrV61udxOEEBYwwr2oqMjmlgSu7GSPQPMqJsNdCBHdnHAP1UCEchZfxIS7Uuo6pdQipdQppVS2Uuo/SqledrdLCOE8kViOCfX07IgJd2AJEA/cCnQAtgJfKKXq29oqIYQIkhXX3UREuCul6gDNgWla6+1a6z3AGKA60MrWxgkhRBCsuqAyUiaD/gikAw8ppTYD+cDjwHFgU3kbKaUSgASfp8IyRSa3IDcchxHCcTp2hJLrh0KkKseO/Rfw7wSpsV7Yb3+bzIsvQijHr/Xrw8aNIdsdYO2V8hER7lprrZTqDSwCcoBi3MHeR2udVcGmY4GXwtBEj5z8HJ5c+mQ4DymEYxw7BpmZodyjAhqY3urMGffDvb0zWb0Eiq3hrpSaSOXh2wn36PxN4ATQHfgJGI675t5Ja320nG2nAjN8vk8EDgfT5ork5OfQ54M+7Du9z6pDCOFo9UN+Bkxz7Ngx01Mbk5OTqV69BqAJVcCH8r2FY20ru0fus4CPKnnNAdwnUfsByVrr7JLnR5aM5ocA0y61odY6H3cJB7B2epQR7Gkn0njzrjcZ8q8hlh1LCKcKddkiP7+A66//ud/3RXW5XBQXFzNz5nsMGTKEoqLigG6ybaVwLVpoa7hrrU8Bpyp7nVLKmMVfdgX+YhxwUtg32FcMXkG8y+7PTCGEE4VzNdpISaH/AmeA95VSL+MuyzwKNMU9RdKUjB8zuPzo5SFpWG5BLk8ufZJ9p/fx5l1vEu+KJ/2UexU346sQIjAFBQXkX5EPBf69XisNGg7kH2Dz0c0UFRVZMnIP5Hc87UQaI5eMpFntZky7bRp7Tu/xa7uMHzMCaqOKlMn+SqmOwCtAR6AKsBN4WWu91MQ+koCzjAGqWdJMIYQIrTyMwnNNn7J0pSIm3EPBCPcXv3yRgTcODGpfJ8+dZOA/BnK+8Dwv9XiJ/i37e36WfiqdBxc+yN8G/o2UOilBtrpiszfP5q2NbzGi4wjbbtbrOyKZ1XcWNarWCHsbyv4F1bpe67C3AaQ/DKHqj4KCAu4ZeA/Hjvo3v1IphdaaiRMn8stf/pJ3Nr7D25vfDnl/mPkd/3z350xaNYnq8dVZeN9C6l5e19SxFm5dyCt3vgImw91zS6pYeABJgH5n7Ts6GJlnM3Xi7xM1E9FzN8+96OebjmzSTERvOrIpqONUZvK3kzUT0ZO/nWzpcSqy/vB6nTQ1SXed01Vn52Xb0obsvGzddU5XnTQ1Sa8/vN6WNmgt/WEIZX/k5eXpa665RuOe9lLpw+VyaUC/9957lvaHv7/jczfP1UxEJ/4+UWeezQzoWO+sfcd4f0naRN5FSs3dMY5kH6HlGy3JKchhbv+5PHLTI7a0Y8qqKYxfOZ7JvSYzrsc4W9rghFsVlj2Z3blh57C3AaQ/DE7pj8+zPmfhyoVM6jnJtv6Yt2UeQz8fSmLVRHY/sZurkq4K6/El3E2QYPeSIPGS/nBzSn/QAxZmLWRyr8mMvXmsLU2wO9jBAdMII4UEu5cEiZf0h5tT+kN313ArDKw1MGZH7AYJdz9IsHtJkHhJf7g5pT/oAbqXhq+hf63+lb/eAk4JdpBwr5QEu5cEiZf0h5tT+oMewK2gVipYZU8TnBTsIOFeIQl2LwkSL+kPN6f0hxHsfA1qtT0LhTkt2EHCvVwS7F4SJF7SH25O6Q/fYJcRe2kS7pcgwe4lQeIl/eHmlP5wQrAvyljkyGAHCfeLSLB7SZB4SX+4OaU/nBDsAC9/+7Ijgx0k3EuRYPeSIPGS/nBzSn84IdjTTqQBUL1KdUcGO0i4e0iwe0mQeEl/uDmlP5wQ7KmZqYxcMhKAhfctdGSwg4Q7IMHuS4LES/rDzSn94ZRg772gN81rNwegbg1zi4CFU8yHuwS7lwSJl/SHm1P6w0nB3rpea17v+7o9jTAhpsNdgt1LgsRL+sPNKf2R1S7LUcG+7H+X2bKMslkxG+4S7F4SJF7SH25O6Y+p/5nK2fZnHRXsdvRHIGIy3LN+ypJgL+GE/3GdEiTSH25O6o9JqydRc3PNmA72rJ+yAtvQzOLvkf6g5GYdCeMSyr3RRij4s5C/3NjBTW604SX94WX0x0tfvRTwzTq01rqwsDDgNlTUH+G6IU/m2Ux3XgVwsw7bAzecDyPcGWNdsGtdecdLkLg5LUikP5zXH8HciUnrwMO9sv4IR7h77vg2Bgn3St9sSbgP+WhIEP/klauo4yVI3JwYJHaR/vAq2x92hLs//WF1uPveynPIR0MCCveYrLnf3PhmW44rNV03J9V0pT+isz+UCmx1SCf0R9nJHoHmVUyGux0kSNyiMUgCJf3hFer+0O6/1E1xQn+EchafhHsYSJC4RWuQBEL6w0v6wy3U07Ml3C0m/+O6SZB4SX94VdYf2nu+zDJO6A8rrruJD0G7Is7q71Zbuv+DWQchD15e/jKLMhYxrP0wetTvwarvwj9ZN/1kOqOWj6JpclPGdx7PloNbwt6G8wXnGb1iNPvP7GfGHTPIy82z5d9i/rb5zNk8R/ojgvqjsLCQorpFkO/fPotVMWjYnbWbVd+toqioiLi4uHJfH2h/ZPyYAXmwcOtCNh3Y5F/jypH1UxbjV44nvzCfITcOoTCvkHf/+67n54HmlbL6U9FJlFINgcN2t0MIIQLQSGud6e+LYy3cFXAVkGNTExJxf7g0srENoSbvKTJE23uKtvcDFb+nROCINhHYMVWWKfmH8fuTL9R8pmflaK2z7WpHKMl7igzR9p6i7f1Ape/J9HuUE6pCCBGFJNyFECIKSbiHVz4wCb/P/UcEeU+RIdreU7S9Hwjxe4qpE6pCCBErZOQuhBBRSMJdCCGikIS7EEJEIQl3IYSIQhLuYaKUGqmU2q+UylNKbVJKdbe7TYFSSo1VSm1QSuUopU4opf6llLre7naFUsl71Eqp1+xuSzCUUg2VUn9TSv2olDqvlNqqlOpgd7sCpZSKV0pNKfld+kkp9b1SaoJSKmKyTCnVQym1WCl1pOT/sbvL/FwppSaW/PwnpdQ3SqlWZo8TMf8gkUwpdT/wGvAKcBOwGliqlGpsa8MC1xN4A+gC9MZ9pfO/lVI1bG1ViCilOgGPAdvtbkswlFLJwH+AC0Bf4Abgd0CAd1x2hOeBXwNPAinAc8CzwFN2NsqkGsA23O/hUp4DRpX8vBNwDFihlDK1XKVMhQwDpdR6YLPWeoTPc+nAv7TWY+1rWWgopeoCJ4CeWmub7lMfGkqpy4HNwEhgHLBVa/2Mva0KjFJqGnCz1jpi/0osSyn1BXBcaz3M57lPgfNa68H2tSwwSikN3KO1/lfJ9wo4Arymtf5DyXMJwHHgea312/7uW0buFlNKVQU6AP8u86N/A13D3yJL1Cz5etrWVoTGG8ASrfX/2d2QEOgPbFRKfVJSPtuilHrU7kYFaQ1wm1LqOgClVDugG/Clra0KnaZAfXzyQmudD3yLybyIqYXDbFIHiMP9yevrOO5OjGglI40ZwBqtdZrd7QmGUup/gPa4/xSOBtcCI3D3z++BzsBflFL5Wuv5trYscH/APZjYrZQqwv279aLW+u/2NitkjEy4VF5cY2ZHEu7hU7b+pS7xXCSaBbTFPXqKWEqpq4E/A7/QWufZ3Z4QcQEbtdYvlHy/peTE3AggUsP9fuBB4FfATuBG4DWl1BGt9fu2tiy0gs4LCXfrnQKKuHiUXo+LP50jilLqddx/+vfQWkf6TVA64O6TTT5Lr8YBPZRSTwIJWusiuxoXoKPArjLPpQODbGhLqPwJmKa1/qjk+x1KqWuAsUA0hPuxkq/1cfefwXReSM3dYlrrAmAT7lklvnoDa8PfouCVTNWaBQwEbtVa77e7TSHwFdAG90jQeGwEPgBujMBgB/dMmbJTVK8DDtrQllCpDhSXea6I6Mmy/bgD3pMXJeftemIyL2TkHh4zgAVKqY3Af3FPs2sM/NXWVgXuDdx/Fg8AcpRSxl8lZ7XWP9nXrMBprXOAUucMlFK5wI8RfC5hJrBWKfUC8A/cNffHSh6RajHwolLqB9xlmZtwTxuca2urTCiZkdXc56mmSqkbgdNa6x9Krq14QSm1B9gDvACcBz40dSDj7uLysPaBe2rdAdzLeW7CXcqwvV0BvhddzuNhu9sW4vf5De4paba3JYj30A/YAeThLsk8anebgnw/ibivGTkI/ATsA6YAVe1um4n3cEs5vz/vlfxcARNxl2XycM+UaW32ODLPXQgholC01KmEEEL4kHAXQogoJOEuhBBRSMJdCCGikIS7EEJEIQl3IYSIQhLuQggRhSTchRAiCkm4CxEApVScUmptyY0ifJ+vqZQ6pJSaYlfbhAC5E5MQAVNKtQC2Ao9prT8oeW4+0A7opN2LxglhCwl3IYKglHoa9zogrXHf5OMToLPWequd7RJCwl2IIJTciepr3MvOtgFe11pLSUbYTsJdiCAppVriXnFxB9Bea11oc5OEkBOqQoTAUNzrbTcFGtncFiEAGbkLERSl1M+BVUBf4Dnct+a7XcsvlrCZjNyFCJBS6jLc9+18W2v9f8Bw3CdVH7e1YUIg4S5EMKbh/h16HkBr/QPwO+BPSqkm9jVLCCnLCBEQpVRP3DfVvkVrvabMz5bjvj+xlGeEbSTchRAiCklZRgghopCEuxBCRCEJdyGEiEIS7kIIEYUk3IUQIgpJuAshRBSScBdCiCgk4S6EEFFIwl0IIaKQhLsQQkQhCXchhIhCEu5CCBGF/j9p4SRSoE/fswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# If there is a warning that reads \"The specified user volume\n",
    "# is larger than the simulation domain and has been truncated\",\n",
    "# It has to do with some numerical errors between python and meep.\n",
    "# Ignore.\n",
    "\n",
    "f = plt.figure(dpi=100)\n",
    "sim.plot2D(ax=f.gca())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks pretty good. Simulations at the high enough resolution required to avoid spurious reflections in the bend are very slow! This can be sped up quite a bit by running the code in parallel from the terminal. Later, we will put this notebook's code into a script and run it in parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4. Simulate FDTD and Animate results\n",
    "\n",
    "More detailed meep documentation available [here](https://meep.readthedocs.io/en/latest/Python_Tutorials/Basics/#transmittance-spectrum-of-a-waveguide-bend)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set to true to compute animation (may take a lot of memory)\n",
    "compute_animation = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stop condition: decay to 0.1% of peak value in the last 15.5 time units.\n",
      "on time step 200 (time=1), 0.0383541 s/step\n",
      "on time step 477 (time=2.385), 0.0144558 s/step\n",
      "on time step 755 (time=3.775), 0.0144323 s/step\n",
      "on time step 1029 (time=5.145), 0.0146304 s/step\n",
      "on time step 1306 (time=6.53), 0.0144882 s/step\n",
      "on time step 1583 (time=7.915), 0.0144634 s/step\n",
      "on time step 1857 (time=9.285), 0.0146198 s/step\n",
      "on time step 2135 (time=10.675), 0.0144386 s/step\n",
      "on time step 2409 (time=12.045), 0.0146068 s/step\n",
      "on time step 2687 (time=13.435), 0.0144364 s/step\n",
      "on time step 2965 (time=14.825), 0.0144058 s/step\n",
      "field decay(t = 15.505): 1.0968531133594091e-76 / 1.0968531133594091e-76 = 1.0\n",
      "on time step 3239 (time=16.195), 0.0146098 s/step\n",
      "on time step 3517 (time=17.585), 0.0143887 s/step\n",
      "on time step 3796 (time=18.98), 0.0143683 s/step\n",
      "on time step 4071 (time=20.355), 0.0145537 s/step\n",
      "on time step 4350 (time=21.75), 0.0143623 s/step\n",
      "on time step 4625 (time=23.125), 0.0145803 s/step\n",
      "on time step 4904 (time=24.52), 0.0143763 s/step\n",
      "on time step 5182 (time=25.91), 0.0143901 s/step\n",
      "on time step 5457 (time=27.285), 0.0145944 s/step\n",
      "on time step 5735 (time=28.675), 0.014404 s/step\n",
      "on time step 6010 (time=30.05), 0.0145722 s/step\n",
      "field decay(t = 31.01): 1.9557107984193176e-10 / 1.9557107984193176e-10 = 1.0\n",
      "on time step 6289 (time=31.445), 0.014379 s/step\n",
      "on time step 6567 (time=32.835), 0.0143953 s/step\n",
      "on time step 6842 (time=34.21), 0.014556 s/step\n",
      "on time step 7121 (time=35.605), 0.0143872 s/step\n",
      "on time step 7399 (time=36.995), 0.0144125 s/step\n",
      "on time step 7674 (time=38.37), 0.0145564 s/step\n",
      "on time step 7953 (time=39.765), 0.0143832 s/step\n",
      "on time step 8228 (time=41.14), 0.0145736 s/step\n",
      "on time step 8507 (time=42.535), 0.0143846 s/step\n",
      "on time step 8785 (time=43.925), 0.0143967 s/step\n",
      "on time step 9060 (time=45.3), 0.0145485 s/step\n",
      "field decay(t = 46.515): 8.409707458954242e-07 / 8.409707458954242e-07 = 1.0\n",
      "on time step 9338 (time=46.69), 0.0143913 s/step\n",
      "on time step 9613 (time=48.065), 0.0145733 s/step\n",
      "on time step 9892 (time=49.46), 0.0143734 s/step\n",
      "on time step 10170 (time=50.85), 0.0144098 s/step\n",
      "on time step 10445 (time=52.225), 0.0145753 s/step\n",
      "on time step 10724 (time=53.62), 0.0143822 s/step\n",
      "on time step 11000 (time=55), 0.0145901 s/step\n",
      "on time step 11279 (time=56.395), 0.0143693 s/step\n",
      "on time step 11557 (time=57.785), 0.0143891 s/step\n",
      "on time step 11832 (time=59.16), 0.0145821 s/step\n",
      "on time step 12110 (time=60.55), 0.0143941 s/step\n",
      "on time step 12389 (time=61.945), 0.0143824 s/step\n",
      "field decay(t = 62.02): 0.0007836348644637144 / 0.0007836348644637144 = 1.0\n",
      "on time step 12664 (time=63.32), 0.0145938 s/step\n",
      "on time step 12942 (time=64.71), 0.0144102 s/step\n",
      "on time step 13217 (time=66.085), 0.0145961 s/step\n",
      "on time step 13495 (time=67.475), 0.0144085 s/step\n",
      "on time step 13774 (time=68.87), 0.0143874 s/step\n",
      "on time step 14049 (time=70.245), 0.0145881 s/step\n",
      "on time step 14327 (time=71.635), 0.0143903 s/step\n",
      "on time step 14602 (time=73.01), 0.0145596 s/step\n",
      "on time step 14881 (time=74.405), 0.0143841 s/step\n",
      "on time step 15160 (time=75.8), 0.0143825 s/step\n",
      "on time step 15435 (time=77.175), 0.0145681 s/step\n",
      "field decay(t = 77.525): 1.4724243279845177 / 1.4724243279845177 = 1.0\n",
      "on time step 15714 (time=78.57), 0.0143819 s/step\n",
      "on time step 15993 (time=79.965), 0.0143726 s/step\n",
      "on time step 16269 (time=81.345), 0.014539 s/step\n",
      "on time step 16548 (time=82.74), 0.0143661 s/step\n",
      "on time step 16823 (time=84.115), 0.0145743 s/step\n",
      "on time step 17102 (time=85.51), 0.0143864 s/step\n",
      "on time step 17380 (time=86.9), 0.014416 s/step\n",
      "on time step 17655 (time=88.275), 0.0145817 s/step\n",
      "on time step 17934 (time=89.67), 0.0143798 s/step\n",
      "on time step 18209 (time=91.045), 0.0145599 s/step\n",
      "on time step 18488 (time=92.44), 0.0143826 s/step\n",
      "field decay(t = 93.03): 3.4417968622809 / 3.4417968622809 = 1.0\n",
      "on time step 18767 (time=93.835), 0.014362 s/step\n",
      "on time step 19042 (time=95.21), 0.0145779 s/step\n",
      "on time step 19321 (time=96.605), 0.0143754 s/step\n",
      "on time step 19600 (time=98), 0.0145558 s/step\n",
      "on time step 19879 (time=99.395), 0.0143843 s/step\n",
      "on time step 20158 (time=100.79), 0.0143639 s/step\n",
      "on time step 20433 (time=102.165), 0.0145757 s/step\n",
      "on time step 20712 (time=103.56), 0.0143783 s/step\n",
      "on time step 20991 (time=104.955), 0.0143581 s/step\n",
      "on time step 21266 (time=106.33), 0.0145521 s/step\n",
      "on time step 21544 (time=107.72), 0.0143957 s/step\n",
      "field decay(t = 108.535): 0.8564689296158491 / 3.4417968622809 = 0.2488435441969291\n",
      "on time step 21819 (time=109.095), 0.0145462 s/step\n",
      "on time step 22097 (time=110.485), 0.0144197 s/step\n",
      "on time step 22375 (time=111.875), 0.014394 s/step\n",
      "on time step 22649 (time=113.245), 0.0146112 s/step\n",
      "on time step 22928 (time=114.64), 0.0143857 s/step\n",
      "on time step 23203 (time=116.015), 0.0145816 s/step\n",
      "on time step 23481 (time=117.405), 0.0144068 s/step\n",
      "on time step 23760 (time=118.8), 0.014357 s/step\n",
      "on time step 24035 (time=120.175), 0.0145713 s/step\n",
      "on time step 24314 (time=121.57), 0.014373 s/step\n",
      "on time step 24593 (time=122.965), 0.0143595 s/step\n",
      "field decay(t = 124.04): 0.00013649421971979937 / 3.4417968622809 = 3.965783722324153e-05\n",
      "Normalizing field data...\n",
      "run 0 finished at t = 124.04 (24808 timesteps)\n",
      "Generating MP4...\n"
     ]
    }
   ],
   "source": [
    "# Setup and run the simulation\n",
    "\n",
    "# The following line defines a stopping condition depending on the square\n",
    "# of the amplitude of the Ez field at the port 2.\n",
    "print(f\"Stop condition: decay to 0.1% of peak value in the last {2.0/df:.1f} time units.\")\n",
    "stop_condition = mp.stop_when_fields_decayed(2.0/df,mp.Ez,p2.center,1e-3)\n",
    "if compute_animation:\n",
    "    f = plt.figure(dpi=100)\n",
    "    animate = mp.Animate2D(sim,mp.Ez,f=f,normalize=True)\n",
    "    sim.run(mp.at_every(1,animate), until_after_sources=stop_condition)\n",
    "    plt.close()\n",
    "    # Save video as mp4\n",
    "    animate.to_mp4(10, 'media/bend.mp4')\n",
    "else:\n",
    "    sim.run(until_after_sources=stop_condition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize results\n",
    "\n",
    "Things to check:\n",
    "- Was the simulation time long enough for the pulse to travel through port2 in its entirety? Given the automatic stop condition, this should be the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video src=\"media/bend.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import Video, display\n",
    "display(Video('media/bend.mp4'))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "from IPython.display import Video, display\n",
    "if compute_animation:\n",
    "    # Save video as gif\n",
    "    animate.to_mp4(10, 'media/bend.mp4')\n",
    "    # Play video\n",
    "    display(animate.to_jshtml(10))\n",
    "else:\n",
    "    # Load pre-computed video from .gif\n",
    "    display(Video('media/bend.mp4'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5. Compute loss and reflection of the bend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MPB solved for omega_1(2.01141,0,0) = 0.618698 after 34 iters\n",
      "MPB solved for omega_1(1.87532,0,0) = 0.580676 after 10 iters\n",
      "MPB solved for omega_1(1.87521,0,0) = 0.580645 after 4 iters\n",
      "MPB solved for omega_1(1.87521,0,0) = 0.580645 after 1 iters\n",
      "Dominant planewave for band 1: (1.875210,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.03004,0,0) = 0.623906 after 34 iters\n",
      "MPB solved for omega_1(1.89458,0,0) = 0.586052 after 10 iters\n",
      "MPB solved for omega_1(1.89447,0,0) = 0.586022 after 4 iters\n",
      "MPB solved for omega_1(1.89447,0,0) = 0.586022 after 1 iters\n",
      "Dominant planewave for band 1: (1.894467,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.04866,0,0) = 0.629115 after 24 iters\n",
      "MPB solved for omega_1(1.91383,0,0) = 0.591428 after 10 iters\n",
      "MPB solved for omega_1(1.91372,0,0) = 0.591398 after 4 iters\n",
      "MPB solved for omega_1(1.91372,0,0) = 0.591398 after 1 iters\n",
      "Dominant planewave for band 1: (1.913720,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.06729,0,0) = 0.634325 after 27 iters\n",
      "MPB solved for omega_1(1.93307,0,0) = 0.596804 after 10 iters\n",
      "MPB solved for omega_1(1.93297,0,0) = 0.596774 after 4 iters\n",
      "MPB solved for omega_1(1.93297,0,0) = 0.596774 after 1 iters\n",
      "Dominant planewave for band 1: (1.932968,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.08591,0,0) = 0.639536 after 31 iters\n",
      "MPB solved for omega_1(1.95232,0,0) = 0.60218 after 10 iters\n",
      "MPB solved for omega_1(1.95221,0,0) = 0.602151 after 4 iters\n",
      "MPB solved for omega_1(1.95221,0,0) = 0.602151 after 1 iters\n",
      "Dominant planewave for band 1: (1.952212,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.10453,0,0) = 0.644749 after 31 iters\n",
      "MPB solved for omega_1(1.97156,0,0) = 0.607556 after 10 iters\n",
      "MPB solved for omega_1(1.97145,0,0) = 0.607527 after 4 iters\n",
      "MPB solved for omega_1(1.97145,0,0) = 0.607527 after 1 iters\n",
      "Dominant planewave for band 1: (1.971452,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.12316,0,0) = 0.649962 after 26 iters\n",
      "MPB solved for omega_1(1.99079,0,0) = 0.612932 after 10 iters\n",
      "MPB solved for omega_1(1.99069,0,0) = 0.612903 after 4 iters\n",
      "MPB solved for omega_1(1.99069,0,0) = 0.612903 after 1 iters\n",
      "Dominant planewave for band 1: (1.990687,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.14178,0,0) = 0.655177 after 28 iters\n",
      "MPB solved for omega_1(2.01002,0,0) = 0.618308 after 10 iters\n",
      "MPB solved for omega_1(2.00992,0,0) = 0.61828 after 4 iters\n",
      "MPB solved for omega_1(2.00992,0,0) = 0.61828 after 1 iters\n",
      "Dominant planewave for band 1: (2.009918,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.16041,0,0) = 0.660393 after 33 iters\n",
      "MPB solved for omega_1(2.02924,0,0) = 0.623684 after 10 iters\n",
      "MPB solved for omega_1(2.02915,0,0) = 0.623656 after 4 iters\n",
      "MPB solved for omega_1(2.02915,0,0) = 0.623656 after 1 iters\n",
      "Dominant planewave for band 1: (2.029145,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.17903,0,0) = 0.66561 after 36 iters\n",
      "MPB solved for omega_1(2.04847,0,0) = 0.62906 after 10 iters\n",
      "MPB solved for omega_1(2.04837,0,0) = 0.629032 after 3 iters\n",
      "MPB solved for omega_1(2.04837,0,0) = 0.629032 after 1 iters\n",
      "Dominant planewave for band 1: (2.048367,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.19766,0,0) = 0.670828 after 27 iters\n",
      "MPB solved for omega_1(2.06768,0,0) = 0.634436 after 10 iters\n",
      "MPB solved for omega_1(2.06759,0,0) = 0.634409 after 3 iters\n",
      "MPB solved for omega_1(2.06759,0,0) = 0.634409 after 1 iters\n",
      "Dominant planewave for band 1: (2.067585,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.21628,0,0) = 0.676047 after 33 iters\n",
      "MPB solved for omega_1(2.08689,0,0) = 0.639812 after 10 iters\n",
      "MPB solved for omega_1(2.0868,0,0) = 0.639785 after 3 iters\n",
      "MPB solved for omega_1(2.0868,0,0) = 0.639785 after 1 iters\n",
      "Dominant planewave for band 1: (2.086799,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.2349,0,0) = 0.681268 after 32 iters\n",
      "MPB solved for omega_1(2.1061,0,0) = 0.645188 after 10 iters\n",
      "MPB solved for omega_1(2.10601,0,0) = 0.645161 after 3 iters\n",
      "MPB solved for omega_1(2.10601,0,0) = 0.645161 after 1 iters\n",
      "Dominant planewave for band 1: (2.106008,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.25353,0,0) = 0.686489 after 34 iters\n",
      "MPB solved for omega_1(2.12531,0,0) = 0.650564 after 10 iters\n",
      "MPB solved for omega_1(2.12521,0,0) = 0.650538 after 3 iters\n",
      "MPB solved for omega_1(2.12521,0,0) = 0.650538 after 1 iters\n",
      "Dominant planewave for band 1: (2.125213,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.27215,0,0) = 0.691711 after 27 iters\n",
      "MPB solved for omega_1(2.14451,0,0) = 0.65594 after 10 iters\n",
      "MPB solved for omega_1(2.14441,0,0) = 0.655914 after 3 iters\n",
      "MPB solved for omega_1(2.14441,0,0) = 0.655914 after 1 iters\n",
      "Dominant planewave for band 1: (2.144414,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.29078,0,0) = 0.696935 after 31 iters\n",
      "MPB solved for omega_1(2.1637,0,0) = 0.661316 after 10 iters\n",
      "MPB solved for omega_1(2.16361,0,0) = 0.66129 after 3 iters\n",
      "MPB solved for omega_1(2.16361,0,0) = 0.66129 after 1 iters\n",
      "Dominant planewave for band 1: (2.163611,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.3094,0,0) = 0.70216 after 27 iters\n",
      "MPB solved for omega_1(2.18289,0,0) = 0.666692 after 10 iters\n",
      "MPB solved for omega_1(2.1828,0,0) = 0.666667 after 3 iters\n",
      "MPB solved for omega_1(2.1828,0,0) = 0.666667 after 1 iters\n",
      "Dominant planewave for band 1: (2.182803,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.32803,0,0) = 0.707385 after 36 iters\n",
      "MPB solved for omega_1(2.20208,0,0) = 0.672068 after 10 iters\n",
      "MPB solved for omega_1(2.20199,0,0) = 0.672043 after 3 iters\n",
      "MPB solved for omega_1(2.20199,0,0) = 0.672043 after 1 iters\n",
      "Dominant planewave for band 1: (2.201991,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.34665,0,0) = 0.712612 after 33 iters\n",
      "MPB solved for omega_1(2.22126,0,0) = 0.677444 after 10 iters\n",
      "MPB solved for omega_1(2.22118,0,0) = 0.677419 after 3 iters\n",
      "MPB solved for omega_1(2.22118,0,0) = 0.677419 after 1 iters\n",
      "Dominant planewave for band 1: (2.221175,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.36527,0,0) = 0.71784 after 29 iters\n",
      "MPB solved for omega_1(2.24044,0,0) = 0.68282 after 10 iters\n",
      "MPB solved for omega_1(2.24036,0,0) = 0.682796 after 3 iters\n",
      "MPB solved for omega_1(2.24036,0,0) = 0.682796 after 1 iters\n",
      "Dominant planewave for band 1: (2.240355,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.3839,0,0) = 0.723069 after 29 iters\n",
      "MPB solved for omega_1(2.25962,0,0) = 0.688196 after 10 iters\n",
      "MPB solved for omega_1(2.25953,0,0) = 0.688172 after 3 iters\n",
      "MPB solved for omega_1(2.25953,0,0) = 0.688172 after 1 iters\n",
      "Dominant planewave for band 1: (2.259531,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.40252,0,0) = 0.728299 after 41 iters\n",
      "MPB solved for omega_1(2.27879,0,0) = 0.693572 after 10 iters\n",
      "MPB solved for omega_1(2.2787,0,0) = 0.693548 after 3 iters\n",
      "MPB solved for omega_1(2.2787,0,0) = 0.693548 after 1 iters\n",
      "Dominant planewave for band 1: (2.278703,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.42115,0,0) = 0.73353 after 39 iters\n",
      "MPB solved for omega_1(2.29795,0,0) = 0.698948 after 10 iters\n",
      "MPB solved for omega_1(2.29787,0,0) = 0.698925 after 3 iters\n",
      "MPB solved for omega_1(2.29787,0,0) = 0.698925 after 1 iters\n",
      "Dominant planewave for band 1: (2.297870,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.43977,0,0) = 0.738762 after 30 iters\n",
      "MPB solved for omega_1(2.31711,0,0) = 0.704324 after 10 iters\n",
      "MPB solved for omega_1(2.31703,0,0) = 0.704301 after 3 iters\n",
      "MPB solved for omega_1(2.31703,0,0) = 0.704301 after 1 iters\n",
      "Dominant planewave for band 1: (2.317033,-0.000000,0.000000)\n",
      "MPB solved for omega_1(2.45839,0,0) = 0.743995 after 30 iters\n",
      "MPB solved for omega_1(2.33627,0,0) = 0.7097 after 10 iters\n",
      "MPB solved for omega_1(2.33619,0,0) = 0.709677 after 3 iters\n",
      "MPB solved for omega_1(2.33619,0,0) = 0.709677 after 1 iters\n",
      "Dominant planewave for band 1: (2.336193,-0.000000,0.000000)\n",
      "MPB solved for omega_1(0,2.01141,0) = 0.618698 after 27 iters\n",
      "MPB solved for omega_1(0,1.87532,0) = 0.580676 after 7 iters\n",
      "MPB solved for omega_1(0,1.87521,0) = 0.580645 after 3 iters\n",
      "MPB solved for omega_1(0,1.87521,0) = 0.580645 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,1.875210,0.000000)\n",
      "MPB solved for omega_1(0,2.03004,0) = 0.623906 after 28 iters\n",
      "MPB solved for omega_1(0,1.89458,0) = 0.586052 after 7 iters\n",
      "MPB solved for omega_1(0,1.89447,0) = 0.586022 after 3 iters\n",
      "MPB solved for omega_1(0,1.89447,0) = 0.586022 after 1 iters\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dominant planewave for band 1: (0.000000,1.894467,0.000000)\n",
      "MPB solved for omega_1(0,2.04866,0) = 0.629115 after 27 iters\n",
      "MPB solved for omega_1(0,1.91383,0) = 0.591428 after 7 iters\n",
      "MPB solved for omega_1(0,1.91372,0) = 0.591398 after 3 iters\n",
      "MPB solved for omega_1(0,1.91372,0) = 0.591398 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,1.913720,0.000000)\n",
      "MPB solved for omega_1(0,2.06729,0) = 0.634325 after 24 iters\n",
      "MPB solved for omega_1(0,1.93307,0) = 0.596804 after 7 iters\n",
      "MPB solved for omega_1(0,1.93297,0) = 0.596774 after 3 iters\n",
      "MPB solved for omega_1(0,1.93297,0) = 0.596774 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,1.932968,0.000000)\n",
      "MPB solved for omega_1(0,2.08591,0) = 0.639536 after 29 iters\n",
      "MPB solved for omega_1(0,1.95232,0) = 0.60218 after 7 iters\n",
      "MPB solved for omega_1(0,1.95221,0) = 0.602151 after 3 iters\n",
      "MPB solved for omega_1(0,1.95221,0) = 0.602151 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,1.952212,0.000000)\n",
      "MPB solved for omega_1(0,2.10453,0) = 0.644749 after 30 iters\n",
      "MPB solved for omega_1(0,1.97156,0) = 0.607556 after 7 iters\n",
      "MPB solved for omega_1(0,1.97145,0) = 0.607527 after 3 iters\n",
      "MPB solved for omega_1(0,1.97145,0) = 0.607527 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,1.971452,0.000000)\n",
      "MPB solved for omega_1(0,2.12316,0) = 0.649962 after 24 iters\n",
      "MPB solved for omega_1(0,1.99079,0) = 0.612932 after 7 iters\n",
      "MPB solved for omega_1(0,1.99069,0) = 0.612903 after 3 iters\n",
      "MPB solved for omega_1(0,1.99069,0) = 0.612903 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,1.990687,0.000000)\n",
      "MPB solved for omega_1(0,2.14178,0) = 0.655177 after 25 iters\n",
      "MPB solved for omega_1(0,2.01002,0) = 0.618308 after 7 iters\n",
      "MPB solved for omega_1(0,2.00992,0) = 0.61828 after 3 iters\n",
      "MPB solved for omega_1(0,2.00992,0) = 0.61828 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.009918,0.000000)\n",
      "MPB solved for omega_1(0,2.16041,0) = 0.660393 after 26 iters\n",
      "MPB solved for omega_1(0,2.02924,0) = 0.623684 after 7 iters\n",
      "MPB solved for omega_1(0,2.02915,0) = 0.623656 after 3 iters\n",
      "MPB solved for omega_1(0,2.02915,0) = 0.623656 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.029145,0.000000)\n",
      "MPB solved for omega_1(0,2.17903,0) = 0.66561 after 29 iters\n",
      "MPB solved for omega_1(0,2.04847,0) = 0.62906 after 7 iters\n",
      "MPB solved for omega_1(0,2.04837,0) = 0.629032 after 3 iters\n",
      "MPB solved for omega_1(0,2.04837,0) = 0.629032 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.048367,0.000000)\n",
      "MPB solved for omega_1(0,2.19766,0) = 0.670828 after 31 iters\n",
      "MPB solved for omega_1(0,2.06768,0) = 0.634436 after 7 iters\n",
      "MPB solved for omega_1(0,2.06759,0) = 0.634409 after 3 iters\n",
      "MPB solved for omega_1(0,2.06759,0) = 0.634409 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.067585,0.000000)\n",
      "MPB solved for omega_1(0,2.21628,0) = 0.676047 after 29 iters\n",
      "MPB solved for omega_1(0,2.08689,0) = 0.639812 after 7 iters\n",
      "MPB solved for omega_1(0,2.0868,0) = 0.639785 after 3 iters\n",
      "MPB solved for omega_1(0,2.0868,0) = 0.639785 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.086799,0.000000)\n",
      "MPB solved for omega_1(0,2.2349,0) = 0.681268 after 33 iters\n",
      "MPB solved for omega_1(0,2.1061,0) = 0.645188 after 7 iters\n",
      "MPB solved for omega_1(0,2.10601,0) = 0.645161 after 3 iters\n",
      "MPB solved for omega_1(0,2.10601,0) = 0.645161 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.106008,0.000000)\n",
      "MPB solved for omega_1(0,2.25353,0) = 0.686489 after 36 iters\n",
      "MPB solved for omega_1(0,2.12531,0) = 0.650564 after 7 iters\n",
      "MPB solved for omega_1(0,2.12521,0) = 0.650538 after 3 iters\n",
      "MPB solved for omega_1(0,2.12521,0) = 0.650538 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.125213,0.000000)\n",
      "MPB solved for omega_1(0,2.27215,0) = 0.691711 after 27 iters\n",
      "MPB solved for omega_1(0,2.14451,0) = 0.65594 after 7 iters\n",
      "MPB solved for omega_1(0,2.14441,0) = 0.655914 after 3 iters\n",
      "MPB solved for omega_1(0,2.14441,0) = 0.655914 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.144414,0.000000)\n",
      "MPB solved for omega_1(0,2.29078,0) = 0.696935 after 32 iters\n",
      "MPB solved for omega_1(0,2.1637,0) = 0.661316 after 7 iters\n",
      "MPB solved for omega_1(0,2.16361,0) = 0.66129 after 3 iters\n",
      "MPB solved for omega_1(0,2.16361,0) = 0.66129 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.163611,0.000000)\n",
      "MPB solved for omega_1(0,2.3094,0) = 0.70216 after 37 iters\n",
      "MPB solved for omega_1(0,2.18289,0) = 0.666692 after 7 iters\n",
      "MPB solved for omega_1(0,2.1828,0) = 0.666667 after 3 iters\n",
      "MPB solved for omega_1(0,2.1828,0) = 0.666667 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.182803,0.000000)\n",
      "MPB solved for omega_1(0,2.32803,0) = 0.707385 after 33 iters\n",
      "MPB solved for omega_1(0,2.20208,0) = 0.672068 after 7 iters\n",
      "MPB solved for omega_1(0,2.20199,0) = 0.672043 after 3 iters\n",
      "MPB solved for omega_1(0,2.20199,0) = 0.672043 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.201991,0.000000)\n",
      "MPB solved for omega_1(0,2.34665,0) = 0.712612 after 32 iters\n",
      "MPB solved for omega_1(0,2.22126,0) = 0.677444 after 7 iters\n",
      "MPB solved for omega_1(0,2.22118,0) = 0.677419 after 3 iters\n",
      "MPB solved for omega_1(0,2.22118,0) = 0.677419 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.221175,0.000000)\n",
      "MPB solved for omega_1(0,2.36527,0) = 0.71784 after 28 iters\n",
      "MPB solved for omega_1(0,2.24044,0) = 0.68282 after 7 iters\n",
      "MPB solved for omega_1(0,2.24036,0) = 0.682796 after 3 iters\n",
      "MPB solved for omega_1(0,2.24036,0) = 0.682796 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.240355,0.000000)\n",
      "MPB solved for omega_1(0,2.3839,0) = 0.723069 after 25 iters\n",
      "MPB solved for omega_1(0,2.25962,0) = 0.688196 after 7 iters\n",
      "MPB solved for omega_1(0,2.25953,0) = 0.688172 after 3 iters\n",
      "MPB solved for omega_1(0,2.25953,0) = 0.688172 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.259531,0.000000)\n",
      "MPB solved for omega_1(0,2.40252,0) = 0.728299 after 31 iters\n",
      "MPB solved for omega_1(0,2.27879,0) = 0.693572 after 6 iters\n",
      "MPB solved for omega_1(0,2.2787,0) = 0.693548 after 3 iters\n",
      "MPB solved for omega_1(0,2.2787,0) = 0.693548 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.278703,0.000000)\n",
      "MPB solved for omega_1(0,2.42115,0) = 0.73353 after 37 iters\n",
      "MPB solved for omega_1(0,2.29795,0) = 0.698948 after 6 iters\n",
      "MPB solved for omega_1(0,2.29787,0) = 0.698925 after 3 iters\n",
      "MPB solved for omega_1(0,2.29787,0) = 0.698925 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.297870,0.000000)\n",
      "MPB solved for omega_1(0,2.43977,0) = 0.738762 after 36 iters\n",
      "MPB solved for omega_1(0,2.31711,0) = 0.704324 after 6 iters\n",
      "MPB solved for omega_1(0,2.31703,0) = 0.704301 after 3 iters\n",
      "MPB solved for omega_1(0,2.31703,0) = 0.704301 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.317033,0.000000)\n",
      "MPB solved for omega_1(0,2.45839,0) = 0.743995 after 28 iters\n",
      "MPB solved for omega_1(0,2.33627,0) = 0.7097 after 6 iters\n",
      "MPB solved for omega_1(0,2.33619,0) = 0.709677 after 3 iters\n",
      "MPB solved for omega_1(0,2.33619,0) = 0.709677 after 1 iters\n",
      "Dominant planewave for band 1: (0.000000,2.336193,0.000000)\n",
      "----------------------------------\n",
      "Parameters: radius=8.0\n",
      "Frequencies: [0.58064516 0.58602151 0.59139785 0.59677419 0.60215054 0.60752688\n",
      " 0.61290323 0.61827957 0.62365591 0.62903226 0.6344086  0.63978495\n",
      " 0.64516129 0.65053763 0.65591398 0.66129032 0.66666667 0.67204301\n",
      " 0.67741935 0.6827957  0.68817204 0.69354839 0.69892473 0.70430108\n",
      " 0.70967742]\n",
      "Transmitted fraction: [0.96811397 0.96457314 0.96193938 0.95956887 0.95759774 0.95577341\n",
      " 0.95425478 0.95320403 0.95259306 0.95241691 0.95263553 0.95321373\n",
      " 0.95415149 0.95549316 0.95714643 0.95902923 0.96114761 0.9633693\n",
      " 0.96559774 0.96780174 0.96998597 0.97210452 0.97392884 0.97520382\n",
      " 0.97587294]\n",
      "Reflected fraction: [5.13819077e-06 1.01448304e-06 1.44162393e-06 3.19704150e-06\n",
      " 5.57674063e-06 3.77329573e-06 1.60050053e-06 3.59522431e-06\n",
      " 3.02927753e-06 6.43372767e-07 9.81981707e-08 5.16602092e-07\n",
      " 1.56090295e-06 1.90144715e-06 2.31550685e-06 9.05986589e-07\n",
      " 2.35855658e-06 7.24506092e-06 1.22082455e-06 5.41079527e-06\n",
      " 1.17105358e-05 4.14269422e-06 2.15620830e-05 2.37287919e-05\n",
      " 2.37353074e-04]\n"
     ]
    }
   ],
   "source": [
    "# Every mode monitor measures the power flowing through it in either the forward or backward direction\n",
    "eig_mode1 = sim.get_eigenmode_coefficients(mode1, [1], eig_parity=mp.NO_PARITY)\n",
    "eig_mode2 = sim.get_eigenmode_coefficients(mode2, [1], eig_parity=mp.NO_PARITY)\n",
    "\n",
    "# First, we need to figure out which direction the \"dominant planewave\" k-vector is\n",
    "# We can pick the first frequency (0) for that, assuming that for all simulated frequencies,\n",
    "# The dominant k-vector will point in the same direction.\n",
    "k1 = eig_mode1.kdom[0]\n",
    "k2 = eig_mode2.kdom[0]\n",
    "\n",
    "# eig_mode.alpha[0,0,0] corresponds to the forward direction, whereas\n",
    "# eig_mode.alpha[0,0,1] corresponds to the backward direction\n",
    "\n",
    "# For port 1, we are interested in the +x direction, so if k1.x is positive, select 0, otherwise 1\n",
    "idx = (k1.x < 0) * 1\n",
    "p1_thru_coeff = eig_mode1.alpha[0,:,idx]\n",
    "p1_reflected_coeff = eig_mode1.alpha[0,:,1-idx]\n",
    "\n",
    "# For port 2, we are interestred in the -y direction\n",
    "idx = (k2.y > 0) * 1\n",
    "p2_thru_coeff = eig_mode2.alpha[0,:,idx]\n",
    "p2_reflected_coeff = eig_mode2.alpha[0,:,1-idx]\n",
    "\n",
    "# transmittance\n",
    "p2_trans = abs(p2_thru_coeff/p1_thru_coeff)**2\n",
    "p2_reflected = abs(p1_reflected_coeff/p1_thru_coeff)**2\n",
    "\n",
    "print(\"----------------------------------\")\n",
    "print(f\"Parameters: radius={ring_radius:.1f}\")\n",
    "print(f\"Frequencies: {mode1_freqs}\")\n",
    "print(f\"Transmitted fraction: {p2_trans}\")\n",
    "print(f\"Reflected fraction: {p2_reflected}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "S21 = p2_thru_coeff/p1_thru_coeff\n",
    "S11 = p1_reflected_coeff/p1_thru_coeff\n",
    "\n",
    "S21_mag = np.abs(S21)\n",
    "S21_phase = np.unwrap(np.angle(S21))\n",
    "S11_mag = np.abs(S11)\n",
    "S11_phase = np.unwrap(np.angle(S11))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAI4CAYAAABUVDNLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3iUdbrG8e+THgwltNBDBym2AAZ1j8SKZW1rw7JYsa7urq66ukXPNnfddlTUFRs2svayixUBRUEgFAHpCEiR3kIJSeY5f8yAIQQIJJl3Jrk/1zVX5m0zdwI8vPnNr5i7IyIi0ZcQdAARkbpKBVhEJCAqwCIiAVEBFhEJiAqwiEhAVIBFRAKiAiwiEhAVYKnzzOwEM/vCzDaZ2Xoz+9zM+prZWWY2zsw2mtl3ZjbMzOqXue7iyHXbzGxMgN+CxCnTQAypy8ysAbAUuAl4BUgBfgB8B/QC1gOfAqnAy8ASd78xcu0pQGOgO3CSuw+Idn6JbyrAUqeZWR/gY3dvVIlzLwAecPfe5fZfB1yhAiwHS00QUtfNA0rNbLiZnWFmmfs593+AWVHKJXWACrDUae6+GTgBcGAYsMbM3jGzrLLnmdmpwGDgN9FPKbWVCrDUee4+292vcvc2hNt9WwH/3HXczHIJt/9e6O7zAooptZAKsEgZ7j4HeI5wIcbMjgbeAa5x91EBRpNaSAVY6jQz625md5hZm8h2W2AQMMHMegHvAz9x93cruDbRzNKAJCDBzNLMLDma+SW+qQBLXbcFOBb40sy2AhOAmcAdkUcz4GkzK4w8yn4IdyWwHXiccNe17YTbkUUqRd3QREQCojtgEZGABFqAzWygmc01swVmdk8Fx83MHo4c/8rMjqnstSIisS6wAmxmicBQ4AygBzDIzHqUO+0MoEvkMYRwW1tlrxURiWlB3gH3Axa4+yJ33wnkA+eWO+dc4HkPmwA0MrOWlbxWRCSmJQX43q2Bb8tsLyP8afSBzmldyWsBMLMhhO+eSU9Pz2nbtu0ex0OhEAkJsd8UrpzVJx4yQnzkjIeMEHzOefPmrXX3ZuX3B1mArYJ95btk7Oucylwb3un+JPAkQJ8+fXzy5Ml7HB8zZgwDBgw4UNbAKWf1iYeMEB854yEjBJ/TzJZUtD/IArwMKHs72gZYUclzUipxrYhITAvyd4dJQBcz62BmKcClhId8lvUO8ONIb4hcYJO7r6zktSIiMS2wO2B3LzGzW4EPgETgGXefZWY3Ro4/AYwEzgQWANuAq/d3bQDfhojIIQuyCQJ3H0m4yJbd90SZ5w7cUtlrRUTiSex/fCkiUkupAIuIBEQFWEQkICrAIiIBUQEWEQmICrCISEBUgEVEAqICLCISEBVgEZGAqACLiAREBVhEJCAqwCIiAVEBFhEJiAqwiEhAVIBFRAKiAiwiEhAVYBGRgKgAi4gERAVYRCQgKsAiIgFRARYRCYgKsIhIQFSARUQCogIsIhIQFWARkYCoAIuIBEQFWEQkICrAIiIBUQEWEQmICrCISEBUgEVEAqICLCISEBVgEZGAqACLiAREBVhEJCCBFGAza2xmH5nZ/MjXzH2cN9DM5prZAjO7p8z+h8xsjpl9ZWZvmlmj6KUXEakeQd0B3wOMcvcuwKjI9h7MLBEYCpwB9AAGmVmPyOGPgF7ufgQwD/hlVFKLiFSjoArwucDwyPPhwHkVnNMPWODui9x9J5AfuQ53/9DdSyLnTQDa1HBeEZFqZ+4e/Tc12+jujcpsb3D3zHLnXAgMdPfrIttXAse6+63lznsX+Le7v7iP9xoCDAHIysrKyc/P3+N4YWEhGRkZ1fBd1SzlrD7xkBHiI2c8ZITgc+bl5RW4e5/y+5Nq6g3N7GOgRQWH7qvsS1Swb4//LczsPqAEeGlfL+LuTwJPAvTp08cHDBiwx/ExY8ZQfl8sUs7qEw8ZIT5yxkNGiN2cNVaA3f2UfR0zs1Vm1tLdV5pZS2B1BactA9qW2W4DrCjzGoOBs4GTPYjbeBGRKgqqDfgdYHDk+WDg7QrOmQR0MbMOZpYCXBq5DjMbCNwNnOPu26KQV0Sk2gVVgB8ETjWz+cCpkW3MrJWZjQSIfMh2K/ABMBt4xd1nRa5/FKgPfGRm08zsiWh/AyIiVVVjTRD74+7rgJMr2L8COLPM9khgZAXnda7RgCIiUaCRcCIiAVEBFhEJiAqwiEhAVIBFRAKiAiwiEhAVYBGRgKgAi4gERAVYRCQgKsAiIgFRARYRCYgKsIhIQFSARUQCogIsIhIQFWARkYCoAIuIBEQFWEQkICrAIiIBUQEWEQmICrCISEBUgEVEAqICLCISEBVgEZGAqACLiAREBVhEJCAqwCIiAVEBFhEJiAqwiEhAVIBFRAKiAiwiEhAVYBGRgKgAi4gERAVYRCQgKsAiIgFJOtAJZtYcOB5oBWwHZgKT3T1Uw9lERGq1fRZgM8sD7gEaA1OB1UAacB7QycxeA/7m7pujEVREpLbZ3x3wmcD17r60/AEzSwLOBk4FXj/YNzWzxsC/gfbAYuBid99QwXkDgf8DEoGn3P3BcsfvBB4Cmrn72oPNISISpH22Abv7LyoqvpFjJe7+lrsfdPGNuAcY5e5dgFGR7T2YWSIwFDgD6AEMMrMeZY63JfwfQIUZRURi3X4/hDOzE83siMjzi83sUTP7mZmlVvF9zwWGR54PJ9ysUV4/YIG7L3L3nUB+5Lpd/gHcBXgVs4iIBMLcK65fZjYUOAJIBeYBGcD7wHFAortffshvarbR3RuV2d7g7pnlzrkQGOju10W2rwSOdfdbzewc4GR3v93MFgN99tUEYWZDgCEAWVlZOfn5+XscLywsJCMj41C/lahRzuoTDxkhPnLGQ0YIPmdeXl6Bu/fZ64C7V/gAvo58TQPWES66AAbM2Nd1Za7/mHCPifKPc4GN5c7dUMH1FxFu9921fSXwCFAP+BJoGNm/GGh6oDzuTk5Ojpc3evTovfbFIuWsPvGQ0T0+csZDRvfgcxLuObZXTdrfh3A7IgV6h5ktcffSyLabWfGBKr67n7KvY2a2ysxauvtKM2tJuIdFecuAtmW22wArgE5AB2C6me3aP8XM+rn7dwfKJSISK/ZXgJub2c8J3/Huek5ku1kV3/cdYDDwYOTr2xWcMwnoYmYdgOXApcBl7j4LaL7rpAM1QYiIxKr9fQg3DKhPuO131/Nd209V8X0fBE41s/mEezI8CGBmrcxsJIR7WgC3Ah8As4FXIsVXRKRW2OcdsLs/UFNv6u7rgJMr2L+CcP/jXdsjgZEHeK321Z1PRCQa9jcS7uH9Xejut1V/HBGRumN/TRAFkUcacAwwP/I4Ciit+WgiIrXb/poghgOY2VVAnrsXR7afAD6MSjoRkVqsMtNRtiL84dsuGZF9IiJSBQecjpJwD4WpZjY6sn0icH+NJRIRqSMOWIDd/Vkzew84NrLrHg14EBGpuv31gmjv7osBIgX37XLHDWjt7stqNKGISC21vzvgh8wsgXDhLQDWEO4R0RnII9yP97eEhwyLiMhB2l8viIsi8+9eDlwDtAS2ER6VNhL4g7vviEpKEZFaaL9twO7+NXBflLKIiNQpWhVZRCQgKsAiIgFRARYRCcgBC7CFXWFmv4lstzOzfjUfTUSkdqvMHfBjQH9gUGR7C+HVikVEpAoqMxT5WHc/xsymArj7BjNLqeFcIiK1XmXugIvNLJHI8u9m1gwI1WgqEZE6oDIF+GHgTcLrwv0BGAf8sUZTiYjUAZWZjOclMysgPPTYgPPcfXaNJxMRqeX2NxlP4zKbq4ERZY+5+/qaDCYiUtvt7w64gHC7rwHtgA2R542ApUCHGk8nIlKL7bMN2N07uHtHwsvC/9Ddm7p7E+Bs4I1oBQzSxG/WMXT0AgqWbAg6iojUQpXphtbX3W/cteHu75nZ72owU0woWLKBS5+cQMghOdF49up+nNC5adCxRKQWqUwviLVm9isza29m2WZ2H7CupoMFbcKidbiHnxeXOlc9M5Gf5k9l7Lw1lIY82HAiUitU5g54EOGJ19+MbH/K96Piaq3cjk1ITU6guCREYmICA7o245M5q3lr2gqa10/l3KNaccExbdi2s5QJi9aR27EJOdmZQccWkThSmW5o64Hbo5AlpuRkZ/LSdbl7FNeiklI+mb2aN6Yu59nPFzPss2+wyPmpyQm8dF2uirCIVNoBC3BkNeS9fud295NqJFEMycnO3KOgpiYlckbvlpzRuyXrt+7krtem8/Hs1QDsKA4xavYqFWARqbTKNEHcWeZ5GvAjoKRm4sSPxoelcNOAzoxbsJai4hAOvDBhCUe0acjAXi2DjicicaAyTRAF5XZ9bmZjayhPXCnbTNGmUTrDxi3ixhencMHRrfntOT1pmJ4cdEQRiWGVaYIoOyIuAcgBWtRYojhTtpnizCNa8sgnCxg6egHjF63joQuP5IQu6romIhWrTBNE2RFxJcA3wLU1GSpeJScm8PNTu3Jy9+b87JVpXPH0lwzun83pPVsw9duN6ikhInuoTAE+vPzy82aWWkN5aoUj2zZi5G0/4M/vz+HZzxfz/PglmEFKknpKiMj3KjMQ44sK9o2v7iC1TVpyIr/9YU8u6dMGB0IOO0tCTFhU68ewiEgl7W82tBZAayDdzI6G3V1eGwD1opCtVri4bzvemraCopIQIYeURK2DKiJh+2uCOB24CmgD/L3M/i3AvTWYqVbJyc7k5etzGT1nNSNnrOAvH8yhXZN6nN5Tn2OK1HX7mw1tuLvnAVe5e16ZxznuXqXZ0MyssZl9ZGbzI18rbBQ1s4FmNtfMFpjZPeWO/SRybJaZ/aUqeWpaTnYmd57ejTdvOYGerRpy80tTeHf6iqBjiUjA9lmAzeyKyNP2Zvbz8o8qvu89wCh37wKMimyXf/9EwqsvnwH0AAaZWY/IsTzgXOAId+8J/LWKeaKiYXoyL153LDntMrk9fypvTFkWdCQRCdD+GiQPi3zNAOpX8KiKc4HhkefDgfMqOKcfsMDdF7n7TiA/ch3ATcCD7l4E4O6rq5gnajJSk3jumr7kdmzCHa9OJ3/i0qAjiUhAzD36Uyua2UZ3b1Rme4O7Z5Y750JgoLtfF9m+EjjW3W81s2nA28BAYAdwp7tP2sd7DQGGAGRlZeXk5+fvcbywsJCMjIzq++YqaWep8+jUIr5aW8oVh6dwSvb+R80FlfNgxUPOeMgI8ZEzHjJC8Dnz8vIK3L1P+f2VGQnXDLgeaF/2fHe/5gDXfUzFI+buO9B77nqJCvbt+t8iCcgEcoG+wCtm1tEr+N/E3Z8EngTo06ePDxgwYI/jY8aMofy+aDnxxFJufXkqL369isRGLWneIG2fgzWCzHkw4iFnPGSE+MgZDxkhdnNWZiDG28BnwMdAaWVf2N1P2dcxM1tlZi3dfaWZtSS86Gd5y4C2ZbbbACvKHHsjUnAnmlkIaAqsqWy+WJCalMhjlx/D4KcnMnz8EgxNaylSl1SmU2o9d7/b3V9x99d3Par4vu8AgyPPBxMu8uVNArqYWQczSwEujVwH8BZwEoCZdQVSgLVVzBSI5MQE+nduAoRv74uKNVhDpK6oTAH+j5mdWc3v+yBwqpnNB06NbGNmrcxsJIC7lwC3El4UdDbwirvPilz/DNDRzGYS/nBucEXND/HiuE5NSUsO/1E4ENKSRyJ1QmWaIG4H7jWzIqCYcNusu3uDQ31Td18HnFzB/hXAmWW2RwIjKzhvJ3BF+f3xate0luPmr+G9md/xyCcL6NO+Mf07NQk6mojUoAPeAbt7fXdPcPd0d28Q2T7k4isVy8nO5PZTupI/JJfsJvUY8vxkvl6xOehYIlKDDliAzeyYCh6dzKwyd89ykBrVS2H4Nf3ISEti8LMT+Xb9tqAjiUgNqUwb8GPABGBY5DGBcLvrPDM7rQaz1VmtGqUz/Jp+FBWXMviZiawrLAo6kojUgMoU4MXA0e6e4+45wFHATOAUIKbnYIhnXbPq88xVfVm+cTvXDJ9MUYk+mBOpbSpTgLuX6X2Au39NuCAvqrlYAtCnfWMeGXQ0M5ZtZOi0IopLQ0FHEpFqVJkCPNfMHjezEyOPxwg3P6QS7hUhNei0ni34w/m9+WptKfe8PoM47m0nIuVU5oO0q4CbgZ8S7oI2jvBS9cVAXo0lk90G9WvHxK/m8PqUZYTc6dw8Q+vLidQClVmWfjvwt8ijvMJqTyQVOqdTMltTGvPm1OUasixSS1SmG1oXM3vNzL42s0W7HtEIJ98zM45o0xAIj5bT+nIi8a8ybcDPAo8TXpI+D3geeKEmQ0nF+ndqSlpS+I8s5NAoff9TWIpIbKtMAU5391GE5w5e4u73E5kIR6IrJzuTl67P5da8zrRqmMZfPpjLgtVqBRKJV5UpwDvMLAGYb2a3mtn5QPMaziX7sGt9ufwh/UlONAY/M5HVm3cEHUtEDkFlCvBPCS9DfxuQA1zJ91NJSkDaNanHs1f1Y8O2nVz93CQKi0qCjiQiB6kyk/FMcvdCd1/m7le7+wXuPiEa4WT/erdpyNDLj2HOd1u4+aUpGqghEmf22Q3NzN7Z1zEAdz+n+uPIwcrr1pw/nd+bu17/inten8FfLzoCs4pWcxKRWLO/fsD9gW+BEcCXVLxGm8SAi/u2ZcWm7fzz4/m0apTGHad1CzqSSK008Zt1TFq8odoGQu2vALcgvFrFIOAy4L/AiLLzQkjsuP3kLqzcuINHPllAcUmI+unJGi0nUo0+mPUdN75QgBmkJFXPQKh9FmB3LwXeB96PzPswCBhjZv/r7o9U6V2l2pkZvz+/F/NWb+GJTxdptJxINSouDfHAO7NwwB2KIwOhqvpva78fwplZqpldALwI3AI8DLxRpXeUGpOcmMCJXZoBGi0nUp3+8v4cVmzaQXKikWiQnJRAbseqLxm2vw/hhgO9gPeAB9x9ZpXfTWrcD7o24/GxCykqCRFyyG5SL+hIInHto69XMeyzb7gyN5vzjm7NhEXrotIGfCWwFegK3Fbmk/UqL8opNScnO5OXr8/lvZkr+fekb/nL+3M5tkMTmtVPDTqaSNxZtmEbd746nZ6tGnDfWYeTlpxYrU16+2yCiCzEWX/XIpxlHlqUM8blZGfyq7N68Pw1/Vi9ZQfXDp/EVg3UEDkoO0tC3PLyVEIh57HLjyEtObHa36MyI+EkTh3dLpOhlx3DzOWbuOVlDdQQORh/fn8O07/dyJ8vPILsJofVyHuoANdyJx+exR/O782YuWu49w2tqCFSGR/M+o6nx33D4P7ZnNm7ZY29j5aWrwMG9WvHyk07eHjUfFo2TOPnGqghsk/frg+3+x7RpiH3nnV4jb6XCnAd8bNTuvDdpu08/MkCWjRM57Jj2wUdSSTm7CwJcevLUwB4dNAxpCZVf7tvWSrAdYSZ8Yfze7NmSxG/emsGzeunckqPrKBjicSUP46czfRlm3jiihzaRaELp9qA65DkxASGXn4MvVs35NYRUxgxcSlDRy+gYMmGoKOJBO69GSt57ovFXH18ewb2ahGV91QBrmPqpSTx9FV9aZSezC/fmMHfPpzL5U9NUBGWOm3pum3c9dpXHNmmIb88o2bbfctSAa6Dmmak8sMjWwHhteWKNWRZ6rCiklJueXkKZvDoZceQkhS9sqgCXEcN7NVy91+0EHBk20bBBhIJyB/+O5sZyzfx14uOpG3j6A7dVwGuo3KyMxlxfS4/OqY1Bvxr7EJ2lmightQt//1qJc+PX8K1J3TgtJ7RafctS70g6rCc7ExysjM5tkOTyIoaX/G3i4/UihpS6xUs2cD7M1fy4oSlHNW2EXcP7B5IDhVg4eK+bflu8w7+/tE8WjRM466A/jKKREPBkg1cPmwCOyK/8d14YseotvuWpQIsAPzkpM58t3kHj41ZSFaDNAYf1z7oSCI14vMFa3YX3wSDhWu2BpYlkLJvZo3N7CMzmx/5WuH8bmY20MzmmtkCM7unzP6jzGyCmU0zs8lm1i966WsnM+N/z+nJKYdncf+7s3h/5sqgI4lUu50lIT6dtxYIz6ubUk0Tqx+qoD6EuwcY5e5dgFGR7T2YWSIwFDgD6AEMMrMekcN/ITxJ/FHAbyLbUkVJiQk8Muhojm7biNvypzHxm/VBRxKpNuHpJacweckGrv9BB+48vVvgS3YFVYDPBYZHng8HzqvgnH7AAndf5O47gfzIdRBecWfXnMQNgRU1mLVOSU9J5OnBfWmTmc51wycxf9WWoCOJVFlJyPnJiCl89PUqHjinJ/ed1YNb8joHvl6iBTE9oZltdPdGZbY3uHtmuXMuBAa6+3WR7SuBY939VjM7HPiA8G8RCcBx7r5kH+81BBgCkJWVlZOfn7/H8cLCQjIyMqrvm6sh0c65ZluI33+5g0SDKw9PZnmh071xIp0z9z85STz8POMhI8RHznjIWBJyHi3YyrR1xuXdUzi1fXLUM+Tl5RW4e5/y+2vsQzgz+5jw0vbl3VfZl6hg367/LW4Cfubur5vZxcDTwCkVvYi7Pwk8CdCnTx8fMGDAHsfHjBlD+X2xKIichx+5iYueGM/DU3dGluIuPeCvbPHw84yHjBAfOWM9Y0lpiNvzpzFt3TZ+ddbhXPeDjkFH2kONNUG4+ynu3quCx9vAKjNrCRD5urqCl1gGtC2z3YbvmxoG8/3qzK8Sbq6QatardUN+eGRLHA1ZlvhTUhriZ69M578zVnJJt5SYK74QXBvwO4SLKJGvb1dwziSgi5l1MLMU4NLIdRAuxCdGnp8EzK/BrHXaJX3bkZwY/mXEgX4dGgcbSKQSSkPOHa9O593pK/jlGd05o0P0mx0qI6gC/CBwqpnNB06NbGNmrcxsJIC7lwC3Em7rnQ284u6zItdfD/zNzKYDfyTSxivVLyc7k/wh/Tmpe3NCDm9MWa5ljSSmlYacO1+dztvTVnDXwG7ccGKnoCPtUyADMdx9HXByBftXAGeW2R4JjKzgvHFATk1mlO/lZGfyzFV9+fP7c3h8zEKa10/lZ6d2DTqWyF5KQ84vXpvOm1OX84vTu3HzgM5BR9ovjYSTSrvr9G6s3VLE/42aT9P6qVyZmx10JJHdQiHn7te/4o0py/n5qV25JS+2iy+oAMtBMDP+dEFv1m/dyW/enkmTw1JqdMVYkcqavHg9D7z7NTOWb+Knp3ThtpO7BB2pUjQdpRyUpMQEHr3sGI5u24if5k/ji4Vrg44kddzkxeu55MkJzFi+iaQE4wddmgUdqdJUgOWgpack8sxVfWnXpB5Dni9g5vJNQUeSOqq4NMRv355FaSj8wbC7x1VXSRVgOSSN6qXw/DX9qJ+WxFXPTmLpum1BR5I6ZkdxKTe9WMCslZtJSjASDZIDnlznYKkNWA5Zq0bpvHBtPy58YjxXPvMlvzm7B+8t3En9DhsCH2MvtduWHcVcN3wyExev53fn9qRHq4ZMWLSO3I5N4urvngqwVEnn5vV5enBfBj05nuuenwwO/1k8IfBZpqT2WldYxFXPTmL2ys3885KjOPeo1gBx+fdNTRBSZTnZmfzwyFa4h0fLaciy1JQVG7dz0b/GM2/VFob9uM/u4huvVIClWlx2bPYeQ5aPaadVlqV6LVxTyIWPf8GazUW8cO2x5HVvHnSkKlMBlmqxa8hyn6xEQg7DPvtGqyxLtZm5fBMXPzGenaUhRgzJrTVzkqgAS7XJyc7k1qPT+P15vfhkzmpuz59KSamKsFTNhEXruPTJCaQlJ/LqjcfRq3XDoCNVGxVgqXZX5Gbz67N78N7M77jj1em7+2iKHKxRs1cx+JmJtGiYxms39adD08OCjlSt1AtCasS1J3RgR3EpD30wl7SkRP50QW8SEiqaY1+kYm9NXc4dr06nV6sGPHt1PxoflhJ0pGqnAiw15pa8zhQVl/LwJwtITU7ggXN6YqYiLPtXsGQDj41ewKg5qzmuUxOe/HEfMlJrZ6mqnd+VxIyfndqVHSUhnvx0EWnJifzyjO4qwrJPBZF5HUpCToLBT07qXGuLL6gASw0zM355Rnd2FJeGi3BSAj8/rVvQsSQGbdlRzN1vzKAk8pmBAVOWbqR/p6bBBqtBKsBS48yM+3/Yk6LiEA9/soC1W3fSulF63A0blZozf9UWbnixgMVrt5KUYLh73M3rcChUgCUqEhKMP17QmxWbtvPyl0sxIDU5QUOWhfdmrOTOV6eTnpLIy9fnkpyYEJfzOhwKFWCJmsQEo1+Hxnw2fy0OFBWHGL9wba3/RyYVKykN8dCHc/nX2EUc3a4Rj1+eQ4uGaUB8zutwKNQPWKLquE5NSUtOwAgPWZ64eD3FGqxR56wrLGLwsxP519hFXJHbjvwhubuLb12iO2CJqpzsTF66LpcJi9ayfOMOXv5yKdcNn8zjVxxDvRT9dawLvlq2kRtfKGDt1p08dOERXNSnbdCRAqO/8RJ1OdmZu3/FPKJ1Q+59cwaDhn3Js1f1rZWd7eV7r0z6ll+9PZNmGam8fuNx9G5Te4YVHwo1QUigLu3XjieuyGHOys1c+MQXLN+4PehIUgOKSkr55RszuOv1r+jXvjHv/uSEOl98QQVYYsBpPVvwwrXHsmZLET967AvmrdoSdCSpJgVLNvDge7M56+FxjJi4lJsHdGL4NbVzWPGhUAGWmNCvQ2NevbE/IXcufPwLJi9eH3QkqaKCJRsY9OQEnhi7iAWrC7nr9G7cNbA7iZoTZDcVYIkZ3Vs04PWbjqNpRiqXP/Ul/xq7kKGjF1CwZEPQ0eQglZSG+OsHc9gZ6eGSYOFeL7InfQgnMaVt43q8emN/Lv7XeP703hwN2IhD323awW35U5n4zXoSLdzhsC6MajsUKsASc5pkpPLDI1rxz1Hzdw/Y+HyBBmzEg9FzV3PHK9PZUVzK3y8+kuwmh9WZUW2HQgVYYtIPujbjiU8XUlQcwoH/zljJJX3bktWg7nXWjwfFpSH+GhnV1r1FfR697Bg6N88A6s6otkOhAiwx6fsBG+soDTlPjF3I2Y+MY+hlx9Sa9cBqi2UbtvGTEVOZunQjlx3bjt+c3YO05MSgY8UFFWCJWWUHbAzs1YIbXijgsmETuPfMw7n6+PaaVzgGTFlVwu0Pj6M05Dwy6MYKxVEAACAASURBVGh+eGSroCPFFfWCkLjQNas+b996PHndm/O///man/57Gtt2lgQdq84qKinlgXdn8fDUIto1rsd/bztBxfcQ6A5Y4kaDtGT+dUUOj49dyF8/nMvc77bwxBU5tK9lCzXGsoIlG3h/5kpGzV7NorVbOTU7iUev709qkpocDoUKsMSVhATjlrzO9GrdkNvzp/LDR8fxf5cexUnds4KOVusVLF7PpcMmUFwa7tF718Bu9GCZim8VqAlC4tKJXZvx7q0n0K5xPa55bjJ3vTadRz+Zr0EbNWTpum3c8er03cU3wcA1sqLKAinAZtbYzD4ys/mRrxX2UzGzZ8xstZnNPJTrpXZr27ger990HHndmvHK5GX89cN5XDZsgopwNSouDfH4mIWc9s+xrNpcRFKCkWiQooEV1SKoO+B7gFHu3gUYFdmuyHPAwCpcL7VcWnIifdpnsqs/RFFJiD+/P5vtO0sDzVUbTPt2Iz98ZBx/fn8OJ3Ztxid3nsi/b+jPz0/rppGJ1SSoNuBzgQGR58OBMcDd5U9y90/NrP2hXi91Q27HpqQmL6C4JAQYE7/ZwGn/HMvvz+vNiV2bBR0v7hQWlfDXD+YyfPximtdP5YkrchjYqwUALRumq/BWI/MAGnLMbKO7NyqzvcHd99UM0R74j7v3OsTrhwBDALKysnLy8/P3OF5YWEhGRkYVvpvoUM79W7ChlDnrS+neOJFSh2dnFfHdVie3ZSKDuqfSMPX7PsP6We7blFUlvDh7Jxt2OCe1S+LCrimkJ+27v7V+lpWTl5dX4O599jrg7jXyAD4GZlbwOBfYWO7cDft5nfbAzHL7Kn192UdOTo6XN3r06L32xSLlPDg7ikv8Hx/N9S73jvTev33fX/5yiZeWhtw9djIeSDRzrty43W94frJn3/0fP/0fY71gyfpKXaefZeUAk72CmlRjTRDufsq+jpnZKjNr6e4rzawlsPogX76q10stl5qUyE9P6crZR7Tivjdn8Ms3ZvDGlGVcmZvNpwt3Ur/DBv0qDUxevJ4nxi7k8wXrCLlz18BuXP+DjiQnqoNUNAT1U34HGBx5Phh4O8rXSx3RuXkG+UNy+cuFR/D1is3clj+N1+YX1/neEjuKS/nze7O56InxfDx7NUUl4dnLbh7QWcU3ioL6ST8InGpm84FTI9uYWSszG7nrJDMbAYwHupnZMjO7dn/Xi1TEzLi4T1uuOr797n1FJSF+8/ZMZi7fFFywAKzZUsQ/PprH8Q9+wuNjF+2eJN2Axeu2BRmtTgqkF4S7rwNOrmD/CuDMMtuDDuZ6kf05qXsWT4/7hqLiEAkJxoLVhZz9yDj6ts/k6uM7cFqPLJJq6d3fvFVbePqzb3hz2nJ2loQ45fDmnNClKQ++N4fikpAmTA+IhiJLnbFrissRH09i0Cl96dw8g1cnf8tzXyzm5pem0LpROlf2z+bSvm1pVC/+F410dz6bv5anxn3Dp/PWkJacwMV92nD18R3o1CzcI6B360aaMD1AKsBSp+RkZ7KlU8ruYnPdDzpy9fEdGDV7Fc9+vpgH35vDPz+ex/lHt+HYDo1ZvnF7XBWngiUbGDd/DUUlpXw8ezXzVhXSrH4qvzi9G5f1a0dmudWIy075KdGnAix1XmKCcVrPFpzWswWzV27muc8X81rBt4yYuBSApATj/nN6cknftjH7AdW6wiJe+nIp/zdqPqWhcMtudpN6/O2iIzn7yJaaMCdGqQCLlHF4ywb8+cIjaFY/laGjF+BAScj51Vsz+ePI2fTr0JjjOzXl+M5N6d6iPgkBLbG+taiEiYvX88WCtYxbsI7ZKzfvcTzB4OI+bfhRTptA8knlqACLVCCve3OeGreI4pIQSYkJ3H5yZ1ZuKuLzBWsZM3c2AI0PS6F/pyYc36kpDdOTWLxuK7kdm1b7r/QFSzbwxcK1NK6XwtrCnXy+YC1Tv91AcamTkphATnYmvzi9G40PS+GBd2eV+VCtabXmkOqnAixSgbJr0pVvA165aTufL1jHFwvW8vnCtfz3q5W7jxnzOKJNQzo3r0/zBqk0r59K8/ppezz/euXm3a/bq3UDNm4rZv3WnWzYupP128JfN2wrZsa8Ip5e+CWfL1hLqMyMAb1bN+TaEzpyfOcm9MluTHrK980LXbPq60O1OKICLLIP+/qAqmXDdC7MacOFOW1wd3733695dtxiHHBg1eYdrNpcxNrCIkpChz7XSnoSpCRv2l18EwxuyevMHad1O+jMEptUgEWqwMw4q3crXv5y6e5f/YdenkNOdiahkLN+207WbCli9ZYiVm/ewTvTV/DZ/LXha4ETOjdlYO8WNK6XQqN6KTQ+LIXMw5JplJ7CF+M+pX6HI7n8qQm7X3tAt+bBfsNSrVSARapoX80VCQlG04xUmmakcnjL8Lkdm2UwafH63QX1p6d23e8d6/6aQiT+qQCLVIPK/up/KAVVzQq1lwqwSJSpoMousdmrXESkDlABFhEJSCBLEgXFzNYAS8rtbgqsDSDOwVLO6hMPGSE+csZDRgg+Z7a777VAYZ0qwBUxs8le0VpNMUY5q088ZIT4yBkPGSF2c6oJQkQkICrAIiIBUQGGJ4MOUEnKWX3iISPER854yAgxmrPOtwGLiARFd8BSZ5hZbzNbYmY37eecdDMba2ZVnsHczFLM7FMz04AnqZAKsNQZ7j4DuBT48X5OuwZ4w91Lq+H9dgKjgEuq+lpSO6kAS12zGui5n+OXA28DmFl7M5u564CZ3Wlm90f2zzGzp8xsppm9ZGanmNnnZjbfzPqVeb23Iq8pshcVYKlrHgRSzSy7/AEzSwE6uvviSrxOZ+D/gCOA7sBlwAnAncC9Zc6bCfStYmappVSApc4ws4HAYcB/qfguuCmwsZIv9427z3D3EDALGOXhT7RnAO13nRRpythpZvWrkl1qJxVgqRPMLA34C3Az4SLZq4LTtgNpZbZL2PPfSNljRWWeh8psh9h7lsFUYMfBp5baTgVY6opfAc9HmhcqLMDuvgFIjBRrgFVAczNrYmapwNkH+6Zm1gRY4+7Fh5xcai0VYKn1zKwbcCrwz8iufd0BA3xIuC2XSNH8X+BL4D/AnEN4+zxg5CFcJ3WABmKIlGFmRwM/d/crq+n13gB+6e5zq+P1pHbRHbBIGe4+FRhdXQMxgLdUfGVfdAcsIhIQ3QGLiAREBVhEJCAqwCIiAVEBFhEJiAqwiEhAVIBFRAKiAiwiEhAVYBGRgKgAi4gERAVYRCQgdWqxwIEDB/ratWuDjiEidUxBQcEH7j6w/P46VYABJk+evMf2mDFjGDBgQDBhDoJyVp94yAjxkTMeMkLwOc2swv11qglCd78iEpCmFe2sUwVYRCSWqACLiAREBVhEJCAqwCIiAVEBFhEJiAqwiEhAVID3o2DJBoaOXkDBkg1BRxGRWqjODcSorIIlGxj05HhKQk5KUgIvXZdLTnZm0LFEpBbRHfA+jF+4lp2lTshhZ0mICYvWBR1JRGoZFeB96N+pKalJ4R9PyKFBmn5ZEJHqpQK8DznZmbx8fS43D+hE60Zp/Om9OUz8Zn3QsUSkFlEB3o+c7EzuGtidN28+nhYN07jq2YlMWqwiLCLVQwW4Epo3SCP/+lxaNEjjqmcmUrBERVhEqk4FuJKaN0hjxJBcmjdIY/Azk9Q1TUSqTAX4IGQ1SGPE9bk0zUhh8DMTmbJURVhEDp0K8EFq0TB8J9wkI4XBT08kf9JSDdYQkUOiAnwIWjZMZ8T1udRLSeSe12fwtw/ncvlTE1SEReSgxHwBNrM0M5toZtPNbJaZPRDZf7+ZLTezaZHHmdHM1apROucd0xpAgzVE5JDEfAEGioCT3P1I4ChgoJnlRo79w92PijxGRjvYaT1a7DFYo/FhKdGOICJxLOYLsIcVRjaTIw8PMNJuuwZr3HhiR1o0SOUP/52tZggRqTRzj4latl9mlggUAJ2Boe5+t5ndD1wFbAYmA3e4+17Vz8yGAEMAsrKycvLz8/c4XlhYSEZGRpUzrt8R4s8Td7CpyLmjTxpdMhOr/JplVVfOmhYPOeMhI8RHznjICMHnzMvLK3D3PnsdcPe4eQCNgNFALyALSCR8F/8H4JkDXZ+Tk+PljR49eq99h2rlxu0+4KHR3uPX7/mkb9ZV2+u6V2/OmhQPOeMho3t85IyHjO7B5wQmewU1KeabIMpy943AGGCgu69y91J3DwHDgH6BhiPcRS1/SC5ZDdL48TMTNXeEiOxXzBdgM2tmZo0iz9OBU4A5ZtayzGnnAzODyFdeVoNwEd41d8SX6hkhIvsQ8wUYaAmMNrOvgEnAR+7+H+AvZjYjsj8P+FmQIctqHinCLRumcdWzk3h+/GIN1hCRvcT8JLfu/hVwdAX7rwwgTqU1r59G/pD+nP/Y5/zm7VkkGFpZQ0T2EA93wHGrWf1UzjmyFaDBGiKyNxXgGnby4Vl7DNbY9VxERNWghu0arHFrXmc6NK3HX96fy8dfrwo6lojEABXgKMjJzuTO07vx1s0ncHjL+tz4YgHvz1wZdCwRCZgKcBQ1rJfMC9cdyxFtGnLLy1N5d/qKoCOJSIBUgKOsQVoyz197LDntMrk9fypvTl0WdCQRCYgKcAAyUpN47pq+HNuhCT9/ZTqvTv426EgiEoCY7wdcW9VLSeKZq/oy5IXJ/OK1r1i4ppD6acnkdmyifsIidYQKcIDSUxIZ9uM+XDZsAk+MXYQBqckarCFSV6gJImBpyYkM6NYMCE9yXFSswRoidYUKcAw4vnMz0pLDfxQOLFu/LdhAIhIVaoKIATnZmbx0XS5fLFzL5MXrGTHpWxpnpHDnad0ws6DjiUgNUQGOETnZmeRkZxIKOfe9NZOhoxeyfWeIX599uIqwSC2lAhxjEhKMP57fi9SkBJ75/Bt2lJTy+3N7BR1LRGqACnAMMjN++8MepKck8viYhewoLuWsprG/dp+IHBwV4BhlZtx1ejfSkxP5+0fz+LZFIv9zYojkRH1uKlJb6F9zDDMzbju5C/ee2Z1J35Vy04sFjF+4VqtriNQSugOOA0P+pxNLv1nEi7NX88mc1YBW1xCpDXQHHCdOyU5mYM8WhDw8sXuxVtcQiXsqwHHk+v/pSHJiuEtaCOjZqkGwgUSkSmK+AJtZmplNNLPpZjbLzB6I7G9sZh+Z2fzI11r/u3hOdib5Q/pzwdGtSUww/jRyDqs37wg6logcopgvwEARcJK7HwkcBQw0s1zgHmCUu3cBRkW2a72c7Ez+fslRDL+6H99u2MbF/xrPsg0auiwSj2K+AHtYYWQzOfJw4FxgeGT/cOC8AOIF5vjOTXnh2mNZt3UnFz8xnkVrCg98kYjElJgvwABmlmhm04DVwEfu/iWQ5e4rASJfmweZMQjhJolcikpCXPyvCcxeuTnoSCJyEMw9fkZYmVkj4E3gJ8A4d29U5tgGd9+rHdjMhgBDALKysnLy8/P3OF5YWEhGRkaN5q4O+8u5ojDEQ5N2sDPk3JGTRsdGiVFO9714+HnGQ0aIj5zxkBGCz5mXl1fg7n32OuDucfUAfgvcCcwFWkb2tQTmHujanJwcL2/06NF77YtFB8q5dN1WP+HPo7zHr9/z4V98449+Mt8nL14fnXBlxMPPMx4yusdHznjI6B58TmCyV1CTYr4JwsyaRe58MbN04BRgDvAOMDhy2mDg7WASxoa2jevx6g3H0eiwFH7z9iz++uFcLn9qgkbMicSwmC/AhO9uR5vZV8Akwm3A/wEeBE41s/nAqZHtOq1FwzTOP6o1AO6wU4M1RGJazA9FdvevgKMr2L8OODn6iWJbXvfmPPXZInaUhAg5bNleHHQkEdmHmC/AcnBysjN56fpcxs1fw5h5a3ji00U0rJfCTQM6BR1NRMpRAa6Fdq2ucdOAztzx6nT+/P4cNm0v5u6BWuJIJJaoANdiKUkJ/POSo2iQlsQTYxeyeUcxvzu3F4kJKsIisUAFuJZLTDB+f14vGqYn89iYhWzeXszfLz6KlKR4+PxVpHZTAa4DzIy7BnanYXoyf3pvDoVFJTx+eQ7pKcEN2BARFeA65YYTO9EwPZlfvjmDHz/zJT85qQszlm8it2MTTewuEgAV4Drm0n7tqJ+WzG35Uxj8zETMtLqGSFDUEFgHnXVES84/ug1OeHUNDdgQCYYKcB01qF+73R/EhRxaNUwPOJFI3aMCXEflZGcy4vpcrj6uPZn1krn/3VkULFkfdCyROkUFuA7Lyc7kt+f05J1bT6DxYSlc/tSXjJ67OuhYInWGCrCEZ1K7sT+dm2dw/fDJvDV1edCRROoEFWABoGlGKiOuz6VP+0x++u9pPPv5N0FHEqn1VIBlt/ppyTx3dT9O75nFA+9+zd8+nLtrEnwRqQHqByx7SEtOZOhlx3DfmzN55JMFrNu6k/OPbs3Eb9ZrwIZINVMBlr0kJSbw4I960zgjhcfHLCR/4lJAAzZEqpuaIKRCZsbdA7uT160ZIQ/3FS7WgA2RaqUCLPt160ldSE4MT18ZcujWon7AiURqDxVg2a+c7Ezyh/TnwpzWJCcm8MC7s1i8dmvQsURqBRVgOaCc7Ez+etFRvHJjf7YWlfKjx7/gq2Ubg44lEvdUgKXSjmrbiNdu7E9aciKXPjmBz+avCTqSSFyL+QJsZm3NbLSZzTazWWZ2e2T//Wa23MymRR5nBp21LujYLIM3bj6Odo3rcc1zk3h7mkbNiRyqmC/AQAlwh7sfDuQCt5hZj8ixf7j7UZHHyOAi1i1ZDdL49w39OaZdJrfnT+OpzxYFHUkkLsV8AXb3le4+JfJ8CzAbaB1sKmmYnszwa/pxRq8W/P6/s/nTyNmEQho1J3IwLJ6GmppZe+BToBfwc+AqYDMwmfBd8oYKrhkCDAHIysrKyc/P3+N4YWEhGRkZNRm7WsRqzpA7L369k0++LaF300Sy65VwVMt0OmfG7npzsfqzLC8ecsZDRgg+Z15eXoG799nrgLvHxQPIAAqACyLbWUAi4bv4PwDPHOg1cnJyvLzRo0fvtS8WxXLOUCjkd7823bPv/o9n3/0f73bfSJ+8eH3QsfYpln+WZcVDznjI6B58TmCyV1CTYr4JAsDMkoHXgZfc/Q0Ad1/l7qXuHgKGAf2CzFiXmRltG9fDIts7SkKMmr0q0Ewi8SDmC7CZGfA0MNvd/15mf8syp50PzIx2NvlebscmpCYn7C7CrxV8y6I1hYFmEol1MV+AgeOBK4GTynU5+4uZzTCzr4A84GeBpqzjcrIzeem6XH7UJZk//6g3pSH40eNfMGXpXs3yIhIR87Ohufs42H1jVZa6ncWYnOxMtnRKYUDfdhzboQmDn53IZcMm8MigYzi1R1bQ8URiTjzcAUscat/0MF6/6Ti6ZdXnhhcm89KXS4KOJBJzVIClxjTNSGXEkFxO7NqM+96cqRU2RMpRAZYaVS8liWE/7sOlfdvyyCcL+MVrX1FcGgo6lkhMiPk2YIl/SYkJ/OmC3rRomMY/P57PgtWFnNi1Gf/TtZlW15A6TXfAEhVmxk9P6crNAzoy7duN/N+o+Vw2bAIFS9RLQuquqN4Bm1kf4AdAK2A74b67H7v7+mjmkOAclpqMAQ4UlYT471crdBcsdVZU7oDN7CozmwL8EkgH5gKrgROAj8xsuJm1i0YWCdauARsJkY6F+ZO+ZfxCrTMndVO07oAPA4539+0VHTSzo4AuwNIo5ZGA7BqwMWHROjo1O4y/fTiPHz/zJQ9deCTnHa1J7qRuiUoBdvehBzg+LRo5JDbkZGfubnbo36kpN7wwmZ/+exrLN27n5gGdCI8+F6n9olKAzezh/R1399uikUNiz655he9+7Sse+mAuyzZs43fn9iIpUZ8PS+0Xrb/lBZFHGnAMMD/yOAoojVIGiVGpSYn845KjuDWvMyMmfst1z0+msKgk6FgiNS5aTRDDIfxhHJDn7sWR7SeAD6ORQWKbmXHn6d1onZnOr96aySX/Gs8zV/Ulq0Fa0NFEaky0f89rBdQvs50R2ScCwKB+7XhqcB8Wr93KWQ9/xgPvzlJfYam1ol2AHwSmmtlzZvYcMAX4Y5QzSIzL69ac+8/pydrCnTz7+WIGPakBG1I7RbUAu/uzwLHAm5FH/13NEyJlrd5StLuv8M7SEMM+1crLUvsE8VFzEbAS2AB0NbP/CSCDxLjcjk1ISUog0SDB4P1Z3/HQB3O08rLUKtEeinwdcDvQBpgG5ALjgZOimUNiX9kBG33bZ/Lm1OUMHb2QJeu28deLjiQtOXZXXRaprGjPhnY70BeY4O55ZtYdeCDKGSROlB2w0bd9Y9o3OYw/vTeHlZt28OSVOTTJSA04oUjVRLsJYoe77wAws1R3nwN0i3IGiUNmxg0nduKxy49h5vJNXPD4FyzUop8S56JdgJeZWSPgLcKT8LwNrIhyBoljZ/ZuyYghuRTuKOGCx77gy0WayEfiV7R7QZzv7hvd/X7g14SXmz8vmhkk/h3TLpO3bjmephkpXPH0l7w5dVnQkUQOSdQKsJklmNnMXdvuPtbd33H3nQe4rq2ZjTaz2WY2y8xuj+xvbGYfmdn8yFdNKluHtG1cjzduOp4+2Y352b+nc9dr0xk6er76C0tciVoBdvcQMP0Q5v0tAe5w98MJ95q4xcx6APcAo9y9CzAqsi11SMN64Yl88ro145XJy3jog3lcrlU2JI5Euw24JTDLzEaZ2Tu7Hvu7wN1XuvuUyPMtwGygNXAusGsQx3DUlFEnpSQl0Kd9JrsmsNxREuKTOasCzSRSWRbNZcLN7MSK9rv72Epe3x74FOgFLHX3RmWObXD3vZohzGwIMAQgKysrJz8/f4/jhYWFZGRkVPI7CI5y7tuCDaX8ZdIOikPhpY4apcAv+qXTOqPi+wv9LKtPPGSE4HPm5eUVuHufvQ64e40/iBT6qpxDeOKeAuCCyPbGcsc3HOg9cnJyvLzRo0fvtS8WKef+TV683h/9ZL6PmLjEc373kff6zfs+du7qCs/Vz7L6xENG9+BzApO9gpoUrSaI0Wb2k/Ltv2aWYmYnmdlwYPC+LjazZOB14CV3fyOye5WZtYwcb0l4jTmpo3KyM7klrzOX9m3H27ceT+vMdK5+bhIvTFgSdDSRfYpWAR5IeOL1EWa2wsy+NrNFhCdlHwT8w92fq+hCC69P8zQw293/XubQO3xftAcDb9dUeIkvrRul89pNx3Fi12b8+q2Z3P/OLEo1h4TEoGhNyL4DeAx4LHI32xTY7u4bK3H58cCVwAwz27V23L2Ep7Z8xcyuJbyY50XVn1ziVUZqEsN+3Ic/jpzN0+O+Ycm6rTw86GjqpyUHHU1kt2jPBYGHV8NYeRDnjwP2tUrjydUSSmqlxATj12f3oEPTw/jtO7O48PHxPH3V3p+DiARFKx9KrXdFbjbPXd2XFZu2c9bD4xg+q0h9hSUmqABLnfCDLs34/bm92Ly9mNHflnDpk+NVhCVwUS/AZpZtZqdEnqebWf0DXSNSHZZt3I5FGrOKS52/aoJ3CVhUC7CZXQ+8BvwrsqsN4ZnRRGrcrlU2EoBEM8YvWs/NL01h286SoKNJHRXtO+BbCPdq2Azg7vOB5lHOIHXUrlU2LuiSzCs35PKrsw7nw6+/48LHx7Ni4/ag40kdFO0CXORlZj8zsyTCo0dFoiInO5OzO6WQ074x1/2gI08P7svS9ds459HPmbpUbcISXdEuwGPN7F4g3cxOBV4F3o1yBpHd8ro3582bj6NeSiKXPDmBt6ctDzqS1CHRLsD3AGuAGcANwEjgV1HOILKHLln1eeuW4zm6bSNuz5+m1ZclaqK9IkbI3Ye5+0WEZyj7MjJRhUigGh+WwgvXHsugfm0ZOnohN71UwNYifTgnNSvavSDGmFkDM2tMeFn6Z83s7we6TiQaUpIS+OP5vfnN2T346OtVnPnwZ/zpvdnqLyw1JtpNEA3dfTNwAfCsu+cAp0Q5g8g+mRnXnNCBe888nCXrtvGvsYsY9KRW2ZCaEe0CnBSZOvJi4D9Rfm+RSisqCZEQGbSxszTEsE8XBRtIaqVoF+D/BT4AFrj7JDPrSHhKSpGYsmvQRqJBgsH7s77j/ndmUVIaCjqa1CJRnQ3N3V8l3PVs1/Yi4EfRzCBSGbsGbUxYtI6+7TP5YNYqnh73DQtWF/LoZUfTqF5K0BGlFohqATazNOBaoCeQtmu/u18TzRwilZGTnUlOdniZwX4dmtCtRX1+9eZMzh36OcN+3IeuWZrGRKom2k0QLwAtgNOBsYTngtgS5Qwih+TiPm0ZMSSXrUWlXPDYF3z8tVZflqqJdgHu7O6/Bra6+3DgLKB3lDOIHLKc7Eze/cnxdGh6GNe/MJmhoxegruxyqKJdgIsjXzeaWS+gIdA+yhlEqqRlw3ReuaE/Zx/Rioc+mMtt+dPYvrM06FgSh6K9JNGTZpYJ/JrwopoZwG+inEGkytJTEnn40qM4vGV9HvpgLrOWb+S0ni04tUeL3e3GIgcS7V4QT0WejgU6RvO9RaqbmXHzgM4kmvGn9+bwxNhFPDNuMSOG5KoIS6VEeyhyqpldZmb32v+3d+/hUVX3Gse/v9xAiIAkEgWEEKAogoIJCEWQtFVB9KF4KAURaU8rVsXLo7XWU289vVl7rIKKioDSHiU93pACFSlGLkpQUC6Ru1wUUQIFlKgEQtb5YyZ0CAkEMpm998z7eZ55yOy9Z/KygB8ra/Zay+y+ysdxXjPFzErMrDji2ANm9qmZLQ8/Lq//9CLVK69wR0zaeGr+R94GksCI9Rjwa8BgoBz4KuJxLM8BA6o5/ohzrlv4MTuqKUVOQNVJG3NX7+DuV1ZRVq5xYTm2WI8Bt3bOVVdMa+ScW2Bm2fUTR6TuIidt9GzXnMK1JUx46yPW79jHkyMvoEWThsd/E0lIse4Bv2Nm0brtbKyZrQwPUWjATTyV2/Y0bsrvMOZ8vQAAE4VJREFUQI/s5vxiwNk8fnV3Vm//kisfX6SdNqRGFot7GM1sFaGth1KAjsAmoAwwwDnnzjvO67OBmc65LuHnWcCu8Hv+Bjizptl0ZjaG0NrDZGVl5RYUFBxxvrS0lPT09JP9rcWMckZPrDJ+sq+C8e/vZ89+x+hz0+jbOvWEXq+2jB6vc+bn5y9zzuUddcI5V+8PoO2xHrV4fTZQfKLnqj5yc3NdVYWFhUcd8yPljJ5YZtxdWuaufmaxa3vXTHff9FXuQPmhWr9WbRk9XucElrpqalKshiB2AEOAOwl9oPapc25r5eNE3yy8pGWlIUBxTdeKeOm0xmlM/XFPruvbjqmLtzJy0hJ2lZZ5HUt8IlYfwk0lNAtuITAQ6AzcWpsXmtk0oD+QaWbbgPuB/mbWjdAQxBZC+8uJ+FJKchK/GtSZzi2b8MuXVzHgkQVc3vVMBndvpfuFE1ysCnBn51xXADObDLxb2xc650ZUc3hytIKJxMqQ7q0pL3fc+fJK/lK0lWnvfUzBmN4qwgksVkMQlWtA4JzTToeSsEpKyw5P2jh4yPG7Was5qEXeE1asCvD5ZvZl+LEPOK/yazP7MkYZRDwXOWkjOcl4/+O9XP1MESX79nsdTTwQkyEI51xyLL6PiN9FTtrolZPBtj1fc9fLK7li/CKevCZXwxEJJtYTMUQSXuWkjdy2pzG4WytevbEPDVOTGT5xMf9btFXrCycQFWARj51zZhP+PvYi+nTI5J7pxdz18kr2H9Q6Eokg1mtBiEg1mjZKZfLoHoz753rGv7mRtZ/vY3R7fTgX79QDFvGJ5CTj9ks7MXFULpt2fsUD73zDOx/t8jqW1CMVYBGfufTcM3htbB9OTTOumbSE+18r5onCDSzbqkV94o2GIER8qP3p6dzb+xSmbmrI1MWh2foNUzby/HXabSOeqAcs4lOnpBh9O2YSnrfB/vIKZq3a7mkmiS4VYBEf65WTSYPUpMOz514o+pjZqz7zNpREjQqwiI9VTty4I/zh3Dktm3Dj8+/z25mawhwPNAYs4nO5bU87PO7bv1MLfjdrNZMWbWblti94fGR3WpyqLY+CSj1gkQBJS0ni14O78OgPu7Hq0y8YNH4R723Z7XUsOUkqwCIB9P3urXj1pm/TOC2ZEROLmLxos6YwB5AKsEhAnX1GE2bcfBHfObsFv5m5mrHTPqC0TKu9BonGgEUCrEnDVJ4elctT8zfxpzlrWf7JXi47N4tBXVvqfuEAUA9YJODMjBv6t+e+Kzrz6Z5vmLJoC8MnLtbMuQBQARaJE18dOHTUbhtl5VpVzc9UgEXiRHW7bQx7uohte772OprUQGPAInGi6m4bJV/u5xcvreSKxxbxyA+7kd+phdcRpQrf94DNbIqZlZhZccSx5mY218w2hH/Vpw0iHLnbxsCuZzLj5os4o0lDfvzsezz8xjoOVehWNT/xfQEGngMGVDn2S2Cec64jMC/8XESqaJfZmOk39WFYXmsee3MjoyYvYee+Mq9jSZjvC7BzbgFQdarPYGBq+OupwPdjGkokQBqmJvPQ0PN5aOh5LNu6h0HjF/LuZs2e8wMLwuwZM8sGZjrnuoSf73XONYs4v8c5V+0whJmNAcYAZGVl5RYUFBxxvrS0lPT09HpKHj3KGT1ByAj1k/PjLw/xxPIydn7j+MG30hiQnYKZHf+FNUjktjwR+fn5y5xzeUedcM75/gFkA8URz/dWOb+nNu+Tm5vrqiosLDzqmB8pZ/QEIaNz9Zfzi28OuJ/9dalre9dM94Mn33YPz1nrlm7ZfVLvlehtWVvAUldNTfL9EEQNdpjZmQDhX0s8ziMSGE0apjJh5AX8uE82727Zw/g3NzLimSJN3PBAUAvwDGB0+OvRwGseZhEJHDMjM73B4YkbB8ormPDWRi3oE2O+L8BmNg1YDHQys21m9hPgQeASM9sAXBJ+LiInIHLiRpLBvDUl3PTC+3y5/6DX0RKG7ydiOOdG1HDquzENIhJnIiduXNiuOcu27uGhOeso/nQRE0ZeQJdWTb2OGPd83wMWkfpTOXEjL7s511/cnr+N6cXBQxVcNeEd/rp4i4Yk6pkKsIgclpfdnFm39KVPhwzufe1Dxk77gH0akqg3KsAicoTmjdOYPLoHdw04m9eLP+eKxxZR/OkXXseKSyrAInKUpKTQGsMFY3pRdjA8JFG0VUMSUaYCLCI16pHdnFm3XETv9hncO72YayYt4ZG563TPcJSoAIvIMWWkN+DZH/Vg5IVtePujfzFuniZuRIsKsIgcV1KS0bLZKUdO3CjUxI268v19wCLiD5UTNw6WV+CAeWtL2LM7me4XHqBZozSv4wWSesAiUiuVEzduv7QTL17fm3sGncOKnYcYNH4Ry7ZqecuToR6wiNRabtvTDm93n5vdnKTdm3lunTHs6SJuv+Rb3HBxe5KSTn55y0SjHrCInLScpsnMvOUiBnY5gz/NWcfoZ9/VjhsnQAVYROqkScNUHhvRnT9c1ZV3N+9m4LiFLNqwy+tYgaACLCJ1ZmaM6NmG18b2oVmjVEZNWcL/zFlH+aEKr6P5msaARSRqzj6jCTPG9uGBGR/yeOFG5q3ZwcWdTueSzmccHjuWf1MPWESiqlFaCg8NPZ/bvteRNZ/v46n5mxg+cbEmblRDBVhE6kVqctLhiRsHDzl+M3M1+w8e8jaUz6gAi0i9iNxxIznJWP7JXgY//jbrd+zzOppvaAxYROpF5I4bvXIyKC0r547/W86Vjy3i3is6M/LCNpgl9j3DKsAiUm8iJ24A/OPWftzx4grumV7Mwg07+eN/nJfQ05g1BCEiMXP6qQ147kc9+NXl5/Dm2hIGjlvIkk3/8jqWZ1SARSSmkpKM6/rl8MoNfWiQksSIZ4r489z1CXnPcKALsJltMbNVZrbczJZ6nUdEaq9r66bMvKUvQ7q3Zvy8DQyfWMS2PV97HSum4mEMON85p3mPIgGU3iCFh4edT9+OmdwzvZiB4xYypm8OSUlGr5yMuJ+8EQ8FWEQC7vvdW9G9TTN+MnUpD89dD0DDlCSev65XXBdhC/KK9ma2GdgDOOBp59zEaq4ZA4wByMrKyi0oKDjifGlpKenp6TFIWzfKGT1ByAjByBntjDM2HuCVjQcPP/9OmxSu7dygzu/rdVvm5+cvc87lVT0e9B5wH+fcdjNrAcw1s7XOuQWRF4SL8kSAvLw8179//yPe4K233qLqMT9SzugJQkYIRs5oZzy13R5mby3iQHkFFQ7mf1JOj84dGNMvh+Q6rDPs17YMdAF2zm0P/1piZq8CPYEFx36ViPhV5OSNLi2b8Leln/DH19cyf30Jfx7WjZbNTvE6YlQF9i4IM2tsZqdWfg1cChR7m0pE6iq37WnclN+Bizu14ImrL+ChoeexctsXDBy3kNmrPvM6XlQFtgADWcAiM1sBvAvMcs697nEmEYkiM2NY3lnMuqUv2RmNuPH597nzxRV8VVbudbSoCOwQhHNuE3C+1zlEpP61y2zMSzd8m0f/uZ4Jb33Ee1t28+jw7nQ7q5nX0eokyD1gEUkgqclJ3HnZ2Uy7rhcHyisY+uQ7PFG4kUMVwb2TSwVYRAKlV04G/7i1H5eFNwId8UwRcz78nCcKNwZu0ffADkGISOJq2iiVx0d0J79TC+55dRXX/3UZBjRITeL5nwZn8oZ6wCISSGbG0NzWjOrdFgjNxtp/sIL560q8DXYCVIBFJNAGdDmThqlJVE7TeGHJxyz+KBhLXKoAi0igVU7e+PllnfjDVV059ZRUrp5UxO9nr6Gs3N970GkMWEQCL3LnjcHdWvLbWWuYuGATC9bvZNzw7h6nq5l6wCISVxqlpfD7IV2ZdG0eu0rLuPLxRczZcpAKH96upgIsInHpe52zeP22fvTrmMm0tQcYNWUJn33xjdexjqACLCJxKzO9Ac9cm8ePzk3j/a17ueyRBcxcud3rWIepAItIXDMz+p+Vyuxb+5JzejpjX/iAa6cs4c9z13k+cUMFWEQSQrvMxrz0s94My2vNgvW7GD9vIyMmFnlahFWARSRhpCQn0TajMZVrux84VMEfZq9h/0FvbldTARaRhNIrJ4O0lCSSDZKTjKVb93DFY4tYuW1vzLPoPmARSSiRu270ysngq7JyfvHSSoZMeIebv9OBm/I7kJocm76pCrCIJJzIiRsAc27rx/0zinn0nxsoXFvCw8O60aFF/W/iqSEIEUl4TRul8ujw7jxx9QV8vPtrBo1fyLNvb673yRsqwCIiYYPOO5M5t/WjT4dMfv331VwzeQmf7q2/yRsqwCIiEVo0acjk0Xk8eFVXVnyylwGPLODlZdtwLvq9YY0Bi4hUYWYM79mGb7fP5I4Xl3PHiyt4Y/XnDO9xFqs/20evnIyoLPoe6AJsZgOAcUAyMMk596DHkUQkjrTJaETBmN5MWriJP81Zx5wPd5BkkJYSnZ03AjsEYWbJwBPAQKAzMMLMOnubSkTiTXKScf3F7bmmV2jnjQoHB8srKNpU90XfA1uAgZ7ARufcJufcAaAAGOxxJhGJU1ee35KGqaEJHKkpSfTKyajzewZ5CKIV8EnE823AhR5lEZE4V3UCRzTGgK0+PtmLBTP7AXCZc+6n4eejgJ7OuZurXDcGGAOQlZWVW1BQcMT7lJaWkp5e/zdc15VyRk8QMkIwcgYhI3ifMz8/f5lzLu+oE865QD6A3sCciOd3A3cf6zW5ubmuqsLCwqOO+ZFyRk8QMjoXjJxByOic9zmBpa6amhTkMeD3gI5m1s7M0oDhwAyPM4mI1Fpgx4Cdc+VmNhaYQ+g2tCnOuQ89jiUiUmuBLcAAzrnZwGyvc4iInIwgD0GIiASaCrCIiEdUgEVEPBLY+4BPhpntBLZWOZwJ7PIgzolSzugJQkYIRs4gZATvc7Z1zp1e9WBCFeDqmNlSV90N0j6jnNEThIwQjJxByAj+zakhCBERj6gAi4h4RAUYJnodoJaUM3qCkBGCkTMIGcGnORN+DFhExCvqAYuIeEQFWETEI3FbgM1sipmVmFnxca7rYWaHzGxoxLEBZrbOzDaa2S99nHOLma0ys+VmttTLnGbW38y+CGdZbmb3RZyLSXvWMaNv2jIi63Iz+9DM5kcc90VbHiejb9rSzO6M+PMuDv8bah4+F7N/5zWqbo3KeHgA/YALgOJjXJMMvEloQZ+hEcc+AnKANGAF0NlvOcPHtwCZfmhPoD8ws4bsMWnPk83ow7ZsBqwG2oSft/BhW1ab0W9tWeXaK4E3Y92Wx3rEbQ/YObcA2H2cy24GXgZKIo7FdK+5OuSMqVrmrE7M2rMOGWOqFjmvBl5xzn0cvr7yz91PbVlTxpg6wT/zEcC08Ne+2FMybgvw8ZhZK2AI8FSVU9XtNdcqVrmqOkZOAAe8YWbLwlsvea23ma0ws3+Y2bnhY75qT6rPCP5qy28Bp5nZW+E814aP+6kta8oI/mpLAMysETCAUEcGfNKWgV4PuI4eBe5yzh0ys8jjVs21Xt6rV1NOgD7Oue1m1gKYa2Zrwz0CL7xPaL57qZldDkwHOuKv9qwpI/irLVOAXOC7wCnAYjMrwl9tWW1G59x6/NWWla4E3nbOVfaWfdGWiVyA84CCcFHLBC43s3JC/xOeFXFda2B77OMdVm1O59x059x2CP34Z2avEvqxypO/6M65LyO+nm1mE8wsEx+1Z00ZnXO7/NSWhNpsl3PuK+ArM1sAnI+P2pKaM673WVtWGs6/hx/AJ22ZsEMQzrl2zrls51w28BJwo3NuOj7ba66mnGbW2MxOBTCzxsClwDHvpKhPZnaGhf+XMLOehP5u/QsftWdNGf3WlsBrQF8zSwn/6HwhsAYftWVNGX3YlphZU+BiQpkr+aIt47YHbGbTCH3qnWlm24D7gVQA51x146mEz8V0r7mTzQlkAa+G60kK8IJz7nUPcw4Fbgj/FPENMNyFPm6OWXuebEYz81VbOufWmNnrwEqgApjknCsOv9YXbVlTRjPLwUdtGb5sCPBGuLdO+Jwv9pTUVGQREY8k7BCEiIjXVIBFRDyiAiwi4hEVYBERj6gAi4h4RAVYRMQjKsAiIh5RAZaEYWZdzWyrmd1wjGtOMbP5ZpYche+XZmYLzCxuJzxJ3agAS8Jwzq0iNOX02mNc9p+Ellk8FIXvdwCYB/ywru8l8UkFWBJNCXDuMc6PJLxmgJllW8ROC2b2czN7IHx8rZlNstAuC8+b2ffM7G0z2xBeZ6LS9PB7ihxFBVgSzYNAAzNrW/VEeFGWHOfcllq8TwdgHHAecDahBcovAn4O/FfEdcVAjzpmljilAiwJw8wGAI2BWVTfC84E9tby7TY751Y55yqAD4F54cWHVgHZlReFhzIOVK4QJhJJBVgSgpk1BB4CbiRUJLtUc9k3QMOI5+Uc+W8k8lxZxNcVEc8rOHqVwQbA/hNPLfFOBVgSxT3AX8LDC9UWYOfcHiA5XKwBdgAtzCzDzBoAV5zoNzWzDGCnc+7gSSeXuKUCLHHPzDoBlxDa3glq7gEDvEFoLJdw0fxvYAkwE1h7Et8+n9Bu1iJH0XrAIhHMrDtwu3NuVJTe7xXgbufcumi8n8QX9YBFIjjnPgAKozURA5iu4is1UQ9YRMQj6gGLiHhEBVhExCMqwCIiHlEBFhHxiAqwiIhHVIBFRDzy/zqChPPaX65UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot S21\n",
    "f, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(5, 8))\n",
    "ax1.plot(1/mode1_freqs, 10 * np.log10(S21_mag), '.-')\n",
    "ax1.set_title(\"S21\")\n",
    "ax1.set_xlabel(r\"$\\lambda$ (um)\")\n",
    "ax1.set_ylabel(\"Magnitude (dB)\")\n",
    "ax1.set_ylim(None, 0)\n",
    "ax1.grid()\n",
    "\n",
    "ax2.plot(1/mode1_freqs, S21_phase, '.-')\n",
    "ax2.set_xlabel(r\"$\\lambda$ (um)\")\n",
    "ax2.set_ylabel(\"Phase (rad)\")\n",
    "ax2.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAI4CAYAAABUVDNLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxU1fnH8c8z2SEsgUDYs7DJKpgAAWWJilrcd0WtWhW1/qyttT+1tnVvrVtb6779xI1Yd0URRFlE1kQQ2QIJSdj3sARISDLP748ZaIQQJsnM3CHzvF+veTlz7517v2B4cubMueeIqmKMMSb4XE4HMMaYcGUF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2IQVETlFRGaLyC4R2SEi34vIIBFpLyKficgGEVERSTnsfZd537dPRKY7Et40OlaATdgQkebARODfQCugI/AgUA64ga+Ai4/y9h3AP4HHAp/UhAuxO+FMuBCRDGCqqras5ZhIoAJIVdWiGvbfCFytqqMCldOED2sBm3CyEqgSkfEi8gsRSXA6kAlvVoBN2FDV3cApgAKvAFu9/b5JziYz4coKsAkrqrpcVa9T1U5AX6ADnr5dY4LOCrAJW6q6AngDTyE2JuisAJuwISIniMjvRaST93Vn4Epgrvd1LBDjPTzG+/rgeyO8ryMBl4jEikhUcP8EprGxAmzCyR5gCDBPRPbiKbxLgN979+8HSr3PV3hfH3SN9/ULwHDv81eCkNk0YjYMzRhjHGItYGOMcUjIFmAROUtE8kQkX0TucTqPMcb4W0h2QYhIBJ5B86OBdcAC4EpVXeZoMGOM8aNQbQEPBvJVdbWqHgCygfMdzmSMMX4V6XSAo+gIrK32eh2eb68PEZFxwDiAuLi49M6dOx9xErfbjcsVqr9jfs6yBoZlDQzLWjcrV67cpqptDt8eqgVYatj2s74SVX0ZeBkgIyNDc3JyjnjD9OnTGTVqVCDy+Z1lDQzLGhiWtW5EpLim7aH6K2wdUL1J2wnY4FAWY4wJiFAtwAuA7iKSKiLRwBXAZw5nMsYYvwrJLghVrRSR/wEmAxHA66q61OFYxhjjVyFZgAFU9UvgS6dzGGNMoIRqF4QxxjR6VoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhVoCNMcYhIVeAReQBEVkvIou8jzFOZzLGmECIdDrAUfxDVZ90OoQxxgRSyLWAjTEmXIRqAf4fEVksIq+LSILTYYwxJhBEVYN/UZGpQLsadt0HzAW2AQo8DLRX1V/VcI5xwDiApKSk9Ozs7CNOVlpaSnx8vB+TB45lDQzLGhiWtW6ysrJyVTXjiB2qGrIPIAVYcqzj0tPTtSbTpk2rcXsosqyBYVkDw7LWDZCjNdSukOuCEJH21V5eCCxxKosxxgRSKI6CeFxEBuDpgigCbnY2jjHGBEbIFWBVvcbpDMYYEwwh1wVhjDHhwgqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xAqwMcY4xJECLCKXishSEXGLSMZh++4VkXwRyRORM53IZ4wxwRDp0HWXABcBL1XfKCK9gSuAPkAHYKqI9FDVquBHNMaYwHKkBayqy1U1r4Zd5wPZqlquqoVAPjA4uOmMMSY4RFWdu7jIdOAuVc3xvn4WmKuqb3tfvwZMUtUPanjvOGAcQFJSUnp2dvYR5y8tLSU+Pj5wfwA/sqyBYVkDw7LWTVZWVq6qZhy+PWBdECIyFWhXw677VPXTo72thm01/oZQ1ZeBlwEyMjJ01KhRRxwzffp0atoeiixrYFjWwLCs/hGwAqyqp9fjbeuAztVedwI2+CeRMcaEllAbhvYZcIWIxIhIKtAdmO9wJmOMCQinhqFdKCLrgKHAFyIyGUBVlwL/AZYBXwG32QgIY0xj5cgwNFX9GPj4KPseBR4NbiJjjAm+UOuCMMaYsGEF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHGIF2BhjHHLMW5FFpC1wMp4VKvbjWc0iR1XdAc5mjDGN2lELsIhkAfcArYCFwBYgFrgA6CoiHwBPqeruYAQ1xpjGprYW8BjgJlVdc/gOEYkEzgFGAx8GKJsxxjRqRy3AqvqHWvZVAp8EJJExxoSJWvuARWQkUKKqi0XkMmAEUAA8r6rlwQhojDGNVW19wM8B/YEYEVkJxOOZJH0Y8DpwVVASGmNMI1VbCzhLVXuLSCywHmirqlUi8hKwODjxjDGm8aptHHAZgKqWAcUHlwZSzzr2FUHIZowxjVptLeC2InInnqXiDz7H+7pNwJMZY0wjV1sBfgVoVsNzgFcDlsgYY8JEbcPQHgxmEGOMCTe1jYJ4prY3qupv/B/HGGPCR21fwuV6H7HAScAq72MAUBX4aMYY07jV1gUxHkBErsMzJK3C+/pFYEpQ0hljTCPmy3SUHfj5F3Dx3m31JiKXishSEXGLSEa17Skisl9EFnkfLzbkOsYYE8qOOR0l8BiwUESmeV+PBB5o4HWXABcBL9Wwr0BVBzTw/MYYE/KOWYBV9f9EZBIwxLvpHlXd1JCLqupyABFpyGmMMea4Jp4b22rYIZKiqkVHfaOnenZU1XX1vrjIdOAuVc05eE1gKbAS2A38SVW/O8p7xwHjAJKSktKzs7OPOKa0tJT4+Pj6xgsqyxoYljUwLGvdZGVl5apqxhE7VLXGB/A+nrl+fwn0AdoCXYBTgYeB2cDoWt4/FU9Xw+GP86sdMx3IqPY6BmjtfZ4OrAWaH+0aBx/p6elak2nTptW4PRRZ1sCwrIFhWesGzypCR9Su2kZBXCoivfHMevYroD2wD1gOfAk8qp55Io72/tOP+WvhyPeUA+Xe57kiUgD0AHLqei5jjAl1tfYBq+oy4L4gZUFE2gA71DPrWhrQHVgdrOsbY0wwObIqsohcKCLrgKHAFyIy2btrBLBYRH4EPgBuUdUdTmQ0xphA82UYmt+p6sfAxzVs/xBbY84YEyYcaQEbY4zxoQCLx9Ui8hfv6y4iMjjw0YwxpnHzpQX8PJ6+2iu9r/cAzwUskTHGhAlf+oCHqOpJIrIQQFVLRCQ6wLmMMabR86UFXCEiEYDCoaFi7oCmMsaYMOBLAX4Gz4iFtiLyKDAL+GtAUxljTBjwZTKed0QkFzgNz4KcF6h3Mh1jjDH1V9uSRK2qvdwCTKi+z26QMMaYhqmtBZyLp99X8EzCU+J93hJYA6QGPJ0xxjRiR+0DVtVUVU0DJgPnqmqiqrYGzgE+ClZAY4xprHz5Em6Qqn558IWqTsKzKoYxxpgG8GUc8DYR+RPwNp4uiauB7QFNZYwxYcCXFvCVQBs8Q9E+wTMx+5W1vsMYY8wx+TIMbQdwRxCyBNX8wu3Mzt/O8B5tSE9OcDqOMSYMHbMAe1dDPmLhOFU9NSCJgmBOwTaufGUeArw4s4B3bsy0ImyMCTpf+oDvqvY8FrgYqAxMnOD4Yc1OwPNbpaLSzdzV260AG2OCzpcuiNzDNn0vIjMClCcoMtNaE+kSKt1KZISLzLTWTkcyxoQhX+YDblXtkSgiZwLtgpAtYNKTE3jo/L4A3JbVzVq/xhhH+NIFUf2OuEqgELghkKGC4bKMTjw2aTnrSvY5HcUYE6Z8KcC9Dl9+XkRiApQnaCIjXJzSPZEZK7eiqoiI05GMMWHGl3HAs2vYNsffQZwwskcbNu8uZ+XmUqejGGPCUG2zobUDOgJxIjIQTxcEQHOgSRCyBdyIHm0AmLlyKz3bNXM4jTEm3NTWBXEmcB3QCXi62vY9wB8DmClo2reIo3vbeGas3MpNI9KcjmOMCTNHLcCqOh4YLyIXq+qHQcwUVCN7tOHNucXsP1BFXHSE03GMMWHkqH3AInK192mKiNx5+KMhFxWRJ0RkhYgsFpGPRaRltX33iki+iOR5h7wF1IgebThQ6WZuoc0vZIwJrtq+hGvq/W880KyGR0N8DfRV1f7ASuBeABHpDVwB9AHOAp73LggaMINTWxET6WJG3tZAXsYYY45QWxfES97/Pujvi6rqlGov5wKXeJ+fD2SrajlQKCL5wGACOOoiNiqCzLTWzFxlBdgYE1yiesQ8Oz8/wLMM/U1ACtUKtqr+yi8BRD4H3lPVt0XkWWCuqr7t3fcaMElVP6jhfeOAcQBJSUnp2dnZR5y7tLSU+Pj4Y2aYXFTBhBUHeGJEHG2a+DIyz/98zRoKLGtgWNbACIWsWVlZuaqacfh2X27E+BT4DpgKVPl6QRGZSs23LN+nqp96j7kPz9117xx8Ww3H1/gbQlVfBl4GyMjI0FGjRh1xzPTp06lp++E6bdnDhBUzOdC6K6OGJB/z+EDwNWsosKyBYVkDI5Sz+lKAm6jq3XU9saqeXtt+EbkWz/pyp+l/m+HrgM7VDusEbKjrteuqa5t4OraMY+bKrVzlUAE2xoQfXz5vTxSRMf68qIicBdwNnKeq1Sdj+Ay4QkRiRCQV6A7M9+e1j5KHET0SmZ2/nYoqd6AvZ4wxgG8F+A48RXi/iOwWkT0isruB130Wz0iKr0VkkYi8CKCqS4H/AMuAr4DbVNXnbo+GGNG9DXvKK1nonSvYGGMCzZf5gP1+j66qdqtl36PAo/6+5rEM65ZIhEuYuXIrg1NbBfvyxpgw5Mt8wCfV8OgqIr70Hx83WsRFMbBzSxuOZowJGl+6IJ7HM1b3Fe9jLpANrBSRMwKYLehG9GjDT+t3sb203Okoxpgw4EsBLgIGqmq6qqYDA4AlwOnA4wHMFnQjerRBFWblb3M6ijEmDPhSgE/wfjkGgKouw1OQVwculjP6dWxBQpMoZqy0bghjTOD50o+bJyIv4Ol2ALgcT/dDDFARsGQOiHAJp3Rvw8yV23C7FZfLVskwxgSOLy3g64B84LfA74DV3m0VQFaggjllRPdEtpWWs3xTQ0faGWNM7XwZhrYfeMr7OFyjW8tn5KFVMrbRp0MLh9MYYxozX4ahdReRD0RkmYisPvgIRjgntG0eywntmjFj5RanoxhjGjlfuiD+D3gBz6Q5WcCbwFuBDOW0kT3akFtcwt7ySqejGGMaMV8KcJyqfoNn6spiVX0AODWwsZw1skcbKqqUOQW2SoYxJnB8KcBlIuICVonI/4jIhUDbAOdyVHpKAnFRETYczRgTUL4U4N/iWYb+N0A6cA1wbSBDOS0mMoKhXW2VDGNMYPkyCmKB92kpcH1g44SOkT3a8O2KLRRt20tKYtNjv8EYY+roqAVYRD6r7Y2qep7/44SOEd7haI98sYxbR3UjPTnB4UTGmMamthbwUGAtMAGYR83LBTVaO0rLEWDq8i1My9vKNZldGN69DamJTencqglREZ7em9ziEuau3k5mWmsr0saYOqmtALcDRgNXAmOBL4AJ1eeFaMzmFu449LzKrbwxu5g3ZhcDnluWOyfEkdA0msXrduF2KzFRLt65MdOKsDHGZ7UtS1+FZ1WKr7zzPlwJTBeRh1T138EK6JTMtNbERLmoqHQTFenixavTaR4XReHWvRRu20vh9r0sKNxBlduznN2BSjdzV2+3AmyM8VmtX8J5C+/ZeIpvCvAM8FHgYzkvPTmBd27MPKJ74aQu/y2wucUlXPXqXMoqPOvIZaa1diSrMeb4VNuXcOOBvsAk4EFVXRK0VCEiPTmh1hbtwSL95JQ85hRsp22zmCCmM8Yc72obB3wN0APPopyzvQty+mtRzkYjPTmBJy89EZfAewvWOh3HGHMcOWoBVlWXqjbzPppXezRT1ebBDBnqOraMY1TPtryXs9aWtTfG+MyXO+GMD8YO7sLWPeV8s3yz01GMMccJK8B+MqpnG9o1j+Xd+dYNYYzxjRVgP4mMcHH5oM58t2ora3fsq/d5cotLeG5aPrnFJX5MZ4wJRY4UYBF5QkRWiMhiEflYRFp6t6eIyH4RWeR9vOhEvvq6fFBnBJgwf0293j+/cDuXvzSHp6bkcdWrc60IG9PIOdUC/hroq6r9gZXAvdX2FajqAO/jFmfi1U+HlnFk9WzLf3LW1evLuOemFVDpVtwKFd4bO4wxjZcjBVhVp6jqweUm5gKdnMgRCGOHdGFbaTlTl9Xty7hNe93MLtjGwYWY3Qr9O9madMY0ZqKqzgYQ+Rx4T1XfFpEUYCmeVvFu4E+q+t1R3jcOGAeQlJSUnp2dfcQxpaWlxMfHByh5zdyq3DVjP+2bCn8YFOfze/46Zy/r9wnj+kXzwxY336+v5ITWLu5MjyXSFVrzIDnx91pfljUwLGvdZGVl5apqxhE7VDUgD2AqsKSGx/nVjrkP+Jj//iKIAVp7n6fjmY2t+bGulZ6erjWZNm1ajdsD7R9f52ny3RO1eNten45/e26RJt89UbPnFx/a9kHOWk2+e6LeMeEHdbvdgYpaL079vdaHZQ0My1o3QI7WULsC1gWhqqerat8aHp8CiMi1wDnAVd6AqGq5qm73Ps8FCvDcjXdcuXxQZ1wCExYc+8u4jbv289iXK+jVysVlGZ0Pbb84vRN3ndGDTxZt4InJeYGMa4xxiFOjIM4C7gbOU9V91ba3EZEI7/M0oDuw2omMDdG+RRynnpDE+zlrOVB59C/jVJU/f7KECreb6/vGIPLzrobbsroxdkgXnp9ewFtziwOa2Ya/GRN8x1ySKECexdPd8LW36MxVz4iHEcBDIlIJVAG3qOqOo58mdI0d0pmpyzczdflmxvRrX+MxExdvZOryLdw3phdt3Ue2lkWEh87rw+ZdZdz/6RKSmsVwRp92fs86v3A7V74yD1UlOtLmNTYmWJwaBdFNVTvrYcPNVPVDVe2jqieq6kmq+rkT+fxhZI+2dGwZx7vzau6GKNl7gAc+W8qJnVpw/ckpRz1PZISLf48dSL+OLfhN9kJ+WOP/Fuq789ZQZcPfjAk6uxMuQCJcwuWDOjMrfxtF2/Yesf/hicvYtb+Cv1/Sn8iI2v83NImO5LXrBpHUPJYbx+dQWMP5GqJo+3/v3BMRm9fYmCCxAhxAl2V0JsIlZB82TeX0vC18tHA9vx7VlRPa+TaxXGJ8DOOvHwzAdf83n29XbPZLn+26kn38uG4nl6Z3pEdSPAi0ahrdoHMaY98p+MYKcAC1axHLqSe05YPc/34ZV1peyX0fL6Fb23huO7Vbnc6XktiU167NYOPO/dwwPscvtyz/J2cdAHec3oO3bxhCXFQEf/zop4PDBI2pswnz13Dpi7N5YnIeV77ceG6pzy3a4fdfKlaAA8xzZ9wBpizbBMATX61gw679/P3i/sRERtT5fAO7JDCmX3tUaXCfbWWVm/dz1jK8exs6JTShbfNY/jimF3NWb+d9b2E2xhdVbmXSTxu5+IXZ3PvRT3iXSuRAlZu/fLKE7aXlzgZsoNziEi5/eS5PTM7jqlf890vFCnCAjejeho4t45gwfw05RTt4c24x1w5NadAog2uGphy6ZTkq0lXvPtsZK7eycVcZYwf/d/zx5RmdGZzaike+WMaWPWX1zmjCQ2l5Ja/PKmTUk9O49Z0f2LqnnBtOTiE20kWEQKRLWLF5N6c9PYP3FqzB7T4+P1l9u2Izld7sFVX++6LaqWFoYSPCJVwxqDNPfb2Sn9YtoHXTaP5wZs8GnTM9OYF7x/Ti0S+Wc8MpqfUu5hPmryUxPobTeiUd2uZyCX+7qB+/+Nd3PPj5Mp4be1KDsprGaft+N3/9cjkT5q1hT3klg1ISuG9Mb0b3TiLCJYzp3+HQgrbNYyO57+Ml3P3hT3yYu55HL+xL96RmTv8R6mTlpj0AuKRhjZ7DWQEOgl7tPV+07S6rJDrCxYpNexo8zvbGU1L5MHcdXy3ZxO9H98RVx/kiNu8uY1reFm4ankbUYaMwuraJ5/asbjz19UouHLCZ03snHeUsJlzkFpcwd/V2EuOjmZW/nS8W70ekkDH92nPDKakM6NzyZ8cfvqBt9rhMPshdx6NfLmfMM99x84iu/M+p3YiNqns3XLAtWb+LqSu2cE6/9vTq0Pxnq6Q3lBXgIMjbvAcBFKhyez6+NPR/oIhw66iu3JG9iKnLN9f5Bo33c9ZS5VauGNS5xv03j+zKxMUb+fOnSxiS1opmsVENymsC52Bx9Gdh+Nn5i3ZwxStzqajyfARvEhXBGclR/PmK4XRs6duEUy6XcNmgzpzaqy1//WI5z07L5/PFG3jkgr4M797G75n9RVV58POlJDSJ5tGL+tEizr//DqwPOAgy01oTE+XpE/Pnx5ez+7WnU0IcL84oqNOoBbdbyV6wlmFdW5OS2LTGY6IjXTx2cT827S7jSZuLIiRt2VPGSzMKuPylOTw5OTCT+KsqT0zJO1R8BbhpRBpXnBDtc/GtLjE+hqcvH8C7Nw7BJcI1r83njuyFbN0Tml/Sfb54IwuKSrjrjJ5+L75gLeCgSE9O4J0bM/3eSomMcHHT8DTu/2wpC4pKGJzayqf3zcrfxrqS/fzvWSfUetzALglcOzSF8XOKOG9AR7s92UE79x3gp/W7WLxuFz+u3clP63excdfPvyQtr/DPp6uDVJUnp+Qxd/UOIsTzGS4q0sWIHm3YU7ihQece1i2RSXcM5/npBbw4vYBpK7YwdnAX4mMjGdo1MSR+1vYfqOJvXy6nT4fmXH6UT4oNZQU4SA7vE/OXyzI6869vVvHC9HwGpw726T3ZC9aQ0CSKM/scu2/3rjN7MmXpJu79aDETbx9OdKR9aAq0OQXb+HTRBmIiXezYV8HidTsprna3YmpiUwantqJ/p5bERrp4aOIyyivdKLDVTyNXVJUnJufx/PQCrhzcmYtP6sS8wh2HGhDTCxt+jdioCO4c3YPzTuzAbyb8wIszPfNuxUTm8+5Nzs9H8sKMAjbuKuNfVwwkIkBzclsBPs7FRUdw3bAUnv56JSs27T7mnXVb95QzZelmrh2W4tM45PiYSB6+oC83jM/hxRkF/Oa07v6KbmqQU7SDsa/O42CPUmLTaDJSWnH5oM6c2KklfTu2OOKj8Antm/N9/jbmFGzjjdnFtGoaw+2ndjtidj1fqSqPT87jhekFjB3ShUfO74vLJWSk+PYJq666tY3n7P7tWb5xDwqUV7qZvHSTowV4Xck+XppRwDn92/v8ybI+rDnTCPxyaDJNoiN4acaxZ+788Id1VLqVKwf7/pHqtF5JnNO/Pc9+m0/+ltKGRDXHMH520aHiGyFw/SmpvHhNOr8e1Y2TuyXW2A+ZnpzAb07rzls3DOGikzry9Ncr+ftXefW6m1FV+ftXnuJ7VbXiG2iZaYnERLkOjW9/b8FaFq/bGfDrHs1fv1yOCPxxTK+AXscKcCPQskk0Vwzqwmc/bmBdyb6jHqeqvLdgLYNSEujWtm7jMO8/tw9x0RHc/u4PfF5woNHcXhpKKt1KTnEJAvX6wjYywsWTl5zIVUO68OKMAh74bGmdbnxQVR77agUvzijg6swuPByk4gv//Z7k92f05JkrBhAfE8nlL83lm+V1W1vRH+YUbOfLnzZx68hudKjHF411YQW4kbhxeCoCvPrd0Tvn5q7eQeG2vVwxqEudz9+mWQxXD+nC8k17+HBVhV+/cZ+Vv42HJi4jt+i4nPrZb2atr2TjrjL+96ye3HlGz3rNy+xyCY9c0Jebhqcyfk4xd3+4mCofirCq8tikFbw0Y3XQi+9B6ckJ3JbVjfMGdOTj24bRtW1Tbnozh3fmBXYxguoqq9w8+PlSOraM4+aRaQG/nvUBNxIdWsZx3oAOZC9Yw29O617jjGbZC9bQPDaSs/vXPEH8sTSJ+W+f8YFK/3zjPnf1Nq55zdPn+fqsQk7plsi5J7ZnWNdEOrdq0qBzH0/KKqr4NL+CgV1acsvIrvXuvwXPGPE/julFk+hI/vXNKsoq3Tx92YlH3HBzkKryt0kreHnmaq7JTOah8/s06Pr+0LZZLO+NG8r/vPsD9328hHUl+/nDGXW/4aiuJixYy4pNe3hu7ElBuUnECnAjcsvIrnz0w3rGzy7id6N/vpReyd4DTPppE1cO7lzvH6zMtERiIvM937grR9z9VB9//XIF1bsqF64pYVb+NgC6tGrCyd1aM6xrIsO6tqZ1fEzAbzpwyltziikpV547s6dfip+I8LvRPYiLjuCxSSvYf6CKZ8cOPOL/vary1y+X88p3hfxyaDIPnud88T2oaUwkr/wygz9/upQXphewYed+Hr+kfpNY+WLnvgM8PSWPIamtGNPP/yvP1MQKcCPSI6kZp/dqy/g5Rdw8Mo0m0f/93/vRwvUcqHJzxeC6dz8clJ6cwLs3ZfKvz+bz3fpK/pPjuZmjvv9g31uwhsXrdhHpElQ9Y0zf/NVgmsdFMSt/G9/nb2fijxuZMN8zn3Jyqyas27kft1uJiXTxTggMVfKHPWUVPD89nz6tXQzrmujXc98ysitNoiP4y6dLuenNHF6+JoO4aE8BU1Ue/WI5r84q5NqhyTwQQsX3oMgIF3+9sC+dEuJ4YnIem3aV8fI1GbRo4v+bIv45dRW79lcE9e/B+oAbmVtGdmXnvgreqzYJvKqSPX8NJ3ZueWheivpKT07gV/1iuHN0Dz5dtIF3jrLk0rH8uHYnf/50Kad0S2TCuMz/9nmmtKJ7UjOuPzmVV6/NYOFfRvPxr4fxhzN7ouqZ9lCBsko3t7/7A09OzmN2/jbKKqoa9Ody0ivfFVKyr4JLegRmIvxfDk3hiUv6833+Nq59fT57yipQVR7xFt/rhqWEZPE9SES4Lasb/7x8AD+sKeGSF2fX+mVzfeRt2sNbc4sZO6RLg/+N1IW1gBuZjJRWZCQn8Op3hVydmUxUhIsf1pSwakspf8QhwycAACAASURBVL+4n9+uc1tWNxYUl/DQ58sY0NkzPtVX20vLufXtXNrEx/DMlQNp1TSaQUcZYxoZ4WJglwQGdkkgM601V70ylwNVblwixMdG8sKMAp6dlk90pIuM5ASGdW3N0K6J9O/UgsXrdjGx4ADNUktCtqW8vbSc175bzZh+7UhtsSdg17k0w9P19Lv3FnHhc9/TLC6KhWt2ct2wFO4/t3fIFt/qLhjYkbbNY7j5rVwufH42/3fdoDr93B2NqvLQxKU0jY7gztENm6mwrqwAN0K3jurKDeNz+PzHDVx0UifenbeWptERnNO/g9+u4XIJ/7x8AGc/8x23vpPLxNuH+3SvfGWVm9snLGTb3gN8dOuwOi1/lJ6cwDs3/fyW7j1lFSwo2sHs/O18X7CdJ6esBFYSG+niQJWnr3pi0dyQXen5uWkF7K+o4s7RPVm3LCeg1zr3xA6s37mfxyatADxTpZ7bv/1xUXwPGtY1kQ9vHcb1/7eAy16aw+9G9+BApbtB3wlMWbaZ7/O388C5vYO+HJd1QTRCWT3b0iMpnpdmrGbXvgq++GkD5w3oSNMY//6+bdU0mmfHDmTjzjL+8P6PPg38f2JKHrMLtvPoBX3r1Xo5OFTp4D+2ZrFRnHpCEn86pzeT7hjOD38ezfNXncQJ7ZvhVs8MdAdCdKXn9Tv38/bcYi5J70S3tvFBuWaVWw/d7IAqcwuPv6F/PZKa8dGvh5HUPIZHv1jumYionqtUlFVU8cgXy+iRFM/VmckBSFs7K8CNkMsl3DyiK3mb9/D79xdRVuGu051vdZGe3Ip7fnECU5Zt5rVZtU8Q8MXijYfGmV6aEZg8rZpGM6Zfe/58Th9ivPNWqEJGCLZ+/zV1JeBZjy9YMtNaEx3p/5n5gi2peSznndgR4NB3Ag9PXMqPa3fW6Q7A12YVsnbHfv5yTp9jrk4eCFaAG6nzBnSgQ4tYpi7fQttmMVR4FwUNhBtOSeXMPkk8NmkFucU1t6hWbd7DHz74kYFdWvKXc/oELMtBB0dsDG0fgQLTV24N+DXrIn9LKR/kruPqzOR6TetYXwfvOKvvjR6hZESPNsR6b1+OEGHpht2c/9z3nPnPmbz63Wq2HWMduk27ynhuWj5n9E7ilO7+HX3iK0cKsIg8LCKLRWSRiEwRkQ7V9t0rIvkikiciZzqRrzGIinBxZl/PWMate8q56rV5Abt9WER4/JIT6dAyjtveWXjEAoy7yyq4+a1cmkRH8MJV6UGbUS09OYGbT4zlikGdeXFGAfNCqBviqSl5xEVFcFtW16Bf+/BunONV9duX/3PLUHL/PJq/XtiPpjGRPPLFcjL/+g03vZnDD5srqag6sgHy969WUOlW/nR2bwfSezjVAn5CVfur6gBgIvAXABHpDVwB9AHOAp4XkdBfsyREJTSJOrQSR0NWT/ZFi7gonr/qJHbsO8Dv/vPjoTkI3G7l9//5keId+3hu7Em0axEbsAxH8+dzetOlVRPu/M+P7C6rCPr1D7d43U4mLdnEDcPTaB0f43Sc41r1XybNY6MYO6QLH//6ZKbeOYIbTkll0dqdPLOwnKF/+4ZHJi4jz7u227vz1vDxwvWc2789XVo7d8elIwVYVXdXe9kUT40AOB/IVtVyVS0E8gHfJrk1Rzi5W5uArMRxNH07tuD+c3szc+VWnp2WD3jmVP162WbuG9OLIQ71NzaNieQflw9g0+4y7v90qSMZqntich4JTaK4aXiq01EarW5tm3HvmF7MuedUfntSDOnJCbwxu4gz/zmT05+azn0f/wR4vpdwcmIpx4ahicijwC+BXUCWd3NHYG61w9Z5t5l6CNRKHLUZO7gLCwp38PTXK5lTsI05q3dw/oAOXH9ySsCvXZuTuiRw+6nd+OfUVWSd0JbzTvTfkLy6mF2wje9WbeO+Mb1snb0giIxwMaBtJL8dlcH20nI+WbSBF6fnH2rxHVxi3qnuGKnPnKE+nVhkKlDTDdX3qeqn1Y67F4hV1ftF5Dlgjqq+7d33GvClqn5Yw/nHAeMAkpKS0rOzs4+4UGlpKfHxwRne01CNKeuybZU8nuPpBxbgDxkx9E505nd99axVbuXReWVs2uvm4ZPjaB0X3A+AqsrDc8soKVP+PiKO6Iifj79tTD8DoeTwrKtKKnl8QTlVboh0wf8OiqVbQmB7OrOysnJVNeOIHarq6ANIBpZ4n98L3Ftt32Rg6LHOkZ6erjWZNm1ajdtDUWPK+uy3qzT1nomafPdETb1noj777argBKvB4VkLt5Zqrz9P0itemqNVVe6gZpm8ZKMm3z1RJ8wrrnF/Y/oZCCU1Zc0p2qHPfrtKc4p2BCUDkKM11C6nRkFUX9fmPGCF9/lnwBUiEiMiqUB3YH6w85mGqT7WNDrExpqmJDbl/nN7M2f1dl7/3g8Lm/moyu1Z4DItsSmXpHcK2nVNzUJlJIhTfcCPiUhPwA0UA7cAqOpSEfkPsAyoBG5T1eN3lpUw5UTfc11cltGZb5Zv4fGv8ji5W2JQJl/5dNF6Vm4u5dmxAx0Z8G9Ck1OjIC5W1b7qGYp2rqqur7bvUVXtqqo9VXWSE/lMw4VKC6MmIsJjF/enRZMofpu9KOAzqc1bvZ0HP19KamITxvSt32T4pnGyX8UmLLVqGs0Tl/Qnb/MenpicF7Dr5BaXcNWr89i1v5L1JftZuNa5hSZN6LECbMLWqJ5t+eXQZF6bVcisVdsCco05Bduo9N6UUuXWkJwUyDjHCrAJa/f+ohdd2zTl9gk/8NSUPL8Pyj+4fI7rOJ/8xgSGFWAT1uKiI7hlZFdK9lXw72/z/braM3hWfE5oEsVvT+9x3E9+Y/zPCrAJe1v2lHPwlgh/zh1csLWUGSu3ct2wVH5zWncrvuYIVoBN2MtMa31o7uCDr/3hrTnFREUIY4fUfyFU07hZATZh7+BSR0PTWuFWaBHX8OHxe8oq+CB3Hef070CbZjbjmamZFWBj8BTh565KJy4qghemr27w+T7MXUdpeSXXDUtpeDjTaFkBNsarVdNorhzchU8WrWftjvove+52K2/OKWZA55ac2LmlHxOaxsYKsDHV3DQiFZfAK9/VvxU8c9VWVm/b6/gUnCb0WQE2ppr2LeK4aGAn3luwlq17al9T7GjGzy6iTbMYfmG3HZtjsAJszGFuGdWViip3vWZLK9y2l2l5W7lqSJegrX1njl/2E2LMYVITmzKmX3vemlPMrv11W0PuzTlFNvTM+MwKsDE1+PWobpSWV/LWnCKf31NaXsn7Oes4u1972jYL/uKj5vhjBdiYGvTu0JxTT2jL698Xse9ApU/v+egHz9Cza23omfGRFWBjjuLXo7qyY+8BsuevPeaxbrfyxuwiTuzUgoFd7JZj4xsrwMYcRUZKKwantuKV71ZzoNJd67Gz8rexeuterrOhZ6YOrAAbU4vbsrqxcVcZnyxcX+txb8wuIjE+mjH9bOiZ8Z0VYGNqMaJ7In07NueFGQVUeSdWP1zRtr1My9vC2CHJh+b/NcYXVoCNqYWIcNuobhRu28ukJRtrPObNOcVEiHCVDT0zdWQF2JhjOLNPO9LaNOW5aQWo/rwVvLe8kvdz1jKmX3uSmtvQM1M3VoCNOQaXS7h1ZFeWb9zN9JVbf7bvox/WsceGnpl6sgJsjA8uGNiRji3jeH5a/qFtqsr4OcX079SCk7rYrGem7qwAG+ODqAgX40aksaCohPmFOwD4Pn87+VtKuXZoCiJyjDMYcyQrwMb46PJBnUmMj+Y5byv4jdmFJMZHc86JNvTM1I8jBVhEHhaRxSKySESmiEgH7/YUEdnv3b5IRF50Ip8xNYmNiuD6k1OZsXIrX/60kW9WbOHKwV1s6JmpN6dawE+oan9VHQBMBP5SbV+Bqg7wPm5xKJ8xNbpmaDLNYiL5zYSFCNC/o/X9mvpzpACr6u5qL5sCNY9wNybENI+N4ow+SVS6FbfC7dk/kFtc4nQsc5ySw8c1Bu3CIo8CvwR2AVmqulVEUoClwEpgN/AnVf3uKO8fB4wDSEpKSs/Ozj7imNLSUuLj4wOS398sa2AEIuuHKw/w+WrPPMEu4KLuUZzTNbrB5w33v9dACYWsWVlZuaqaccQOVQ3IA5gKLKnhcf5hx90LPOh9HgO09j5PB9YCzY91rfT0dK3JtGnTatweiixrYAQia07RDu35py817Z6J2vNPX2pO0Q6/nDfc/14DJRSyAjlaQ+2KDFTFV9XTfTz0XeAL4H5VLQfKve/PFZECoAeQE5iUxtRdenIC79yYydzV28lMa016sk0/aeonYAW4NiLSXVVXeV+eB6zwbm8D7FDVKhFJA7oD9V+e1pgASU9OsMJrGsyRAgw8JiI9ATdQDBwc7TACeEhEKoEq4BZV3eFQRmOMCShHCrCqXnyU7R8CHwY5jjHGOMKxURD+JCJb8bSkD5cIbAtynPqyrIFhWQPDstZNsqq2OXxjoyjARyMiOVrT0I8QZFkDw7IGhmX1D5sLwhhjHGIF2BhjHNLYC/DLTgeoA8saGJY1MCyrHzTqPmBjjAlljb0FbMKUiPQTkWIRubWWY+JEZIaINHg+SRGJFpGZIuLU2HpzHLICbBolVf0JuALPhE9H8yvgI1Wt8sP1DgDfAJc39FwmfFgBNo3ZFqBPLfuvAj6FQ4sBLDm4Q0TuEpEHvNtXiMirIrJERN4RkdNF5HsRWSUig6ud7xPvOY3xiRVg05g9BsSISPLhO0QkGkhT1SIfztMN+BfQHzgBGAucAtwF/LHacUuAQQ3MbMKIFWDTKInIWXgm+/+CmlvBicBOH09XqKo/qaobz3zV33inGPwJSDl4kLcr44CINGtIdhM+rACbRkdEYoHHgV/jKZJ9azhsPxBb7XUlP//3UH1febXn7mqv3Rw5n0oMUFb31CYcWQE2jdGfgDe93Qs1FmBVLQEivMUaYDPQVkRai0gMcE5dLyoirYGtqlpR7+QmrFgBNo2Kd5rT0cA/vZuO1gIGmIKnLxdv0XwImIdnodgV9bh8FvBlPd5nwpTdiGHClogMBO5U1Wv8dL6PgHtVNc8f5zONn7WATdhS1YXANH/diAF8YsXX1IW1gI0xxiHWAjbGGIdYATbGGIdYATbGGIdYATbGGIdYATbGGIdYATbGGIdYATbGGIdYATbGGIdYATbGGIdYATbGGIc0igUEzzrrLN22bZvTMYwxpka5ubmTVfWsw7c3igIMkJOTc8S26dOnM2rUqOCHqQfLGhiWNTAsa92ISI3bG0UXhLV+jTEhLrGmjY2iABtjzPHICrAxxjjECrAxxjjECrAxxjjECrAxxjjECrAxxjgkbAtwbnEJz03LJ7e4xOkoxpgw1WhuxKiL+YXbuerVeVS5lehIF+/cmEl6coLTsYwxYSYsW8CzC7ZTUaW4FSoq3cxdvd3pSMaYMBSWBXh49zZEujy3BkZGuMhMa+1wImNMOArLApyenMBr12YQFSEMSWtl3Q/GGEeEZQEGGNmzLdcOTeH7/O2s3bHP6TjGmDAUtgUY4IbhqbgEXvlutdNRjDFhKKwLcPsWcVw4sCPvLVjLttJyp+MYY8JMWBdggJtHduVAlZs3vi9yOooxJsyEfQHu2iaes/q0Y/ycIvaUVTgdxxgTRsK+AAPcMrIre8oqeXfeGqejGGPCiBVg4MTOLTm5W2tenVVIWUWV03GMMWHCCrDXr0d1Y+uecj5euN7pKMaYMGEF2GtY19b079SCl2YUUOVWp+MYY8KAFWAvEeHWkV0p2r6PSUs2Oh3HGBMGrABXc2afdqS1acoL0wtQtVawMSawrABX43IJt4zoytINu/lulS11b4wJLCvAhzl/YAfaNY/l+en5TkcxxjRyVoAPExMZwY3DU5m7egc/rLHVMowxgWMFuAZXDu5Ci7goXpxe4HQUY0wj5lgBFpHOIjJNRJaLyFIRucO7/QERWS8ii7yPMcHO1jQmkmuHpTBl2WYe+nyprRtnjAkIJ1vAlcDvVbUXkAncJiK9vfv+oaoDvI8vnQg3sEtLAF7/voirXp1rRdgY43eOFWBV3aiqP3if7wGWAx2dynO4ZRt2H3pu68YZYwJBQmG8q4ikADOBvsCdwHXAbiAHTyv5iOaniIwDxgEkJSWlZ2dnH3He0tJS4uPj65Upv6SKvy8oo8INLoE/Do6lW0JEvc7li4ZkDTbLGhiWNTBCIWtWVlauqmYcsUNVHX0A8UAucJH3dRIQgad1/ijw+rHOkZ6erjWZNm1ajdt9lVO0Q09/arr2/vMk3Vte0aBzHUtDswaTZQ0MyxoYoZAVyNEaapejoyBEJAr4EHhHVT8CUNXNqlqlqm7gFWCwU/nSkxP420X92Hugig9z1zkVwxjTSDk5CkKA14Dlqvp0te3tqx12IbAk2NmqS09OYEDnlrw2qxC3TdJjjPEjJ1vAJwPXAKceNuTscRH5SUQWA1nA7xzMiIhw4/BUirbvY+ryzU5GMcY0MpFOXVhVZwFSwy5Hhp3V5qw+7ejYMo5XZxVyRp92TscxxjQSdiecDyIjXFx/cgrzC3eweN1Op+MYYxoJK8A+unxQZ+JjInltVqHTUYwxjYQVYB81i43iikGdmbh4Ixt27nc6jjGmEbACXAfXnZyCqjJ+dpHTUYwxjYAV4DrolNCEX/Rrz7vz11BaXul0HGPMcc4KcB3dNDyNPWWVvJ+z1ukoxpjjnBXgOhrQuSUZyQm8/n2hrZ5sjGkQK8D1cOPwVNbu2M+UpZucjmKMOY5ZAa6H0b3b0aVVE161IWnGmAawAlwPES7hVyenkFtcYuvGGWPqzQpwPV2a0ZlmsZG89p21go0x9WMFuJ6axkQydkgXJi3ZyNod+5yOY4w5DlkBboDrhqUgwO/eW2Rrxhlj6swKcANs2FmGAjnFJYx9xRbuNMbUjRXgBqi+UOcBW7jTGFNHVoAbIDOtNdGRnr9CBXq3b+ZsIGPMccUKcAOkJyfwzo2Z/OoUT1/w7AJrARtjfOfYihiNRXpyAunJCezaV8mbc4q54ZQ02rWIdTqWMeY4YC1gP/nt6d1xq/Lvb1c5HcUYc5ywAuwnnVs14fJBnXlvwVrWbLdxwcaYY7MC7Ee3n9qdCJfwz29WOh3FGHMcsALsR0nNY7l2WAofL1zPqs17nI5jjAlxVoD97JaRXWkaHcnTX1sr2BhTOyvAftaqaTQ3nJLKpCWb+GndLqfjGGNCmBXgALhxeCotm0Tx5JQ8p6MYY0KYYwVYRDqLyDQRWS4iS0XkDu/2ViLytYis8v43wamM9dUsNopbRnZlxsqtLCja4XQcY0yIcrIFXAn8XlV7AZnAbSLSG7gH+EZVuwPfeF8fd64dmkKbZjE88VUeqrZ2nDHmSI4VYFXdqKo/eJ/vAZYDHYHzgfHew8YDFziTsGHioiO4/dRuzC/awcxV25yOY4wJQRIKrTMRSQFmAn2BNarastq+ElU9ohtCRMYB4wCSkpLSs7OzjzhvaWkp8fHxAUp9bJVu5e6Z+2keLfxlaCwictRjnc5aF5Y1MCxrYIRC1qysrFxVzThih6o6+gDigVzgIu/rnYftLznWOdLT07Um06ZNq3F7ML23YI0m3z1RJ/20sdbjQiGrryxrYFjWwAiFrECO1lC7HB0FISJRwIfAO6r6kXfzZhFp793fHtjiVD5/uGhgR9LaNOWRL5bx7LerbNJ2Y8whTo6CEOA1YLmqPl1t12fAtd7n1wKfBjubP0VGuLhgQEfWleznqSkruepVWznDGOPhZAv4ZOAa4FQRWeR9jAEeA0aLyCpgtPf1cc3l7fpVoMJWzjDGeDk2H7CqzgKO9q3UacHMEmhDuyYSFbGKiirF5RIy01o7HckYEwLsTrggSE9OYMJNmbRuGk275rEM7Nzy2G8yxjR6VoCDJCOlFX85tzdrS/YzZdlmp+MYY0KAFeAgOrtfe1JaN+HZaavs7jhjjBXgYIqMcPHrUd1Ysn4301dudTqOMcZhVoCD7IKBHenYMo5/f2OtYGPCnRXgIIuOdHHLyDR+WLOTOTYczZiwZgXYAZdmdKZtsxie/Tbf6SjGGAdZAXZAbFQE40akMbtgO7nFNl+wMeHKCrBDxg7pQqum0dYKNiaMWQF2SJPoSG44JZVpeVtZst7WjjMmHFkBdtAvhybTPDbSWsHGhCkrwA5qFhvFdSen8tXSTazf43Y6jjEmyKwAO+z6YSk0jY7g89UHnI5ijAkyK8AOS2gazdVDk5m3sYrVW0udjmOMCSIrwCHgxlPSiHTBC9MLnI5ijAkiK8AhoE2zGEZ1juTjhetZu2Of03GMMUFiBThE/CI1ChF45ItlPDct35YtMiYMOLYihvm5VrEuRvVoy+Slm/l62WaiI128c2Mm6ckJDTpvbnEJc1dvJzOtdYPPZYzxLyvAISQ1sSkAbv3v2nENKZrzC7dz9avzqXS7/VbQjTH+Y10QIeTMvu2IivAsk6dARgOKZXllFfd9vIQDVW7cCgdsMVBjQo4V4BCSnpxA9rihnNUnCbfCm3OKqXLXfc7g/QeqGPdmLqu2lBLhXZLZrRDpOtoaqMYYJzSoC0JEMoDhQAdgP7AEmKqqNsVXPaUnJ5B+TQavfreaR75YTnxMJI9d3A8R34rnnrIKbhifw4KiHfz94n50a9uMr5Zs5Jvlm/nbpBVs33uAP5zZk6gI+91rjNPqVYBF5DrgN0AhkAvkAbHAKcDdIrIE+LOqrvFTzrBz4/A0du+v4Jlv82kWG8l9Z/c6ZhHeue8A174+n6UbdvOvKwZy3okdAE9Rv+vMnjz6xXJenrmanKId/HvsSXRsGReMP4ox5ijq2wJuCpysqvtr2ikiA4DugBXgBvjd6B7sLqvk1VmFtIiL4vbTuh/12K17yrnmtXms3rqXF65OZ3TvpJ/tj4mM4KHz+zIktTV3f7iYs5/5jqcvO5FTT0g6yhmNMYFWr8+hqvrc0Yqvd/8iVf2m/rEMgIjwl3N6c9HAjjz19Ure+L6wxuM27NzPZS/NoXj7Pl6/btARxbe6s/u3Z+Ltp9ChRRy/eiOHv01aTkWVTQRkjBPq2wXxTG37VfU3Pp7ndeAcYIuq9vVuewC4CTi4bPAfVfXL+uRsDFwu4fFL+lNaXskDny+jWWwUF6d3OrS/aNternp1Hrv3V/DWDYPJSGl1zHOmJDblo18P45EvlvHSjNXkFpVw04g08reU2nhhY4Kovt/E5HofscBJwCrvYwBQVYfzvAGcVcP2f6jqAO8jbIvvQZERLp65ciAnd2vN/364mMlLNwGwcvMeLn1pDvsOVDJhXKZPxfeg2KgIHrmgH89cOZAl63dx81u5PDk5j6tenWt34RkTJPXtghivquPx9PNmqeq/VfXfwGl4irCv55kJ2IgJH8RGRfDyNRn069iC299dyEOfL+W8Z2dRWeXmvZuH0rdji3qd97wTO3B1ZjLgGXtcVuFmRt4WPyY3xhyNqNZ9nOmhN4vkAUMPDjsTkQRgrqr2rMM5UoCJh3VBXAfsBnKA36vqEU0yERkHjANISkpKz87OPuLcpaWlxMfH1+nP5BRfs5YeUB6cs4+t3h74KBfcPSiWbgkR9b52fkkVf19QRoW3K7hZNNzSP5Y+iTWfszH+vYYCyxoYoZA1KysrV1UzDt/e0FuRHwMWisg07+uRwAMNPOcLwMN4GmQPA08Bvzr8IFV9GXgZICMjQ0eNGnXEiaZPn05N20NRXbIuc6/ghekFKJ4bLMpbJjNqVLd6X3sUMPAkz5wRifHRvDxzNU/k7OWazGTuHXMCTaJ//mPiS9ZQmYOisf4MOM2y+keDCrCq/p+ITAKGeDfdo6qbGnjOzQefi8grwMSGnK8xOq1XEq9/X0hFpZuoSBeZaa0bfM705IRDhfL8AR15cnIer31fyMxVW3nq0hN97l+urHKTPX8t93++lCq3Eh3hYsI4m4PCmJr4YzKecmAjni/keohID2/fbr2ISHtV3eh9eSGeu+tMNenJCbxzY2bAWpixURH86ZzenN47ibve/5FLX5rDuOFp/G50D2KjjuyW2LnvADNWbuXbFVuYnreVXfsrDu07UOXmb18u5+0bh9T4XmPCWUNvRb4RuAPoBCwCMoE5wKk+vn8Cnk/AiSKyDrgfGOW9kUOBIuDmhmRsrKq3WAMlM601X/12BI9+sZyXZq5mWt4WbhqextyCA+xssZ5Nu8v4dvkWcop34FZo3TSa03slkZrYhGe/zfeMLxYhp7iE856dxT8uH0CfDvX7stCYxqihLeA7gEF4vnjLEpETgAd9fbOqXlnD5tcamMn4UXxMJH+7qB9n9EnizvcW8YcPFgPw4apFAPRu35zbsrpx6gltObFTS1zeCX+Gdk081ELfU1bB/36wmAue+57fnt6DW0Z2PTRJkDHhrKEFuExVy0QEEYlR1RUi4vMICHP8yOrZlqszk/n3t/kACHDzyDTu+UWvGo8/vIU++bcj+NMnS3hich7TVmzh6csG0KV1k2BENyZkNXRKrHUi0hL4BPhaRD4FNjQ8lglFo3q2JTbKhQuIiXIxunc7n9+b0DSaZ8cO5J+XDyBv8x7O+tdMJsxfQ0OGQRpzvGvoKIgLvU8f8A5FawF81eBUJiQd/PJvwtQFXHn6oDr3QYsIFwzsyODUVtz1/o/c+9FPTF22mbFDurBi0x7Hh6wZE2z1LsAi4gIWH7yBQlVn+C2VCVnpyQns6RrdoELZoWUcb98whDdmF/G3Scv5ZsUWBE+r2pZNMuGk3l0QquoGfhSRLn7MY8KEyyX86pRUrql2G7Qtm2TCTUP7gNsDS0XkGxH57ODDH8FMeDi7fwdiIz0/hm6Fts1iHE5kTPA0dBSEz0POjKlJenIC79yUyTfLN/OfnLU8PjmPzLTWdG5lIyRM41ff+YBFPY7a73vwmPpHM+Hi4JC1Cwd25OIXZnP9Gwv48JZhtGgS5XQ0YwKqvl0Q1fWlmwAAGohJREFU00Tk9sP7f0UkWkROFZHxwLUNj2fCSfekZrz8ywzWbN/HuLdyKK+sy9TSxhx/6luAz8Iz8foEEdkgIstEZDWeSdmvxDOh+ht+ymjCSGZaa564tD/zCndw1/uLcbvtQ5RpvOrVBaGqZcDzwPMiEgUkAvtVdac/w5nwdP6AjmzYWcbfv1pBx5Zx3POLE5yOZExANHg2NFWtwDMbmjF+c8vINNbv3MeLMwromBB3aLiaMY2JP6ajNMbvRIQHzu3Dxp1l3P/pEto3j+X0WlZ7NuZ41NBxwMYETGSEi3+PHUjfji24fcJCflxrPVymcWlwARaRZBE53fs8TkSaNTyWMR5NoiN57dpBJDaL5v/bu/P4qOpzj+OfZ2ayScKWQEBAIKxlh7AqVaJSUbFaREW9aguKVvG29dri9bbq9bZXW8t1qStYBEXgauuCYAVRQEUoEJAlsgTCvu8QhITJPPePmXBHSMgyMzknyfN+vfJK5szJyTc/wpOT35zz/EZNXsq2g985HcmYqImoAIvIPcDfgNdCm5oT7IxmTNQ0Sklg0s/64g8oI8Yv4s+z15O99Zx1Wo2pdiI9A34AuITgCsaoai7QONJQxpytTaNkHv5RB3YdPcWL8zZy+4TFVoRNtRdpAS5Q1cLiByLiI9hXxZioO3ryNMXraJzyB/h60wFH8xgTqUgL8AIReRRIEpHBwLvAR5HHMuZc/TNSSYjznCnCq7YfsYbuplqLtAA/AuwHVhNcPPNj4LeRhjKmJMUN4R++qgM39mrGp2v38ezcXKdjGVNpka6IEQAmABNEpCHQ3BrwmFgqbtyjqng9wguf5dI4JYF/sRs1TDUU6bL084Efh47zDbBfRBao6kNRyGZMqUSE//5JVw7mF/LYh2tIS05gSJfyr1FnjBtEOgVRT1WPAcOAN1Q1E7gy8ljGlM3n9fDibb3o3qI+/zp9BUs2H3I6kjEVEmkB9olIU+BmYGYU8hhTIUnxXibe1YfmDZK4e/JS1u857nQkY8ot0gL8JDAb2KiqS0Ukg2BLSmOqTIM68bw5si9J8V7umriEnUdOOh3JmHKJqACr6ruq2k1V7w89zlPVG8v7+SIyUUT2iciasG0NReRTEckNvbclck2Zmje4gMkj+3Ki0M9dE5dw5LvCsj/JGIdFeityoog8ICIvh4rpRBGZWIFDTCLY3D3cI8BnqtoO+Cz02JgydWxSlwl39mbboe8YOWkpJwttRQ3jbpFOQbwFNAGuAhYQ7AVR7kk4Vf0COPuVk+uByaGPJwM3RJjR1CL9M1J5/pYerNh+hDv++k9mbCy0W5aNa0kkl+2KyApV7Skiq1S1W2h1jNmqenkFjtEKmKmqXUKPj6hq/bDnD6vqOdMQIjIaGA2Qnp6eOX369HOOnZ+fT3JyckW/LUdY1uiavq6AT7b4ASXOI4ztk0jbBl6nY51XdRjXYpa1YrKysrJVtffZ2yNtyH469P6IiHQB9gCtIjxmuajqeGA8QO/evXXQoEHn7DN//nxK2u5GljW6cnQjs7esRxFOByA/pQWDBrV3OtZ5VYdxLWZZoyPSKYjxoRfJfgfMAL4F/hThMfeGLm0j9H5fhMcztdDZfSNmrtrNwfwCRzMZc7ZIr4J4XVUPq+oCVc1Q1caq+mqEmWbw/0va3wV8GOHxTC1U3DfixnZxjB3Sge2HvmP4q4usobtxlUhvRU4AbiQ47XDmWKr6ZDk/fxowCEgTkR3A48DTwDsiMgrYBtwUSUZTe2W2bMDxNvEMGtSWvq0bMmryMoa9spA3ftqXrs3rOR3PmIinID4keNWCHzgR9lYuqnqrqjZV1ThVba6qf1XVg6p6haq2C723+0tNxDJbNuRv911Mgs/LiPGL+GLDfqcjGRPxi3DNVfXs63iNcaW2jZN57/6LuWviEkZOWsqfhndjWK/mTscytVikZ8Bfi0jXqCQxpgqk103knfsG0KdVQx56ZyWvzN9kTd2NYyp1BiwiqwkuPeQDfiYieUABIICqarfoRTQmuuomxjFpZB8efncVf/xkHXuPneJ3Qzvh9UjZn2xMFFV2CmJoVFMYU8USfF6ev6UH6SkJvP7VZtbvOUb/jFQGtmtEZktrP2KqRmUL8F7gPqAtweWI/qqq/qilMqYKeDzCb4d2wh8IMOnrrSzKO8TL8zcx9Z7+VoRNlajsHPBkoDfB4ns1MC5qiYypYo1SEs/csFHgDzDjm52O5jG1R2XPgDupalcAEfkrsCR6kYypWsV3zRX6AwQUpi3ZTt/WqVzbranT0UwNV9kCXNwDAlX1i9iLF6b6Kr5rbnHeQTo0SeGV+Zt4YOpyvt3dhn8b3AGPvThnYqSyBbi7iBwLfSxAUuhx8VUQdaOSzpgqUrzaMsAP26Xx2Ac5vDRvE+t2H+fZET2omxjncEJTE1VqDlhVvapaN/SWoqq+sI+t+JpqLcHn5ekbu/Lk9Z2Zv2E/P3lpIXn7852OZWqgSG/EMKZGEhHuHNCKKaP6cehEIde/tJD5660xn4kuK8DGnMeANqnMGDOQZvWTGDlpKa8tsDvnTPRE2gvCmBqvRcMLeO/+i/n1u6t46h/r+Hb3MUb0acHybUfon5Fq1wybSrMCbEw5XBDv48XbetJpfl2emb2eGSt3IUC8z8Pbd9uNG6ZybArCmHISER7Iasuwns1QhYBCoT/A4ryDTkcz1ZQVYGMq6Pb+LUnwBf/rBBRydh2lwF/kcCpTHdkUhDEVlNmyAVPv6c/CjftZv+c4s1bvYfuhRbx4W09aptZxOp6pRqwAG1MJ4Tdu/DhnD79+dyXXvvAVT9/YlaHdLnQ4nakubArCmAhd1bkJH//ih7RLT2bM1BU8+v5qTp22KQlTNivAxkRB8wYX8M69A7j3sgym/nMbN7y0kE1295wpgxVgY6Ikzuvh36/+AW/8rA/7jhdw3V++4r3lO5yOZVzMCrAxUZbVoTEf/+sP6dKsHg+9s5KfvbGE5+ZuIHvrYaejGZexAmxMDDSpl8jUu/txU2Yz5q3fz3Nzc7l1wmIrwuZ7rAAbEyM+r4dWackUtxMu9Ad4ZvY6u2bYnGEF2JgY6p+RSrzPg1fAK8LivENc/+JCcnYddTqacQHXXgcsIluA40AR4FfV3s4mMqbiwlfb6J+RypHvCnnkvdVc/+JCHry8HfdntSHOa+dBtZVrC3BIlqoecDqEMZEIv2kDYM4vG/DERzk8O3cDc9fuZdzN3WmfnuJgQuMU+9VrTBVrUCee50f05JXbe7HryEmGvvAVr8zfRMD6DNc64tbm0iKyGTgMKPCaqo4/6/nRwGiA9PT0zOnTp59zjPz8fJKTk6sgbeQsa2y4PeuxAuXNbwtYtreIVinKvd0voGmy+8+L3D6u4dyQNSsrK7ukaVQ3F+ALVXWXiDQGPgUeVNUvStq3d+/eumzZsnO2z58/n0GDBsU2aJRY1tioDllVlRkrd/Ho37/Brx5u63cRqXXiGdAmzbV9hqvDuBZzQ1YRKbEAu3YOWFV3hd7vE5H3gb5AiQXYmOpMRLi+RzN0z3qmbE7kjYVbAIj3bWTaPdbsvSZz5d86IlJHRFKKPwZ+BKxxNpUxsVU/0UNWx0aELhum0B/gPz/K4dCJQkdzmdhxZQEG0oGvRGQlsASYpaqfOJzJmJjrn5FGQlzoumGPsGbnUQY9M4/JX2/BXxRwOp6JMldOQahqHtDd6RzGVLWzrxuum+jjiY9yeHxGDtOWbOOJH3emf0aq0zFNlLiyABtTm5193fCUUf2YnbOH/5q5lhHjFzO0W1MeveYHXFg/ycGUJhrcOgVhjAkREYZ0acrchy7jF1e049Nv93LFuAW8+HmuNX6v5uwM2JhqIiney68Gt2d4ZnP+MGstf56zgXeW7eD2fhfhDwTon+Hey9ZMyawAG1PNtGh4Aa/ekclXuQcY+94qnvrHOgC8nlwevLwtwzOb06x+EiJSxpGM06wAG1NNDWyXxog+LfifORtQoCigPDc3l+fm5pJeN4HeLRvSKzSf3PnCutb0x4WsABtTjV3cJo2X4jZy2h8gzufh9zd04bvCIrK3HmbZlsPMWr0bgMQ4D92a16d3ywbUS4rj6MnTXNGxMZmtGjr8HdRuVoCNqcbOvmyteA74zgGtANhz9BTLtwWLcfa2w7y2YBNFoe4DL8/fRNN6ibRplEyz+kk0a5BEs/pJXFg/ieYNkmhSL/HMWXP21sPnfA0TOSvAxlRzZ1+2Fq5JvUSu6dqUa7o2BeD5uRt4/rNcAgoC1E2MI7/Az+fr97H/eMH3PtcjkF43kXpJceTuzUdR4n0e3r67f6y/pVrDCrAxtcjAdo14ZcGmM1MW/z2s65nifep0EbuPnmLn4ZPsOnKSHUdOsvPwSZZsPkhRqGnXqdMBPlmzm0vqOPld1BxWgI2pRUqbsgBIjPPSOq0OrdO+X12ztx7m9gmLKfAHUGDiws1sa+GjR99C6l8QX8XfQc1iBdiYWuZ8Uxal7f/2PcGi3SE9hblr9/K/S7ez+Jn5PHh5W+4c0Ip4n11hURlWgI0xZQov2ld2Sqdz/AHm7KvD72et5a3FWxk7pCNXd2li1x5XkP3aMsZUWIsUD2+N6sfkkX1J9Hm5/+3lDH91Ecu3HXY6WrViZ8DGmEq7rH0jLmmTyrvZOxg3ZwPDXv6aod2acm3XpuQdOGGXrZXBCrAxJiI+r4db+17Edd0vZPyCTby6YBMzV+1GgASfh7dtVY9S2RSEMSYqkhN8PPSjDowc2BoIrqZ7yh9g5qpdzgZzMSvAxpioGtypCYlxHjyh1+OmLN7K5K+34NYFgJ1kUxDGmKgKv9a4Y5MUpizeyuMzcpi/fh9/Gt6dRikJTkd0DSvAxpioC79s7fKOjXlr8Vb+MGstQ577gmdu6sblHdMdTugONgVhjIkpEeHOAa346MGBNEpJYOSkZTz24RpbzQMrwMaYKtI+PYUPHriEUQNb8+airVz3l6/4dtcxp2M5ygqwMabKJMZ5+d3QTrw5si9HTp7mhpcW8vqXeQQCtfMFOivAxpgqd2n7Rsz+5aVc1qERv5+1lrveWMLcb/fy0ryNZG+tPXfT2YtwxhhHNKwTz/g7Mpm2ZDtPzFjDl7kHgjdvxAV7DteGmzfsDNgY4xgR4bZ+F/EvoRU8FCg4HWBx3kFHc1UV1xZgERkiIutFZKOIPOJ0HmNM7FzbtSmJoZaWCuSfOu1soCriygIsIl7gJeBqoBNwq4h0cjaVMSZWinsO/+rKdvS8qD6vLMjjzUVbnI4Vc64swEBfYKOq5qlqITAduN7hTMaYGMps2YBfXNme6aP7M7hTOo99mMPrX+Y5HSumxI33Z4vIcGCIqt4denwH0E9Vx4TtMxoYDZCenp45ffr0c46Tn59PcnJy1YSOkGWNDcsaG7HO6g8or60qYOmeIm5sF8d1bSq/9JEbxjUrKytbVXuf84Squu4NuAl4PezxHcBfSts/MzNTSzJv3rwSt7uRZY0NyxobVZH1tL9IfzFtubYcO1PHzVmvgUCgUsdxw7gCy7SE2uXWy9B2AC3CHjcHrKedMbWIz+th3M09iPN6eOGzXAr9AcYO6VCjlj1yawFeCrQTkdbATmAEcJuzkYwxVc3rEf54YzfifR5eXbCJQn+A3w39QY0pwq4swKrqF5ExwGzAC0xU1RyHYxljHODxCL+/oQvxPg8TF26msKiIJ3/cBY+n+hdhVxZgAFX9GPjY6RzGGOeJCI8N7USCz3vmTPipYd3wVvMi7NoCbIwx4USEsUM6EO8LzgmfLlKeGd4Nn7dqrqadt34f/8w7yOBOTaJ2m7QVYGNMtSEiPDS4PQk+D8/MXs/eo6cY0DaVi9ukxbR3xLIthxg1aSkBhUkLt0RtoVG33ohhjDGleiCrLT+9uCVf5x1k3JwN3P764ph2UXvx840Ud8w8XRS9XhVWgI0x1VKjlESKZ4Bj2cBn+bbDfJm7H4+AVyDO56F/RmpUjm1TEMaYaql/RioJcR4KTgdQIBavxx397jQPTl1Bk3pJPDWsK6t3HqV/RqrNARtjarfi1Ze/yt3PjJW7eHneJq7q3ISMRtG57VhVGfv3Vew9dop37htAr4sacGn7RlE5djGbgjDGVFvFDXzeHNWPOJ+He9/K5kSBPyrHnrJ4K5/k7OE3QzrQ66LYvMBnBdgYU+01q5/EX27tyab9+fzmb6uKe8hUWs6uo/zXzLVkdWjE3QMzopTyXFaAjTE1wiVt0/jNkI7MWr2bCRG0scwv8DNm6goa1Ilj3M09YnrHnRVgY0yNce+lGVzdpQlP/2MdX288UOHPV1V++/5qth48wfMjetKwTuXbYJaHFWBjTI0hIjxzU3cyGiUzZtoKdh45WaHPfzd7Bx98s4tfXtk+apeanY8VYGNMjZKc4OO1OzIp9Ae4f0o2hUXlmw/O3Xucxz5cw8VtUnkgq22MUwZZATbG1DhtGiUz7uburNxxlClrC8vc/2RhEWOmrqBOvI/nbulRZU1+rAAbY2qkqzo3YUxWW77Y4Wfakm3n3ffJmTms33ucZ2/pQeO6iVWU0AqwMaYG+9Xg9nRJ8/L4hzms2FZyr4gZK3cxbcl2fj6oTdRvtCiLFWBjTI3l9Qj3dUsgvV4CP5+ynP3HC773/JYDJ3j0vdVktmzAQ4PbV3k+K8DGmBotOV545fZMDn9XyIPTluMvCgBQ4C9izLTleD3CC7f2JK6K+gqHswJsjKnxujSrx1PDurI47xB//GQdAE99vI41O4/x55u606x+kiO5rBmPMaZWGNarOSu3H2HCl5vZcuA7Pl27l2u7NmFwp3THMtkZsDGm1viPazvRsUkKn67dC8Bn6/bFtJF7WawAG2NqjXifh8s7Nj7z+LQ/do3cy8MKsDGmVrniB+kkxnmivrpFZdgcsDGmVilu5L4472BUV7eoDCvAxphaJ7NlA0cLbzGbgjDGGIe4rgCLyBMislNEvgm9XeN0JmOMiQW3TkE8q6p/djqEMcbEkuvOgI0xpraQSBevizYReQL4KXAMWAb8m6qec6W0iIwGRgOkp6dnTp8+/Zxj5efnk5wcnSWqY82yxoZljQ3LWjFZWVnZqtr7nCdUtcrfgLnAmhLergfSAS/Bs/M/ABPLOl5mZqaWZN68eSVudyPLGhuWNTYsa8UAy7SE2uXIHLCqXlme/URkAjAzxnGMMcYRrnsRTkSaquru0MOfEDwzPq/s7OwDIrK1hKfSgIovjeoMyxobljU2LGvFtCxpoxvngN8CegAKbAHuDSvIFT3WMi1p3sWFLGtsWNbYsKzR4bozYFW9w+kMxhhTFewyNGOMcUhNL8DjnQ5QAZY1NixrbFjWKHDdHLAxxtQWNf0M2BhjXMsKsDHGOKRaFmARmSgi+0TkvNcIi0gfESkSkeFh24aIyHoR2Sgij7g86xYRWR3qCrfM6awiMkhEjoZ1qnss7DlXjWsZWV01rmF5vxGRHBFZELbdVeNaRlZXjauI/Drs339N6P9Xw9BzVTqupSrp9ji3vwGXAr2ANefZxwt8DnwMDA/btgnIAOKBlUAnN2YNbd8CpLllXIFBwMxS8rtqXEvL6tJxrQ98C1wUetzYxeNaYlY3jutZ+14HfO7UuJb2Vi3PgFX1C+BQGbs9CPwd2Be2rS+wUVXzVLUQmE6w/0TMRJC1ypUza0ncOq6uUI6stwHvqeq20P7FPwduHNfSsla5Cv4M3ApMC31c5eNammpZgMsiIs0I3sb86llPNQO2hz3eEdrmmPNkheDdgHNEJDvU/c0NBojIShH5h4h0Dm1z3biGlJQV3Deu7YEGIjI/lOnO0HY3jmtpWcF94wqAiFwADCF4kgMuGlfX3QkXJc8BY1W1SETCt0sJ+zp9HV5pWQEuUdVdItIY+FRE1oV+6ztlOdBSVfNDK5V8ALTDneNaWlZw37j6gEzgCiAJWCQii3HnuJaYVVU34L5xLXYdsFBVi8+WXTOuNbUA9wamhwpaGnCNiPgJ/qZrEbZfc2BX1cf7nhKzquoHqroLgn/micj7BP90cuwHWlWPhX38sYi8LCJpuHBcS8uqqgfcNq4Ex++Aqp4ATojIF0B3XDiulJ51gwvHtdgI/n/6AVw0rjVyCkJVW6tqK1VtBfwNuF9VPwCWAu1EpLWIxBP8h5nhYNRSs4pIHRFJARCROsCPKEdnuFgSkSYS+k0hIn0J/vwcxIXjWlpWN44r8CHwQxHxhf5c7gesxYXjSilZXTquiEg94DKCuYu5Zlyr5RmwiEwj+Cp3mojsAB4H4gBUtaS5VELP+UVkDDCb4CuhE1U1x41ZCTamfz9UQ3zAVFX9xOGsw4Gfh/6aOAmM0ODLym4c1xKziojrxlVV14rIJ8AqIAC8rqprQp/rqnEtLauIZOCycQ3t9hNgTuiMndBzVf7zWhq7FdkYYxxSI6cgjDGmOrACbIwxDrECbIwxDrECbIwxDrECbIwxDrECbIwxDrECbIwxDrECbGokEekqIltF5Ofn2SdJRBaIiDcKXy9eRL4QkWp5c5NxhhVgUyOp6mqCt5jeeZ7dRhJsrVgUha9XCHwG3BLpsUztYQXY1GT7gM7nef52Qj0CRKSVhK2sICIPi8gToe3rROR1Ca6q8LaIXCkiC0UkN9RnotgHoWMaUy5WgE1N9jSQICItz34i1IQlQ1W3lOM4bYHngW5AR4JNyQcCDwOPhu23BugTYWZTi1gBNjWSiAwB6gCzKPksOA04Us7DbVbV1aoaAHKAz0JNiFYDrYp3Ck1lFBZ3BTOmLFaATY0jIonAn4D7CRbJLiXsdhJIDHvs5/v/H8KfKwj7OBD2OMC5HQUTgFMVT21qIyvApib6LfBmaHqhxAKsqocBb6hYA+wFGotIqogkAEMr+kVFJBXYr6qnK53c1CpWgE2NIiIdgMEEl3qC0s+AAeYQnMslVDSfBP4JzATWVeLLZxFc2dqYcrF+wKbWEpGewEOqekeUjvce8O+quj4axzM1n50Bm1pLVVcA86J1IwbwgRVfUxF2BmyMMQ6xM2BjjHGIFWBjjHGIFWBjjHGIFWBjjHGIFWBjjHGIFWBjjHHI/wGdOpiAz3aMngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot S11\n",
    "f, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(5, 8))\n",
    "ax1.plot(1/mode1_freqs, 10 * np.log10(S11_mag), '.-')\n",
    "ax1.set_title(\"S11\")\n",
    "ax1.set_xlabel(r\"$\\lambda$ (um)\")\n",
    "ax1.set_ylabel(\"Magnitude (dB)\")\n",
    "ax1.set_ylim(None, 0)\n",
    "ax1.grid()\n",
    "\n",
    "ax2.plot(1/mode1_freqs, S11_phase, '.-')\n",
    "ax2.set_xlabel(r\"$\\lambda$ (um)\")\n",
    "ax2.set_ylabel(\"Phase (rad)\")\n",
    "ax2.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Milestones\n",
    "\n",
    "Goal: Compute the transmission profile for bend radii between 1.5um and 10um. \n",
    "\n",
    "- Q: Is the reflection significant for any radius? What explain the loss?\n",
    "- Q: What is the formula total size of the simulation region? How many pixels are there?\n",
    "- Q: If each pixel can host 3-dimensional E-field and H-field vectors with 64bit complex float stored in each dimension, how many megabytes of data needs to be stored at each time step? Is it feasible to save all this information throughout the FDTD simulation?\n",
    "- Bonus: Collect the simulation runtime for each radius. How does it change with different radii?\n",
    "- Bonus: At what resolution does the accuracy of the simulation start degrading? In other words, if halving the resolution only results in a 1% relative difference in the most important target metric, it is still a good resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write to csv file\n",
    "import csv\n",
    "with open(f'sparams.r{ring_radius:.1f}um.csv', mode='w') as sparams_file:\n",
    "    sparam_writer = csv.writer(sparams_file, delimiter=',')\n",
    "    sparam_writer.writerow(['f(Hz)','real(S11)','imag(S11)','real(S21)','imag(S21)'])\n",
    "    for i in range(len(mode1_freqs)):\n",
    "        sparam_writer.writerow([mode1_freqs[i] * 3e14,\n",
    "                                np.real(S11[i]),np.imag(S11[i]),\n",
    "                                np.real(S21[i]),np.imag(S21[i])\n",
    "                               ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guide for automation\n",
    "\n",
    "We have transcribed the code in every cell into a file called `bent_waveguide.py`, and we have added an argument parser to help you.\n",
    "\n",
    "Simply call `python bent_waveguide -r 4.5`, for example, to run the computation with a radius of 4.5 um, outputting to a csv file called `sparams.r4.5um.csv` for further processing.\n",
    "\n",
    "If MPI is activated, and you have enough processor cores, run it with `mpirun -np 4 python bent_waveguide -r 4.5` for a 4x speedup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "all_tools [~/home/ELE559/anaconda3/envs/all_tools/]",
   "language": "python",
   "name": "conda_all_tools"
  },
  "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": 5
}