{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resonator benchmark (COMSOL)\n",
    "\n",
    "See on [github](https://github.com/flexcompute/tidy3d-notebooks/blob/main/HighQ_Si.ipynb), run on [colab](https://colab.research.google.com/github/flexcompute/tidy3d-notebooks/blob/main/HighQ_Si.ipynb), or just follow along with the output below.\n",
    "\n",
    "In this example, we reproduce the findings of Zhang et al. (2018), which is linked [here](https://www.osapublishing.org/ol/abstract.cfm?uri=ol-43-8-1842).\n",
    "\n",
    "This notebook was originally developed and written by Romil Audhkhasi (USC). \n",
    "\n",
    "The paper investigates the resonances of silicon structures by measuring their transmission spectrum under varying geometric parameters.\n",
    "\n",
    "The paper uses a finite element solver (COMSOL), which matches the result from Tidy3D.\n",
    "\n",
    "<img src=\"img/Si_struct.png\" alt=\"diagram\" width=\"500\"/>\n",
    "\n",
    "(Citation: Opt. Lett. 43, 1842-1845 (2018).  With permission from the Optical Society)\n",
    "\n",
    "To do this calculation, we use a broadband pulse and frequency monitor to measure the flux on the opposite side of the structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:42.649437Z",
     "iopub.status.busy": "2021-11-12T19:03:42.648542Z",
     "iopub.status.idle": "2021-11-12T19:03:43.853353Z",
     "shell.execute_reply": "2021-11-12T19:03:43.853026Z"
    }
   },
   "outputs": [],
   "source": [
    "# get the most recent version of tidy3d\n",
    "!pip install -q --upgrade tidy3d\n",
    "\n",
    "# make sure notebook plots inline\n",
    "%matplotlib inline\n",
    "\n",
    "# standard python imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# tidy3D import\n",
    "import tidy3d as td\n",
    "from tidy3d import web"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set Up Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.856324Z",
     "iopub.status.busy": "2021-11-12T19:03:43.855909Z",
     "iopub.status.idle": "2021-11-12T19:03:43.865321Z",
     "shell.execute_reply": "2021-11-12T19:03:43.864991Z"
    }
   },
   "outputs": [],
   "source": [
    "nm = 1e-3\n",
    "\n",
    "# define the frequencies we want to measure\n",
    "Nfreq = 1000\n",
    "wavelengths = nm * np.linspace(1050, 1400, Nfreq)\n",
    "freqs = td.constants.C_0 / wavelengths\n",
    "\n",
    "# define the frequency center and width of our pulse\n",
    "freq0 = freqs[len(freqs)//2]\n",
    "freqw = freqs[0] - freqs[-1]\n",
    "\n",
    "# Define material properties\n",
    "n_SiO2 = 1.46\n",
    "n_Si = 3.52\n",
    "SiO2 = td.Medium(epsilon=n_SiO2**2)\n",
    "Si = td.Medium(epsilon=n_Si**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.868380Z",
     "iopub.status.busy": "2021-11-12T19:03:43.867965Z",
     "iopub.status.idle": "2021-11-12T19:03:43.876850Z",
     "shell.execute_reply": "2021-11-12T19:03:43.876549Z"
    }
   },
   "outputs": [],
   "source": [
    "# space between resonators and source\n",
    "spc = 1.5\n",
    "\n",
    "# geometric parameters\n",
    "Px = Py = P = 650 * nm  # periodicity in x and y\n",
    "t = 260 * nm            # thickness of silcon\n",
    "g = 80 * nm             # gap size\n",
    "L = 480 * nm            # length in x\n",
    "w_sum = 400 * nm        # sum of lengths in y\n",
    "\n",
    "# resolution (should be commensurate with periodicity)\n",
    "dl = P / 32\n",
    "dz = 20 * nm\n",
    "\n",
    "# computes widths in y (w1 and w2) given the difference in lengths in y and the sum of lengths\n",
    "def calc_ws(delta):\n",
    "    \"\"\" delta is a tunable parameter used to break symmetry.\n",
    "        w_sum = w1 + w2\n",
    "        delta = w1 - w2\n",
    "        w_sum + delta = 2 * w1\n",
    "    \"\"\"\n",
    "    w1 = (w_sum + delta) / 2\n",
    "    w2 = w_sum - w1\n",
    "    return w1, w2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.880213Z",
     "iopub.status.busy": "2021-11-12T19:03:43.879785Z",
     "iopub.status.idle": "2021-11-12T19:03:43.889861Z",
     "shell.execute_reply": "2021-11-12T19:03:43.889518Z"
    }
   },
   "outputs": [],
   "source": [
    "# total size in z and [x,y,z]\n",
    "Lz = spc + t + t + spc\n",
    "sim_size = [Px, Py, Lz]\n",
    "\n",
    "# sio2 substrate\n",
    "substrate = td.Box(\n",
    "    center=[0, 0, -Lz/2],\n",
    "    size=[td.inf, td.inf, 2*(spc+t)],\n",
    "    material=SiO2,\n",
    "    name='substrate'\n",
    ")\n",
    "\n",
    "# creates a list of structures given a value of 'delta'\n",
    "def geometry(delta):\n",
    "    w1, w2 = calc_ws(delta)\n",
    "    center_y = (w1 - w2) / 2.0\n",
    "\n",
    "    cell1 = td.Box(\n",
    "        center=[0, center_y + (g + w1)/2., t/2.],\n",
    "        size=[L, w1, t],\n",
    "        material=Si,\n",
    "        name='cell1'\n",
    "    )\n",
    "\n",
    "    cell2 = td.Box(\n",
    "        center=[0, center_y - (g + w2)/2., t/2.],\n",
    "        size=[L, w2, t],\n",
    "        material=Si,\n",
    "        name='cell2'\n",
    "    )\n",
    "\n",
    "    return [substrate, cell1, cell2]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.892517Z",
     "iopub.status.busy": "2021-11-12T19:03:43.892065Z",
     "iopub.status.idle": "2021-11-12T19:03:43.900701Z",
     "shell.execute_reply": "2021-11-12T19:03:43.900973Z"
    }
   },
   "outputs": [],
   "source": [
    "# time dependence of source\n",
    "gaussian = td.GaussianPulse(freq0, freqw)\n",
    "\n",
    "# plane wave source\n",
    "source = td.PlaneWave(\n",
    "            source_time=gaussian,\n",
    "            injection_axis='-z',\n",
    "            position= Lz/2 - spc + 2*dl,\n",
    "            polarization='x')\n",
    "\n",
    "# Simulation run time.  Note you need to run a long time to calculate high Q resonances.\n",
    "run_time = 7e-12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.903466Z",
     "iopub.status.busy": "2021-11-12T19:03:43.903138Z",
     "iopub.status.idle": "2021-11-12T19:03:43.911931Z",
     "shell.execute_reply": "2021-11-12T19:03:43.911566Z"
    }
   },
   "outputs": [],
   "source": [
    "# monitor fields on other side of structure (substrate side) at range of frequencies\n",
    "monitor = td.FreqMonitor(\n",
    "        center=[0., 0., -Lz/2 + spc - 2 * dl],\n",
    "        size=[td.inf, td.inf, 0],\n",
    "        freqs=freqs,\n",
    "        store=['flux'],\n",
    "        name='transmitted_fields')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Case Studies\n",
    "\n",
    "Here we define the three simulations to run\n",
    "\n",
    "- With no resonators (normalization)\n",
    "- With symmetric (delta = 0) resonators\n",
    "- With asymmetric (delta != 0) resonators\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.915324Z",
     "iopub.status.busy": "2021-11-12T19:03:43.914989Z",
     "iopub.status.idle": "2021-11-12T19:03:43.936336Z",
     "shell.execute_reply": "2021-11-12T19:03:43.935931Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing simulation...\n",
      "Mesh step (micron): [2.03e-02, 2.03e-02, 2.00e-02].\n",
      "Simulation domain in number of grid points: [32, 32, 206].\n",
      "Total number of computational grid points: 2.11e+05.\n",
      "Total number of time steps: 19987.\n",
      "Estimated data size (GB) of monitor transmitted_fields: 0.0000.\n",
      "Initializing simulation...\n",
      "Mesh step (micron): [2.03e-02, 2.03e-02, 2.00e-02].\n",
      "Simulation domain in number of grid points: [32, 32, 206].\n",
      "Total number of computational grid points: 2.11e+05.\n",
      "Total number of time steps: 199868.\n",
      "Estimated data size (GB) of monitor transmitted_fields: 0.0000.\n",
      "Initializing simulation...\n",
      "Mesh step (micron): [2.03e-02, 2.03e-02, 2.00e-02].\n",
      "Simulation domain in number of grid points: [32, 32, 206].\n",
      "Total number of computational grid points: 2.11e+05.\n",
      "Total number of time steps: 199868.\n",
      "Estimated data size (GB) of monitor transmitted_fields: 0.0000.\n"
     ]
    }
   ],
   "source": [
    "# normalizing run (no Si) to get baseline transmission vs freq\n",
    "# can be run for shorter time as there are no resonances\n",
    "sim_empty = td.Simulation(size=sim_size,\n",
    "                    mesh_step=[dl, dl, dz],\n",
    "                    structures=[substrate],\n",
    "                    sources=[source],\n",
    "                    monitors=[monitor],\n",
    "                    run_time=run_time/10,\n",
    "                    pml_layers=[0,0,15])\n",
    "\n",
    "# run with delta = 0\n",
    "sim_d0 = td.Simulation(size=sim_size,\n",
    "                    mesh_step=[dl, dl, dz],\n",
    "                    structures=geometry(0),\n",
    "                    sources=[source],\n",
    "                    monitors=[monitor],\n",
    "                    run_time=run_time,\n",
    "                    pml_layers=[0,0,15])\n",
    "\n",
    "# run with delta = 20nm\n",
    "sim_d20 = td.Simulation(size=sim_size,\n",
    "                    mesh_step=[dl, dl, dz],\n",
    "                    structures=geometry(20 * nm),\n",
    "                    sources=[source],\n",
    "                    monitors=[monitor],\n",
    "                    run_time=run_time,\n",
    "                    pml_layers=[0,0,15])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:43.938883Z",
     "iopub.status.busy": "2021-11-12T19:03:43.938547Z",
     "iopub.status.idle": "2021-11-12T19:03:44.247867Z",
     "shell.execute_reply": "2021-11-12T19:03:44.248207Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAv4AAAHwCAYAAADeh2tUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9G0lEQVR4nO3de5hkZX3u/e/tMICRUVAMZwEjvls0iGbE+BojBoxgjBijEQ8IUTcx0ZityVYMHklM8LCTmC1GiBpAVAQUGRVFxCNvxDAoAiMiA0EZGA4i4CgKM/B7/1irpaao7unu6a6q7vX9XFdfvQ5PrfXr6u6n7nrqqVWpKiRJkiQtbvcbdQGSJEmS5p/BX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDhiL4J/kbUlOGXUdkhYX+xZJku41FsF/HCXZP8maIZ3rtUluSPLTJB9OstUUbQ9I8v0kdyT5SpLde/Zt1d7+p+3xXjcPtT44yZlJfp7kh0le1Lf/Re32nyf5dJIHz3UN0kI2rL4ljb9Pcl2S25N8Ncmjp2i/R9un3NH2MQf27Z92PzWN2qbqx96T5Mok69o2L53teSRJGzP4j1iSZwBHAQcAuwMPB94+SdvtgU8BbwYeDKwEPtHT5G3AXu1xnga8PslBs6jpbUneNsnu44C7gB2AFwP/NhEm2u/HA4e1++8A3j/T80uaE88HXgY8haa/+CbwkSnafxz4DvAQ4GjgjCQPhZn1U5syjX7s58AfAg8CDgfem+T/nc25JEkbmzL4J/nfST7Zt+1fk7w3yZOS/Kzn65dJrpnkOHskqSRHJrk+ydokfzPFeU9vR5ZuT/L13lGqJCcmOS7J59oRoW8l+Y2e/f8jyblJfpLkiiR/MsV5/jTJ5e1xrk7yZ+32BwCfB3bu+fl27rvtlkkuTvKX7fqSJP9fkrdMdZ8OcDjwoapaVVW3An8HHDFJ2+cCq6rq9Kr6JU3Qf2yS/9FzrL+rqlur6nLg33uPleRZbc23JfnPJPvMpND2fvlj4M1V9bOqOh9YQRP0oXki8Jmq+npV/Yzmgf25SZa1t68kj+g53olJ/r5d3j/JmiSvT3JT+zfynCTPTPKD9vf5tzOpV+MryW+0v9PHt+s7J7m5/Tvoet/yhCQ3JlnSs+25Sb67yTt2Y3sC51fV1VV1N3AKsPck9T4SeDzw1qr6RVV9EriU5v8dNtFPzeS+YRP9WFW9taq+X1X3VNW3gG8AT2rPc59XS5Jck/bViTSDFqcnOaW97y9N8sgkb2z7lWuT/P4M7kNJWlQ2NeJ/CnBQkm0BkmwBHAqcXFXfrKptqmobYDvgWzQjRlN5Gs2I9O8Db0jfS8k9Pt+2+3Xg28BH+/YfSjPatB2wGnhHW98DgHOBj7W3PRR4f5KBD3bATcCzgAcCfwr8c5LHV9XPgYOB6yd+xqq6vveGVXUX8BLgmCSPohkNW9JTy4vagD3Z18PaQz0a6H1A/y6wQ5KHDKh3o7ZtnVcBj06yHbDTgGNNjMY/Dvgw8Gc0I3rHAysys5frHwlsqKofDDrHgPquonl14JHTPP6OwNbALsBbaJ64vAT4LZpRyzcn2XMG9WpMtX8bbwBOSfJrwH8AJ1XVV+1b6kLglvZnmXAYcHJby3T7llOB32iD71Ka8P6FSep9NHB1Va3r2Tbp/zY9/dQs7ptJ+7H+hknuDzwBWDXJsQb5Q5pXNrajeQXjHJrHul2AY2j6PknqpCmDf1WtBb5O85IxwEHAj6vqor6m/wqso3l5eCpvr6qfV9WlNA/0L5zkvB+uqnVVdSf3jgY9qKfJmVX1X1W1geaBe992+7OAa6rqP6pqQ1V9B/hkT/395/lcVV1Vja8BX6QJmNNSVZcBfw98Gvgb4LB2ZI2q+lhVbTvF14/aw2wD3N5z2InlZQNO2d92ov2ydh/c91gTxzkSOL6qvlVVd1fVScCdwG9P9+dtz/HTSc6/qfqmYz3wjqpaTxNatgfe2/4trAK+Bzx2BvVqjFXVv9OE62/RPGkd1H90sm8BTqJ50kua98k8gyZYz6RvWQucD1wB/KKt9bWTnG9T/7tT9VMzum+mca5eH6B5knDOJMca5BtVdU77OzwdeChwbE+/ssfEYJYkdc105vj/6gGo/b7RHNH2Jez9gRdV1T3ttt6X6R/W0/zanuUfAhu9xN3edkmSY5NcleSnwDXtru17mt3Qs3wH94be3YEn9o5+0Uw/2XHQD5bk4CQXtC9P3wY8s+8803FSe96zq+rKGd4W4Gc0o4ITJpbXTaPtRPt17T6477EmjrM78Nd9981utL+DJJ/t2X4UcFRP289O4/zT2b8pt0w8caIJKgA39uz/Bff+rrU4/DvwGOD/tmH8Vzret5wC/GE7mv4nNGF27QxuD82rZk+g+T/fmuaVjC+3r7D0m+n/dm8/Nel9k+Rhvb+zaZ4LgCTvpvnb+JOqqun8wK3+PuPHA/oV+xFJnTSd4P9pYJ8kj6EZ2fnVS+NJnkIz1/OQqvrVSHDPS9jb9Iw+QfMANOFhwEYvcbdeBBwCHEjz5q49Jk43jVqvBb7WN/q1TVX9eX/DdorLJ4H3ADtU1bbA2T3nme4DzfuBzwLPSPI7Pcd/cV9I6f+aCC2r2HgU+7HAjVV1y4BzbdS2DQW/QTNf9laaEb7+Y028RH4tzWh6733za1X1cYCqetbEduBYmhGyiXbPao/xA2CLJHtNco7++h4ObNXeDpog1Rs6BoYmdUOSbYB/AT4EvC09V4Dqet9SVdfRvBn3uTTTfH414DKDvmVf4BNVtaYdiT+RZvrLoCk4q4CHp30/TmvS/2027qcmvW+q6ke9v7NBx+rtx3q2vZ1mStTv9/7+ad74+2s97ZbQjOhLkqZhk8G/mjdfnUHzMvN/TTzYJtkNOA14aW0853sqb07ya2neUPenbHwlhwnLaKag3ELTwf/DNI8NTQB/ZJLDkixtv56QZg5+vy1pQunNwIYkB7PxnNobgYf0TQPYSJLDaOafHwG8BjipDTNU1Uf7Qkr/10RoORl4eZK925ef3wScOMkpzwQek+SPk2xNM6J3SVV9v+dYb0qyXZo3yv3PnmP9O/DKJE9M4wFJ/qDvgX5K1czF/RTN+xoekOTJNEFqIpR8lGaU8intg/kxwKfq3nnDFwMvakdeDwKeOt1za1F6L7Cyql4BfI5mWod9y71OBl4P/CbN/x0wo77lQuD5SXZIcr+2v1pKM71qI+39fDHw1iRbJ/kjYB+aJzATtUzWT83kvoFN9GNJ3kjzJO3AAQMgPwC2bvuupW0ds76sqCR1zXQv53kSzYNP7zSfA2gu2XhGz0jTpt6A9TWaB53zgPdU1RcHtDmZ5qX662jmdF8wzRppA+bv07y57Hqal+3fyYAHhrbta2gCxq00DzQrevZ/n+YNhVe3L1/3X3njYTSjlS+t5go3H6O5LN0/T7fe9jxfAN4FfAX4Ec3P/tae86xK8uK27c00V9l4R1vzE9ufdcJbad4k90Oa+/rd7fGpqpU0TwTe1952NZNfPWgqfwHcn+bNix8H/rya+fe0319J8wTgJpqg9Rc9t/0rmjfe3UYzFeDTszi/FoEkh9C8Z2hixPx1wOPbv/VO9y09zqSZRnNmVd0x3Vp7vJNmfvzFNP9zrwX+uKpuA0jygSQf6Gl/KLC8rflY4HltnzNlPzWT+6Ztv6l+7B9oXrVZ3fP7/9v2trfT9CkfpPk9/hwYyuetSNJikOlMnWxD7veBHftedp3eSZI9gP8GlrZvuJKkzbbY+5YkVwF/VlVfGnUtkqSFb5Mj/knuRzMad+psQr8kaeaS/DHN+wG+POpaJEmLwxZT7Wznad9I87LujD8BVpI0c0m+SvMm3MMmrmgkSdLmmtZUH0mSJEkL23Tf3CtJkiRpATP4S5IkSR0w5Rz/xeL+S+5f/8+uW466jJG785d38/2bfv7jqvIDb6RZsj9p2J9I0sLTieD/oPv/Ov/1zc+MuoyRe/Pzzuf7N73qh6OuQ1rI7E8a9ieStPCM1VSfJLsl+UqS77UfXPVXA9okyb8mWZ3kkiSPH0Wtksab/YkkSRsbtxH/DcBfV9W3kywDLkpyblV9r6fNwcBe7dcTgX9rv0tSL/sTSZJ6jNWIf1Wtrapvt8vrgMuBXfqaHQKcXI0LgG2T7DTkUiWNOfsTSZI2NlbBv1eSPYDHAd/q27ULcG3P+hru+2BOkiOTrEyycv369fNWp6TxZ38iSdKYBv8k2wCfBP5XVf10NseoqhOqanlVLV+6dOncFihpwbA/kSSpMXbBP8lSmgfpj1bVpwY0uQ7YrWd913abJG3E/kSSpHuNVfBPEuBDwOVV9U+TNFsBvLS9GsdvA7dX1dqhFSlpQbA/kSRpY+N2VZ8nA4cBlya5uN32t8DDAKrqA8DZwDOB1cAdwJ8Ov0xJC4D9iSRJPcYq+FfV+UA20aaAVw2nIkkLlf2JJEkbG6upPpIkSZLmh8FfkiRJ6gCDvyRJktQBBn9JkiSpAwz+kiRJUgcY/CVJkqQOMPhLkiRJHWDwlyRJkjrA4C9JkiR1gMFfkiRJ6gCDvyRJktQBBn9JkiSpAwz+kiRJUgcY/CVJkqQOMPhLkiRJHWDwlyRJkjrA4C9JkiR1gMFfkiRJ6gCDvyRJktQBBn9JkiSpA7YYdQHDkNyPLXfeZ9RljFy2vHjUJUgLnv1Jw/5EkhYeR/wlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQPGKvgn+XCSm5JcNsn+/ZPcnuTi9ustw65R0sJgfyJJ0sa2GHUBfU4E3gecPEWbb1TVs4ZTjqQF7ETsTyRJ+pWxGvGvqq8DPxl1HZIWPvsTSZI2NlbBf5qelOS7ST6f5NGTNUpyZJKVSVauX79+mPVJWjjsTyRJnbHQgv+3gd2r6rHA/wU+PVnDqjqhqpZX1fKlS5cOqz5JC4f9iSSpUxZU8K+qn1bVz9rls4GlSbYfcVmSFiD7E0lS1yyo4J9kxyRpl/ejqf+W0VYlaSGyP5Ekdc1YXdUnyceB/YHtk6wB3gosBaiqDwDPA/48yQbgF8ChVVUjKlfSGLM/kSRpY2MV/KvqhZvY/z6ay/NJ0pTsTyRJ2tiCmuojSZIkaXYM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSdIkklyT5MBR1yHNBYO/JEnSApLkxCR/P8/n2CrJh5P8NMkNSV63ifavbdv9tL3dVj379kjylSR3JPn+fDyRSrJvkovac1yUZN+efUnyziS3tF/vTJK5rmEhMPhLkiSp39uAvYDdgacBr09y0KCGSZ4BHAUc0LZ/OPD2niYfB74DPAQ4GjgjyUNnWlCSrybZf8D2LYGzgFOA7YCTgLPa7QBHAs8BHgvsA/wh8GczPf9iYPCXJEkLSpLfSPKTJI9v13dOcnOS/ZM8P8lFfe1fl+SsSY71tiRnJPlEknVJvp3ksZO03S/JN5PclmRtkvf1hEuSVJJXJrmybXNc78hykpcluTzJrUnOSbL7FD/j6e0I+u1Jvp7k0e32I4EX0wTxnyX5zIDbTuyb+Fqf5MQp79T7Ohz4u6q6taouB/4dOGKKth+qqlVVdSvwdxNtkzwSeDzw1qr6RVV9ErgU+ON2//2SHJXkqnY0/rQkD55hrfsDWwD/UlV3VtW/AgF+r6e+/1NVa6rqOuD/9NS3f5I1vQfrnd7V/n2cnuSU9u/j0iSPTPLGJDcluTbJ78+w3pEx+EuSpAWlqq4C3gCckuTXgP8ATqqqrwIrgD2TPKrnJocBJ09xyEOA04EHAx8DPp1k6YB2dwOvBbYHnkQzwv0XfW2eBTyBZmT5T4BnACQ5BPhb4LnAQ4Fv0IyET+bzNCPuvw58G/ho+7Of0C6/q6q2qao/7L9hVU3s2wZ4FHAz8Im2jve3T0oGfV3SttkO2An4bs9hvws8epJaHz2g7Q5JHtLuu7qq1k1yrL+kGY1/KrAzcCtw3BT3y2Tnv6SqqmfbJT3nGFTfZD/LIH8IfITm1YTvAOfQZOhdgGOA42dY78gY/CVJ0oJTVf8OrAa+RRNSj26330kTcl8C0I6U7wF8dorDXVRVZ1TVeuCfgK2B3x5wzouq6oKq2lBV19AEvqf2NTu2qm6rqh8BXwH2bbe/EvjHqrq8qjYA/wDsO9mof1V9uKrWtT/P24DHJnnQFD/DfSS5P/Bp4L1V9fn2uH9RVdtO8rVPe9Nt2u+39xzudmDZJKfaZkBb2vb9+/qP9Urg6HY0fuJnfV6SLWbwo27qHIPq22YG8/y/UVXntL+302meuB3b/r2cCuyRZNsZ1DsyBn9JkrRQ/TvwGOD/tqFxwknAi9pgdxhwWlXdmeTFPdNfPt/T/tqJhaq6B1hDM/q8kXaKx2fbKTg/pQnv2/c1u6Fn+Q7uDdG7A++dGF0HfkIzHWWXAedZkuTYdvrLT4Fr2l3959qUDwFXVNU7Z3i7n7XfH9iz7YHAugFtJ9r3t6Vt37+v/1i7A2f23C+X07yysgNA7ysSwO8An+3ZdtQk5+8/x6D6ftb3CsFUbuxZ/gXw46q6u2cd7v09jzWDvyRJWnCSbAP8C024fVvvvPCqugC4C3gK8CKaaRpU1UcnpsBU1cE9h9ut57j3A3YFrh9w2n8Dvg/sVVUPpJm6M91R42uBP+sbYb9/Vf3ngLYvopl+dCDwIJpXLOg51yYDaxuKHwm8vG/7B/rm//d+rQJo5+mvpXkz7ITHAqsmOd2qAW1vrKpb2n0PT7Ksb//Esa4FDu67X7Zu5+LTux04H3hWz7Zje86/T98I/j495xhU38S+nwO/1nP/LKEZ0V+UDP6SJGkhei+wsqpeAXwO+EDf/pOB9wHrq+r8TRzrt5I8t51e8r+AO4ELBrRbBvwU+FmS/wH8+Qzq/QDwxp436T4oyfMnabusreEWmlD6D337b6S5cs5ASQ4GXgP8UVX9ondfVb2y58lP/1fvvPeTgTcl2a79Wf8ncOIkpzwZeHmSvdspL2+aaFtVPwAuBt6aZOskf0QTyj/Zc7+8Y2LKU5KHtu+HmImv0rxK8Jo0lyF9dbv9yz31vS7JLkl2Bv6652f5AbB1kj9o39fxJuBXlyJdbAz+kiRpQWmD4UHcG7xfBzw+yYt7mn2EZhrQKdM45FnAC2jeWHoY8Nx2/na/v6EZjV9HM83oE9OtuarOBN4JnNpO37kMOHiS5icDPwSuA77HfZ+EfAjYu53u8ukBt38Bzaj15T2j+f1PjDblrcBVbR1fA95dVV8ASPKw9pgPa3+2LwDvonlPw4/a27y151iHAstp7t9jgedV1c3tvvfSvCH7i0nWtT/rE2dSaFXdRfMG4ZcCtwEvA57TbofmvRifobma0GU0TxSPb297O80btD9Ic3//nGaq16KU6U9vWrh2WrZnrV3336MuY+SO2v9k3vm1wy+qquWjrkVaqOxPGvYnGnftG1tvAh5fVVdO0e5twCOq6iXDqk0aFUf8JUnSYvTnwIVThX6pa2ZyqSRJkqSxl+QamjfCPme0lUjjxeAvSZIWlaraYwZt3zZ/lUjjxak+kiRJUgc44i9JkmZkqzygtttm0V7qXFpw7rpjA7fec92Pq2rKf0yDvyRJmpHttnkoq9deNuoyJLXefvDpvOf8l/1wU+2c6iNJ0iKW5KAkVyRZ3X6aa//+Vya5NMnFSc5Psvco6pQ0/wz+kiQtUkmWAMfRfFDU3sALBwT7j1XVb1bVvjQfwvRPw61S0rAY/CVJWrz2A1ZX1dXtp5ieChzS26Cqftqz+gBg8X+yp9RRzvGXJGnx2gW4tmd9DfDE/kZJXgW8DtgS+L3hlCZp2MZqxD/Jh5PclGTgO4bS+Nd2nuIlSR4/7BolLQz2J9L0VdVxVfUbwBuANw1qk+TIJCuTrFy/fv1wC5Q0J8Yq+AMnAgdNsf9gYK/260jg34ZQk6SF6UTsT6TrgN161ndtt03mVCb5tNuqOqGqllfV8qVLl85dhZKGZqyCf1V9HfjJFE0OAU6uxgXAtkl2Gk51khYS+xMJgAuBvZLsmWRL4FBgRW+DJHv1rP4BcOUQ65M0RGMV/Kdh0FzFXQY19CVJSZtgf6JFr6o2AK8GzgEuB06rqlVJjkny7LbZq5OsSnIxzTz/w0dTraT5tmjf3FtVJwAnAOy0bE+vUCBp1uxPtJBV1dnA2X3b3tKz/FdDL0rSSCy0Ef+ZzlWUpMnYn0iSOmWhBf8VwEvbq3H8NnB7Va0ddVGSFiT7E0lSp4zVVJ8kHwf2B7ZPsgZ4K7AUoKo+QPNS5TOB1cAdwJ+OplJJ487+RJKkjY1V8K+qF25ifwGvGlI5khYw+xNJkja20Kb6SJIkSZoFg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDthi1AVIkiRNZtmyZaMuQZqVdevWjbqE+3DEX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIWsSQHJbkiyeokRw3Y/7ok30tySZLzkuw+ijolzT+DvyRJi1SSJcBxwMHA3sALk+zd1+w7wPKq2gc4A3jXcKuUNCwGf0mSFq/9gNVVdXVV3QWcChzS26CqvlJVd7SrFwC7DrlGSUNi8JckafHaBbi2Z31Nu20yLwc+P2hHkiOTrEyycv369XNYoqRh2WLUBUiSpNFL8hJgOfDUQfur6gTgBIAdl+1RQyxN0hwx+EuStHhdB+zWs75ru20jSQ4EjgaeWlV3Dqk2SUPmVB9JkhavC4G9kuyZZEvgUGBFb4MkjwOOB55dVTeNoEZJQ2LwlyRpkaqqDcCrgXOAy4HTqmpVkmOSPLtt9m5gG+D0JBcnWTHJ4SQtcE71kSRpEauqs4Gz+7a9pWf5wKEXJWkkHPGXJEmSOsDgL0mSJHWAwV+SJEnqAIO/JEmS1AEGf0mSJKkDDP6SJElSBxj8JUmSpA4w+EuSJEkdYPCXJEmSOsDgL0mSJHWAwV+SJEnqgLEL/kkOSnJFktVJjhqw/4gkNye5uP16xSjqlDTe7EskSdrYFqMuoFeSJcBxwNOBNcCFSVZU1ff6mn6iql499AIlLQj2JZIk3de4jfjvB6yuqqur6i7gVOCQEdckaeGxL5Ekqc+4Bf9dgGt71te02/r9cZJLkpyRZLdBB0pyZJKVSVauX79+PmqVNL7mrC8B+xNJ0uIwbsF/Oj4D7FFV+wDnAicNalRVJ1TV8qpavnTp0qEWKGlBmFZfAvYnkqTFYdyC/3VA76jbru22X6mqW6rqznb1g8BvDak2SQuHfYkkSX3GLfhfCOyVZM8kWwKHAit6GyTZqWf12cDlQ6xP0sJgXyJJUp+xuqpPVW1I8mrgHGAJ8OGqWpXkGGBlVa0AXpPk2cAG4CfAESMrWNJYsi+RJOm+xir4A1TV2cDZfdve0rP8RuCNw65L0sJiXyJJ0sbGbaqPJEmSpHlg8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSVrEkhyU5Iokq5McNWD/7yb5dpINSZ43iholDYfBX5KkRSrJEuA44GBgb+CFSfbua/Yjmk+u/thwq5M0bGP3yb2SJGnO7AesrqqrAZKcChwCfG+iQVVd0+67ZxQFShoeR/wlSVq8dgGu7Vlf026T1EEGf0mStElJjkyyMsnK9evXj7ocSbNg8JckafG6DtitZ33XdtuMVdUJVbW8qpYvXbp0ToqTNFwGf0mSFq8Lgb2S7JlkS+BQYMWIa5I0Ip14c+/ddTf/dfV/jbqMkVv3y5+OugRpwbM/adifLAxVtSHJq4FzgCXAh6tqVZJjgJVVtSLJE4Azge2AP0zy9qp69AjLljRPOhH8JUnqqqo6Gzi7b9tbepYvpJkCJGmRc6qPJEmS1AGdGPFfkiXs9/D9Rl3GyH1q6++PugRpwbM/adifSNLC44i/JEmS1AEGf0mSJKkDOjHVZ1OSzPkxq2pszytp/szH//Wo2J9I0uLiiL8kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDthitjdMsh2wM/AL4JqqumfOqpLUKfYn6rIkTwJeAjwF2Inm/+Ay4HPAKVV1+wjLG7l169aNugRp0ZhR8E/yIOBVwAuBLYGbga2BHZJcALy/qr4y51VKWnTsTyRI8nngeuAs4B3ATTT/B48EngacleSfqmrF6KqUtFjMdMT/DOBk4ClVdVvvjiS/BRyW5OFV9aE5qk/S4mV/IsFhVfXjvm0/A77dfv2fJNsPvyxJi9GMgn9VPX2KfRcBF212RZI6wf5Egv7Qn+SB9Dw2V9VPBjwxkKRZ2Zw5/vsAe7BxB/WpOahJUsfYn6jrkvwZ8Hbgl0C1mwt4+MiKkrTozCr4J/kwsA+wCph4E14BPlBLmhH7EwmAvwEe4+i+pPk02xH/366qvee0EkldZX8iwVXAHaMuQtLiNtvg/80ke1fV9+a0GkldZH8iwRuB/0zyLeDOiY1V9ZrRlSRpsZlt8D+Z5sH6BpoOKkBV1T5zVpmkrrA/keB44MvApdw75U2S5tRsg/+HgMOwg5K0+exPJFhaVa8bdRGSFrfZBv+b/TARSXPE/kSCzyc5EvgMG0/1+cnoSpK02Nxvlrf7TpKPJXlhkudOfM1FQUkOSnJFktVJjhqwf6skn2j3fyvJHnNxXkkjY38iNZ9g/UbgP2k+w+IiYOVcHNj/A0kTZjvif3+aEYnf79m22ZffS7IEOA54OrAGuDDJir43/b0cuLWqHpHkUOCdwAs257ySRsr+RJ1XVXvOx3H9P5DUa1bBv6r+dK4Lae0HrK6qqwGSnAocAvR2UIcAb2uXzwDelyRVVUhacOxPJEjy0kHbq+rkzTy0/weSfmW2H+D1H9z7yYK/UlUv28x6dgGu7VlfAzxxsjZVtSHJ7cBDgP6PPT8SOBLgwVvtzM9+9rNJT7pu3brNLPu+pjrfqM579913z/n5pM1lfzK+7E+G6gk9y1sDBwDfprnq1eaYt/8DSQvPbKf6fLZneWvgj4DrN7+cuVNVJwAnAOy4bA9HLaTxZX+izquqv+xdT7ItcOpoqhnM/wNp4ZvtVJ9P9q4n+Thw/hzUcx2wW8/6ru22QW3WJNkCeBBwyxycW9II2J9IA/0cmIt5//4fSPqV2V7Vp99ewK/PwXEuBPZKsmeSLYFDgf7L/K0ADm+Xnwd82XmI0qJif6LOSfKZJCvar88CVwBnzsGh/T+Q9CuzneO/jo3n5N4AvGFzi2nnFr4aOAdYAny4qlYlOQZY2V7r+0PAR5KsBn5C04lJWqDsTyQA3tOzvAH4YVWt2dyD+n8gqddsp/osm+tCeo59NnB237a39Cz/Enj+fJ1f0nDZn6jLJq6eU1Vf21Sb2Z7D/wNJE2Y01WdTH+qRxq6bVZGkTrA/kQD4SpK/TPKw3o1Jtkzye0lO4t5pOJK0WWY64v/uJPcDzqL5VMGbaa7C8QjgaTSXH3srzeXCJGkq9icSHAS8DPh4kj2B22j+D5YAXwT+paq+M7ryJC0mMwr+VfX8JHsDL6bpqHYC7gAup3kZ8R3tS4aSNCX7E+lX02zeD7w/yVJge+AXVXXbSAuTtCjNeI5/+zHfR89DLZI6xv5EuldVrQfWjroOSYvXXF3OU5IkSdIYM/hLkiRJHWDwlyRpxNor+2w36jokLW6zCv5JzkvyzL5tJ8xNSZK6xP5EAmAH4MIkpyU5KElGXZCkxWe2I/57Am9I8taebcvnoB5J3WN/os6rqjcBe9F8iu4RwJVJ/iHJb4y0MEmLyqw+uZfmOsMHAP+a5DPAS+asIkldcxv2JxJVVUluAG4ANgDbAWckObeqXj/a6kZn2bJ5+3BvaV6tW7du1CXcx2yDf6pqA/AXSY4AzqfpoCRppuxP1HlJ/gp4KfBj4IPA/66q9e2H3F0JdDb4S5o7sw3+H5hYqKoTk1wKvGpuSpLUMfYnEjwYeG5V/bB3Y1Xdk+RZI6pJ0iIzq+BfVcf3rV9E88mbkjQj9icSVNVbp9h3+TBrkbR4eTlPSZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR2wxagLGIq6m7vXfXXUVYze3beNugJp4bM/adifSNKC44i/JEmS1AEGf0mSFqEkD05ybpIr2+/bTdLuC0luS/LZYdcoabi6MdUnS1iybP9RVzF6S24edQXSwmd/0rA/WQiOAs6rqmOTHNWuv2FAu3cDvwb82TCLkzR8jvhLkrQ4HQKc1C6fBDxnUKOqOg9YN6SaJI2QwV+SpMVph6pa2y7fAOwwymIkjV43pvpIkrQIJfkSsOOAXUf3rlRVJanNPNeRwJEAD95q5805lKQRMfhLkrRAVdWBk+1LcmOSnapqbZKdgJs281wnACcA7Lhsj816EiFpNJzqI0nS4rQCOLxdPhw4a4S1SBoDBn9JkhanY4GnJ7kSOLBdJ8nyJB+caJTkG8DpwAFJ1iR5xkiqlTTvnOojSdIiVFW3AAcM2L4SeEXP+lOGWZek0XHEX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdcDYBP8kD05ybpIr2+/bTdLu7iQXt18rhl2npPFnfyJJ0n2NTfAHjgLOq6q9gPPa9UF+UVX7tl/PHl55khYQ+xNJkvqMU/A/BDipXT4JeM7oSpG0wNmfSJLUZ5yC/w5VtbZdvgHYYZJ2WydZmeSCJM+Z7GBJjmzbrVy/fv1c1yppvNmfSJLUZ4thnizJl4AdB+w6unelqipJTXKY3avquiQPB76c5NKquqq/UVWdAJwAsOOyPSY7lqQFyv5EkqSZGWrwr6oDJ9uX5MYkO1XV2iQ7ATdNcozr2u9XJ/kq8DjgPg/UkhY3+xNJkmZmnKb6rAAOb5cPB87qb5BkuyRbtcvbA08Gvje0CiUtFPYnkiT1Gafgfyzw9CRXAge26yRZnuSDbZtHASuTfBf4CnBsVflALamf/YkkSX2GOtVnKlV1C3DAgO0rgVe0y/8J/OaQS5O0wNifSJJ0X+M04i9JkiRpnhj8JUmSpA4w+EuSJEkdYPCXJEmSOsDgL0mSJHWAwV+SJEnqAIO/JEmS1AEGf0mSJKkDDP6SJElSBxj8JUmSpA4w+EuSJEkdYPCXJGkRSvLgJOcmubL9vt2ANvsm+WaSVUkuSfKCUdQqaTgM/pIkLU5HAedV1V7Aee16vzuAl1bVo4GDgH9Jsu3wSpQ0TAZ/SZIWp0OAk9rlk4Dn9Deoqh9U1ZXt8vXATcBDh1WgpOEy+EuStDjtUFVr2+UbgB2mapxkP2BL4KpJ9h+ZZGWSlevXr5/bSiUNxRajLkCSJM1Oki8BOw7YdXTvSlVVkpriODsBHwEOr6p7BrWpqhOAEwB2XLbHpMeSNL4M/pIkLVBVdeBk+5LcmGSnqlrbBvubJmn3QOBzwNFVdcE8lSppDDjVR5KkxWkFcHi7fDhwVn+DJFsCZwInV9UZQ6xN0ggY/CVJWpyOBZ6e5ErgwHadJMuTfLBt8yfA7wJHJLm4/dp3JNVKmndO9ZEkaRGqqluAAwZsXwm8ol0+BThlyKVJGhFH/CVJkqQOMPhLkiRJHWDwlyRJkjrA4C9JkiR1gMFfkiRJ6gCDvyRJktQBBn9JkiSpAwz+kiRJUgcY/CVJkqQO8JN7JUnS2Fq3bt2oS5AWDUf8JUmSpA4w+EuSJEkdYPCXJEmSOsDgL0mSJHWAwV+SJEnqAIO/JEmS1AEGf0mSJKkDDP6SJElSBxj8JUmSpA4w+EuSJEkdYPCXJEmSOsDgL0mSJHWAwV+SJEnqAIO/JEmS1AFjE/yTPD/JqiT3JFk+RbuDklyRZHWSo4ZZo6SFwf5EkqT7GpvgD1wGPBf4+mQNkiwBjgMOBvYGXphk7+GUJ2kBsT+RJKnPFqMuYEJVXQ6QZKpm+wGrq+rqtu2pwCHA9+a9QEkLhv2JJEn3NU4j/tOxC3Btz/qadtt9JDkyycokK9evXz+U4iQtKPYnWtSSPDjJuUmubL9vN6DN7km+neTidnrcK0dRq6ThGGrwT/KlJJcN+Dpkrs9VVSdU1fKqWr506dK5PrykEbM/kTbpKOC8qtoLOK9d77cWeFJV7Qs8ETgqyc7DK1HSMA11qk9VHbiZh7gO2K1nfdd2m6SOsT+RNukQYP92+STgq8AbehtU1V09q1ux8GYCSJqBhfYPfiGwV5I9k2wJHAqsGHFNkhYm+xMtdjtU1dp2+QZgh0GNkuyW5BKaqW/vrKrrh1WgpOEam+Cf5I+SrAGeBHwuyTnt9p2TnA1QVRuAVwPnAJcDp1XVqlHVLGk82Z+oK6Y75a2qCqhBx6iqa6tqH+ARwOFJJnuC4HtdpAVunK7qcyZw5oDt1wPP7Fk/Gzh7iKVJWmDsT9QVU015S3Jjkp2qam2SnYCbNnGs65NcBjwFOGPA/hOAEwB2XLbHwCcRksbb2Iz4S5KkObUCOLxdPhw4q79Bkl2T3L9d3g74HeCKoVUoaagM/pIkLU7HAk9PciVwYLtOkuVJPti2eRTwrSTfBb4GvKeqLh1JtZLm3dhM9ZEkSXOnqm4BDhiwfSXwinb5XGCfmR77rjs28PaDT9/sGiXNjR9eMuVMvl8x+EuSpBm59Z7rfvye81/2w77N2wM/HkU9Uxi3msatHhi/msatHhi/miarZ/dN3dDgL0mSZqSqHtq/LcnKqlo+inomM241jVs9MH41jVs9MH41bU49zvGXJEmSOsDgL0mSJHWAwV+SJM2FE0ZdwADjVtO41QPjV9O41QPjV9Os6zH4S5KkzdZ+wNdYGbeaxq0eGL+axq0eGL+aNqceg78kSZLUAQZ/SZIkqQMM/pIkacaSPDjJuUmubL9vN6DN7km+neTiJKuSvHIMato3yTfbei5J8oJR1tO2+0KS25J8dh5rOSjJFUlWJzlqwP6tknyi3f+tJHvMVy3TrOd327+dDUmeN5+1TLOe1yX5Xvs3c16STV4zfwg1vTLJpe3/1/lJ9t7UMQ3+kiRpNo4CzquqvYDz2vV+a4EnVdW+wBOBo5LsPOKa7gBeWlWPBg4C/iXJtiOsB+DdwGHzVANJlgDHAQcDewMvHBASXw7cWlWPAP4ZeOeI6/kRcATwsfmqY4b1fAdYXlX7AGcA7xqDmj5WVb/Z/n+9C/inTR3X4C9JkmbjEOCkdvkk4Dn9Darqrqq6s13divnPHdOp6QdVdWW7fD1wE3CfDyQbVj1tHecB6+apBoD9gNVVdXVV3QWc2tbWq7fWM4ADkmRU9VTVNVV1CXDPPNUw03q+UlV3tKsXALuOQU0/7Vl9AFCbOqjBX5IkzcYOVbW2Xb4B2GFQoyS7JbkEuBZ4Zxu2R1pTT237AVsCV41DPfNoF5r7f8KadtvANlW1AbgdeMgI6xmmmdbzcuDz81rRNGtK8qokV9GM+L9mUwfdYs7KkyRJi0qSLwE7Dth1dO9KVVWSgaONVXUtsE87xefTSc6oqhtHWVN7nJ2AjwCHV9WsR5Xnqh4tDEleAiwHnjrqWgCq6jjguCQvAt4EHD5Ve4O/JEkaqKoOnGxfkhuT7FRVa9sQfdMmjnV9ksuAp9BMJRlZTUkeCHwOOLqqLphtLXNVzxBcB+zWs75ru21QmzVJtgAeBNwywnqGaVr1JDmQ5gndU3umsI20ph6nAv+2qYM61UeSJM3GCu4dXTwcOKu/QZJdk9y/Xd4O+B3gihHXtCVwJnByVc36Cchc1TMkFwJ7Jdmz/fkPbWvr1Vvr84AvV9V8vUIxnXqGaZP1JHkccDzw7KoaxhO46dS0V8/qHwBXbuqgBn9JkjQbxwJPT3IlcGC7TpLlST7YtnkU8K0k3wW+Brynqi4dcU1/AvwucER7GcSLk+w7wnpI8g3gdJo31K5J8oy5LKKds/9q4BzgcuC0qlqV5Jgkz26bfQh4SJLVwOuY/ApEQ6knyROSrAGeDxyfZNUo66G58tI2wOnt38y8PlGZZk2vTnNZ2otpfmdTTvMByPw9mRsfOy7bo1avvWzUZYzc2w8+nfec/7KLqmr5qGuRFir7k4b9iSQtPI74S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZKkRSzJ/ZN8LcmSOTjWQ5N8YS7q0vAZ/CVJkha3lwGfqqq7N/dAVXUzsDbJkze/LA2bwV+SJGkBaj/k6pIkWyd5QPthTo8Z0PTFtJ8anGT/JJ/tOcb7khzRLl+T5B/bD6hameTxSc5JclWSV/Yc79PtMbXAbDHqAiRJkjRzVXVh+wmyfw/cHzilqjb6hMEkWwIPr6prpnnYH1XVvkn+GTgReDKwNXAZ8IG2zcr2nFpgDP6SJEkL1zHAhcAvgdcM2L89cNsMjrei/X4psE1VrQPWJbkzybZVdRtwE7DzrCvWyDjVR5IkaeF6CLANsIxmZL7fL/q2b2Dj/Nd/mzvb7/f0LE+sTwwYb90eVwuMwV+SJGnhOh54M/BR4J39O6vqVmBJkomA/0Ng7yRbJdkWOGAW53wkzdQfLTBO9ZEkSVqAkrwUWF9VH2sv1fmfSX6vqr7c1/SLwO8AX6qqa5OcRhPc/xv4zixO/TTgc5tTu0bD4C9JkrQAVdXJwMnt8t3AEydpehzwWuBLbdvXA68fcLw9epZPpHlz7332Ac8GDtmM0jUiTvWRJElaxKrq28BX5uoDvIB/aqcQaYFxxF+SJGmRq6oPz9Fxbqa5jr8WIEf8JUmSpA4w+EuSJEkdYPCXJEmSOsDgL0mSJHWAwV+SJEnqAIO/JEmS1AEGf0mSJKkDDP6SJElSB4xN8E/y/CSrktyTZPkU7a5JcmmSi5OsHGaNksaffYkkSYON0yf3XgY8Fzh+Gm2fVlU/nud6JC1M9iWSJA0wNsG/qi4HSDLqUiQtYPYlkiQNNjZTfWaggC8muSjJkZM1SnJkkpVJVq5fv36I5UlaIKbVl4D9iSRpcRjqiH+SLwE7Dth1dFWdNc3D/E5VXZfk14Fzk3y/qr7e36iqTgBOANhx2R4166IljZ1h9iVgfyJJWhyGGvyr6sA5OMZ17febkpwJ7AcMfLCWtDjZl0iSNHMLaqpPkgckWTaxDPw+zRv5JGna7EskSV00NsE/yR8lWQM8CfhcknPa7TsnObtttgNwfpLvAv8FfK6qvjCaiiWNI/sSSZIGG6er+pwJnDlg+/XAM9vlq4HHDrk0SQuIfYkkSYONzYi/JEmSpPlj8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQO2GHUBQ3HPPdx9w2WjrmL07vr5qCuQFj77k4b9iSQtOI74S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIHGPwlSZKkDjD4S5IkSR1g8JckSZI6wOAvSZIkdYDBX5IkSeoAg78kSZLUAWMT/JO8O8n3k1yS5Mwk207S7qAkVyRZneSoIZcpaQGwP5Ek6b7GJvgD5wKPqap9gB8Ab+xvkGQJcBxwMLA38MIkew+1SkkLgf2JJEl9xib4V9UXq2pDu3oBsOuAZvsBq6vq6qq6CzgVOGRYNUpaGOxPJEm6r7EJ/n1eBnx+wPZdgGt71te02+4jyZFJViZZuX7DhkFNJHWD/YkkScAWwzxZki8BOw7YdXRVndW2ORrYAHx0c85VVScAJwDs+ICH1eYcS9L4sT+RJGlmhhr8q+rAqfYnOQJ4FnBAVQ16cL0O2K1nfdd2m6SOsT+RJGlmxmaqT5KDgNcDz66qOyZpdiGwV5I9k2wJHAqsGFaNkhYG+xNJku5rbII/8D5gGXBukouTfAAgyc5JzgZo36z3auAc4HLgtKpaNaqCJY0t+xNJkvoMdarPVKrqEZNsvx54Zs/62cDZMzr23Xfxiy++efMKXATuvnWfUZcgDYX9yfyzP5GkhSeDp74uLkluBn44ye7tgR8PoYxhnWdT59q9qh46pDqkRcf+ZCP2J5K0gHQi+E8lycqqWr5YzjPsc0m6l/2JJGmcjdMcf0mSJEnzxOAvSZIkdYDBv/1QnkV0nmGfS9K97E8kSWOr83P8JUmSpC5wxF+SJEnqgE4E/yQPTnJukivb79sNaLNvkm8mWZXkkiQv6Nl3YpL/bj8I6OIk+/bd9qAkVyRZneSoAcfeKskn2v3fSrJHz743ttuvSPKMTfwcmzrP65J8r63/vCS79+y7u6d+P51UmoX57kvaNvYnkqR50YmpPkneBfykqo5tH+C2q6o39LV5JFBVdWWSnYGLgEdV1W1JTgQ+W1VnDDj2EuAHwNOBNcCFwAur6ns9bf4C2KeqXpnkUOCPquoFSfYGPg7sB+wMfAl4ZFXdPcvzPA34VlXdkeTPgf2r6gXtvp9V1Tazuf8kNeazL2lva38iSZo3nRjxBw4BTmqXTwKe09+gqn5QVVe2y9cDNwHT+WCa/YDVVXV1Vd0FnNqeb7LznwEckCTt9lOr6s6q+m9gdXu8WZ2nqr5SVXe0qxcAu06jfknTN599CdifSJLmUVeC/w5VtbZdvgHYYarGSfYDtgSu6tn8jvYl739OslXP9l2Aa3vW17TbGNSmqjYAtwMPmeZtZ3KeXi8HPt+zvnWSlUkuSPKcKW4naXLz2ZeA/YkkaR5tMeoC5kqSLwE7Dth1dO9KVVWSSec3JdkJ+AhweFXd025+I82D/JY0l7Z7A3DMXNQ9H5K8BFgOPLVn8+5VdV2ShwNfTnJpVV01+AhSd9mXbMz+RJIWj0UT/KvqwMn2JbkxyU5VtbZ9ML5pknYPBD4HHF1VF/Qce2KE784k/wH8Tc/NrgN261nftd3GgDZrkmwBPAi4ZZq3ncl5SHIgTUB5alXd2fMzXNd+vzrJV4HHsfEopCRG2peA/YkkaR51ZarPCuDwdvlw4Kz+Bkm2BM4ETu5/4137AE87j/Y5wGU9uy8E9kqyZ3uMQ9vzTXb+5wFfruZd1SuAQ9urdOwJ7AX81yQ/wybPk+RxwPHAs6vqpp7t201MKUiyPfBk4HtImqn57EvA/kSSNI8WzYj/JhwLnJbk5cAPgT8BSLIceGVVvaLd9rvAQ5Ic0d7uiKq6GPhokocCAS4GXjlx4KrakOTVwDnAEuDDVbUqyTHAyqpaAXwI+EiS1cBPaB5kadudRvOguQF41aArcMzgPO8GtgFOb3IFP6qqZwOPAo5Pcg/Nk71je6/eIWna5q0vAfsTSdL86sTlPCVJkqSu68pUH0mSJKnTDP6SJElSBxj8JUmSpA4w+EuSJEkdYPCXJEmSOsDgP8eS3D/J15IsmYNjPTTJF+aiLkmSJHWbwX/uvQz41GTXz56JqroZWJvkyZtflqSFxoEESdJcMvhPU5JjkvyvnvV3JPmrAU1fTPtpnkn2T/LZntu8b+IDfZJck+Qfk1ycZGWSxyc5J8lVSXo/1OfT7TEldY8DCZKkOWPwn74PAy8FSHI/mk/LPKW3QZItgYdX1TXTPOaPqmpf4BvAicDzgN8G3t7TZiXwlM2oW9KYSfKEJJck2TrJA5KsSvKYAU0dSJAkzZktRl3AQlFV1yS5JcnjgB2A71TVLX3Ntgdum8FhV7TfLwW2qap1wLokdybZtqpuA24Cdt686iWNk6q6MMkK4O+B+wOnVNVlvW1mO5CQ5J9pBhKeDGwNXAZ8oG2zsj2nJKmDDP4z80HgCGBHmlcA+v2C5oF2wgY2flVl642bc2f7/Z6e5Yn1id/N1u1xJS0uxwAXAr8EXjNgvwMJkqQ55VSfmTkTOAh4AnBO/86quhVYkmQi4P8Q2DvJVkm2BQ6YxTkfSTNiJ2lxeQiwDbCM+w4KgAMJkqQ5ZvCfgaq6C/gKcNoUb7b7IvA7bftrgdNogvtpwHdmcdqnAZ+bxe0kjbfjgTcDHwXe2b/TgQRJ0lxzqs8MtG/q/W3g+VM0Ow54LfAlgKp6PfD6/kZVtUfP8ok0c3Lvsw94NnDIrIuWNHaSvBRYX1Ufay/V+Z9Jfq+qvtzXdGIg4UtVdW2SiYGE/8aBBEnSDKWqRl3DgpBkb+CzwJlV9debaPsy4KTNvQRfkocCT66qT2/OcSQtTEkeD7y2qg6bo+N9HTikfTVBktQxBn9JGmMOJEiS5orBX5IkSeoA39wrSZIkdYDBX5IkSeoAg78kSZLUAQZ/SZIkqQMM/pIkSVIH/P8h5ZKEa7LjKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Structure visualization in various planes\n",
    "\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 8))\n",
    "sim_d0.viz_eps_2D(normal='x', position=0, ax=ax1)\n",
    "sim_d0.viz_eps_2D(normal='y', position=g, ax=ax2)\n",
    "sim_d0.viz_eps_2D(normal='z', position=0, ax=ax3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:44.251199Z",
     "iopub.status.busy": "2021-11-12T19:03:44.250851Z",
     "iopub.status.idle": "2021-11-12T19:03:44.260568Z",
     "shell.execute_reply": "2021-11-12T19:03:44.260238Z"
    }
   },
   "outputs": [],
   "source": [
    "# define a function that runs the simulations on our server and computes the flux\n",
    "\n",
    "def get_flux(sim, task_name='Si_resonator'):\n",
    "    project = web.new_project(sim.export(), task_name=task_name)\n",
    "    web.monitor_project(project['taskId'])\n",
    "    print('Downloading results')\n",
    "    web.download_results(project['taskId'], target_folder='output')\n",
    "    sim.load_results('output/monitor_data.hdf5')\n",
    "    with open('output/tidy3d.log', 'r') as f:\n",
    "        print(f.read())\n",
    "    return np.squeeze(sim.data(monitor)['flux'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:03:44.262776Z",
     "iopub.status.busy": "2021-11-12T19:03:44.262447Z",
     "iopub.status.idle": "2021-11-12T19:07:30.726778Z",
     "shell.execute_reply": "2021-11-12T19:07:30.726367Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uploading the json file...\n",
      "Project 'normalization-16b6e1d68c55d7e8' status: success...                                \n",
      "\n",
      "Downloading results\n",
      "Applying source normalization to all frequency monitors using source index 0.\n",
      "Simulation domain Nx, Ny, Nz: [32, 32, 206]\n",
      "Applied symmetries: [0, 0, 0]\n",
      "Number of computational grid points: 2.1094e+05.\n",
      "Using subpixel averaging: True\n",
      "Number of time steps: 19987\n",
      "Automatic shutoff factor: 1.00e-05\n",
      "Time step (s): 3.5023e-17\n",
      "\n",
      "Compute source modes time (s):     0.0517\n",
      "Compute monitor modes time (s):    0.0546\n",
      "\n",
      "Rest of setup time (s):            0.0932\n",
      "\n",
      "Starting solver...\n",
      "- Time step    318 / time 1.11e-14s (  1 % done), field decay: 1.00e+00\n",
      "- Time step    799 / time 2.80e-14s (  4 % done), field decay: 1.18e-05\n",
      "- Time step   1598 / time 5.60e-14s (  8 % done), field decay: 4.34e-13\n",
      "Field decay smaller than shutoff factor, exiting solver.\n",
      "\n",
      "Solver time (s):                   1.6589\n",
      "Post-processing time (s):          0.0047\n",
      "\n",
      "Uploading the json file...\n",
      "Project 'Si-resonator-delta-0-16b6e1e153d4b898' status: success...                                \n",
      "\n",
      "Downloading results\n",
      "Applying source normalization to all frequency monitors using source index 0.\n",
      "Simulation domain Nx, Ny, Nz: [32, 32, 206]\n",
      "Applied symmetries: [0, 0, 0]\n",
      "Number of computational grid points: 2.1094e+05.\n",
      "Using subpixel averaging: True\n",
      "Number of time steps: 199868\n",
      "Automatic shutoff factor: 1.00e-05\n",
      "Time step (s): 3.5023e-17\n",
      "\n",
      "Compute source modes time (s):     0.0600\n",
      "Compute monitor modes time (s):    0.0772\n",
      "\n",
      "Rest of setup time (s):            0.3107\n",
      "\n",
      "Starting solver...\n",
      "- Time step    318 / time 1.11e-14s (  0 % done), field decay: 1.00e+00\n",
      "- Time step   7994 / time 2.80e-13s (  4 % done), field decay: 6.99e-03\n",
      "- Time step  15989 / time 5.60e-13s (  8 % done), field decay: 5.25e-04\n",
      "- Time step  23984 / time 8.40e-13s ( 12 % done), field decay: 5.74e-04\n",
      "- Time step  31978 / time 1.12e-12s ( 16 % done), field decay: 3.27e-05\n",
      "- Time step  39973 / time 1.40e-12s ( 20 % done), field decay: 1.33e-04\n",
      "- Time step  47968 / time 1.68e-12s ( 24 % done), field decay: 1.53e-05\n",
      "- Time step  55963 / time 1.96e-12s ( 28 % done), field decay: 2.24e-05\n",
      "- Time step  63957 / time 2.24e-12s ( 32 % done), field decay: 9.00e-06\n",
      "Field decay smaller than shutoff factor, exiting solver.\n",
      "\n",
      "Solver time (s):                   16.4182\n",
      "Post-processing time (s):          0.0058\n",
      "\n",
      "Uploading the json file...\n",
      "Project 'Si-resonator-delta-20-16b6e1f5210d2e18' status: success...                                \n",
      "\n",
      "Downloading results\n",
      "Applying source normalization to all frequency monitors using source index 0.\n",
      "Simulation domain Nx, Ny, Nz: [32, 32, 206]\n",
      "Applied symmetries: [0, 0, 0]\n",
      "Number of computational grid points: 2.1094e+05.\n",
      "Using subpixel averaging: True\n",
      "Number of time steps: 199868\n",
      "Automatic shutoff factor: 1.00e-05\n",
      "Time step (s): 3.5023e-17\n",
      "\n",
      "Compute source modes time (s):     0.0577\n",
      "Compute monitor modes time (s):    0.0684\n",
      "\n",
      "Rest of setup time (s):            0.3133\n",
      "\n",
      "Starting solver...\n",
      "- Time step    318 / time 1.11e-14s (  0 % done), field decay: 1.00e+00\n",
      "- Time step   7994 / time 2.80e-13s (  4 % done), field decay: 7.52e-03\n",
      "- Time step  15989 / time 5.60e-13s (  8 % done), field decay: 1.05e-02\n",
      "- Time step  23984 / time 8.40e-13s ( 12 % done), field decay: 6.16e-03\n",
      "- Time step  31978 / time 1.12e-12s ( 16 % done), field decay: 5.09e-03\n",
      "- Time step  39973 / time 1.40e-12s ( 20 % done), field decay: 6.53e-03\n",
      "- Time step  47968 / time 1.68e-12s ( 24 % done), field decay: 2.04e-03\n",
      "- Time step  55963 / time 1.96e-12s ( 28 % done), field decay: 1.41e-03\n",
      "- Time step  63957 / time 2.24e-12s ( 32 % done), field decay: 4.56e-03\n",
      "- Time step  71952 / time 2.52e-12s ( 36 % done), field decay: 3.37e-03\n",
      "- Time step  79947 / time 2.80e-12s ( 40 % done), field decay: 1.25e-03\n",
      "- Time step  87941 / time 3.08e-12s ( 44 % done), field decay: 2.30e-03\n",
      "- Time step  95936 / time 3.36e-12s ( 48 % done), field decay: 3.46e-03\n",
      "- Time step 103931 / time 3.64e-12s ( 52 % done), field decay: 2.51e-03\n",
      "- Time step 111926 / time 3.92e-12s ( 56 % done), field decay: 9.81e-04\n",
      "- Time step 119920 / time 4.20e-12s ( 60 % done), field decay: 7.87e-04\n",
      "- Time step 127915 / time 4.48e-12s ( 64 % done), field decay: 1.58e-03\n",
      "- Time step 135910 / time 4.76e-12s ( 68 % done), field decay: 1.44e-03\n",
      "- Time step 143904 / time 5.04e-12s ( 72 % done), field decay: 5.84e-04\n",
      "- Time step 151899 / time 5.32e-12s ( 76 % done), field decay: 8.73e-04\n",
      "- Time step 159894 / time 5.60e-12s ( 80 % done), field decay: 2.06e-03\n",
      "- Time step 167889 / time 5.88e-12s ( 84 % done), field decay: 2.05e-03\n",
      "- Time step 175883 / time 6.16e-12s ( 88 % done), field decay: 1.12e-03\n",
      "- Time step 183878 / time 6.44e-12s ( 92 % done), field decay: 7.67e-04\n",
      "- Time step 191873 / time 6.72e-12s ( 96 % done), field decay: 1.17e-03\n",
      "- Time step 199867 / time 7.00e-12s (100 % done), field decay: 1.03e-03\n",
      "\n",
      "Solver time (s):                   49.5821\n",
      "Post-processing time (s):          0.0057\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# run all simulations, take about 2-3 minutes each with some download time\n",
    "flux_empty = get_flux(sim_empty, task_name='normalization')\n",
    "flux_delta0 = get_flux(sim_d0, task_name='Si-resonator-delta-0')\n",
    "flux_delta20 = get_flux(sim_d20, task_name='Si-resonator-delta-20')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get Results and Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:07:30.731416Z",
     "iopub.status.busy": "2021-11-12T19:07:30.731092Z",
     "iopub.status.idle": "2021-11-12T19:07:31.395243Z",
     "shell.execute_reply": "2021-11-12T19:07:31.395510Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEaCAYAAAAcz1CnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbpklEQVR4nO3deZgddZ3v8ffHhE0Ia4JgEiCBIEaRrWUZM4jgIMQLwavDIijbNfeOMnBHYB4QBS6OI8i94x00AlFRUGQVMQxRcFgEHQNZWAMGQgAJsoSdTCCQ8J0/6tdYOZzTXZ0+dZauz+t5+uk6Vb+q+p7q5HxObb9SRGBmZtX1rnYXYGZm7eUgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOzinMQ2JAlKSRtk4YvkPS1AvMslTS+/OqaT9IoSX+UtE6L1/sVST9Iw++R9KCktVpZgw2OfB+BDVWSApgQEQvbXUsrSPp/wJKIOLuEZU8B/g8wHngDuBc4NiIerdP2e8CDEfGdZtdh5fAegXUEScPbXUMrNfv9pm/gRwI/beZy07K3AS4BTgQ2AMYB04CVDWa5FPifza7DyuMgsEGR9JikkyTdK+llSVdIWjs3/QuSFkp6QdIMSe/NTQtJX5L0MPCwpL0kLZb0j5KelfSUpIMkTZb0UFrGV3Lz7yrpD5JeSm2/K2nNBnX+WNI/peHr0iGg3p+3JB2Vq2mb3DzTJF0v6VVJd0jaOrfMfSUtSO/7e5J+K+l/NFj/mZKulvRTSa8AR+VrSm32krS46LatsRvwUkTk579V0tcl/T7Vf6OkkbnpB0qan7bfrZLe32DZOwKPRsRNkXk1In4eEX/Kvbd8AN0BjJe0ZYPlWYdxEFgzHAzsR/ZN8UPAUQCS9ga+maZvDjwOXF4z70FkH2IT0+vNgLWB0cDpwPeBI4BdgL8GviZpXGq7EvgHYCSwB7AP8MX+io2IAyJivYhYD/hb4GngpgbNDyU7JLIRsBD4RnpvI4GrgVOBTYAFwF/1s+opaZ4Nyb41F1F329axfaqh1meBo4FNgTWBk1L92wKXAf8bGAXMBK5rEKTzgO0kfVvSxySt11fBEbGCbFvt0Fc76xwOAmuG8yLizxHxAnAd2TdIgMOBiyJiXkQsJ/vQ3EPSVrl5vxkRL0TEa+n1m8A3IuJNstAYCfxr+hY6H3iA9AETEXMjYlZErIiIx4ALgY8WLTp9GF4MHBwRTzRo9ouIuDN9uF2ae2+TgfkRcU2adh5ZoPTlDxFxbUS8lXu//Wm0bWttCLxaZ/yPIuKhtL4rc/MfAlwfEb9J2/r/AutQJ8wiYhGwF1k4Xwk8l/Zm+gqEV1NN1gUcBNYM+Q/AZUDvB8R7yfYCAIiIpcDzZB8ovWo/gJ+PiN5jz70fls/kpr/Wu3xJ20r6N0lPp8Mt/0wWHP2StAHwS+CrEfG7Ppr29d7erj2yqy4W07dGYdOXRuuv9SIwYgDz1/5t3kr15f825KbPioiDI2IU2Z7ZnsBpfdQ9Anipj+nWQRwEVqY/A28fJ5a0LtlhlCdzbQZz2dr5wB/JrgxaH/gKoP5mkvQu4GfALRExfTXX/RQwJrdM5V83UPte/xN4d+71ZqtZC2RX8Ww7gPa1fxsBY1n1b1NXRMwGrgE+WG96OhG+DXDPAOqxNnIQWJkuA46WtGO6quWfgTvSYZxmGAG8AiyVtB3wdwXn+wawLnDCINZ9PbB9Opk9HPgSA/8gvxuYLGljSZuRHa9fXXcCG0qq+42+jiuBT0raR9IaZFcELQf+o7ahpEnppP+m6fV2wIHArAbL3hV4LCIebzDdOoyDwEoTEf8OfA34Odk36K3JTr42y0lkJ0NfJTupfEXB+Q4DdgdezF05dPhAVhwRz5GdaP4W2eGuicAcsg/Ton5C9q35MeBGitdfr543gB+TnVgv0n5Bavsd4DngAOCAtJxaL5F98N8naSnwa+AXZO+9nsOBCwZQvrWZbygza4J0uGkxcHhE3NKmGkYBtwM7DeBkdLNr2BT4barh9XbUYAPnIDBbTZI+QXbN/GvAyWSHh8a360PYbHWVdmhI0kXppqD7G0yXpPOU3Wx0r6Sdy6rFrCR7AI/wl0MrBzkErBuVtkcgaU9gKXBJRLzj6gJJk4G/J7seezeya8V3K6UYMzNrqLQ9goi4DXihjyZTyEIiImIW2RUPm5dVj5mZ1dfOjr5Gs+oNNovTuKdqG0qaCkwFWHfddXfZbrvtWlKgmdlQMXfu3OfSDYHv0BU9PqabfqYD9PT0xJw5c9pckZlZd5HU8L6Odt5H8CTZnYy9xlDgrkYzM2uudgbBDODz6eqh3YGXI+Idh4XMzKxcpR0aknQZWY+FI1Mf62cAawBExAVk3d5OJuuudhlZV7lmZtZipQVBRBzWz/QguwHHzMzayH0NmZlVXFdcNdQM06dP5+yzm/5M79WS9fjbfq5jVa5jVa5jVZ1QxxlnnMEhhxzS9OVWJghGjx7NpEmT2l0GndK3k+tYletYletYVafUsdFGG5Wy3K7rdM73EZiZDZykuRHRU2+azxGYmVWcg8DMrOIcBGZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOzinMQmJlVnIPAzKziHARmZhXnIDAzqzgHgZlZxTkIzMwqzkFgZlZxDgIzs4pzEJiZVZyDwMys4hwEZmYV5yAwM6s4B4GZWcU5CMzMKs5BYGZWcQ4CM7OKcxCYmVWcg8DMrOIcBGZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnGlBoGk/SQtkLRQ0il1pm8h6RZJd0m6V9LkMusxM7N3Ki0IJA0DpgH7AxOBwyRNrGn2VeDKiNgJOBT4Xln1mJlZfWXuEewKLIyIRRHxBnA5MKWmTQDrp+ENgD+XWI+ZmdVRZhCMBp7IvV6cxuWdCRwhaTEwE/j7eguSNFXSHElzlixZUkatZmaV1e6TxYcBP46IMcBk4CeS3lFTREyPiJ6I6Bk1alTLizQzG8rKDIIngbG512PSuLxjgSsBIuIPwNrAyBJrMjOzGmUGwWxggqRxktYkOxk8o6bNn4B9ACS9nywIfOzHzKyFSguCiFgBHAfcADxIdnXQfElnSTowNTsR+IKke4DLgKMiIsqqyczM3ml4mQuPiJlkJ4Hz407PDT8AfKTMGszMrG/tPllsZmZt5iAwM6s4B4GZWcU5CMzMKs5BYGZWcQ4CM7OKcxCYmVWcg8DMrOIcBGZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOzinMQmJlVnIPAzKziHARmZhXnIDAzqzgHgZlZxTkIzMwqzkFgZlZxDgIzs4pzEJiZVZyDwMys4hwEZmYV5yAwM6s4B4GZWcU5CMzMKs5BYGZWccP7ayBpW+BkYMt8+4jYu8S6zMysRfoNAuAq4ALg+8DKcssxM7NWKxIEKyLi/NIrMTOztihyjuA6SV+UtLmkjXt/Sq/MzMxaosgewZHp98m5cQGMb345ZmbWav3uEUTEuDo/hUJA0n6SFkhaKOmUBm0OlvSApPmSfjbQN2BmZoNT5KqhNYC/A/ZMo24FLoyIN/uZbxgwDfgbYDEwW9KMiHgg12YCcCrwkYh4UdKmq/UuzMxstRU5R3A+sAvwvfSzSxrXn12BhRGxKCLeAC4HptS0+QIwLSJeBIiIZ4sWbmZmzVHkHMGHI2KH3OubJd1TYL7RwBO514uB3WrabAsg6ffAMODMiPh17YIkTQWmAmyxxRYFVm1mZkUV2SNYKWnr3heSxtO8+wmGAxOAvYDDgO9L2rC2UURMj4ieiOgZNWpUk1ZtZmZQbI/gZOAWSYsAkd1hfHSB+Z4ExuZej0nj8hYDd6TzDY9KeogsGGYXWL6ZmTVBv0EQETelk7rvS6MWRMTyAsueDUyQNI4sAA4FPlvT5lqyPYEfSRpJdqhoUcHazcysCRoGgaS9I+JmSf+9ZtI2koiIa/pacESskHQccAPZ8f+LImK+pLOAORExI03bV9IDZIebTo6I5wf1jszMbED62iP4KHAzcECdaQH0GQQAETETmFkz7vTccABfTj9mZtYGDYMgIs5Iv4ucDzAzsy7V71VDkk6QtL4yP5A0T9K+rSjOzMzKV+Ty0WMi4hVgX2AT4HPA2aVWZWZmLVMkCJR+TwYuiYj5uXFmZtbligTBXEk3kgXBDZJGAG+VW5aZmbVKkRvKjgV2BBZFxLL0LAKfQDYzGyKK7BHsQXYT2UuSjgC+CrxcbllmZtYqRXsfXSZpB+BE4BHgklKrMjOzlikSBCvSjV9TgO9GxDRgRLllmZlZqxQ5R/CqpFOBI4A9Jb0LWKPcsszMrFWK7BEcAiwHjo2Ip8l6ET231KrMzKxlivQ++jTwL7nXf8LnCMzMhoy+eh/9XURMkvQqWSdzb08i6y9u/dKrMzOz0vXV6dyk9Nsnhs3MhrAiJ4uRtBHZ08bebh8R88oqyszMWqffIJD0deAosieH9XYtEcDe5ZVlZmatUmSP4GBg64h4o+xizMys9YpcPno/sGHJdZiZWZsU2SP4JnCXpPvJ7icAICIOLK0qMzNrmSJBcDFwDnAf7n7azGzIKRIEyyLivNIrMTOztigSBLdL+iYwg1UPDfnyUTOzIaBIEOyUfu+eG+fLR83MhogifQ19rBWFmJlZe/R7+aikEyStr8wPJM2TtG8rijMzs/IVuY/gmIh4BdgX2AT4HHB2qVWZmVnLFAkCpd+TgUsiYn5unJmZdbkiQTBX0o1kQXCDpBH4fgIzsyGjyFVDxwI7AosiYpmkTYCjS63KzMxapshVQ29JegaYKKlQt9VmZtY9inRDfQ7Zc4sfAFam0QHcVmJdZmbWIkW+4R8EvC8ilvfX0MzMuk+Rk8WLgDXKLsTMzNqjUKdzwN2SbmLVvoaOL60qMzNrmSJBMCP9mJnZEFTkqqGLV3fhkvYD/hUYBvwgIurekSzp08DVwIcjYs7qrs/MzAauyFVDE8ieUjYRWLt3fESM72e+YcA04G+AxcBsSTMi4oGadiOAE4A7Bly9mZkNWpGTxT8CzgdWAB8DLgF+WmC+XYGFEbEoPfj+cmBKnXZfJ3sC2uuFKjYzs6YqEgTrRMRNgCLi8Yg4E/hkgflGA0/kXi9O494maWdgbERc39eCJE2VNEfSnCVLlhRYtZmZFVUkCJZLehfwsKTjJH0KWG+wK07L/BfgxP7aRsT0iOiJiJ5Ro0YNdtVmZpZTJAhOAN4NHA/sAhwBHFlgvieBsbnXY9K4XiOADwK3SnqM7AloMyT1FFi2mZk1SZ8ni9MJ30Mi4iRgKQPrbG42MEHSOLIAOBT4bO/EiHgZGJlb163ASb5qyMystRruEUgaHhErgUmrs+CIWAEcB9wAPAhcGRHzJZ0l6cDVqtbMzJqurz2CO4GdgbskzQCuAv6zd2JEXNPfwiNiJjCzZtzpDdruVaBeMzNrsiJ3Fq8NPA/sTdbrqNLvfoPAzMw6X19BsKmkLwP385cA6BWlVmVmZi3TVxAMI7tMtN7ziR0EZmZDRF9B8FREnNWySszMrC36uo+g3p6AmZkNMX0FwT4tq8LMzNqmYRBExAutLMTMzNqjSBcTZmY2hDkIzMwqzkFgZlZxDgIzs4pzEJiZVZyDwMys4hwEZmYV5yAwM6s4B4GZWcU5CMzMKs5BYGZWcQ4CM7OKcxCYmVWcg8DMrOIcBGZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOzinMQmJlVnIPAzKziHARmZhXnIDAzqzgHgZlZxTkIzMwqzkFgZlZxpQaBpP0kLZC0UNIpdaZ/WdIDku6VdJOkLcusx8zM3qm0IJA0DJgG7A9MBA6TNLGm2V1AT0R8CLga+FZZ9ZiZWX1l7hHsCiyMiEUR8QZwOTAl3yAibomIZenlLGBMifWYmVkdZQbBaOCJ3OvFaVwjxwK/qjdB0lRJcyTNWbJkSRNLNDOzjjhZLOkIoAc4t970iJgeET0R0TNq1KjWFmdmNsQNL3HZTwJjc6/HpHGrkPRx4DTgoxGxvMR6zMysjjL3CGYDEySNk7QmcCgwI99A0k7AhcCBEfFsibWYmVkDpQVBRKwAjgNuAB4EroyI+ZLOknRganYusB5wlaS7Jc1osDgzMytJmYeGiIiZwMyacafnhj9e5vrNzKx/HXGy2MzM2sdBYGZWcQ4CM7OKcxCYmVWcg8DMrOIcBGZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOzinMQmJlVnIPAzKziHARmZhXnIDAzqzgHgZlZxTkIzMwqzkFgZlZxDgIzs4pzEJiZVZyDwMys4hwEZmYV5yAwM6s4B4GZWcU5CMzMKs5BYGZWcQ4CM7OKcxCYmVWcg8DMrOIcBGZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOziis1CCTtJ2mBpIWSTqkzfS1JV6Tpd0jaqsx6zMzsnUoLAknDgGnA/sBE4DBJE2uaHQu8GBHbAN8GzimrHjMzq6/MPYJdgYURsSgi3gAuB6bUtJkCXJyGrwb2kaQSazIzsxrDS1z2aOCJ3OvFwG6N2kTECkkvA5sAz+UbSZoKTE0vl0pasJo1jaxddofrpnq7qVbornq7qVbornq7qVYYXL1bNppQZhA0TURMB6YPdjmS5kRETxNKaoluqrebaoXuqrebaoXuqrebaoXy6i3z0NCTwNjc6zFpXN02koYDGwDPl1iTmZnVKDMIZgMTJI2TtCZwKDCjps0M4Mg0/Bng5oiIEmsyM7MapR0aSsf8jwNuAIYBF0XEfElnAXMiYgbwQ+AnkhYCL5CFRZkGfXipxbqp3m6qFbqr3m6qFbqr3m6qFUqqV/4CbmZWbb6z2Mys4hwEZmYV1/VBIOkiSc9Kuj83bmNJv5H0cPq9URq/l6SXJd2dfk7PzdNndxgl1vq3kuZLektST037U1M9CyR9opW1DrReSVtJei23bS/ITdtF0n2p3vPKuGmwQa3nSvqjpHsl/ULShrlpnbht69bbodv266nOuyXdKOm9abxSHQvT9J1z8xyZ/k8+LOnIeutqQ70d95mQm3aipJA0Mr0ub9tGRFf/AHsCOwP358Z9CzglDZ8CnJOG9wL+rc4yhgGPAOOBNYF7gIktqvX9wPuAW4Ge3PiJqY61gHGpvmGtqnU16t0q365mOXcCuwMCfgXs36Ja9wWGp+Fzcv8OOnXbNqq3E7ft+rnh44EL0vDkVIdSXXek8RsDi9LvjdLwRi3cto3q7bjPhDR+LNmFNo8DI8vetl2/RxARt5FdcZSX77riYuCgfhZTpDuMQatXa0Q8GBH17pSeAlweEcsj4lFgYaqzJbWuRr11Sdqc7D/hrMj+1V5C/3+PAWtQ640RsSK9nEV2Lwt07rZtVG9dbd62r+Rergv0XnUyBbgkMrOADVOdnwB+ExEvRMSLwG+A/Zpd62rU20jbPhOSbwP/yKp1lrZtuz4IGnhPRDyVhp8G3pObtoekeyT9StIH0rh63WGMbkGdfWlUUyfW2mucpLsk/VbSX6dxo8lq7NWueo8h+zYF3bFt8/VCB25bSd+Q9ARwONB7SKVjt22DeqHDPhMkTQGejIh7aiaVtm2HahC8LX1T6k3VecCWEbED8B3g2nbVNQQ9BWwRETsBXwZ+Jmn9NtcEgKTTgBXApe2upYg69Xbkto2I0yJiLFmdx7W7nv40qLejPhMkvRv4CqsGVemGahA8k3aZenefn4Vs9zAilqbhmcAa6URMke4wWq1RTZ1YK+kwy/NpeC7Z8dVtyWrLH+Joab2SjgL+G3B4+lIAHbxt69Xbqds251Lg02m4Y7dtztv1duBnwtZk563ukfRYWu88SZv1UdOgax2qQZDvuuJI4JcAkjbrvapC0q5k7/95inWH0WozgEOVPbxnHDCB7MRgJ9aKpFHKnkGBpPFk9S5Kh+hekbR72vafJ/09WlDTfmTHWQ+MiGW5SR25bRvV26HbdkLu5RTgj2l4BvD5dIXL7sDLqc4bgH0lbaTsKr5907iWaFRvp30mRMR9EbFpRGwVEVuRHebZOSKepsxt2+yz4K3+AS4j23V+M220Y8m6sr4JeBj4d2Dj1PY4YD7ZFQCzgL/KLWcy8BDZt63TWljrp9LwcuAZ4IZc+9NSPQvIXQ3SiloHWi/ZN6z5wN1ku9sH5JbTA9yf6v0u6Y72FtS6kOzY6d3p54IO37Z16+3QbfvztN57geuA0amtyB5I9QhwH6teWXZMeo8LgaNbvG0b1dtxnwk10x/jL1cNlbZt3cWEmVnFDdVDQ2ZmVpCDwMys4hwEZmYV5yAwM6s4B4GZWcU5CMzMKs5BYGZWcQ4CsxxJS5u8vA0lfTH3eqt6fc83mHed1MncsCbUsaak2ySV9pxy614OArNybQh8sb9GDRwDXBMRKwdbRGRdKd8EHDLYZdnQ4yCwjifpZEnHp+FvS7o5De8t6dI0fK2kucqenjY1jTtb0pdyyzlT0klp+AhJdyp7KtWF9b5112uTvtE/KOn7aV03Slontf+asida/U7SZWldZwNbp2WcmxY9rN78dRxOru8gSddI+qf0zf5Pkj6exl8l6btpvY9LmiTpJ5IekvTD3PKuTcs0W4WDwLrB7UBvH/w9wHqS1kjjbkvjj4mIXdL04yVtAlwBHJxbzsHAFZLeT/bN+CMRsSOwkpoPyH7aTACmRcQHgJeAT0v6MFmfQDsA+6c6IHtC3iMRsWNEnNxo/to3nDo6Gx8Rj+VGbw+8FBF7Aifk6tmerBO6ScCFwA/JOq+bCHxS0lqp3f3Ah2vXZebjhdYN5gK7KOuDfzlZx2s9ZEFwfGpzvKRPpeGxwISImCVpU2XPpx0FvBgRT0g6DtgFmJ06nlyH1FV5zj4N2twGPBoRd+dq2woYCfwyIl4HXpd0XR/vp978tUaShQTwdj/1G5A9uQpgDeAlSWuTHX76/2l8AD+M9GAmSSuBNwAiYqWkNySNiIhX+6jPKsZBYB0vIt6U9ChwFPAfZD1IfgzYBnhQ0l7Ax4E9ImKZpFuBtdPsVwGfATYj20OArBfHiyPi1D5WW7eNpK3IwqjXSrKQGIgi87/GX94DZN/u5+bOF3yI7Bv+B4B5EfFWGr8DcH6qdQzw51i1Z8m1gNcHWK8NcT40ZN3iduAksm/ktwP/C7grfchtQPZtf5mk7cge7N3rCrK+5D9DFgqQnTT9jKRNASRtLGnLmvUVaZP3e+AASWtLWo/s4TIArwIjBvpmI3v27LD0jR+ywz9355p8iCwQtyfrQrl2PGSh0DtMOlz2XES8OdB6bGhzEFi3uB3YHPhDRDxD9q329jTt18BwSQ+SnZyd1TtTRMwn+yB+svdwSUQ8AHwVuFHSvWQP+948v7IibWrazyZ7cMi9ZM8avo/swSHPA7+XdH/uZHFRNwKT0nBtEHyQbI/g7fEpNNZJIQKrhgJke1HXD7AGqwA/j8CsSSStFxFL0/H824CpETFvEMvbGfiHiPhck+q7BjglIh5qxvJs6PA5ArPmmS5pItmx/YsHEwIAETFP0i2Shg32XoJ0FdK1DgGrx3sEZmYV53MEZmYV5yAwM6s4B4GZWcU5CMzMKs5BYGZWcQ4CM7OKcxCYmVXcfwGPzEDP4E/GogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot normalizing run (transmission through substrate alone)\n",
    "plt.plot(wavelengths / nm, np.abs(flux_empty), color='k')\n",
    "plt.ylim([0,1])\n",
    "plt.title('normalizing run (no Si)')\n",
    "plt.xlabel('wavelength ($nm$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normalizing run computes the transmitted flux for an air -> SiO2 interface, which is just below unity due to some reflection.\n",
    "\n",
    "While not technically necessary for this example, since this transmission can be computed analytically, it is often a good idea to run a normalizing run so you can accurately measure the *change* in output when the structure is added.  For example, for multilayer structures, the normalizing run displays frequency dependence, which would make it prudent to include in the calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-12T19:07:31.398409Z",
     "iopub.status.busy": "2021-11-12T19:07:31.398069Z",
     "iopub.status.idle": "2021-11-12T19:07:31.581770Z",
     "shell.execute_reply": "2021-11-12T19:07:31.582025Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEqCAYAAAACibeEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABMD0lEQVR4nO2dd3gU5fbHPyeNkNDbBQEVEZEivQiCXURsF3tDrvATGwqKBRWxF1RQFCvqFVREvSpiBRVQwEqRbkEEifQWSgop7++PdzbZhE3YJLs7s9nzeZ59ZnZ2duZkdzPfOee87zlijEFRFEVRgiXObQMURVGU6EKFQ1EURSkTKhyKoihKmVDhUBRFUcqECoeiKIpSJlQ4FEVRlDIRNuEQkddEZIuILC/hdRGRZ0RktYgsFZFO4bJFURRFCR3h9DheB/qW8voZQAvnMQR4IYy2KIqiKCEibMJhjPkW2FHKLucCk43lB6CWiDQKlz2KoihKaHAzx9EYWO/3PM3ZpiiKoniYBLcNCAYRGYINZ5Gamtr56KOPdtmi0LJwoV127uyuHa6yYgVkZUHz5iz9uxY5OZCYCO3auW2YolQOFi5cuM0YUz8Ux3JTOP4Bmvo9b+JsOwBjzMvAywBdunQxCxYsCL91EUTELivZn1U2uneHn36CIUOo89jt7NwJNWvG+GeiKCFERNaF6lhuhqqmA1c6o6uOBdKNMRtdtEdxE1+xzZ072bvXrmZmumeOoiglEzaPQ0TeBk4E6olIGnAvkAhgjHkR+AzoB6wGMoCrwmWLEgXs2QPA/p37yMmxm7KyXLRHUZQSCZtwGGMuPcjrBrghXOePFvLz3bbAIzjCsW+7VYt69WDbNgpyHYqieAedOe4yvrvrmGf3bgAydu0HoG5du1nDVYriPVQ4XEaFA+t2+TyO9FygUDg0XKUo3kOFw2VUOKAgGw5k7C4qHOpxKIr3UOFwGRUOCsJUAPuy4gEVDkXxMiocLqPCQUGYiipVyMi2wlGrlt2Um+uOSYqilIwKh8uocFDocTRoQEaW/UnWqGE36eejKN5DhcNl9MJIEeHYt9+Ova1Z025Sj0NRvIcKh8uocFDU49hvpxapcCiKd1HhcBkVDgpzHA0asC83CdBQlaJ4maiojluZ0QsjRT0ObHJcPQ5F8S7RJxzZ2W5bEFL273fbAg/gE4769dlHNnFxhmrVbMlgFVZF8R5RF6rKX74SbrwR1q8/+M5RgF4YscKRnAw1apBBCilVTUF9KvU4FMV7RJ1wrIg7hkUT5sPRR8PTT0f9DDEVDmyOo0YNSE4mgxRSq+aT4PjC+vkoiveIOuHYn59AZxbRt+ocNtz8OHTpAj/+6LZZ5UYvjFiPo3p1SEpiH6mkJOerx6EoHibqhOOww+xyxvauNGYD49efB8ceC/fdF5X5DxUOrHDUqAGJiTZUlaweh6J4magTjnr1bDHV4cPt8+F7HqR5tU2suP896338/rur9pUVvTBygHCkJucVCId6HIpSfvLyYMIEqFMntMeNOuEA26P7qadgyxZo1QrW7P0XbVnBc3/2hTZt4NlnC1uRehwVDgpzHImJNlSVlKehKkWpIMuWQc+edixR09QdIT12VAqHj/r1YeVKeOst+3xo5hN0SV1J+k2jYODAqGjmoMJBYY7DF6qqkquhKkUpJzk5MHIkdO4My5bmM/7oF/glrW5IzxHVwuHjsstg+3bo1QsWpregFun88MbvcNJJsGaN2+aVil4YOTDHkZSrHoeilIMlS2zEfswY6Nt6HcviO3LT+tuQZ58N6XkqhXCAjeHNmQOPPGKf9+AHnlvSC9OjJyxc6KptpeHXwyh2KS4ciepxKEpZyM2F0aOha1dYtcrwatun+GjJ4TTvWgeWLoWhQ0N6vkojHADx8XDnnbBokX0+NPMJTkl/n709ToP33nPXuBLYudMuq1Rx1w7X2L/fjobzT44n5ajHoShB8s8/cPLJ8OCDcMoxm1mT2o5Bf95tvYyvv4Yjjgj5OSuVcPjo2NEmzmvUgNnZx9E1YRFrL77DThj0GLt22WVcpfwmgsBX4NCZx2E9jv3qcShKEEyZYscDzZ0LY3t9yGeLGtKkeRVYvNh6GWG6sFTay1X9+rB1Kxx/PPyaeTi9k39i5c0v21iWh0Zc+TyO/Hx37XANX52qGjUwCU6oKkE9DkUpjf37bXTl8suhbo0cfjn6Em6Zdx4yfDjMnw8tW4b1/JVWOACSkuCbb+DllyEtsx5tWMn3d39sldgj4uHzODxiTuTxE46s/CQMcaQk7CfeFslVj0NRivH339CtGzz2GPxfn3UsSz+U9ptnwrRpdp5CBOLelVo4fFx9NXz0EVSrZujJ97zz/Da70QO3s75IjQpHDTKcXhwpCdmI2JyVB74iRfEMX39t52YsXWp469/vMXHm4aQ0b2QHAJ17bsTsiAnhADjnHPj0U1uq+xLe4f5XG8OQIXZqpYv4pprErHD45Th8wpEab0vHJCaqcCgK2OvDhAlw6qkgJo/vuw3nsmkXwVVX2dBUs2YRtSdmhANsvmPJErt+H/cz7r+1MJdf4ap4+Ir7ao6jBhk5NrGREm/VNCFBQ1WKsn8/3HCDnQHes2MmP8f3oPviF+Gll+DVV6Fq1YjbFFPCAdCuHaSl2fURjOOWd7phrrnWtSt3zHscgYQjTj0ORQHYtw9OOAFeeAGuPGMLc9Y1o2HWWvj2WxsxEXHFrpgTDoDGjW2EpFEjeJqbufjV0+Daa125eqtwBBIO9TgU5a+/4Jhj4IcfYMyVy5k0+zAS61SH77+H7t1dtS0mhQOgWjVYuxbq1oX3uIgrJh5vfcEIX8H9y2nFpHj4chypqWTst5M3UuJs/E49DiVWWbDAJsE3bID/XT2D299sb8Ml330HzZu7bV7sCgfY4bobNkDt2oa3uII7n2sMd98dURuiQTgmTrReWljwFTiMi2Nfth2DmyLqcSixy6xZ0KcPpKcbvrn4ec6f2Bf69bMv1K/vtnlAjAsHWPFISxPatDE8xp088GiCLcseIfw733pVOIYMsQIbFvt8daqAjCwrHKmSAVjhUI9DiSU++ghOOQXq1zMs7jOS7pNvsP+AH34Iqalum1dAzAsHQEoKzJsn9DjWcC8PcM9N6fYbDDPGWI/Dl9/yqnD4CLtwWL0gBbuioSollpgwAfr3hxbN85nV8DJafvS4LUD14osU1ODxCCocDrVqwcefCCkphocYxaPn/QwzZ4b1nPv326VvNJ3Xh+SGZdTynj02VMWBwqGhKiVWGD/epljbt83lpxqn0vi79+C112DUKNdGTpWGCocfdevC33/bL+mu/Id47dxpsHx52M63b59dpqTYpdc9jrAIWyCPw9gPRj0OJRYYO9a2wu7eOZf5CSdSa/k8+OADO7nPo6hwFKNuXdi8Gdq33s/grOeZdtJ4uyEM+A7bsKFdqnCAkE+VfJv4UY9DqeyMGQO33gr9Tsthdv4JpCz/Cd5/35a68DAqHAFo0ACmvp9Egzo59N82kS9PebTwdjiEbNxol02a2GXMCodfqCpVMpBcqxbqcSiVmUcftS1ej+uey3vbTqTqigXW0zj7bLdNOygqHCVw9NEw42s7Ia3PiqeZdeojIQ/yb9hgl76hrjGb43A8jn37nKG4jpuhHodSWXnhBbjrLutpzMw+gZQVC+zIqbPOctu0oFDhKIUOHWDZMrt+yvcP8dOgF0N6fF+oqlEju4w5j8OYA0JVKXGZBWqhHodSGRk/Hq6/Hnp0y+N/u04lZaUjGv36uW1a0KhwHIS2beHbb+wVvfvkG0h79sOQHXvLFls6v2ZN+zzmhCMjwx7UXzjis9XjUCotU6faRHinjvl8mXw2VRfNh3feiSrRgDALh4j0FZHfRGS1iIwM8PqhIjJbRBaLyFIR8eSn1/t4YcLT9ta36U392f7lopAcd8sWm0/xdXeMOeHw1alychz79jkFDv2EQz0OpbLwxRdw2WXQ6mjD7AaXkDr3C3j9dfj3v902rcyETThEJB54DjgDaA1cKiKti+02CnjXGNMRuAR4Plz2VJQbhiUw4VFbV6lxn9ZsXV7xkVabN1vh8A3Tjrkch1+BQ7DCUS0hs2CCi4aqlMrC99/DBRfAoYcavm55PTVmvAfPPw9XXOG2aeUinB5HN2C1MWaNMWY/MBUo3qLKADWc9ZrAhjDaU2FuGFmdp2/7h2yS6dE1h307sit0PJ/HES0zx0MubL4Ch45w7N0L1RI0Oa5ULtassTPCwTCr9300+uhFOw732mvdNq3chFM4GgPr/Z6nOdv8uQ+4QkTSgM+AG8NoT0i4aUxjzuy0gT+zmtC1+Y4KXdhiXjiKeRxWOApDVfHxrjdoVJQK8c8/0KMHZGfDnEFvcMSbD8CIEXD77W6bViHcTo5fCrxujGkC9APeEJEDbBKRISKyQEQWbN26NeJGFrUFPl5wCP9usYJVuxrRvOHeIhVug8UYO49DcxwU5Dj27oVqiUVzHCocSrSyezdcfLG9QXzvull0eXag3fD4426bVmHCKRz/AE39njdxtvkzGHgXwBjzPZAM1Ct+IGPMy8aYLsaYLvU9UFZYBN75pSXn1p/P+h3VuOWKLWW+6H/8sY3fa46DIsnx1MT9RTwOzXEo0Uhurs15f/cdvHH3r5z6ZF/bym/SpMI7xSgmnH/Bz0ALEWkmIknY5Pf0Yvv8DZwCICKtsMLhrksRJEkpCbyzuCWtE3/nhfcb0LFdbpnEY+lSuxwwIIZDVenpdlmzJsY4HkfSfvU4lKjGGDsjfPZsuPXqdK54qTcccYSdq1GlitvmhYSwCYcxJhcYCswAVmFHT60QkQdExFeIZQRwtYgsAd4G/mOM1y+fhVRpXI9l83bThPUsWZ7AySeZoC/+GzbYulj/+pcKBzVrkp1tRaJa0v4CN0M9DiUamTDBFi4cPDCHMd8fb2+Epk+H2rXdNi1khNVnMsZ8Zow5yhjT3BjzsLNttDFmurO+0hhznDGmvTGmgzEmvHXMw0Bcty78+fxM6rKNOd8I118fnACsXw+HHOIcI1ZzHOnpkJwMVaqwd6/dVK1KThHhUI9DiSbmzbN577ZtDc+k/wdZsdxO8DvqKLdNCynRH2zzAEnXDWbdZXfRk/m8+CKcdFLpF7y1a2HGDOjVyz6PlhxHWITDmTbvKzGfmlQoHBqqUqKJP/+E006DevVgRt+nSZk2xSbCTz/dbdNCjgpHiEh9+SlmtLoZgG++sRe9nTsD7/vcc9Z7HenMpY+WUFXIL+J+wlGSx6GhKiUa2LsXLrrI/l7/d+M3HPLkLTaBecstbpsWFlQ4QkVqKtU+eous1LocXXUtAHXq2Bmj/mzbBk8+adcPPdQuYzpUVYpwqMehRAP5+XD55bBoEbz16N90f/As6NYNXn7Zk937QoEKRyhp0YIqLz/LL5ktGdB6AQA9e9rfziefwLRp4BtN3KJF4ds0VOUnHMm56nEoUcXzz9vc93VXZXHR8yfa4eUffmjzd5UUFY5Qc9llVLluMJNXduW9uxYXbD77bF/ZAfu78u9IGy03JWERjlq1AL8ch59wqMeheJ1Zs2yv8ON65jNh7Vl2uOSHHxaOfKmkqHCEg7FjoW1bLnilL+sXbOaaawpfSkqyv62kpAPf5vVQVURyHMl5OgFQiQq2b7chqlq14N12DxM3+2sbnure3W3Two4KRzioWtUOwduzhya3XsKLE+zkQGNszZpq1YruHrMex65dBwpH1bwioSpjvB/CU2IPY2xb8E2b4JO7vuOQF0fDddfBlVe6bVpEUOEIF61b2/6Qc+bAvfcG9RavexwhvYDn5NjOTcWFIyW/SKgKNFyleI9HHrHlREbdtJvjxpwD7dvDuHFumxUxVDjCycCBMGiQ7Uo/a1aJu8Wkx+GrU1VMOFKr5hfxOECFQ/EWs2fDqFHQ9/R87l94lg0jvPtupU6GF0eFI9w884ydNTpggB2LG8WE9ALuV24EbHI8ORnik+KtQuXnq8eheI4tW+D8821eY1KrMcTNnwsvvVTpZoYfDBWOcJOaCm+/bUVj0KBS41ExFaoqJhx79zq5Hz+18HkcmiBXvEBeHgwZYlNz796xkAbj74bBg20/2BhDhSMSdOxoO359/LEd9F2MmAxV+YTDGY57gHDk5qrHoXiKsWPho49gxDV7Oe2pfjaP+cwzbpvlCiockWLYMOjXz3b/8p/E4UesexypqRQRDvU4FK/w228wejQcc4zh4V/Pt22P330XUlLcNs0VVDgihQj897+2TeqAAbB/f5GXooFw5jj27HH6OfmEIydHk+OKJ8jPL6xT+PHZE0maM9PWTm/d2l3DXESFI5I0aACvvAK//AL33HPAyzHlcezYYZdOj4KCuYABQlXqcShuctttsG4djBm+kcOevNGWgLjqKrfNchUVjkhzzjlwzTXwxBMFQ3SjxeMIqXBs326rOzo5jgLhSEy0r/uFqtTjUNxizhx4+mk456x8hn5xlv29vvRS9PzThgkVDjcYO9YO37vyysI7b2LQ46hdu6A0cHq6jeJpclzxCrt22RBVrVrwypGPEb9kEUycWFipNIZR4XCD1FSYMsUOCh8yBMHjiuEQ0gv49u22d65DoFCVJscVN3nkEZuKfOT/1lD/mXvscPpzzjn4G2MAFQ636NQJHnoI3n8f5n4LxJjH4SccOTmQmVlyjkM9DiXSvPmmjSZfemEu13xwOjRtCk895bZZnkGFw01GjIATT0QmT3LbkqAIuXDUqQMUqz6iHofiMjk5hY37Xqhxh+0JO2mSE0tVQIXDXeLjYfJkiLcXS7M/x2WDSiekF3A/j6PIyNwAwqEehxJJbrgBtm6FN25fRs1Xx8HNN8MJJ7htlqdQ4XCbpk2R/xtk1wPMKvcCvta2OaHUtYMJR06ODsdVIs60aTb/fcOQHK549xxo2RIefthtszyHCocX6H4sAOaFF2HePJeNORDfCFm/OYsVIzvbVjVUj0PxEFlZ8OCDdv028zisXWtVJIaq3gaLCocHKBgS3rgxXHFF4ZXUI/i6FYZMOHxDkAMJh988Dk2OK5FkwABYtAjefuQvDnt1tK1o2Lu322Z5EhUOD2GeeBLS0myQ1UP4ruXZ2SE64PbtdllMOIrP49DkuBIpfv4Z/vc/6N3LcMl759sqD2PGuG2WZ1Hh8AAFHkeHDraS2ltv2XkeHiHkHodPOJxRVUV0RIfjKhEmN9cWcwD4b69XYfFiePbZgqoGyoGocHgIY4C77oKePW3/4rVr3TYJCEOOo5jHsX27HWCmw3EVNxg2zGrFa2O20vyZYXaS3/nnu22Wp1Hh8ABFyt4kJNjZR8bYoKsHbrfDLRzbtlnnIy4OTY4rEWXjRjuYMSnJMODLK+2PcMKEmK9FdTBUODxEwczxZs3sr3nePHjsMVdt8idkwrF5s102aABYHalXz3lNq+MqEeSmm+xy3t1fkPDVF/Doo3aWuFIqKhweIODNzeWXw6WXwr33wk8/RdymQIRMODZutC5GlSqA9TgOEA7tx6GEmddftwnx22/MoOuzV0L37jZErBwUFQ4PUaRWlYj1Oho3tj2N9+513a6QCcemTdCwYcHTbdv86h1qclyJADt22JYaVavCqO232FK4EydScLeilIoKhwcoMZxaqxa88QasWWMzeC4TUuFo1KjgaZFQVYB+HBqqUkKNr0jDmEG/UX3KS3DHHXDMMe4aFUWocHid44+HO++E116zlXRdIOQex8aNBR6HMSWEqtTjUMLEihVw331w2in53PjVuXDkkTBqlNtmRRUqHB6ixLLq990HXbrA1VfbCYIuERLhMKZIqGrPHlsDK1CoSj0OJRy0a2dvRl7q+gr89hs884yWFSkjKhwe4KAj/xIT7YTA7GwYODDE9c0Pjk/QQjJzfM8e23zDEQ7fyFz1OJRI8NVX9t+nV7dsmk0YYedsnHGG22ZFHSocHqLURk4tWtg7o1mzYNy4iNnkT2ZmCA6yaZNdOsKxbZt9Gkg41ONQQsnu3XDaadCqFXxy6A3W1dXmTOVChcMDBD3XaNAg6N/fzi5fvDisNvnjE7SQCMfGjXbpJMd9wlFaqEo9DiUUjBhhl3eft4qa/3vVJsSPOMJdo6IUFQ4PcdDWsSJ2yGD9+naeR0ZGROzyEZLTFfM4SgxV+fXjUOFQKsratfDKKxAfb7h02sVw2GFWOJRycVDhEJGjRGSiiMwUkVm+RySMixXKVN2gbl3bxnLVKrjttrDZ5E9IPY5//rHLYh6HhqqUcJGfD3372vVFI6YQt2KZDVGlpLhrWBSTEMQ+7wEvAhMBvfcLIwf1OHyceqr1u8eOtYm9s84Kq10+QiIc69ZBtWpQuzZghSMuzilwCNqPQwk5779vB0+dcUo27V66Afr0gX//222zoppgQlW5xpgXjDE/GWMW+h7BHFxE+orIbyKyWkRGlrDPRSKyUkRWiIh3aolHkHLVU3v4YWjf3uY9fOGfMBFSj+Pvv22YwPmjfR1kfe1pC1bU41BCQF4eXHSRXX+r6Z228+T48VrEsIIEIxwfi8j1ItJIROr4Hgd7k4jEA88BZwCtgUtFpHWxfVoAdwLHGWPaAMPL/BdUIoL2OMDWeZoyxQ5vveqqMr65fIQkx7FuHRx6aMHTIuVGwP5DJySox6GEhEmT7HJw/x3Unjwerr8ejj7aXaMqAcEIx0DgNuA7YKHzWBDE+7oBq40xa4wx+4GpwLnF9rkaeM4YsxPAGLMlWMMrE+W++WndGp58Er74wpaCDhMh9TjWrbMeh0ORWeM+HOFQj0OpCNu2weDB0KSJ4fn0y2089N573TarUnBQ4TDGNAvwCGYMW2Ngvd/zNGebP0cBR4nIfBH5QUT6Bm965aNcTsP110O/fjZRvmJFyG3yJzOzgo7N3r22upyfcBSpU+WjmHCox6GUhwED7HLc5YtImvWFrcBQ56DBEiUIghlVlSgiN4nI/5zHUBFJDNH5E4AWwInApcBEEakVwIYhIrJARBZs3bo1RKf2DhUKt4rYOlY1a9oqullZIbPLh79YVGj2+N9/22VpoSooEA4Rm/JQ4Ygtpk2DDRsqdow//4Qvv7Tr5304AFq21JLpISSYUNULQGfgeefR2dl2MP4B/DuiNHG2+ZMGTDfG5Bhj/gJ+xwpJEYwxLxtjuhhjutSvXz+IU0cn5b6b/9e/4L//haVL7eTAMFKhcNW6dXbpeBzGlOJx5OQAtsq1hqpiB2PsHNdevcp/jPx8W9oNYMWdbxL/+yob0k0M1f2uEoxwdDXGDDTGzHIeVwFdg3jfz0ALEWkmIknAJcD0YvtMw3obiEg9bOhqTbDGVxZK8zh27bJNyebOLbiWBqZfPxg61I5PnzkzpPb5C1qFEuS+HuqOcKSn28KJJYWqfKvqccQOvjJsf/1V/mN89pn9v+nULpfWL95kh6+feWZI7FMswQhHnog09z0RkSMIYj6HMSYXGArMAFYB7xpjVojIAyJyjrPbDGC7iKwEZgO3GWO2l/WPqCwU9zimTLHTHe66y1ZXT0oq7CMQkMcftwnzgQMhhCE9YwqLh1bI4/j9dzvp6pBDANtmBGyn3CIkJhYIh3ocsUVF63fm58PZZ9v19zs+ZO9Oxo3T4bchJhjhuA2YLSJzROQbYBYwIpiDG2M+M8YcZYxpbox52Nk22hgz3Vk3xphbjDGtjTHHGGOmlvcPiWYC/abXr7dVRXw47bm54QY46qgSSpxXrWrVZudOuPLKkFbR9U2yrZBw/PqrjTU7czX+/NNubt682H7qccQsFf3JPvGEXfY7YR9NJz8M//d/2qApDAQzquprbN7hJuBGoKUxZna4DYtF/D0O3+jaZ56x2zdvhtWr7bY//oDjjrPu+AG0bw9PP22H6D7+eMjsCplw+I2h91Ufadq02H5+wqEeR2xRkVF7e/fCSGea8ft1h9i5Tg88EBrDlCKUKBwicrKzPA84EzjSeZzpbFNCRHGPwxj45BObIBw6tHB78+b2wn3ccbBggQ1j7dkT4IDXXAMXX2y7ms2dGxIbq1a1y3LnODIzbXLcTzg2bbJRqQNGSBYTDvU4YoeKeBy+NMYj16eR/MEUuPVWO3BECTmleRwnOMuzAzwiUxwpxvDdbc2cCStXwiWXHCgqyckwb57tPwNWWA64IxeBl1+2yYNLLy2sJFgBuyrscfzxhz1Qy5YFmzZvtiG4A0J1GqqKWcorHJmZ9v8CYPDSYbaC9IigIupKOShROIwx9zrLqwI8BkXOxMpP8Qunr9XGlVeW/J6PPrI3VJMn28T5AdSoAe++a0UjBPkOn3D89ls5D/Drr3bp53H4tR4vioaqYpby/EyNsTdSxsB7o5bQYN4H1tuuXj30BipAcBMAh4lIDbG8IiKLRKRPJIyLNXwex9q1dojqwX73Y8bYyXPff19CxZGOHe3w3M8/L8waltOuo46yF/FyR76WLbMHOOqogk0rV5ZQNshvHod6HLFFeYRjxgzbEhYMF3zyH+tpX3NNiC1T/AlmVNUgY8xuoA9QFxgAPBZWq2IMf4/DGJg9u0hEp0Ti4mwEqGVLuPFG+PDDADtde60tD3r33RXKd6SkQNeuNgFZLhYtskOFnWTJnj125FibNgH2VY8jZimPcPhahv/3mh/hl1/gwQdtYlwJG8EIh++y1g+YbIxZ4bdNCTFbt9rpDuefH9z+tWtboUlKgvPOs+8tgq9r4BFHwAUXFJb9KAM+TyglpQLJ8UWLrAfkkJZml4cfHmBfv3kc6nHEFmUVDv/Q6cCZl9tRhZdeGlqjlAMIRjgWishMrHDMEJHqQOgmCCgFGFNYo8evDuBBadQI3n7brl92WYCLe40aNimSlQXnnmt7EpQRkQoIx8aNdghVp04Fm3zC0bh42UtQjyOGKctwXGMKQ52/3fOmnVH62GN+zV2UcBHMJzwYGIktPZIBJAJXhdWqGMM/VOUTDmdyddCcd56d+7dwYQnJ8latrLosWQL/+U+Zbu38PY5yaE5htt9PONY7dZObNAmwvw7HjVnK4nH85z92eUybPI56aQSceCKcfno4zFKKEYxw9AB+M8bsEpErgFFAenjNik38PY6yCgfY4btHH23F45lnAuzQr5+dFPi//9k4cBkQgdTUcnocP/xg7wI7dCjYtGiR7SAb0LMqNhxXPY7YIVjhMMaOKAR4vMPbsGWLHS2ipUUiQrDVcTNEpD221MifwOSwWhVjBPI4Ag5TDeI4vvz3gw+W0J5jxAg7PPe++6yABIHP44iLs57Cxo1lNGz2bOjcucgwsSVLrI74em4UQT2OmCVY4XjtNbtsdVQup0+7ziYFu3ULn2FKEYLtOW6w3fsmGGOeA3SAdBjweRz169tkd3moV88Oc927F9q2tX2TiiACL70Exx5rBcQXRjoIIvDBB3b90UfLYFBGBvz4I5x0UpHNf/0FRx5Zwnt0AmDMEoxw/PCDLUEFMOf4e5GsTHj44fAaphQhGOHYIyJ3AlcAn4pIHDbPoYQIf49j40ab7K4IrVrZhoBQ2AWtCMnJduxu3bo2Wb55c6nH83kc7drZpa/gYlB8952dk3HiiQWbMjNtnaoDquL60H4cMUswwnGCU9Oi1ZH7aTDpCdsfNpjx60rICEY4LgaygcHGmE3Yhkzln02mlIgxtnBh7doVP9YDD0D37rY3weefB9ihYUOYPt3OLO/f/6Ct/UQKI1tlCiN/8YUdXuvXmWfJErsssWipehwxy8GEY9euwsrQ8zrcaH8g2kc84gRTHXeTMWacMWau8/xvY4zmOEKI/4V49247ejYUfP21DXv161fY+6IIHTvCpEl26vmQISX+1/o8jnr17Lyq9GCHRhhj+4CeckqR/MYPP9hl9+4lvE/7ccQspQ3Hzcsr9KAbN9hPnfcnwrBh5RtJolSI0qrjznOWe0Rkt99jj4jsjpyJsYMx9qJcs2Zojpeaauf+AZxVUlnKCy+0ifLJk23FxBL+c33ilp1dhuolK1bYphv//neRzT/+aIfhlvj/rsnxmKU0j+P0023VaIDVxw6wd1i33x4Zw5QilFbksJezrG6MqeH3qG6MCdE9sQLh8zjApjB694ZVq+Ctt0rYafRomxR54YWA4hFIS3YHc+vw9tv2j/OV8nX48cdSvA3Q4bgxTEnCkZNjPWiAV0b/TfL0d+Hmm0MT11XKTFBTLEWktoi0E5FOvke4DYtFjAm9cIAt096sGVxxha1tdQAidgz8rbfa3rQ33XSAWvjEzVdb6qBDcrOybGn3s88uku3futWOqApWONTjiC1KEo7RowvXe8x/0rrlw4ZFxijlABIOtoOIPAj8B1hDYakRA5wcPrNiC99FOTPT3lmFKlTlIznZ3vwfe6wtTpufHyDBLWInB+bnF/ZoHj8eRIpoyPjxcOqpVjhKHcjiK+l+441FNv/4o12WxeNQ4YgdAgnH2WcXhqh6ddxL66+fteHVWrUiaZriRzAex0VAc2PMCcaYk5yHikYY8CWdQ+1xgL1Q+3IKAWeVgxWLJ5+0IYBnn4Xhwws8D5/Q+Nq8/vxzKSfzic/RR9vEuB8//mi9iM6dS3m/1qqKWfyF4+237Qgqn2gATGt0vXobHiAY4VgO1AqzHTGN76LsyxuEQzjADoOtUsW6/b5+3wGNGTvW/mM+8wycfz4mv9DlaNECeva0g7FK5PXX7cnuuecA12b2bDsMNzW1lPdrP46YxV84LrusaHX0RW//Rt3P3rA3NuptuEowwvEosFhEZojIdN8j3IbFIj6PI9ShKh/16sGXX1qBatu2lAuyiG0A9eST9nZv315kQ1rBS+eeawdMBexIu3q1FZ3evW3xLD82bID58+1ArlJRjyNmKW04bsd371RvwyMEIxyTgDHY5k1j/R5KiPDdlIczVOWjd287s3zXLnjvvYMYNWIEfPcdBoF337OzCvfvp317u4vfZHDL9u22ZlBCArz55gHlrWfNssti1UcOJDHR3nrm52tyPMYoKTlet1aurXag3oYnCEY4MowxzxhjZhtjvvE9wm5ZDBIJ4QD49FO7vPTSIO7mu3SB1FTk6JZ2hu4RR9D9R5skWbHCFiUlNxc+/thWLfz1V3jnHTj00AMO9cILNkfS6WBj8hKcMRt5eTocN8YoSTjW9R6g3oaHCEY45orIoyLSQ4fjhofiHke4QlU+mjWz8WMoZW6HH8aI7c/5+efQsiW17h1GbWz1xNEdptvyJeecY4dvzZ8PfQ5sSb9qlS1bdeGFQXT19AlHbq56HDFGIOH4+c3fSP14qnobHiIY4egIHAs8QmGY6slwGhWrRMrjAJu/PuYY2wxn9eqD7y9xAn372llYixbxwMnW6Xxp4zlknX6uDSMsW2Y9lAD4ipdedFEQxvkJhybHY4tAwtHlfc1teI1galWdFOChw3FDSKRGVfmTmAj332/Xzzij9H0PSFh27MjQr/sXPH25+6u2rEhycsD3b9xoPZtWrQ4yf8NHMY9DQ1WxQ8BQleY2PMdBhUNEholIDbG8IiKLROTAWIRSYdLT7bW3vL04ysqZZ9rl6tW2I19pBKqI6+vANmyYnRFeEoMG2WXQTQfV44hZ/G9SGjeGmcfdr96GBwkmVDXIGLMb6APUBQZgR1gpIcI/xxEJb8NHUlLhZMDSRliVNERywAC44w673qABvPLKgft8+aWtrN61q+2LHhQ+4cjJUY8jxvB5HLfcAmmf/MJp8+9Tb8ODBCMcvnvNfsBkY8wKv21KCIm0cICtCNK7txWQTZtK3q+kHhyPPWZFAeDqq+Gqq+zcvXffte/p08fWt/r88zL08dDkeMziE44zz8QO/1Zvw5MEIxwLRWQmVjhmiEh1CmtWKSHAP8fh17YiYjz/vO3wWqyIbQGlTcoC+OqrwvXXX7eezMUXF27zNRsMmkSnwaQTqlKPI3bwCUfcmtWa2/AwwQjHYGAk0NUYkwEkAVeF1aoYJTu7xPxyWPFVvP3555IbAZbmLdSoYYUn0IipVatsmZIyoR5HzFIgHJNfV2/DwwQzqiof2Ay0FpHjgTZo7aqwkJ8fucS4PyI25AR2eG5xDuZxAFStauf9GWP/jrw8u3700eUwqFhyHILrRa1EPwXCMXeOehseJpiy6mOwfcdXAr57PwN8G0a7Ygr/u3k3hAMKi9hOnWqrkhanLH3GRcrYl7w4xTwOZ9W1z0aJHD7hkNRU9TY8zEGFA/g30NIYU0IQQwklvvB+pOnSBdq1g6VLYfNm+Ne/Cl8LxuMIKQE8Dg1XxQZm9Z9Ac+IuuUi9DQ8TTI5jDeDS5Sw28ILHAXZEFAQemlshD6KslOBxKJWXRx+1v7GJD28BsMKheJagihwCv4jISyLyjO8RbsNiFTeFY/hwO+nq+++LbnfN43DmcYB6HJWdu+6yy/c39AAgrqYLwwuVoAkmVDXdeShhwv9u3q1QlY9jj4Vvv7V3+Al+vw63PA6/VSWGiAvmllZxjYMKhzGmtF5vSohxOwF8zjnw/vswY0ZhSZKIexx+8zjU46j8vP76gdtUOLxNMKOqWmC7ALYGCmYZGGOOCKNdMYVXchxQWBF97dqi2932OFQ4Kifr1xfm1vxR4fA2wXw9/wVeAHKBk4DJwJvBHFxE+orIbyKyWkRGlrLf+SJiRCRwTe4Ywm3haNDALocOLdzm5qgqTY5XXvLyAvb7AiJ8o6KUmWCEo6ox5mtAjDHrjDH3AWce7E0iEg88B5yB9VYuFZHWAfarDgwDfiyL4ZUJL+U4/O/0cnIK19XjUEJNaRWV1ePwNsF8PdkiEgf8ISJDRaQ/UC2I93UDVhtj1hhj9gNTgXMD7Pcgtqd5VrBGV2bc9jgABg60y1277FI9DiUcjBpV8msqHN4mmK9nGJAC3AR0Bq4ABgbxvsbAer/nac62ApwWtE2NMZ8GZW0lxUs5DijMc+zYUbjNbY9DhaPy8eqrheszP9hb5DUVDm9T6tfjhJsuNsbsNcakGWOuMsacb4z5oaIndryYccCIIPYdIiILRGTB1tL820qA26EqgDp17HLaNLt0cx6HhqoqJ/410Sb3/5DT+lcr8jtT4fA2JX49IpJgjMkDepXz2P8ATf2eN3G2+agOtAXmiMhabF/z6YES5MaYl40xXYwxXerXr19Oc7yLl3IcALVr2+VIZziDMe57HP75FiW6mTEDJjmD/Gek9GfAaycdsI8Kh7cpbTjuT0AnYLGITAfeA/b5XjTGfHCQY/8MtBCRZljBuAS4zO/96UA933MRmQPcaoxZUMa/oVLhBeHwNWby7w8eUeHwm8fhW1XhqDzcfLNd1mQXfW7vELAmlY6q8jbBzBxPBrYDJ2Or4oqzLFU4jDG5IjIUmAHEA68ZY1aIyAPAAmOMzkZ38P8nSQjmGwkzcXFwxhmlj3oJK5rjqLTs2GF7tAA8XvVeGHZ/wP3U4/A2pV2mGojILcByCgXDR1BRb2PMZ8BnxbaNLmHfE4M5ZmXHC8IB9iZw9erC526FqtTjqFzMmlW4PuSOOiVWwFXh8DalXabiscNuA10yIp0urdT4X5R9w0/dpmZNOxw34olxCCgc6nFEP7t3w6WXQrzk8Xe1NjCs5DE2KhzepjTh2GiMeSBiliiAdzyOOnVg587CC7Ymx5WK0rix/T09yzAOGXFpqf02VDi8TWmXKU1PRQiv5TgAWre2/+S+eHRE8VMLDVVVDjZtgr3OVI1zq8+GYfNL3V+Fw9uU9vWcEjErlAK8Ihxt29qlTzjcznFoqCq6ecCJXXTjR5qOOHh3PxUOb1PiZcoYs6Ok15TQ4sUcR926drlzpwsn96szoqGq6Cc7G775xq5PrT4Ehn1z0PfocFxvo7ruMbzicfhuCH3CEdF/ZBH7QeioqkrBgw/CypUwhUtpNuK8oHqJq8fhbTxymVJ8eEU4UlOtLa54HFAgHDqPI7rZsQMeftiuX1Ljcxj2QlDvU+HwNvr1eAAvJsdF7I2hKx4HqMdRSbj8cru8nueQW24OytsAFQ6v45HLlOLDKzkOKCocEaeYcKjHEZ188YVd3ld9HAxbGPT7VDi8jX49HsCLHgdA9eqwZ49dd8vj0OR49PLf/9rl/Yym/ogrg/Y2QIXD63joMqWAd4Uj4iQk6DyOKCYzE0aPhsNTNnNHwosw7PcyvV+Fw9vo1+MB1OMIgIaqopqXXoK0NBiVcRdVbrmhTN4G6HBcr+Ohy5QC3spxVKvmoseRmKihqiglPb2wdPrgmu/DsLVlPoZ6HN5Gvx4PoB5HADTHEbWMG2eXA3ndKkgZvQ1Q4fA6HrpMKeBd4Yg4jnDExdmLiIaqogNj4LXX7PqzNUbBsOXlOo4Kh7fRr8cDeNXjqFat8E7fLY8DbNRKPY7oYMoUm9t4nYFUHzGkXN4GqHB4HQ9dphTwVo6jenUXT15MONTj8D7798PVV0Pnmn9whXwKw/4s97E0Oe5tVNc9gFc9Dn/hcNPjcEbmKh7nhRfsMNwB6ROIv/Vm2w1MqZR46DKlgLeEo1o1F0/upxYaqvI++/bB8OF2/fraU+HGP1y1Rwkv6nF4AC+WVQfveBwaqvI+n3xilzczjsTbb4YaNdw1SAkrHrq/VcC7whFxEhNt0BwNVXmd/Hy4+244NHkzj6SOg6G/um2SEmZUODyA/928l0aT+IeqXPE4MjIADVV5nXffhT//hLcZRvIDw1yOcSqRQIXDY3jJ40hJcfHkGqqKCnJz4dproWXK35yXOg+uf9Vtk5QI4KH729jFqx6Hf6JeR1UpgXjmGVti5MaMMSTdcbPtAKZUetTj8BheFY6I45fj0FCVN8nIgCeegNT4TIbU/QiuK1sFXCV6UeHwAOpxBCApqUhyXENV3mPaNNi0CSZyI4l33+5ybFOJJCocHsOrwhFxqlSB7GxAPQ4vkpcH995rODp5Lf+p9zVc85zbJikRRIXDA3h1HoerHkeVKkVCVepxeItJk2D1auEd7iDhvlH2+1JiBhUOj6Eeh0NSUoHHkZBgS1ko3iAvD4YPNxxT5Q/Oa7oMBk5x2yQlwnjoMqWAd4XDbY9DQ1Xe4cknYc8eYVT2KBIeus9bdXKUiKDfuAeIhuR4xPHzODRU5R3S02HsWEPrxNVc1Pp3uPBCt01SXECFw2NojsOhShUbE8nLIyEhXj0OjzB5MmzdKrzBUHjoIW/d6SgRQ791D6AeRwB8ydb9+/1H5ioukpkJDz1k6J30A32O3QNnnum2SYpLqMfhMbzUwMbf+3FlHgdAdjZJSVVVODzAmDGwZYvwDiORRx721o9ViSgqHB7A///PS/+Lrtri53H4TelQXGLbNnjsMUO/xK844fhEOOkkt01SXESFQwkKNz0OFQ73efBByM4W7uEe5JFn3DZHcRkVDg/gJS/DM6jH4Rn++gteeMHQP346x15wOHTr5rZJisuocChB4cqoKlCPwwM8+STk5RqejL8dHv7MbXMUD+ChMTyxi3ocAfCFqhyPIzfXdppTIsvatfDyy4ZrzEscccMZ0Ly52yYpHkA9DqVU4uLsBdttj8NZpWrVCNsR49xwA8Tl53JX9Wfhnrlum6N4hLB6HCLSV0R+E5HVIjIywOu3iMhKEVkqIl+LyGHhtMereNnjcG1CYrHkOOhcjkgzezZ89hnckf8oTe4eCHXrhv2cGzfC6tVhP41SQcLmcYhIPPAccBqQBvwsItONMSv9dlsMdDHGZIjIdcDjwMXhskkpO77ue655HE6oCjTPEUny8+Geewz1EnYxvMF7cNNPETlvw4YROY1SQcLpcXQDVhtj1hhj9gNTgXP9dzDGzDbGZDhPfwCahNEez6IeRwD81MLP+VAixMcfw/z5wp25D1Ln0ds0RqgUIZzC0RhY7/c8zdlWEoOBz8Noj1IOfMLh2jwO9TgiTm4u3HZrPkckrOPGdt/C5Ze7bZLiMTyRHBeRK4AuwAklvD4EGAJw6KGHRtCyyOBlj8NXr8oryXEl/EycCH+sjmMaN5H4xCPeqrypeIJwehz/AE39njdxthVBRE4F7gbOMcYEvDQYY142xnQxxnSpX79+WIxVAuN6clw9joiSng5335XPCfFzOef0/dCnj9smKR4knMLxM9BCRJqJSBJwCTDdfwcR6Qi8hBWNLWG0xdOoxxEA9Thc4d57YeeuOB5mFDL+abfNUTxK2ITDGJMLDAVmAKuAd40xK0TkARE5x9ntCaAa8J6I/CIi00s4nOISXkiOq3BEhl9/hQkTDBczleOGd4WWLd02SfEoYc1xGGM+Az4rtm203/qp4Tx/tKAeRwAChKp0Hkd4eexRQ1x+LmPrPQajv3XbHMXDaMkRpVTU44gN5s6FSZOFW8xYGj8xHGrUcNskxcN4YlRVrKMeRwB8HkdWlgpHmMnPh5uH5dE0bhP3dPwcrpzttkmKx1HhUErFNY9DxE46y8xU4QgzTz8NCxfHM5mRpD73eET6F+fk5JCWlkZWVlbYzxVrJCcn06RJExITE8N2DhUOpVRc8zgAUlIgI0NnjoeRtDS4/948TpS5XDEgHrp3j9B506hevTqHH3444mWXO8owxrB9+3bS0tJo1qxZ2M6jOQ4P4OX/G1fnfjnCoR5H+LjrLsOevcKE1JHImMcidt6srCzq1q2rohFiRIS6deuG3ZNT4VBKxQsehwpHeFi0CN54Q7iZp2jz5FURrzCoohEeIvG5aqjKA3j5/8dLHocOxw0dublw9VW5/Eu2M6rrl3C1dvZTgkeFQykV14ocgnocYeTpp2HR0gTejRtG7VefjEhCXKk86K/FA6jHUQKOcMTHQ2IiZGa6aEslYs0aGD0qj7OZzgV3NIe2bd02yTXS09Pp378/nTt35phjjuGVV14p13G++OILWrZsyZFHHsljj0UuV+QW6nEopeJqjqNqVdhiS5ilpsK+fS7YUMkwBq4dkkdCTibPH/Y4cs+XbpvkKu+//z7Vq1dn4cKFAGSW4+4kLy+PG264gS+//JImTZrQtWtXzjnnHFq3bh1qcz2DehweQD2OEnA8DlDhCBVvvAFffh3Po/l30GTivTHfoKlTp0588803dOnShXvvvZcqvrhoGfjpp5848sgjOeKII0hKSuKSSy7ho48+CoO13kE9DqVUvDCqClQ4QsHWrXDLTTn04Geu+79cOO00t02yDB8Ov/wS2mN26GATOaWQnp7O7bffztKlS0lNTeXkk0+mQ4cO9O/fv2Cf3r17s2fPngPe++STT3LqqbbU3j///EPTpoUdJJo0acKPP/4Ykj/Dq6hweAD1OEogJaUgsaHCUXFuvimX3emGVxqNJm7sB26b4zovvfQSp59+OjVr1gSgR48ebNq0qcg+c+fOdcM0z6PCoZSKlzyOvXtdsKGS8Mkn8NbUBEZzP60nj/RWEcODeAbhYvHixQwcOLDI87PPPrvIPsF4HI0bN2b9+sIu2WlpaTRuXFqX7OhHhcMDqMdRAikpkJUF+fmkpMSxa5eLtkQxW7bAlZflcAyruOv/tsKp2s0AoHbt2ixevJi+ffvy6aefsnv3bnr27Flkn2A8jq5du/LHH3/w119/0bhxY6ZOncqUKVPCZbYnUOFQSsV1jwMgM5PU1FT+OaDxsBIM1/1fDvv25DO50UiqjHvHbXM8w2233cbFF1/M1KlTadasGR988AFx5ZjPkpCQwIQJEzj99NPJy8tj0KBBtGnTJgwWewcVDg+gHkcJ+Eb8ZGSQmpqqOY5y8PYUwwcfJ/KgjKbD+/dA9epum+QZmjVrxk8//RSSY/Xr149+/fqF5FjRgAqHUio+j8MYF07u8zgyMjQ5Xg7++AOGDMrlOH5k5AMp0KOH2yYplQQVDg8QDR5Hbq4LJ09Ntcu9e1U4ykhWFlx8bhZJ2ft4u+cEEu58y22TlEqECodSKj6PIy/PhZM7wyTZvZvUVDvAKj9fyyoFw8hbc1m8Kpnp1QfT9L1xLscclcqG/gt6APU4SsAnHLt2FTgfWq/q4HzwvmH8cwlcz3Oc/fZlcMghbpukVDJUOJRScdXjqFXLLtPTC4RDw1Wls3o1DLw8h3Ys4ak7t8KZZ7ptklIJ0VCVB1CPowQCeBwqHCWzZw9c0G8fidnZfHD8eJIenOi2SUolRT0OpVQ84XHs2kW1anZVZ48HxhgYdFkmy/5I5s2Gt9F82ljNayhhQz0OD6AeRwlUrWobcaSnU7u23bRzpwt2RAFjH8nmf59U5fHEu+k3YxgFH5iihAH1OJRScdXjELHhql27qFPHbtqxwwU7PM577+Rxx6gELuB/3PpBT2jXzm2TooZQNHJav349J510Eq1bt6ZNmzaMHz++4LXK2uBJPQ4PoB5HKdSqBenp1K1rn6pwFGXeXMOAy/PpxgJeH7cDOWuI2yZFFaFo5JSQkMDYsWPp1KkTe/bsoXPnzpx22mm0bNmy0jZ4Uo9DKRVXPQ6wwuHncWzf7pIdHmT5cjjztGwOyVvPR9d+QerNKhplJRSNnBo1akSnTp0AqF69Oq1ateKff/4pU4On8847j1GjRnH88cdz6KGH8tVXXwFw4YUXMnToUHr16sVhhx3GvHnzGDBgAEcddRSDBw8u/x9eQdTj8BBezGX6bHJNOJxQVbVqNt2hHodl3To4s/duUrL3MefC52nw3ONum1RuXOrjFLJGTv6sXbuWxYsX0717d2bOnBl0g6dly5bRs2dPvv32Wz788EPeeustTj31VJYtW0aPHj2YMGECjzzyCIMHD2bOnDnUr1+fJk2akJ2dXS6xqygqHB7AF6ry4oxo14Wjdm1IS0ME6tRRjwMgLQ1O7pLO7l35fH3aeA59e4w3fzweJ9SNnPbu3cv555/P008/TY0y9DvJyMggPT2dm2++GYCcnBxq1apFVlYWu3btYvjw4QCICIMHD6ZRo0YAxMfHk5SUFPR5QokKh4fwssfhWo6jUSP48ksA6tZVj2PLFjityw62bkvg614P0umTx7z5wykDLvVxClkjJ7AX+/PPP5/LL7+c8847Dwi+wdPKlSvp3Lkz8c73uHTpUtq2bcuKFSvo1KlTQan3JUuWcN111xUc65BDDkFcSpCqcHgA33fvxf9/13McjRpBejrs20edOqkxLRzr1hrO6LaNdVtTmdH7Ybp++Si4dMdZGQhVIydjDIMHD6ZVq1bccsstBduDbfC0bNkyOnToUPB86dKlnHvuuSxbtoz27dsX2d7OGTG3ZMmSgnU3UP/WA/hKlntROFz3OHx1ljZupG7d2A1VLVqQz7Ft97BhayKf93mK3rPuBxdi25WJ2267jQ8//JD27dszceLEcjdymj9/Pm+88QazZs2iQ4cOdOjQgc8++6xIg6dWrVpx0UUXBWzwVFw4li9fTtu2bYtsz8rKIjMzk9rO/Bx/EXED9Tg8QH6+XXoxTO26x+EnHA0bHsn8+S7Z4SKff5jFhRca6uTt4qtLJtPmrbu8+WOJMkLVyKlXr16YEhrWBNPgady4cUWer1mzBoCxY8cWbEtOTuavv/4qeH7nnXeW19yQoL8+D+ATDvU4AuATjg0bOPxw2LYtdsqOGAPPP7yTs89L4Ki8Vfxw3wzaTLlbRUNxHfU4PIDvbt6L1wPXPQ5nBAkbNnDYYXZ13Tqo5C2d2bMHrjkrjbe/bUK/hBlMnSpUP/9qt81SFEA9Dk/gK+A3YIC7dgTCdY+jdm3bJ/vPPzn8cLtp3TqXbIkQP32bRaemW3nn20Y81HACHy9rRvXz+7htlqIUoB6HB0hNtQOHfALiJVz3OESgdWtYsaLA41i71iVbwsyePfDE9X/x6JtNOIQM5pz3Fr2nXFdpk+DGGNeGk1ZmSsq3hBIVDo9QhvlCEeWss6BbNxg92kUj2raF6dNp2BCSk22zolLJzbXJkC1bbAOPqlWt59K4caESeoicHJj46DbufySRLdnNuKzadJ57uw61zhrutmlhIzk5me3bt1O3bl0VjxBijGH79u0kJyeH9Tze+y9SPEXNmlBClYTI0aYNvPoqcVs307btv1iyJMA+O3fCpEnw+ecwbx6/ZzRmPsfxBy3YREOqkE092UG9OvnUa1yFhu0a0LbPITTo2wmpXy/ifxLY5Pf7L2/nrpF5/LGrAcfHzWX6oFV0nzDAil0lpkmTJqSlpbF161a3Tal0JCcn06RJk7CeI6zCISJ9gfFAPPCKMeaxYq9XASYDnYHtwMXGmLXhtEmJQrp1s8s5c+jW7WImTYKsLOt9sG0bPPAAvPoqWzNSmFDvft6tMolfMxoCkBCfz79qZrM/B7bvrUL+9jj7S1sKvAmN2EC36l9zbJs99Dy9Bt2uakPyYf8K65+zZ7dh+lN/8swE4adtzWnDcj45dQr9XrsAado7rOf2ComJiTRr1sxtM5RyIuGKh4lIPPA7cBqQBvwMXGqMWem3z/VAO2PMtSJyCdDfGHNxacft0qWLWbBgQVhsVjxKXp4dXXXqqXxx5RTOOAM+/Tiffpteg9tvJyM9h3HtXmfM7/9mX2Y8J50E/fvDaadB8+aF0an8fOuYbNsGa3/NYvFnG1n4XRY/ra7D31lWLOLJpUvycnoe9g/Hds2ne7+6HHZG68JuhOUkKyOfz55fy9TXMvjk1+ZkmqocJuu496S5XDmxN/FHHFbBD0lRSkdEFhpjuoTkWGEUjh7AfcaY053ndwIYYx7122eGs8/3IpIAbALqm1KMUuGIUa6/Hl55hewfFnPIyS1pbVby+e6efHT0SO5Kv52/NybRvz88/DC0alX2w29Yl8N3U9Yy/7N0vlteg0W7mpFLIgC12Em75N9p12QH7dsZmrSvxxFtU6jTvDZ1j6qLJFcBEYyBPTty2PzbLjav2sHaX3ay8GfDgj9qsGjH4WSQSn22cNEh87jk8gR63nUicbU8mtxSKh3RIhwXAH2NMf/nPB8AdDfGDPXbZ7mzT5rz/E9nn20lHVeFI0bZvNnmOrZv53UGchWvF7zUvj2MHw8nnBC602VmwuJZO/np4838ujiTJatTWbazMftMapH9hHwSyUEw5BNHDkVrR1Ulg47Jv9L1iO306wcnj+hIQkN3cipKbBNzwiEiQwBfl5q2wPKwGB0Z6gElCmMUoPa7RzTbDmq/27Q0xlQPxYHCmRz/B2jq97yJsy3QPmlOqKomNnVZBGPMy8DLACKyIFSq6QZqv7tEs/3RbDuo/W4jIiEL1YRz5vjPQAsRaSYiScAlwPRi+0wHfAXxLwBmlZbfUBRFUdwnbB6HMSZXRIYCM7DDcV8zxqwQkQeABcaY6cCrwBsishrYgRUXRVEUxcOEdR6HMeYz4LNi20b7rWcBF5bxsC+HwDQ3UfvdJZrtj2bbQe13m5DZH7bkuKIoilI50eq4iqIoSpnwhHCIyGsissUZnuvbVkdEvhSRP5xlbWf7iSKSLiK/OI/Rfu/pKyK/ichqERnpou0XisgKEckXkS7F9r/Tse83ETndTdvLar+IHC4imX6f/Yt+r3UWkWWO/c9IhCrXlWD/EyLyq4gsFZEPRaSW32vR8PkHtD+KPv8HHdt/EZGZInKIs10c21Y7r3fye89A53/9DxEZGOhcHrDf89cev9dGiIgRkXrO89B+9sYY1x/A8UAnYLnftseBkc76SGCMs34i8EmAY8QDfwJHAEnAEqC1S7a3AloCc4AufttbO3ZVAZo59sa7ZXs57D/cf79ix/kJOBYQ4HPgDBft7wMkOOtj/H470fL5l2R/tHz+NfzWbwJedNb7ObaJY+uPzvY6wBpnWdtZr+1B+z1/7XG2N8UOSloH1AvHZ+8Jj8MY8y12VJU/5wKTnPVJwL8PcphuwGpjzBpjzH5gqnOMsBLIdmPMKmPMbwF2PxeYaozJNsb8Bax27HbFdsfWstgfEBFphP1n+8HYX+NkDv59hYQS7J9pjPG1nvoBO4cIoufzL8n+gHjw89/t9zQV8CVSzwUmG8sPQC3H9tOBL40xO4wxO4Evgb7ht77M9peEZ649Dk8Bt1PU7pB+9p4QjhL4lzFmo7O+CfAvWdpDRJaIyOci4msi2hhY77dPmrPNS5RkYzTY7qOZiCwWkW9ExFfKtTHWZh9esn8Q9k4LovPz97cfouTzF5GHRWQ9cDngC+lEzedfgv3g8WuPiJwL/GOMKd58IKSfvZeFowDnLsqnnouAw4wx7YFngWlu2RWDbAQONcZ0BG4BpoiIZ6v0icjdQC7wltu2lIcA9kfN52+MudsY0xRr+9CD7e81SrDf09ceEUkB7qKo0IUFLwvHZseV8rniW8C6kcaYvc76Z0CikwAKpsSJ25RkYzTYjhPi2e6sL8TGdY/C2uofTnHdfhH5D3AWcLlz4wFR9PkHsj+aPn8/3gLOd9aj5vP3o8D+KLj2NMfm7paIyFrHjkUi0pAQf/ZeFg7/ciQDgY8ARKShb8SIiHTD/g3bCa7EidtMBy4RkSoi0gxogU1qRoPtiEh9sX1WEJEjsPavcUKKu0XkWOe7uRLn+3LJzr7YGO85xpgMv5ei4vMvyf4o+vxb+D09F/jVWZ8OXOmM8DkWSHdsnwH0EZHaYkdP9nG2uUJJ9nv92mOMWWaMaWCMOdwYczg27NTJGLOJUH/24c78B/MA3sa64TnOHzsYqAt8DfwBfAXUcfYdCqzAjlz4Aejpd5x+2OZRfwJ3u2h7f2c9G9gMzPDb/27Hvt/wG/nihu1ltR9757UC+AXrtp/td5wu2KrFfwITcCaXumT/amzc9hfn8WKUff4B7Y+iz/99x5alwMdAY2dfAZ5zbFxG0RF7g5y/ezVwlcuff0n2e/7aU+z1tRSOqgrpZ68zxxVFUZQy4eVQlaIoiuJBVDgURVGUMqHCoSiKopQJFQ5FURSlTKhwKIqiKGVChUNRFEUpEyociqIoSplQ4VAUP0Rkb4iPV0tErvd7fnig/gklvLeqU8wwPgR2JInItyIS1nbRSmygwqEo4aUWcP3BdiqBQcAHxpi8ihphbLnvr4GLK3osRVHhUDyPiNwmIjc560+JyCxn/WQRectZnyYiC8V2LhzibHtMRG7wO859InKrs36FiPwktpPbS4Hu6gPt43gMq0RkonOumSJS1dn/HrFd4OaJyNvOuR4DmjvHeMI5dHyg9wfgcvxqTonIByLykOM5/C0ipzrb3xORCc5514lILxF5Q0R+F5FX/Y43zTmmolQIFQ4lGpgL+HpPdAGqiUiis+1bZ/sgY0xn5/WbRKQu8A5wkd9xLgLeEZFW2Dvv44wxHYA8il1QD7JPC+A5Y0wbYBdwvoh0xdaSag+c4dgBtnvln8aYDsaY20p6f/E/2CmWd4QxZq3f5mOAXcaY44FhfvYcgy122At4CXgVWySxNXCmiFRx9lsOdC1+LkUpKxrvVKKBhUBnsb0nsrEF/rpgheMmZ5+bRKS/s94UaGGM+UFEGojtGV0f2GmMWS8iQ4HOwM9OsdOqOGX7/TilhH2+Bf4yxvziZ9vhQD3gI2NMFpAlIh+X8vcEen9x6mFFBSjotVAT290NIBHYJSLJ2HDY0852A7xqnCZoIpIH7AcwxuSJyH4RqW6M2VOKfYpSKiociucxxuSIyF/Af4DvsFVLTwKOBFaJyInAqUAPY0yGiMwBkp23vwdcADTEeiBgK4VOMsbcWcppA+4jIodjxctHHlZUykIw78+k8G8A6z0s9Mt3tMN6EG2ARcaYfGd7e+AFx9YmwAZTtJJpFSCrjPYqShE0VKVEC3OBW7F3/HOBa4HFzkWxJtabyBCRo4Fj/d73DrY/wgVYEQGbJL5ARBoAiEgdETms2PmC2cef+cDZIpIsItWwTZgA9gDVy/rHGtv/Od7xKMCGo37x26UdVkCPwZb5Lr4drIj41nHCd9uMMTlltUdR/FHhUKKFuUAj4HtjzGbsXfNc57UvgAQRWYVNRv/ge5MxZgX2wv2PL3xjjFkJjAJmishS4Evn2Pi976D7FNv/Z2yznKXYHuHLsM1ytgPzRWS5X3I8WGYCvZz14sLRFutxFGx3RKaqIzpQVETAemmfltEGRTkA7cehKCFCRKoZY/Y6+YhvgSHGmEUVOF4n4GZjzIAQ2fcBMNIY83sojqfELprjUJTQ8bKItMbmJiZVRDQAjDGLRGS2iMRXdC6HM0prmoqGEgrU41AURVHKhOY4FEVRlDKhwqEoiqKUCRUORVEUpUyocCiKoihlQoVDURRFKRMqHIqiKEqZUOFQFEVRyoQKh6IoilIm/h+yFiibEr57lQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x324 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot transmission, compare to paper results, look similar\n",
    "fig, ax = plt.subplots(1, 1, figsize=(6, 4.5))\n",
    "plt.plot(wavelengths / nm, flux_delta0 / flux_empty, color='red', label='$\\delta=0$')\n",
    "plt.plot(wavelengths / nm, flux_delta20 / flux_empty, color='blue', label='$\\delta=20~nm$')\n",
    "plt.xlabel('wavelength ($nm$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.xlim([1050, 1400])\n",
    "plt.ylim([0, 1])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results Comparison\n",
    "\n",
    "Compare this plot to published results computed using COMSOL (FEM):\n",
    "\n",
    "<img src=\"img/Si_plot.png\" alt=\"diagram\" width=\"400\"/>\n",
    "\n",
    "(Citation: Opt. Lett. 43, 1842-1845 (2018).  With permission from the Optical Society)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}