{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transmittance Spectrum of a Waveguide Bend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have computed the field patterns for light propagating around a waveguide bend. While this can be visually informative, the results are not quantitatively satisfying. We'd like to know exactly how much power makes it around the bend (transmittance), how much is reflected (reflectance), and how much is radiated away (scattered loss). How can we do this?\n",
    "\n",
    "The basic principles are described in Introduction. The computation involves keeping track of the fields and their Fourier transform in a certain region, and from this computing the flux of electromagnetic energy as a function of $\\omega$. Moreover, we'll get an entire spectrum of the transmittance in a single run, by Fourier-transforming the response to a short pulse. However, in order to normalize the transmitted flux by the incident power to obtain the transmittance, we'll have to do _two_ runs, one with and one without the bend (i.e., a straight waveguide)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using MPI version 3.1, 1 processes\n"
     ]
    }
   ],
   "source": [
    "import meep as mp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "resolution = 10  # pixels/um\n",
    "\n",
    "sx = 16  # size of cell in X direction\n",
    "sy = 32  # size of cell in Y direction\n",
    "cell = mp.Vector3(sx, sy, 0)\n",
    "\n",
    "dpml = 1.0\n",
    "pml_layers = [mp.PML(dpml)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also define a couple of parameters to set the width of the waveguide and the \"padding\" between it and the edge of the cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pad = 4  # padding distance between waveguide and cell edge\n",
    "w = 1  # width of waveguide"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to define the waveguide positions, we will have to use arithmetic to define the horizontal and vertical waveguide centers as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "wvg_xcen = 0.5 * (sx - w - 2 * pad)  # x center of horiz. wvg\n",
    "wvg_ycen = -0.5 * (sy - w - 2 * pad)  # y center of vert. wvg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We proceed to define the geometry. We have to do two simulations with different geometries: the bend, and also a straight waveguide for normalization. We will first set up the straight waveguide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "geometry = [\n",
    "    mp.Block(\n",
    "        size=mp.Vector3(mp.inf, w, mp.inf),\n",
    "        center=mp.Vector3(0, wvg_ycen, 0),\n",
    "        material=mp.Medium(epsilon=12),\n",
    "    )\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The source is a `GaussianSource` instead of a `ContinuousSrc`, parameterized by a center frequency and a frequency width (the width of the Gaussian spectrum)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcen = 0.15  # pulse center frequency\n",
    "df = 0.1  # pulse width (in frequency)\n",
    "sources = [\n",
    "    mp.Source(\n",
    "        mp.GaussianSource(fcen, fwidth=df),\n",
    "        component=mp.Ez,\n",
    "        center=mp.Vector3(-0.5 * sx + dpml, wvg_ycen, 0),\n",
    "        size=mp.Vector3(0, w, 0),\n",
    "    )\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how we're using our parameters like `wvg_ycen` and `w`: if we change the dimensions, everything will shift automatically.\n",
    "\n",
    "Finally, we have to specify where we want Meep to compute the flux spectra, and at what frequencies. This must be done after specifying the `Simulation` object which contains the geometry, sources, resolution, etcetera, because all of the field parameters are initialized when flux planes are created. As described in Introduction, the flux is the integral of the Poynting vector over the specified `FluxRegion`. It only integrates one component of the Poynting vector and the `direction` property specifies which component. In this example, since the `FluxRegion` is a line, the `direction` is its normal by default which therefore does not need to be explicitly defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (0,-11.5,0)\n",
      "          size (1e+20,1,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAIhCAYAAAD6ovlZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df5QlZX3n8fd3+DnTw/xkIMgMNI1KGNQoxECIrhN/0QeXk3h0E4+ziK5ms4mJmKjZJGqWgIkQ1CQqcaOrgGTc/FAx4UjaH8dMokYQBUUnngUZGroRmRFnuqHHcWTm2T/ubWia7pnue+upW1X3/TqnT3mrbj31rfHeD09XPU91pJSQJBVrSa8LkKQmMlwlKQPDVZIyMFwlKQPDVZIyMFwlKQPDVZIyMFwlKQPDVZIyMFwlKQPDVZIyMFwlKYPDe12Anigivg8sA8Z6XYvUxzYAe1JKP9XJzuFTsaonIiY5jGNY0+tKpD72Q1i+dDkPPfRQdLK7PddqGmMNG7/97W+zbtm6Be0w9ZMpht47BMADb3qAgSMHOjrw+OQ4w1uGGd09yuCqQUY2j7B+xfqO2uolz6NaenUeU/umOP7dxwOw/Q3bGThiYd+LnXt28rSnPY2T1p3U8bEN1wpbt2wdxy0/bkHvndo39ej/HjhyoKNwHZsY4/yPnc/o7lGGVg+x9aKtbFi5YdHt9JrnUS1VOY/jBo7ruNPRCW9oCWh9ATZdu4ntu7bX/ovseVRHU86jE4arGvMF8DyqpSnn0SnDtc815QvgeVRLU86jG4ZrH2vKF8DzqJamnEe3+jpcI+KsiPj9iPhkRIxHRIqIecemRcQl0++Z5+fyMuvvRlO+AJ5HtTTlPL7+va933Ua/jxZ4O/BLHez3ZeC7c6zv/v+REjTlC+B5VEtTzuPm8Zu54G8v6Lqdfg/XrwC3A7e0f0aBoxaw3/9JKV2Tr6x8mvIF8DyqpSnncfP4zTzn6uewP+3vuq2+DteU0hUzX0d0NBGjNpryBfA8qqUp5zEdrI8ceITD4jD2013A9vU1137SlC+A51EtTTmPmcF6+JLDueEVN3TdZl/3XLvw/Ih4JnA0MA78c0qp8OutUz+ZetzMq4O+d8b7Zu8ze+rhja+8kTVL1yy47arwPKqlDudxsO/FtFvuu4UXXvdC9qf9HBaH8bn/+jlOXnVy18f2wS0zRMRe4KiU0pzXByLiEuB/zbP7J4BXp5QeXsTxts2z6VTWcRSvX2hLkgp3FWxct5Ft27Z1dL3QywKL813gzcAZwHJajyTbDNwHvAy4rnelSaoSe64zHKrnepD9TgC+BawFfj6ldFOXdWxjHRsf2L7wp1vNfPpPN0/Fkppkvu/FoZ7SNbVviuOHju+q5+o11wKklO6PiKtp9WqHga7CdVqnT7fqdD+pyaa/F2U9pcvLAsW5s708oadVSJpXmaMbDNfirG4vq3OrVNKjxifHSx025mWBAkRr9sFL2y9v7UUN45PjvTisVBvT11jLGo9rz3WBImJdRLw+Io6ZtX458AHgbOD7wCfLrm1sYozhLcNlH1aqlbL/EkJf91wj4iW0Ht4y7cj2+pk3pC5LKX0aGADeD1weEbcA9wPrgDNpjRLYDbw8pbSnjNqnTV9DGt09WuZhpdoZXDVY6gyyvg5XWuF49hzrz571HoAHgSuAc4CnAucC+4G7gWuAP08p3Zet0jnMvDg/uGrQgJUOYmTzSKlTc/s6XNtPtrpmge99CPj9nPUsxuy7nje+8kZ++qqf7nVZUmWV/Vdz+zpc62qu4SRrlq7pdVmSZvCGVs005SlEUtMZrjVisEr1YbjWhMEq1YvhWgMGq1Q/hmvFGaxSPRmuFVb2XGhJxXEoVoWVPRdaUnHsuVaYwSrVl+FaYWXPhZZUHMO1wsqeCy2pOIZrhZU9F1pScQxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkWcYnx7tuw3CVpBnGJsYY3jLcdTuGqyS1jU2MsenaTYzuHu26LcNVkngsWLfv2s7gqsGu2zNcJfW9mcE6tHqIkc0jXbd5eAF1SVJtzQ7WrRdtZc3SNV23a89VUt+aK1g3rNxQSNuGq6S+lDNYwXCV1IdyBysYrpL6TBnBCoarpD5SVrCC4SqpT4xPjpcWrOBQrMYoYi601GTDW4YZ3T1aSrCCPddGKGoutNRkZQYrGK61V+RcaKnJBlcNlhasYLjWWtFzoaUmG9k8UlqwguFaWznmQktNtn7F+lKP5w2tGso1F1pScey51kyZ4/Qkdc5wrRGDVaoPw7UmDFapXvo6XCPirIj4/Yj4ZESMR0SKiLSA/V4dEV+NiIcj4ocRcWNEnJurToNVqp9+v6H1duCXFrNDRPwFcDHwI+CzwNHAi4AXR8TLU0qfKrJAg1Wqp34P168AtwO3tH9GgaPme3NEvJBWsD4I/HxK6c72+p8HtgJXR8TWlNLuIoobnxzn/I+db7BKNdTX4ZpSumLm64g41C6/216+YzpY2+18JSL+N/AG4LXAu4uor+y50JKK09fXXBcjIpYCz2+//Pgcb5led0FRxzRYpfoyXBfuNFqXDHamlOZ6BNWt7eUzijpg2XOhJRWnry8LLNJJ7eWcz/ZLKU1FxG5gdUQck1J66FANRsS2eTadCuXPhZZUHHuuC7e8vdxzkPdMtZfHFHHAsudCSyqOPdceSimdMdf6do92Y8nlSCqQPdeFe7i9XHaQ9wy0l4e8JCCp2QzXhbu3vZzzd/WIGABWAbsWcr1VUrMZrgv3/4AfA+si4sQ5tp/ZXt5eXkmSqspwXaCU0o+AL7Rf/pc53vLy9vKGciqSVGWG6+K8p718W0Q8ZXple/rrrwO7gQ/3ojBJ1dLXowUi4iW0Ht4y7cj2+ptmrLsspfRpgJTS5yPiL2k9X+AbEfG59j4vAgJ4TVHPFZBUb30drsA64Ow51p896z2PSim9MSK+AfwWrVDdB3yeVgj/e65CJdVLX4drSuka4Jqy9pPUP7zmKkkZGK6SlIHhKkkZGK6SlIHhKkkZGK6SlIHhKkkZGK6SlIHhKkkZGK6SlIHhKkkZGK6SlIHhKkmzjE+Od92G4SpJM4xNjDG8ZbjrdgxXSWobmxhj07WbGN092nVbhqsk8Viwbt+1ncFVg123Z7hK6nszg3Vo9RAjm0e6brOv/xKBJM0O1q0XbWXN0jVdt2vPVVLfmitYN6zcUEjbhqukvpQzWMFwldSHcgcrGK6S+kwZwQqGq6Q+UlawguEqqU+MT46XFqzgUKzGKGIutNRkw1uGGd09Wkqwgj3XRihqLrTUZGUGKxiutVfkXGipyQZXDZYWrGC41lrRc6GlJhvZPFJasILhWls55kJLTbZ+xfpSj+cNrRrKNRdaUnHsudZMmeP0JHXOcK0Rg1WqD8O1JgxWqV4M1xowWKX6MVwrzmCV6slwrbCy50JLKo5DsSqs7LnQkopjz7XCDFapvgzXCit7LrSk4hiuFVb2XGhJxTFcK6zsudCSimO4SlIGhqskZWC4SlIGhqskZWC4SlIGhqskZWC4SlIGhqskZWC4SlIGhqskZWC4SlIGhqskZWC4SlIGhmsHImJrRKSD/Az3ukZJveWfeenOJ4CH51h/X9mFSKoWw7U7b04pjfa6CEnV42UBScrAcJWkDLws0J3XRsRa4ABwB/CplNK9Pa5JUgUYrt1526zX74qIy1JKly1k54jYNs+mU7srS1KveVmgM/8GXEgrBJcBpwFvBR4BLo2Ii3tYm6QKsOfagZTSH81adQfwpxHxNeAzwCUR8cGU0o8O0c4Zc61v92g3FlKspEUbnxzvug17rgVKKX0W+BqwCji7x+VI6sDYxBjDW7qfB2S4Fu/O9vKEnlYhadHGJsbYdO0mRnePdt2W4Vq81e3lVE+rkLQo08G6fdd2BlcNdt2e4VqgiFgHPLf98tZe1iJp4WYG69DqIUY2j3TdpuG6SBFxbkT8ckQcNmv9IHA9MAD8U0qp+yvikrKbHaxbL9rK+hXru27X0QKL91TgauD7EXErsBs4GTgLOBrYBvxa78qTtFBzBeuGlRuY2tf9VT3DdfFuBj5AazTAs2ldY50CvgH8A/CBQw3BktR78wVrUQzXRUopfQf4zV7XIalzuYMVvOYqqc+UEaxguErqI2UFKxiukvrE+OR4acEKXnNtjCLmQktNNrxlmNHdo6UEK9hzbYSi5kJLTVZmsILhWntFzoWWmmxw1WBpwQqGa60VPRdaarKRzSOlBSsYrrWVYy601GRFTGldDG9o1dBcw0nWLF3T67IkzWDPtWbKHKcnqXOGa40YrFJ9GK41YbBK9WK41oDBKtWP4VpxBqtUT4ZrhZU9F1pScRyKVWFlz4WWVBx7rhVmsEr1ZbhWWNlzoSUVx3CtsLLnQksqjuFaYWXPhZZUHMNVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCVplvHJ8a7bMFwlaYaxiTGGtwx33Y7hKkltYxNjbLp2E6O7R7tuy3CVJB4L1u27tjO4arDr9gxXSX1vZrAOrR5iZPNI120eXkBdklRbs4N160VbWbN0Tdft2nOV1LfmCtYNKzcU0rbhKqkv5QxWMFwl9aHcwQqGq6Q+U0awguEqqY+UFaxguErqE+OT46UFKzgUqzGKmAstNdnwlmFGd4+WEqxgz7UjEbE0Ii6NiDsiYm9EfC8iPhIRJ/ainqLmQktNVmawQhfhGhF/HhHLiiymDiLiaOALwNuB5cA/AmPAa4DbImKozHqKnAstNdngqsHSghW667leDHwrIl5YVDE18TbgHOArwFNTSr+aUjobeBOwDvhIWYUUPRdaarKRzSOlBSt0F65/C5wCfCYiro6I1QXVVFkRcSTwW+2Xr08pPTy9LaX0HuB24HkRcVbuWnLMhZaabP2K9aUer+NwTSm9ErgAuA+4CPiPiPiVogqrqF8AVgJ3pZRum2P7x9vLC3IWMddwkrI/OJIOrqsbWimlTwOnA1fR+pX4/0bEP/bqxk4Jfqa9vHWe7dPrn5GrgDLH6UnqXNdDsVJKU8BvR8QW4MO0em3Pi4i/BqYOst+l3R67B05qL+cb9zS9/uSFNBYR2+bZdOpcKw1WqT4KG+eaUropIp4F/Bvwc8Cb53lrAAmoY7guby/3zLN9+j8mxxR9YINVqpfCwrU9BOlDwLOB/cD1HKTnKkgpnTHX+naPduP0a4NVqp+uwzUigtYwpEuAZcA3gdellL7ebdsVND06YL7xvQPt5UNFHdBgleqpqxtaEfF04GbgCuAw4K3AzzY0WAHubS/nuzU/vf6eIg5W9lxoScXpuOcaEe8A3gIcAXwR+LWU0h1FFVZR32wvz5xn+/T624s4WNlzoSUVp5ue6x8CPwJ+I6X0vD4IVoAvAxPAqRHxzDm2v7y9vKGIgxmsUn11E643ABtTSn9dVDFVl1LaB7y//fKqiJi+xkpE/C6t8a3/WtRlkbLnQksqTseXBVJKv1RkITXyDuCFwLnAnRHxRVrjWs8GdgL/ragDlT0XWlJxfOTgIqWU9gK/CFxGa7zrL9MK12uAM1NK24s6llNapfryYdkdSCn9CPij9o8kPYE9V0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0nKwHCVpAwMV0maZXxyvOs2DFdJmmFsYozhLcNdt2O4SlLb2MQYm67dxOju0a7bMlwliceCdfuu7QyuGuy6PcNVUt+bGaxDq4cY2TzSdZuHF1CXJNXW7GDdetFW1ixd03W79lwl9a25gnXDyg2FtG24SupLOYMVDFdJfSh3sILhKqnPlBGsYLhK6iNlBSsYrpL6xPjkeGnBCg7Faowi5kJLTTa8ZZjR3aOlBCvYc22EouZCS01WZrCC4Vp7Rc6FlppscNVgacEKhmutFT0XWmqykc0jpQUrGK61lWMutNRk61esL/V43tCqoVxzoSUVx55rzZQ5Tk9S5wzXGjFYpfowXGvCYJXqxXCtAYNVqh/DteIMVqmeDNcKK3sutKTiOBSrwsqeCy2pOPZcK8xglerLcK2wsudCSyqO4VphZc+FllQcw7XCyp4LLak4hqskZWC4SlIGhqskZWC4SlIGhusiRMSmiEgH+bmp1zVKqgZnaHXmLuBL86yXJMO1Q19KKb2610VIqi4vC0hSBoarJGXgZYHOPCUi3gmsBX5A6/rrSErpQG/LklQVhmtnzm3/zPStiHhZSunOhTYSEdvm2XRqx5VJqgQvCyzOBHAlcA6tXuta4AXATcDTgc9GxMrelSepKvqq5xoR1wOnL3K3V6WUvgqQUroNuG3W9i9ExHOAfwGeC/wm8M6FNJxSOmOeOrcBGxdZp6QK6atwBU4BTlvkPssO9YaU0v6IuIJWuJ7HAsNVUnP1VbimlJ6Zsfnpa60nZDyGpJrwmmtxVreXUz2tQlIlGK7FeVl7eWtPq5BUCYbrIkTEGyNiw6x1ERG/DvwOkIAP9KQ4SZXSV9dcC/BG4F0RcStwN3A0rSFYpwAHgDeklL7ew/okVYThujjvBl4MnEFrqNQRwP3A3wDvTSnd0sPaJFWI4boIKaX3Ae/rdR2Sqs9rrpI0y/jkeNdtGK6SNMPYxBjDW4a7bsdwlaS2sYkxNl27idHdo123ZbhKEo8F6/Zd2xlcNdh1e4arpL43M1iHVg8xsnmk6zYdLSCpr80O1q0XbWXN0jVdt2vPVVLfmitYN6zccOgdF8BwldSXcgYrGK6S+lDuYAXDVVKfKSNYwXCV1EfKClYwXCX1ifHJ8dKCFRyK1RhFzIWWmmx4yzCju0dLCVaw59oIRc2FlpqszGAFw7X2ipwLLTXZ4KrB0oIVDNdaK3outNRkI5tHSgtWMFxrK8dcaKnJ1q9YX+rxvKFVQ7nmQksqjj3XmilznJ6kzhmuNWKwSvVhuNaEwSrVi+FaAwarVD+Ga8UZrFI9Ga4VVvZcaEnFcShWhZU9F1pScey5VpjBKtWX4VphZc+FllQcw7XCyp4LLak4hmuFlT0XWlJxDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaQMDFdJysBwlaRZxifHu27DcJWkGcYmxhjeMtx1O4arJLWNTYyx6dpNjO4e7botw1WSeCxYt+/azuCqwa7b69twjYiBiLgwIt4XETdHxI8jIkXEJQvYd31EXB0R34uIvRFxR0T8cUQcXULpkgo2M1iHVg8xsnmk6zYPL6CuunoK8NHF7hQRTwa+AhwLfBv4IvCzwB8BL4iIF6SUflxkoZLymR2sWy/aypqla7put297rsBDwIeB/wGcRSscF+IaWsH63pTS01NKvwqcBlwP/ALwB8WXKimHuYJ1w8oNhbTdt+GaUrorpfS6lNJfp5RuBX5yqH0i4udoBegO4PdmtPUI8BvtNt4QEf38G4FUCzmDFfo4XDv0kvbyhtm/+qeUHqB1iWA18JyyC5O0cLmDFQzXxfqZ9vLWebZPr39GCbVI6kAZwQr9fUOrEye1l/NN35hef/JCGouIbfNsOnUxRUlamLKCFey5Ltby9nLPPNun2stjSqhF0iKMT46XFqxQ455rRFwPnL7I3V6VUvpqjno6kVI6Y6717R7txsW0VcRcaKnJhrcMM7p7tJRghRqHK3AKrSFQi7Gsy2M+fIh2BtrLh7o8zqIUNRdaarIygxVqHK4ppWf24LD3As8C1s+zfXr9PeWUU+xcaKnJBlcNlhas4DXXxfpme3nmPNun199eQi0Lmgt94ADs3Pn4nwMHyqhO6q3pzz5Tx8LUsXzsvM9x4jHlBCsYrov16fbygog4auaGiDgeeC6wC/hy7kIWOhf6wQfhuOMe//Pgg7mrk3rvwQfhlPUDcOVOuHIn557+5FI/+4brIrRvhn0ZOA64Ynp9e0bWXwFH0JoWe8jZXt2YazjJ+hXzXamQ1Au1veZahPaIgxPaL5/UXr4uIqbvDt2fUnrprN1eQ+vBLRdHxPOB/wCeDQwB/w68M2fN843Tm9o3deidJZWmr8OV1s2p2QP+T2z/wBw3plJKd0bEs4BLgWHgpbRudF0G/GnOJ2KVOQBaUnf6OlxTSoMd7jdGqwdbGoNVqhevudaAwSrVT1/3XKvuwIED3LPrHp5/3fMfDdYvXPgFTjzmRA7MHE914ADpgR0c277smh7YwYETNsCSJe1hV0ue0K7DsdR03Xz2DxTwBTFcK+x7D3+P8z92/qF7rA8+yPKThtg5/frKIdixA9atY8kcv5ssWbJkzvVSk3Tz2V9SwBfEcK2w8/7mPO6ZuIdTVp3CyCtGWHfUOvbu3fvEN+7dy+w/3rV3717Yu5fW249+wra5mpGapJvP/t59rTellDo+vuFaYfdM3MPKAyt50X0v4oNXfnDe9w3s2cMls9ZdfvnlTC1bxp49AzBr6+WXX86yZQ7dUrN189nfxz4AduzY0fHxo5tkVh4RsY11bIzNwTEfP4YlDx/8V5S1Bw7w3cnJx6178ooVPLhkCQcOrGVy8ruP27ZixZNZssRpWmq2bj776fDEREyw5MEl7N+/Pzo5vj3XCksfTUz+cPKQ75vr/8SJyUl2z7N1cnIC2lul5uris38EsKq7G1ve1qiyUh9cKKlIhqskZWC4SlIGhqskZWC4SlIGhqskZWC4SlIGhqskZWC4SlIGztCqsoCIQ8+8izmmMAetfVOaa/9YULtSnXXz2U/R/WMBDNcKO2PjGSzZf+hfLlY/8gh85zuPW3f66aez6/DDeeSR1bM3cfrpp3P44buKLFWqnG4++wcOO8C2+7Zx1FFHHfK98/HBLRU0/eCWv//g33P0YbMfJvhER05McN6FFz5u3Weuu459K1cyMXEkF1543uO2XXfdZ1i5cl+hNUtV081nf+/+vfzKf/8VNhy9gXvvvbejX/MM1wqaDteHxx9m4MiBQ++wcyccd9zj17Ufln2QTVKjdfPZn9o3xfL1y9m4biPbtm3rKFy9odUEa9cyNX43694C694CU+N3w9q1va5K6mtec62wlNLCnoQeQTr2WH7Q7uSmY48lRUBKtHZ//H9406Prpebq5rNfxG/0hmuFRSz8rv7M983cb67dW9sLKVGqrG4++0WMpvGygCRlYLhKUgaGqyRlYLhKUgaGqyRl4GiBhlu7tjVwevY6qenWroW7x6c45S8HAbj74lHWrl3ApJyCGK4Nt2SJs7HUnx797A/8AGj97yUL/F19fHK8++N33YIkNcjYxBjDW4a7bsdwlaS2sYkxNl27idHdo123ZbhKEo8F6/Zd2xlcNdh1e4arpL43M1iHVg8xsnmk6za9oSWpr80O1q0XbWXN0jVdt2vPVVLfmitYN6zcUEjbhqukvpQzWMFwldSHcgcrGK6S+kwZwQqGq6Q+UlawguEqqU+MT46XFqzgUKzGKGIutNRkw1uGGd09Wkqwgj3XRihqLrTUZGUGKxiutVfkXGipyQZXDZYWrGC41lrRc6GlJhvZPFJasILhWls55kJLTbZ+xfpSj+cNrRrKNRdaUnHsudZMmeP0JHXOcK0Rg1WqD8O1JgxWqV4M1xowWKX6MVwrzmCV6slwrbCy50JLKo5DsSqs7LnQkorTtz3XiBiIiAsj4n0RcXNE/DgiUkRccoj90iF+ji6qRoNVqq9+7rk+Bfhoh/tOAR+fZ9v+Dtt8grLnQksqTj+H60PAh4Fb2j8vAS5d4L4/SCm9OlNdjyp7LrSk4vRtuKaU7gJeN/06Il7cw3LmVPZcaEnF6dtrrpKUU9/2XLs0EBFvBU4C9gC3AZ9MKT1c5EGm9k119N7F7Cc1WaffiyK+Q4ZrZ44F3jFr3Xsi4qKU0qcX2khEbJtn06kAx7/7+I6K63Q/qcnK/l54WWDxPgoMAycCy4FnAdcBa4FPRsSze1ibpIqobc81Iq4HTl/kbq9KKX21m+OmlC6ateobwKsiYgz4Q1o92vMW2NYZc61v92g3bn/Ddo4bOG5BdU3tm3r0v8wPvOkBBo4cWNB+3RifHH90osPgqkFGNo/U8iac51EtRZ5Hp9+LHVM7GLpqqKNjTqttuAKnAKctcp9lOQpp+zPgfwKbIuLIlNK+bhscOGKgo5AcOLKz/RZjbGKM8z92fu0nOnge1ZLzPBbzvRjY1/33p7bhmlJ6Zq9rmCmlNBERO4ATaF0iuL/HJWXTlIfJeB7V0pTzmOY114JExBJgRftlY2/XN+UL4HlUS1POYybDtTjDwABwV0ppstfF5NCUL4DnUS1NOY/ZDNdFiIhXzDUaICKeB3yo/fKqcqsqR1O+AJ5HtTTlPOZS22uuRWiPODih/fJJ7eXrImK4/b/vTym9dMYuw8BFEXEHsA34CfBUYPr6798Cf5m36vI15QvgeVRLU85jPn0drrTGqJ48a92J7R+Ae2Zt+zta/2ZnAb9Ia5zrD4F/Bj6SUprvSVm11ZQvgOdRLU05j4Pp63BNKQ0u8v3/TCtI+0JTvgCeR7U05TwOxWuumlNTvgCeR7U05TwWwnDVEzTlC+B5VEtTzmOh+vqyQNXt3LNzwe+d+kkxT8WaPfXwxlfeyJqla2r3pC3Po1p6dR4z298xtWPBM68W892bT6SUum5ExYqISQ7jGNb0uhKpj/0Qli9dzkMPPRSd7G64VlBEfJ/WcxDGel3LIZzaXt7V0yqqyX+b+dXl32YDsCel9FOd7Gy4qmPTz6Od7+le/cx/m/n1y7+NN7QkKQPDVZIyMFwlKQPDVZIyMFwlKQNHC0hSBvZcJSkDw1WSMjBcJSkDw1WSMjBcJSkDw1WSMuOK/Y8AAAWbSURBVDBcJSkDw1WSMjBcVZiI2BQR6SA/N/W6xtwiYmlEXBoRd0TE3oj4XkR8JCJOPPTezRURWw/x2Rg+dCv14p95UQ53AV+aZ31jRcTRwBeAc4D7gX8EBoHXAP85Is5JKW3vXYWV8Ang4TnW31d2IbkZrsrhSymlV/e6iB54G61g/Qrw4pTSwwAR8bvAu4GPAJt6Vl01vDmlNNrrIsrgZQGpABFxJPBb7Zevnw5WgJTSe4DbgedFxFm9qE/lM1ylYvwCsBK4K6V02xzbP95eXlBeSeolLwsoh6dExDuBtcAPaF1/HUkpHehtWVn9THt56zzbp9c/o4Raquy1EbEWOADcAXwqpXRvj2vKwnBVDue2f2b6VkS8LKV0Zy8KKsFJ7eX4PNun159cQi1V9rZZr98VEZellC7rSTUZeVlARZoArqR1U2dt++cFwE3A04HPRsTK3pWX1fL2cs8826fay2NKqKWK/g24kNaf1V4GnAa8FXgEuDQiLu5hbVn4sGw9KiKuB05f5G6vSil99RDtHgb8C/Bc4A9TSu/ssMTKiogPAr8G/ElKaXbvjIh4MnAncGdK6all11dVEfFi4DPAbuBJKaUf9bikwnhZQDOdQqtHsRjLDvWGlNL+iLiCVrieBzQuXHls7OZ8/x4D7eVDJdRSGymlz0bE14CfBc4Gtva2ouIYrnpUSumZGZufvtZ6QsZj9NL0TZn182yfXn9PCbXUzZ20wrVRnw2vuaosq9vLqYO+q76+2V6eOc/26fW3l1BL3TTys2G4qiwvay/nG6pUd1+mdUPv1IiY6zeAl7eXN5RXUvVFxDpal4ugYZ8Nw1WFiYg3RsSGWesiIn4d+B0gAR/oSXGZpZT2Ae9vv7wqIqavsU5Pf30G8K8ppa/3or5eiohzI+KX2zc2Z64fBK6ndT36n1JK8w1jqyVHC6gwETFK69rircDdwNG0hmCdQmvQ+MUppffP20DNtR/cspXWjZn7gS/SGtd6NrAT6MsHt0TEq4Grge/T+mzspvXvchatz8g24PkppR29qjEHw1WFiYjfBl4MnAEcBxzBYyHz3pTSLT0srxQRsRT4A+CVwAbgh8AI8Pam9cwWKiJOB36b1n9kNtC6xjoFfAf4B+ADTRqCNc1wlaQMvOYqSRkYrpKUgeEqSRkYrpKUgeEqSRkYrpKUgeEqSRkYrpKUgeEqSRkYrpKUgeEqSRkYrlIBIuJVEZEi4lsRccQ87zknIvZHxA/azzFVgxmuUgFSSh8FPg88Dfi92dvbgfshWt+5N6WUdpZbocrmU7GkgkTEEPBtIIBnpJTunLHtbcBlwOdTSi/qUYkqkeEqFSgi3gL8GbA1pfSL7XWn0fobWweAp6eU7uphiSqJlwWkYv05cBuwKSJeGxEBfBA4CrjEYO0f9lylgkXEWcDNwCStsL0U+Abw7JTSI72sTeUxXKUMIuJdwJvaL/fT+vtZX+thSSqZ4SplEBFPAsZp3dz6SErptT0uSSXzmquUxx/TClaA8yLimF4Wo/IZrlLBIuI/Aa+l9ZdvPwWcCPxJT4tS6bwsIBUoIo6iNezqNODlwJdo/QnplbSuuzb+z4urxZ6rVKy30QrWf0opfSKl9ACtGVtLgA9FxOE9rU6lsecqFSQingbcCuwFNqaUxtvrA/hX4LnA76WUruxdlSqL4SoVICKWAF8GzgHekFJ636ztp9Ma6/oIcEZKabT0IlUqLwtIxfhNWsF6M3DV7I0ppe8AlwPLgL8qtzT1gj1XqUsRsR74D2ApcGZK6VvzvO8o4HbgqcArUkp/V16VKpvhKkkZeFlAkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjL4/13Ci0eEITXMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim = mp.Simulation(\n",
    "    cell_size=cell,\n",
    "    boundary_layers=pml_layers,\n",
    "    geometry=geometry,\n",
    "    sources=sources,\n",
    "    resolution=resolution,\n",
    ")\n",
    "\n",
    "nfreq = 100  # number of frequencies at which to compute flux\n",
    "\n",
    "# reflected flux\n",
    "refl_fr = mp.FluxRegion(\n",
    "    center=mp.Vector3(-0.5 * sx + dpml + 0.5, wvg_ycen, 0), size=mp.Vector3(0, 2 * w, 0)\n",
    ")\n",
    "refl = sim.add_flux(fcen, df, nfreq, refl_fr)\n",
    "\n",
    "# transmitted flux\n",
    "tran_fr = mp.FluxRegion(\n",
    "    center=mp.Vector3(0.5 * sx - dpml, wvg_ycen, 0), size=mp.Vector3(0, 2 * w, 0)\n",
    ")\n",
    "tran = sim.add_flux(fcen, df, nfreq, tran_fr)\n",
    "\n",
    "plt.figure(dpi=150)\n",
    "sim.plot2D()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the fluxes through a line segment twice the width of the waveguide, located at the beginning or end of the waveguide. Note that the flux lines are separated by length dpml from the boundary of the cell, so that they do not lie within the absorbing PML regions. Again, there are two cases: the transmitted flux is either computed at the right or the bottom of the cell, depending on whether the waveguide is straight or bent.\n",
    "\n",
    "The fluxes will be computed for nfreq=100 frequencies centered on fcen, from fcen-df/2 to fcen+df/2. That is, we only compute fluxes for frequencies within our pulse bandwidth. This is important because, far outside the pulse bandwidth, the spectral power is so low that numerical errors make the computed fluxes useless.\n",
    "\n",
    "As described in Introduction, computing the reflection spectra requires some care because we need to separate the incident and reflected fields. We do this by first saving the Fourier-transformed fields from the normalization run. And then, before we start the second run, we load these fields, negated. The latter subtracts the Fourier-transformed incident fields from the Fourier transforms of the scattered fields. Logically, we might subtract these after the run, but it turns out to be more convenient to subtract the incident fields first and then accumulate the Fourier transform. All of this is accomplished with two commands which use the raw simulation data: get_flux_data and load_minus_flux_data. We run the first simulation as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "field decay(t = 50.050000000000004): 4.825189380557789e-09 / 4.825189380557789e-09 = 1.0\n",
      "field decay(t = 100.05000000000001): 0.02880180987942578 / 0.02880180987942578 = 1.0\n",
      "field decay(t = 150.1): 0.0268934650933857 / 0.02880180987942578 = 0.9337421921042788\n",
      "field decay(t = 200.15): 2.3158397338096397e-13 / 0.02880180987942578 = 8.040604890819488e-12\n",
      "run 0 finished at t = 200.15 (4003 timesteps)\n"
     ]
    }
   ],
   "source": [
    "pt = mp.Vector3(0.5 * sx - dpml - 0.5, wvg_ycen)\n",
    "\n",
    "sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3))\n",
    "\n",
    "# for normalization run, save flux fields data for reflection plane\n",
    "straight_refl_data = sim.get_flux_data(refl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to keep running after the source has turned off because we must give the pulse time to propagate completely across the cell. Moreover, the time required is a bit tricky to predict when you have complex structures, because there might be resonant phenomena that allow the source to bounce around for a long time. Therefore, it is convenient to specify the run time in a different way: instead of using a fixed time, we require that |Ez|2 at the end of the waveguide must have decayed by a given amount (1/1000) from its peak value.\n",
    "\n",
    "The stop_when_fields_decayed routine takes four arguments: dT, component, pt, and decay_by. What it does is, after the sources have turned off, it keeps running for an additional dT time units every time the given |component|2 at the given point has not decayed by at least decay_by from its peak value for all times within the previous dT. In this case, dT=50, the component is Ez, the point is at the center of the flux plane at the end of the waveguide, and decay_by=0.001. So, it keeps running for an additional 50 time units until the square amplitude has decayed by 1/1000 from its peak. This should be sufficient to ensure that the Fourier transforms have converged.\n",
    "\n",
    "Finally, we save the incident flux using get_fluxes which will be used later to compute the reflectance and the transmittance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save incident power for transmission plane\n",
    "straight_tran_flux = mp.get_fluxes(tran)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to run the second simulation which involves the waveguide bend. We reset the structure and fields using reset_meep() and redefine the geometry, Simulation, and flux objects. At the end of the simulation, we save the reflected and transmitted fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "     block, center = (-2,-11.5,0)\n",
      "          size (12,1,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "     block, center = (3.5,2,0)\n",
      "          size (1,28,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "field decay(t = 50.050000000000004): 1.697652269444903e-10 / 1.697652269444903e-10 = 1.0\n",
      "field decay(t = 100.05000000000001): 4.6910710639105265e-07 / 4.6910710639105265e-07 = 1.0\n",
      "field decay(t = 150.1): 2.992872733686264e-07 / 4.6910710639105265e-07 = 0.6379934758846679\n",
      "field decay(t = 200.15): 0.0039278135652722765 / 0.0039278135652722765 = 1.0\n",
      "field decay(t = 250.20000000000002): 0.00015009081939073738 / 0.0039278135652722765 = 0.038212307406279115\n",
      "field decay(t = 300.2): 8.806226395655623e-11 / 0.0039278135652722765 = 2.2420174097660296e-08\n",
      "run 0 finished at t = 300.2 (6004 timesteps)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAIhCAYAAAD6ovlZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df5RdZX3v8fd3JgkkE8gvAkUSGUKVAmIRSkGq1ygqs/Da1qX31muuole7elut2qKtF7UiaIX6o7cq9VavP7DE2x8qWqodrcumVauIoKLUJUgcnAEkAZkZMiEmZJ77xzkDwzCT+XH2s8/e57xfa83aOXuf/ZzvTuZ88sxznmdPpJSQJBWrp90FSFInMlwlKQPDVZIyMFwlKQPDVZIyMFwlKQPDVZIyMFwlKQPDVZIyMFwlKQPDVZIyMFwlKYNl7S5AjxYRPwVWAcPtrkXqYpuBvSmlX1jKyeFdsaonIsbp5QjWt7sSqYv9DFavXM39998fSzndnms1DbOeU77//e+zcdXGBZ0wcWCCLe/dAsDdF91N34q+Jb3wyPgIA9sHGBodon9tP4PbBtl05KYltdVOXsf8Jicn6enpYfc9u7ni8isY/MIgeyf2wjxRkiYTy1cs57TTTuMNf/wGzjrrLKY6aRGzn9yuf4+J/RMc8+5jANj56p30LV/Y+2L33t084QlP4LEbH7vk1zZcK2zjqo0cvfroBT13Yv/EQ3/uW9G3pHAdHhvmgk9cwNDoEFvWbWHHhTvYvGbzottpN69jYabCdU/vHsbuGePOoTvZv3//gs9/zFGPoXeyl74VfYcM16r8exzdd/SSOx1L4QdaAhpvgK1XbWXnfTtrH0hex+JEBMuXL2f58uULPqenp4dly5bN2VOd0in/HkthuKpj3gBeR2vmC8rpenrmj45O+fdYKsO1y3XKG8DraF2RH253yr9HKwzXLtYpbwCvo1o65Tpa1dXhGhFnRsQbIuLTETESESki5vzvOyIumXrOHF+Xl1l/KzrlDeB1VEunXMcNd97QchvdPlvgzcBvLOG8rwE/mmV/6/8iJeiUN4DXUS3DY8M8/eNPr/11XDdyHc/9m+e23E63h+vXgZuA65tfQ8BhCzjv/6aUPpavrHw66Y3sdVRHJwXrUz76FA6mgy231dXhmlK6YvrjxXxaWked9Eb2Oqqj04L1wckH6Y1eDtJawHb1mGs36aQ3stdRHXfvu7sjrmN6sC7rWca1L7y25Ta7uufagmdExOnA4cAI8E8ppcLHWycOTDxi5dUhnzvteTPPmbn08PMv+jzrV65fcNtV4XUUa2qF1sSBCQ7EAdLyBAtcR5CWJfau2surbngVP933U/rX9vO5//a5yv17HOp9MeX6O67nmX/9TA6mg/RGL//83/+Z49ce3/Jre+OWaSJiH3BYSmnW8YGIuAR4yxynfwp4aUppzyJe7+Y5Dp3IRg7jlQttSVLhroRTNp7CzTffvKTxQocFFudHwOuAU4HVNG5Jtg24A3g+8NftK01SldhznWa+nushzjsW+B6wAXhySukbLdZxMxs55e6dC7+71fS7/7RyVyx1j6lhgV27d/HGN76Rz372s427Yi1A77JezjjjDN7+9rdz7pPPnfeuWO0y1/tivrt0Teyf4Jgtx7TUc3XMtQAppbsi4qM0erUDQEvhOmWpd7da6nnqLlPh2re8j+VpOXEg4MDCzo0U9B7sZWXvynnvilUVU++Lsu7S5bBAcW5tbo9taxWS5lTmLA17rsVZ19xW56NSdb3JSbj33rmP9fTAPfcE+/YdQUpHASsX1G5KyzhwYC333beM3bthanRxto7rhg2N12m3kfERLvjEBaVNGzNcCxCNn4We13x4YztqGBkfacfLquLuvReOnvN+61OJtxF4f/NrYQ4ehBtugF//9ak9cw8H7NoFGxf2CzWymhpjLWs+bgX+P6mHiNgYEa+MiCNm7F8NfAA4G/gp8OmyaxseG2Zg+0DZLyvVStm/CaGre64R8RwaN2+ZsqK5f/oHUpellD4H9NH4r/3yiLgeuIvGf/ln0JglMAq8IKW0sI9bCzI1hjQ0OlTmy0q107+2v9QVZF0drjTC8exZ9p894zkA9wJXAOcAjwfOBQ4CPwY+Bvx5SumObJXOYvrgfP/afgNWOoTBbYOlLs3t6nBt3tnqYwt87v3AG3LWsxgzP/X8/Is+zy9d+UvtLksVs2FDY8xzNg/99tfdu3nrW9/KP/7jP7J378I+j+3tXcbpp5/OW97yFs4555xDTsXasGHJ5Req7N/+29XhWlezTSdZv3J9u8tSBfX0zP1h0tRsgZQShx9+PxH3sNDJLhHLWL58lHXrHmTjxkPPFuhWfqBVM51yNyWp0xmuNWKwSvVhuNaEwSrVi+FaAwarVD+Ga8UZrFI9Ga4VNjI+YrBKNeVUrAorey20pOLYc60wg1WqL8O1wspeCy2pOIZrhZW9FlpScQzXCit7LbSk4hiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJM4yMj7TchuEqSdMMjw0zsH2g5XYMV0lqGh4bZutVWxkaHWq5LcNVkng4WHfet5P+tf0tt2e4Sup604N1y7otDG4bbLnNZQXUJUm1NTNYd1y4g/Ur17fcrj1XSV1rtmDdvGZzIW0brpK6Us5gBcNVUhfKHaxguErqMmUEKxiukrpIWcEKhqukLjEyPlJasIJTsTpGEWuhpU42sH2AodGhUoIV7Ll2hKLWQkudrMxgBcO19opcCy11sv61/aUFKxiutVb0Wmipkw1uGywtWMFwra0ca6GlTrbpyE2lvp4faNVQrrXQkopjz7VmypynJ2npDNcaMVil+jBca8Jgleqlq8M1Is6MiDdExKcjYiQiUkSkBZz30oj4ZkTsiYifRcTnI+LcXHUarFL9dPsHWm8GfmMxJ0TE/wZeAzwAfBE4HHgW8OyIeEFK6TNFFmiwSvXU7eH6deAm4Prm1xBw2FxPjohn0gjWe4Enp5Rube5/MrAD+GhE7EgpjRZR3Mj4CBd84gKDVaqhrg7XlNIV0x9HxHyn/GFz+7apYG228/WI+D/Aq4GXA+8uor6y10JLKk5Xj7kuRkSsBJ7RfPjJWZ4yte+5Rb2mwSrVl+G6cCfRGDLYnVKa7RZUNza3TyzqBcteCy2pOF09LLBIj21uZ723X0ppIiJGgXURcURK6f75GoyIm+c4dCKUvxZaUnHsuS7c6uZ27yGeM9HcHlHEC5a9FlpScey5tlFK6dTZ9jd7tKeUXI6kAtlzXbg9ze2qQzynr7mdd0hAUmczXBfuJ83trD+rR0QfsBa4byHjrZI6m+G6cD8Efg5sjIjjZjl+RnN7U3klSaoqw3WBUkoPAF9uPvwvszzlBc3tteVUJKnKDNfFeU9z+6aIeNzUzuby198BRoEPt6MwSdXS1bMFIuI5NG7eMmVFc/83pu27LKX0OYCU0pci4i9o3F/gOxHxz81zngUE8LKi7isgqd66OlyBjcDZs+w/e8ZzHpJSem1EfAd4FY1Q3Q98iUYI/3uuQiXVS1eHa0rpY8DHyjpPUvdwzFWSMjBcJSkDw1WSMjBcJSkDw1WSMjBcJSkDw1WSMjBcJSkDw1WSMjBcJSkDw1WSMjBcJSkDw1WSZhgZH2m5DcNVkqYZHhtmYPtAy+0YrpLUNDw2zNartjI0OtRyW4arJPFwsO68byf9a/tbbs9wldT1pgfrlnVbGNw22HKbXf2bCCRpZrDuuHAH61eub7lde66SutZswbp5zeZC2jZcJXWlnMEKhqukLpQ7WMFwldRlyghWMFwldZGyghUMV0ldYmR8pLRgBadidYwi1kJLnWxg+wBDo0OlBCvYc+0IRa2FljpZmcEKhmvtFbkWWupk/Wv7SwtWMFxrrei10FInG9w2WFqwguFaWznWQkudbNORm0p9PT/QqqFca6ElFceea82UOU9P0tIZrjVisEr1YbjWhMEq1YvhWgMGq1Q/hmvFGaxSPRmuFVb2WmhJxXEqVoWVvRZaUnHsuVaYwSrVl+FaYWWvhZZUHMO1wspeCy2pOIZrhZW9FlpScQxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcF2CiNgREekQXwPtrlFSe/lrXlrzKWDPLPvvKLsQSdViuLbmdSmloXYXIal6HBaQpAwMV0nKwGGB1rw8IjYAk8AtwGdSSj9pc02SKsBwbc2bZjx+V0RcllK6bCEnR8TNcxw6sbWyJLWbwwJL82/Ai2mE4CrgJOCNwIPApRHxmjbWJqkC7LkuQUrpT2bsugX404j4FvAF4JKI+GBK6YF52jl1tv3NHu0phRQradFGxkdabsOea4FSSl8EvgWsBc5uczmSlmB4bJiB7a2vAzJci3drc3tsW6uQtGjDY8NsvWorQ6NDLbdluBZvXXM70dYqJC3KVLDuvG8n/Wv7W27PcC1QRGwEntp8eGM7a5G0cNODdcu6LQxuG2y5TcN1kSLi3Ij4zYjonbG/H7gG6AP+IaXU+oi4pOxmBuuOC3ew6chNLbfrbIHFezzwUeCnEXEjMAocD5wJHA7cDPx2+8qTtFCzBevmNZuZ2N/6qJ7hunjXAR+gMRvgLBpjrBPAd4C/Bz4w3xQsSe03V7AWxXBdpJTSD4Dfa3cdkpYud7CCY66SukwZwQqGq6QuUlawguEqqUuMjI+UFqzgmGvHKGIttNTJBrYPMDQ6VEqwgj3XjlDUWmipk5UZrGC41l6Ra6GlTta/tr+0YAXDtdaKXgstdbLBbYOlBSsYrrWVYy201MmKWNK6GH6gVUOzTSdZv3J9u8uSNI0915opc56epKUzXGvEYJXqw3CtCYNVqhfDtQYMVql+DNeKM1ilejJcK6zstdCSiuNUrAorey20pOLYc60wg1WqL8O1wspeCy2pOIZrhZW9FlpScQzXCit7LbSk4hiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJM4yMj7TchuEqSdMMjw0zsH2g5XYMV0lqGh4bZutVWxkaHWq5LcNVkng4WHfet5P+tf0tt2e4Sup604N1y7otDG4bbLnNZQXUJUm1NTNYd1y4g/Ur17fcrj1XSV1rtmDdvGZzIW0brpK6Us5gBcNVUhfKHaxguErqMmUEKxiukrpIWcEKhqukLjEyPlJasIJTsTpGEWuhpU42sH2AodGhUoIV7LkuSUSsjIhLI+KWiNgXEXdGxEci4rh21FPUWmipk5UZrNBCuEbEn0fEqiKLqYOIOBz4MvBmYDXwWWAYeBnw7YjYUmY9Ra6FljpZ/9r+0oIVWuu5vgb4XkQ8s6hiauJNwDnA14HHp5R+K6V0NnARsBH4SFmFFL0WWupkg9sGSwtWaC1c/wY4AfhCRHw0ItYVVFNlRcQK4FXNh69MKe2ZOpZSeg9wE/C0iDgzdy051kJLnWzTkZtKfb0lh2tK6UXAc4E7gAuB/4iI/1pUYRX1a8Aa4LaU0rdnOf7J5va5OYuYbTpJ2d84kg6tpQ+0UkqfA04GrqTxI/H/i4jPtuuDnRL8cnN74xzHp/Y/MVcBZc7Tk7R0LU/FSilNAL8fEduBD9PotT0tIv4KmDjEeZe2+tpt8Njmdq55T1P7j19IYxFx8xyHTpxtp8Eq1Udh81xTSt+IiCcB/wb8KvC6OZ4aQALqGK6rm9u9cxyf+s/kiKJf2GCV6qWwcG1OQfoQcBZwELiGQ/RcBSmlU2fb3+zRnjL12GCV6qflcI2IoDEN6RJgFfBd4BUppRtabbuCpmYHzDW/t6+5vb+oFzRYpXpq6QOtiDgNuA64AugF3gj8SocGK8BPmtu5Ppqf2n97ES9W9lpoScVZcs81It4GvB5YDnwF+O2U0i1FFVZR321uz5jj+NT+m4p4sbLXQksqTis914uBB4DfTSk9rQuCFeBrwBhwYkScPsvxFzS31xbxYgarVF+thOu1wCkppb8qqpiqSyntB97ffHhlREyNsRIRf0hjfuu/FjUsUvZaaEnFWfKwQErpN4ospEbeBjwTOBe4NSK+QmNe69nAbuB/FPVCZa+FllQcbzm4SCmlfcDTgctozHf9TRrh+jHgjJTSzqJeyyWtUn15s+wlSCk9APxJ80uSHsWeqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyRlYLhKUgaGqyTNMDI+0nIbhqskTTM8NszA9oGW2zFcJalpeGyYrVdtZWh0qOW2DFdJ4uFg3XnfTvrX9rfcnuEqqetND9Yt67YwuG2w5TaXFVCXJNXWzGDdceEO1q9c33K79lwlda3ZgnXzms2FtG24SupKOYMVDFdJXSh3sILhKqnLlBGsYLhK6iJlBSsYrpK6xMj4SGnBCk7F6hhFrIWWOtnA9gGGRodKCVaw59oRiloLLXWyMoMVDNfaK3IttNTJ+tf2lxasYLjWWtFroaVONrhtsLRgBcO1tnKshZY62aYjN5X6en6gVUO51kJLKo4915opc56epKUzXGvEYJXqw3CtCYNVqhfDtQYMVql+DNeKM1ilejJcK6zstdCSiuNUrAorey20pOLYc60wg1WqL8O1wspeCy2pOIZrhZW9FlpScQzXCit7LbSk4hiukpSB4SpJGRiukpSB4SpJGRiuixARWyMiHeLrG+2uUVI1uEJraW4DvjrHfkkyXJfoqymll7a7CEnV5bCAJGVguEpSBg4LLM3jIuIdwAbgHhrjr4Mppcn2liWpKgzXpTm3+TXd9yLi+SmlWxfaSETcPMehE5dcmaRKcFhgccaAdwLn0Oi1bgDOA74BnAZ8MSLWtK88SVXRVT3XiLgGOHmRp70kpfRNgJTSt4Fvzzj+5Yh4CvAvwFOB3wPesZCGU0qnzlHnzcApi6xTUoV0VbgCJwAnLfKcVfM9IaV0MCKuoBGu57PAcJXUuboqXFNKp2dsfmqs9diMryGpJhxzLc665nairVVIqgTDtTjPb25vbGsVkirBcF2EiHhtRGyesS8i4neAPwAS8IG2FCepUrpqzLUArwXeFRE3Aj8GDqcxBesEYBJ4dUrphjbWJ6kiDNfFeTfwbOBUGlOllgN3AVcD700pXd/G2iRViOG6CCml9wHva3cdkqrPMVdJmmFkfKTlNgxXSZpmeGyYge0DLbdjuEpS0/DYMFuv2srQ6FDLbRmuksTDwbrzvp30r+1vuT3DVVLXmx6sW9ZtYXDbYMttOltAUlebGaw7LtzB+pXrW27XnqukrjVbsG5es3n+ExfAcJXUlXIGKxiukrpQ7mAFw1VSlykjWMFwldRFygpWMFwldYmR8ZHSghWcitUxilgLLXWyge0DDI0OlRKsYM+1IxS1FlrqZGUGKxiutVfkWmipk/Wv7S8tWMFwrbWi10JLnWxw22BpwQqGa23lWAstdbJNR24q9fX8QKuGcq2FllQce641U+Y8PUlLZ7jWiMEq1YfhWhMGq1QvhmsNGKxS/RiuFWewSvVkuFZY2WuhJRXHqVgVVvZaaEnFsedaYQarVF+Ga4WVvRZaUnEM1worey20pOIYrhVW9lpoScUxXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkmYYGR9puQ3DVZKmGR4bZmD7QMvtGK6S1DQ8NszWq7YyNDrUcluGqyTxcLDuvG8n/Wv7W26va8M1Ivoi4sUR8b6IuC4ifh4RKSIuWcC5myLioxFxZ0Tsi4hbIuKtEXF4CaVLKtj0YN2ybguD2wZbbnNZAXXV1eOAjy/2pIj4ReDrwFHA94GvAL8C/AlwXkScl1L6eZGFSspnZrDuuHAH61eub7ndru25AvcDHwb+J3AmjXBciI/RCNb3ppROSyn9FnAScA3wa8D/Kr5USTnMFqyb12wupO2uDdeU0m0ppVeklP4qpXQjcGC+cyLiV2kE6C7gj6a19SDwu802Xh0R3fwTgVQLOYMVujhcl+g5ze21M3/0TyndTWOIYB3wlLILk7RwuYMVDNfF+uXm9sY5jk/tf2IJtUhagjKCFbr7A62leGxzO9fyjan9xy+ksYi4eY5DJy6mKEkLU1awgj3XxVrd3O6d4/hEc3tECbVIWoSR8ZHSghVq3HONiGuAkxd52ktSSt/MUc9SpJROnW1/s0d7ymLaKmIttNTJBrYPMDQ6VEqwQo3DFTiBxhSoxVjV4mvumaedvub2/hZfZ1GKWgstdbIygxVqHK4ppdPb8LI/AZ4EbJrj+NT+28spp9i10FIn61/bX1qwgmOui/Xd5vaMOY5P7b+phFoWtBZ6chJ2737k1+RkGdVJ1TK4bbC0YAXDdbE+19w+NyIOm34gIo4BngrcB3wtdyELXQt9771w9NGP/Lr33tzVSdWz6ci5fuDMw3BdhOaHYV8DjgaumNrfXJH1l8ByGsti513t1YrZppOU/Y0j6dBqO+ZahOaMg2ObDx/T3L4iIqY+HborpfS8Gae9jMaNW14TEc8A/gM4C9gC/Dvwjpw1zzVPb2L/xPwnSypNV4crjQ+nZk74P675BbN8MJVSujUingRcCgwAz6PxQddlwJ/mvCNWmROgJbWmq8M1pdS/xPOGafRgS2OwKpeUUrtL6EiOudaAwarcDNjidXXPteomJycZ+tkQ51193kPB+uUXf5njjjiOyenzqSYnSXfv4qjmsGu6exeTx26Gnp7mtKueR7XrdKzulVIiIjh48CA9PT1MTk6ybNkyent7211aRzFcK+zOPXdywScumL/Heu+9rH7sFnZPPX7nFti1CzZupGeWn016enpm3a/u0tP8Jli1ahUTExM88MADba6osxiuFXb+1edz+9jtnLD2BAZfOMjGwzayb9++Rz9x3z5m/vKuffv2wb59NJ5++KOOzdaMusvBgwfp6+vjjjvu4Ec/+hEHDhwgIhwiKIjhWmG3j93OmoNreNYdz+KD7/zgnM/r27uXS2bsu/zyy5lYtYq9e/tgxtHLL7+cVaucutWtJicn6enpobe3l/HxcX74wx9yyy23AI69FslwrbBle5ZxYPsB/u6Bvztkj2LD5OSjwvXqq6/m3p4eJic3MDNcr776anp6XKbV7ZYvX87o6CgHDmRd89K1DNcKW3vtWu656x72znn72IbZ/hHHxscZnePo+PgYNI9KysOPNSqsZ8J/HqmufPdWWETM+mdJ1We4Vtj0MVY/aJDqxXCVpAwMV0nKwHCVpAwMV0nKwHCtMGcIqKoiwu/PeRiuFZYmnSGg/JYSkgcOHHjkndn0KK7QqrDeZQ/fAq6np2fO6Vgxy/6A5pLZ2d449jrU+P6YnJx8xPfVQr4vUkqsWLHiobtqaXaGa4U97hcfx5qVa1i+fPkh7y2w7sEH4Qc/eMS+k08+mfuWLePBB9fNPMTJJ5/MsmX35SpbNbFs2TJ27drFXXfdxTHHHMP69evp7e2d+z/xZvAePHiQ5cuX09/fz8qVKx91XA2Ga4VddNFFcODh8a25vulXjI3Bi1/8iH0XX3wx+9esYWxsxcxDXHzxxaxZsz9X2aqJ3t5e9uzZw549e1i9ejWHH374IX9CAh76PowIjjzySI4//viH9uuRDNcKO++88+hb0Tf/E3fvftSu888/HzZunO0Q559/Phs3FlCgpDk5aNIJNmxgYuTHbHw9bHw9TIz8GDZsaHdVUlez51phKaWF3VMggnTUUdzT7OSmo44iRUBKNE5/5I9s6aH90sO/U2spHA6Ym+FaYYuZSzjzDlpTj2c7vXG8kBLVAQzIPBwWkKQMDFdJysBwlaQMDFdJysBwlaQMnC3Q4TZsgF27Hr1PUl6Ga4fr6cHVWNIijYyPtNyGwwKSNM3w2DAD2wdabsdwlaSm4bFhtl61laHRoZbbMlwliYeDded9O+lf299ye4arpK43PVi3rNvC4LbBltv0Ay1JXW1msO64cAfrV65vuV17rpK61mzBunnN5kLaNlwldaWcwQqGq6QulDtYwXCV1GXKCFYwXCV1kbKCFQxXSV1iZHyktGAFp2J1jCLWQkudbGD7AEOjQ6UEK9hz7QhFrYWWOlmZwQqGa+0VuRZa6mT9a/tLC1YwXGut6LXQUicb3DZYWrCC4VpbOdZCS51s05GbSn09P9CqoVxroSUVx55rzZQ5T0/S0hmuNWKwSvVhuNaEwSrVi+FaAwarVD+Ga8UZrFI9Ga4VVvZaaEnFcSpWhZW9FlpScbq25xoRfRHx4oh4X0RcFxE/j4gUEZfMc16a5+vwomo0WKX66uae6+OAjy/x3Angk3McO7jENh+l7LXQkorTzeF6P/Bh4Prm13OASxd47j0ppZdmqushZa+FllScrg3XlNJtwCumHkfEs9tYzqzKXgstqThdO+YqSTl1bc+1RX0R8UbgscBe4NvAp1NKe4p8kYn9E0t67mLOkzrZUt8XRbyHDNelOQp424x974mIC1NKn1toIxFx8xyHTgQ45t3HLKm4pZ4ndbKy3xcOCyzex4EB4DhgNfAk4K+BDcCnI+KsNtYmqSJq23ONiGuAkxd52ktSSt9s5XVTShfO2PUd4CURMQxcTKNHe/4C2zp1tv3NHu0pO1+9k6P7jl5QXRP7Jx76n/nui+6mb0Xfgs5rxcj4yEMLHfrX9jO4bbCWH8J5HdVS5HUs9X2xa2IXW67csqTXnFLbcAVOAE5a5DmrchTS9GfAHwNbI2JFSml/qw32Le9bUkj2rVjaeYsxPDbMBZ+4oPYLHbyOasl5HYt5X/Ttb/39UzcRN6UAAAd9SURBVNtwTSmd3u4apkspjUXELuBYGkMEd7W5pGw65WYyXke1dMp1THHMtSAR0QMc2XzYsR/Xd8obwOuolk65jukM1+IMAH3AbSml8XYXk0OnvAG8jmrplOuYyXBdhIh44WyzASLiacCHmg+vLLeqcnTKG8DrqJZOuY7Z1HbMtQjNGQfHNh8+prl9RUQMNP98V0rpedNOGQAujIhbgJuBA8Djganx378B/iJv1eXrlDeA11EtnXIdc+nqcKUxR/X4GfuOa34B3D7j2N/S+Ds7E3g6jXmuPwP+CfhISmmuO2XVVqe8AbyOaumU6ziUrg7XlFL/Ip//TzSCtCt0yhvA66iWTrmO+Tjmqll1yhvA66iWTrmOhTBc9Sid8gbwOqqlU65jobp6WKDqdu/dveDnThwo5q5YM5cefv5Fn2f9yvW1u9OW11Et7bqO6e3vmti14JVXi3nvzSVSSi03omJFxDi9HMH6dlcidbGfweqVq7n//vtjKacbrhUUET+lcR+E4XbXMo8Tm9vb2lpFNfl3M7e6/N1sBvamlH5hKScbrlqyqfvRznV3r27m383cuuXvxg+0JCkDw1WSMjBcJSkDw1WSMjBcJSkDZwtIUgb2XCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXCUpA8NVkjIwXFWYiNgaEekQX99od425RcTKiLg0Im6JiH0RcWdEfCQijpv/7M4VETvm+d4YmL+VevHXvCiH24CvzrG/Y0XE4cCXgXOAu4DPAv3Ay4D/HBHnpJR2tq/CSvgUsGeW/XeUXUhuhqty+GpK6aXtLqIN3kQjWL8OPDultAcgIv4QeDfwEWBr26qrhtellIbaXUQZHBaQChARK4BXNR++cipYAVJK7wFuAp4WEWe2oz6Vz3CVivFrwBrgtpTSt2c5/snm9rnllaR2clhAOTwuIt4BbADuoTH+OphSmmxvWVn9cnN74xzHp/Y/sYRaquzlEbEBmARuAT6TUvpJm2vKwnBVDuc2v6b7XkQ8P6V0azsKKsFjm9uROY5P7T++hFqq7E0zHr8rIi5LKV3WlmoyclhARRoD3knjQ50Nza/zgG8ApwFfjIg17Ssvq9XN7d45jk80t0eUUEsV/RvwYhq/VnsVcBLwRuBB4NKIeE0ba8vCm2XrIRFxDXDyIk97SUrpm/O02wv8C/BU4OKU0juWWGJlRcQHgd8G3p5Smtk7IyJ+EbgVuDWl9Piy66uqiHg28AVgFHhMSumBNpdUGIcFNN0JNHoUi7FqvieklA5GxBU0wvV8oOPClYfnbs7199HX3N5fQi21kVL6YkR8C/gV4GxgR3srKo7hqoeklE7P2PzUWOuxGV+jnaY+lNk0x/Gp/beXUEvd3EojXDvqe8MxV5VlXXM7cchn1dd3m9sz5jg+tf+mEmqpm4783jBcVZbnN7dzTVWqu6/R+EDvxIiY7SeAFzS315ZXUvVFxEYaw0XQYd8bhqsKExGvjYjNM/ZFRPwO8AdAAj7QluIySyntB97ffHhlREyNsU4tf30i8K8ppRvaUV87RcS5EfGbzQ82p+/vB66hMR79Dymluaax1ZKzBVSYiBiiMbZ4I/Bj4HAaU7BOoDFp/DUppffP2UDNNW/csoPGBzN3AV+hMa/1bGA30JU3bomIlwIfBX5K43tjlMbfy5k0vkduBp6RUtrVrhpzMFxVmIj4feDZwKnA0cByHg6Z96aUrm9jeaWIiJXA/wJeBGwGfgYMAm/utJ7ZQkXEycDv0/hPZjONMdYJ4AfA3wMf6KQpWFMMV0nKwDFXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJWkDAxXScrAcJUKEBEviYgUEd+LiOVzPOeciDgYEfc072OqDma4SgVIKX0c+BLwBOCPZh5vBu6HaLznLkop7S63QpXNu2JJBYmILcD3gQCemFK6ddqxNwGXAV9KKT2rTSWqRIarVKCIeD3wZ8COlNLTm/tOovE7tiaB01JKt7WxRJXEYQGpWH8OfBvYGhEvj4gAPggcBlxisHYPe65SwSLiTOA6YJxG2F4KfAc4K6X0YDtrU3kMVymDiHgXcFHz4UEavz/rW20sSSUzXKUMIuIxwAiND7c+klJ6eZtLUskcc5XyeCuNYAU4PyKOaGcxKp/hKhUsIv4T8HIav/n2M8BxwNvbWpRK57CAVKCIOIzGtKuTgBcAX6XxK6TX0Bh37fhfL64Ge65Ssd5EI1j/IaX0qZTS3TRWbPUAH4qIZW2tTqWx5yoVJCKeANwI7ANOSSmNNPcH8K/AU4E/Sim9s31VqiyGq1SAiOgBvgacA7w6pfS+GcdPpjHX9UHg1JTSUOlFqlQOC0jF+D0awXodcOXMgymlHwCXA6uAvyy3NLWDPVepRRGxCfgPYCVwRkrpe3M87zDgJuDxwAtTSn9bXpUqm+EqSRk4LCBJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGRiukpSB4SpJGfx/OCANBb1N5oYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim.reset_meep()\n",
    "\n",
    "geometry = [\n",
    "    mp.Block(\n",
    "        mp.Vector3(sx - pad, w, mp.inf),\n",
    "        center=mp.Vector3(-0.5 * pad, wvg_ycen),\n",
    "        material=mp.Medium(epsilon=12),\n",
    "    ),\n",
    "    mp.Block(\n",
    "        mp.Vector3(w, sy - pad, mp.inf),\n",
    "        center=mp.Vector3(wvg_xcen, 0.5 * pad),\n",
    "        material=mp.Medium(epsilon=12),\n",
    "    ),\n",
    "]\n",
    "\n",
    "sim = mp.Simulation(\n",
    "    cell_size=cell,\n",
    "    boundary_layers=pml_layers,\n",
    "    geometry=geometry,\n",
    "    sources=sources,\n",
    "    resolution=resolution,\n",
    ")\n",
    "\n",
    "# reflected flux\n",
    "refl = sim.add_flux(fcen, df, nfreq, refl_fr)\n",
    "\n",
    "tran_fr = mp.FluxRegion(\n",
    "    center=mp.Vector3(wvg_xcen, 0.5 * sy - dpml - 0.5, 0), size=mp.Vector3(2 * w, 0, 0)\n",
    ")\n",
    "tran = sim.add_flux(fcen, df, nfreq, tran_fr)\n",
    "\n",
    "# for normal run, load negated fields to subtract incident from refl. fields\n",
    "sim.load_minus_flux_data(refl, straight_refl_data)\n",
    "\n",
    "pt = mp.Vector3(wvg_xcen, 0.5 * sy - dpml - 0.5)\n",
    "\n",
    "sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3))\n",
    "\n",
    "bend_refl_flux = mp.get_fluxes(refl)\n",
    "bend_tran_flux = mp.get_fluxes(tran)\n",
    "\n",
    "flux_freqs = mp.get_flux_freqs(refl)\n",
    "\n",
    "plt.figure(dpi=150)\n",
    "sim.plot2D()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the flux data, we are ready to compute and plot the reflectance and transmittance. The reflectance is the reflected flux divided by the incident flux. We also have to multiply by -1 because all fluxes in Meep are computed in the positive-coordinate direction by default, and we want the flux in the −x direction. The transmittance is the transmitted flux divided by the incident flux. Finally, the scattered loss is simply 1−transmittance−reflectance. The results are plotted in the accompanying figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxQAAAIpCAYAAAA7EH5WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde1yUZf7/8dcFgpDiIU3ME6amLViaaWmpkKVpB9eyfutqbrIae7C+2mHN1L6ZG7W1lUvb4WvqSm1t267mdlI3S8XUTqZoSnnMc5CRGigIyvX7YxgCYYQZbmYGeD8fDx4T93Xf133NJDCfuT6f6zLWWkRERERERHwREugBiIiIiIhI7aWAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfKaAQkREREREfOZIQGGMucwYM9UY86Yx5oAxxhpjbDX6a26MSTHG7DXGnCx+/IsxppkT4xUREREREWcYa31+3/9TJ8b8B/j5mcettcaHvloCHwNdgN3AeiCu+Gs70M9a+0O1BiwiIiIiIo5wKuXpY+CPwHDgfOBkNfr6C65g4k2gm7X2F9ba7sBfga7AM9Ucq4iIiIiIOMSRGYpynRqTDzT0dobCGHM+cAA4BXSw1maVamsI7AfOBdpYa79zcMgiIiIiIuKDYCvKHoprTB+VDiYArLUngXeAUOD6AIxNRERERETOEGwBRY/ixw0e2t3HL/HDWEREREREpBINAj2AM3Qofjzgod19PKaqHRpjtnpo6grk4UqjEhERERGpC9oDJ6y1rf11w2ALKBoXP57w0H68+DHKgXuFNGzYMKpz586xDvQlIiIiIhJwu3bt4uTJ6qyP5L1gCygcZ62Nq+i4MWZr586dY7du9TSBISIiIiJSu8TFxZGRkeHXDJxgq6HILX48x0N7o+LHHD+MRUREREREKhFsAcW+4sd2Htrdx/f6YSwiIiIiIlKJYAsoNhU/9vLQ7j6+2Q9jERERERGRSgRbQLEMKAIGGGNalW4o3tjuJuA0sCQAYxMRERERkTMEJKAwxtxljPnaGPN46ePW2m+B14Fw4AVjTOmi8SeB84BXtUu2iIiIiEhwcGSVJ2PMDcBDpQ6FFx//pNSxP1pr3yv+75ZAN+D8CrqbDPQFRgJfG2PWA3FAd2AHcK8TYxYRERERkepzatnY84ArKjh+xRnnVMpa+70x5nJgJjACuBnIAp4FHrbWHq3eUEVERERExCnGWhvoMQSEMWZrbGys9qEQERFxgLWW+vqeQqSmGGMwxnh1TfE+FBme9mKrCXV+YzsRERGpGadPnyY7O5ucnBwKCgoCPRyROik8PJyoqChatGhBaGhooIdTIQUUIiIi4rXTp0+zb98+8vPzAz0UkTqtoKCA7Oxsjh8/TocOHYIyqFBAISIiIl7Lzs4mPz+f0NBQoqOjadSoESEhwbYavUjtVlRUxPHjx8nKyiI/P5/s7GxatWpV+YV+poBCREREvJaTkwNAdHQ0TZs2DfBoROqmkJCQkp+vQ4cOkZOTE5QBhT5KEBEREa9Ya0tqJho1ahTg0YjUfe6fs4KCgqBc/EABhYiIiHil9BsapTmJ1LzSP2cKKEREREREpE5RQCEiIiIiIj5TQCEiIiIiIj5TQCEiIiIiIj5TQCEiIiISxHbt2sXNN99My5YtCQkJwRjDqlWr2LNnD8YYEhISAj1Eqee0D4WIiIgEtawsmDcP0tIgJweioiAhAcaPh+joQI+uZhUVFXHrrbeSnp5O3759ufDCCwkJCaF169YBHZcxhpiYGPbs2RPQcUhwUEAhIiIiQSkvDyZNgtRUKCws27Z8OcycCYmJkJICERGBGGHN27NnD+np6QwYMIDVq1eXaxMJBgooREREJOjk5cGwYa5ZCU8KC+Gll2DbNli6FCIj/Tc+fzlw4AAAnTp1CvBIRDxTDYWIiIgEnUmTzh5MlJaWBpMn1+x4qqJ0TcOPP/7IvffeywUXXEBYWBiTSw1w2bJl3HDDDZx33nk0bNiQTp06ce+995KdnV2mP2MM8fHxALz88ssYY7yqmfj000+57bbbOP/88wkPD6ddu3ZMmDCBffv2ebxm2bJlDB8+nOjoaBo2bEj79u258cYbWbRoEQCpqakYYwDYu3dvyZjOHFd6ejpTpkzhsssuK/M8f//733Po0KGzvnZ5eXlMnTqVmJgYGjZsSJcuXXjiiSc8buiWnZ3N9OnTufjii2nUqBFNmjTh4osvZsqUKXz77bcVPseqvP5SdZqhEBERkaCSmelKc/LGggUwa1Zw1FTk5eURHx/P3r17iY+Pp1evXjRv3hyAqVOn8sQTTxAeHk6fPn04//zz2bRpE7Nnz+btt99m7dq1RBc/iTvuuIPMzEz++9//0rlzZ/r37w/ARRddVOkYXnjhBe6++24A+vTpw4ABA9i2bRvz58/n7bffJi0tjZ/97Gdlrrnvvvt45plnCAkJoV+/fnTo0IFDhw6xdu1aDhw4wMiRI+nSpQt33HEHL7/8Mo0aNeLWW28tub70uP70pz+xaNEiLrnkkpJxp6en8+KLL/Kf//yH9evX06ZNm3LjLigoYMiQIWRkZJCQkMDx48dJS0tj6tSp5OTk8Oijj5Y5/6uvvmLIkCEcOHCA1q1bc9111wGwfft2/vznP3PllVcyYsSIkvO9ef3FC9baevkFbI2NjbUiIiLindOnT9uMjAybkZFhT58+Xa69qMjaI0d8/5o+3Vrw/mvGjOrdt6ioeq/LN998YwEL2H79+tkjR46Uaf/Xv/5lAdu9e3e7Y8eOUq9Xkf3f//1fC9hf/OIXZa5ZuXKlBewdd9zh8X7x8fFljn/88cc2NDTUtm3b1q5fv75M27x58yxgr7jiijLH//73v1vAtmnTxm7cuLFM24kTJ+z7779f5hhgY2JiPL4WK1assJmZmWWOnT592j7yyCMWsImJiRU+F/fzOXbsWEnb559/bkNDQ+0555xjc3JySo4XFhbabt26WcBOnjzZnjx5skyfW7ZssTt37iz53pfXP1hU9jNXWmxsrAW2Wn++r/bnzYLpSwGFiIiIbyp7c3PkiG8BQaC/znj/77XSb4o///zzcu09evSwgP3yyy/LtRUVFdmePXva0NBQe/jw4ZLjvgQUP//5zy1g33nnnQrHOXz4cAvYDRs2lBz72c9+ZgH7z3/+s0rPtbKA4mzatm1rW7RoUeaY+7mEhITYr7/+utw1N954owXsypUrS4698cYbFrBxcXH21KlTld7Xl9c/WAR7QKEaChEREREHnX/++fTu3bvMse+++45NmzZx4YUX0r1793LXGGO46qqrOH36NF988YXP9y4qKuLDDz/knHPOKUn/OdOAAQMA+OyzzwA4dOgQX331Fc2aNeP//b//5/O9z5Sdnc2CBQu47777GD9+POPGjWPcuHEUFhaSnZ3NDz/8UO6amJgYunXrVu54165dAcrURHzwwQcATJgwgdDQ0LOOxV+vf32lGgoRERERB3Xo0KHcMfcSrzt27Cgpavbk+++/9/ne33//Pbm5uQCEh4dX6T779+8HXCtJVTa2qnr99ddJSkoqGUtFcnJyOPfcc8sca9euXYXnRkVFAXDy5MmSY+5xd+7cudLx+Ov1r68UUIiIiIijmjaFI0d8v/6ppyA52fvrZsyA++7z/b5Nm/p+bWkRFWyKUVRUBFCmcNiTmJgYn+/tvk/jxo0ZOXLkWc+Ni4vz+T5ns3fvXsaNGwfAX/7yF2644Qbatm1LZPG6vldeeSUff/yxOwW9jJCQmkme8dfrX18poBARERFHGQPNmvl+/V13wZNPlt/M7mzCwlzXVee+Ncn9yXvLli1J9XYJKy+0bNmSiIgIQkJCWLBgQZVmHNq3bw/A7t27sdZWe5ZiyZIlFBQUcP/99zNp0qRy7bt3765W/27uce/atavSc/31+tdXqqEQERGRoNK6NRR/wF1liYnBsWSsJ+3ateOiiy4iIyOD7du319h9GjRoULIPxocfflila9q0acPPfvYzjh49yr///e8qXRMWFsapU6cqbDtSPD1VUfrS6tWrycrKqtI9KnPttdcCMH/+/JIZCE/89frXVwooREREJOikpEDxnm6Vio93nR/sHnroIYqKihg5ciTp6enl2rOzs5k7d2617zN9+nRCQkJITExk1apV5dpzc3P529/+Rl5eXsmxqVOnAnDvvfeyefPmMufn5+ezfPnyMsfatGlDVlYWR48eLde/u4D61Vdf5fjx4yXHDx48yG9/+1ufn9eZbrnlFrp27cqWLVuYMmUKhWdMaW3durXMbIi/Xv/6SClPIiIiEnQiI2HpUtcO2AsWVJz+FBbmmplISYEKyhaCzujRo9m6dSuPPfYYl112GT179qRz585Ya9m1axebN2+mcePG3HnnndW6T//+/Xn++ee56667uPrqq+nevTtdu3YlLCyMPXv2kJ6ezsmTJ7nllltK6hp+9atfsX79ev7617/Sq1cv+vXrR/v27fn2229JT08nJiamzJvw4cOHl5x75ZVXEhERQbdu3fjDH/7A8OHDiYuLY/369XTp0oWrrrqK/Px8Vq5cSc+ePbnyyitZt25dtZ4juGZjFi1axODBg3n66af5xz/+Qb9+/bDWsmPHDrZs2cLixYvp1KkT4L/Xvz7SDIWIiIgEpchImDMH9u93FWkPHgx9+7oek5Ndx+fMqR3BhFtycjJpaWmMHDmSzMxM/vOf/7By5UpOnz7N7373O95++21H7vPb3/6W9evXc8cdd5CTk8O7777Lf//7X3JzcxkzZgzvvvsuTc+oQn/22Wd56623uPbaa8nIyGDRokXs3LmT/v378/DDD5c59/HHH+euu+7i1KlTvPHGG8yfP5/33nsPcK0u9dFHH/G73/2OiIgI3n33Xb766ivuvvtuli9fTlhYmCPPEaB79+5s2rSJ+++/n6ioKJYsWcKKFSswxvDAAw/Qt2/fMuf76/Wvb0xFFfb1gTFma2xsbOzWrVsDPRQREZFapaioiG3btgHQrVu3GluZR0RcvPmZi4uLIyMjI8NaWzPLeFVAvwFERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERERMRnCihERERExC9WrVqFMYZx48YFeijiIAUUIiIiEtyysiA5GYYMgX79XI+PPeY6HkT27NmDMYaEhIRAD6XWmTlzJsYYUlNTK2xPTU3FGMPMmTP9Oi6pmgaBHoCIiIhIhfLyYNIkSE2FwsKybcuXw8yZkJgIKSkQERGIEYqXLr/8cr766iuaNm0a6KGIgxRQiIiISPDJy4NhwyAtzfM5hYXw0kuwbRssXQqRkf4bn/jknHPO4aKLLgr0MMRhSnkSERGR4DNp0tmDidLS0mDy5JodTyVmzpzJBRdcUDycNIwxJV/uegFjDB07dqSgoIBZs2Zx0UUX0bBhQ0aMGAFAfn4+8+fP5+c//zmdOnUiMjKSZs2aMXDgQP75z39WeN9x48ZhjGHVqlWsXr2aQYMGERUVRZMmTbjhhhvIyMgod421ltdee43+/fsTHR1NREQE7du359prr+X555/32P8HH3zAwIEDiYqKolWrVtx5550cO3YMgO+++47f/OY3tG3bloiICC6//HJWrVpV7t4V1VB07NiRRx55BIDExMQyr92qVatISEggMTERgEceeaRMuztFylrL66+/zqhRo+jatSuNGjUiKiqKyy+/nBdeeIGioqIK/5+5+/jyyy8ZPnw4zZs3p1GjRsTHx7Nu3ToP/7fh008/ZdSoUbRt25aGDRty/vnnc8011zB37txy5544cYLHH3+cSy+9lMaNG9O4cWP69u3Lyy+/7LH/2kYzFCIiIhJcMjNdaU7eWLAAZs2C6OgaGVJlevbsyciRI1m0aBHR0dEMHTq0pK1///4l/11UVMSIESNYvXo18fHxXHLJJbRo0QJw1WBMmDCBNm3a0K1bNy6//HIyMzNZt24dH330EV9//bXHGoJ33nmHlJQUevfuzfXXX096ejpLlizh008/ZcuWLbRu3brk3ClTpvDUU0/RsGFDBg4cSMuWLcnMzGTz5s3s3LmTiRMnlut/8eLFPP/88/Tr14+hQ4fyySefMG/ePHbs2MHChQvp168fp0+fZsCAAezZs4dPP/2UoUOH8vnnn3PxxRef9bW79dZb+eCDD9i0aRNXXXUVXbp0KWlr3bo1Q4cO5dSpU6xdu5YePXrQs2fPknb3uSdPnmT06NG0aNGC2NhYevXqRXZ2NuvWrWPixIl89tlnHusz1q9fz8SJE+ncuTPXXXcdX3/9NatXr+aaa67h888/p3v37mXOT0lJ4d5776WoqIjLLruMgQMH8v3337N582b+8Ic/cOedd5ac+9133zF48GA2b95M69atiY+Px1rLunXrGDduHOvXr+evf/3rWV+fWsFaWy+/gK2xsbFWREREvHP69GmbkZFhMzIy7OnTp8ufUFRk7ZEjvn9Nn24teP81Y0b17ltUVK3X5ZtvvrGAjY+Pr7AdsIDt0qWLPXDgQLn277//3i5fvtwWnTGO3bt3244dO9qQkBD7zTfflGm74447LGBDQkLs4sWLS46fOnXKjhw50gL2oYceKjmel5dnGzZsaKOiouzu3bvL9FVYWGhXr17tsf9333235PiPP/5ou3fvbgEbGxtrb7/9dltQUFDSPmPGDAvYX/3qV2X6W7lypQXsHXfcUeb4ww8/bAG7YMGC8i+ctXbBggUWsA8//HCF7YWFhXbx4sVlxmCttd99953t3bu3BWxaWlqF9wRsSkpKmbbJkydbwI4dO7bM8bS0NGuMsVFRUfaDDz4oN4b33nuvzLHrr7/eAnbSpEk2Pz+/5HhmZmbJuJYuXVrhcyqt0p+5UmJjYy2w1frxfbVSnkRERMRZx45B8+a+fyUn+3bfRx+t3n2L03dq2uOPP07btm3LHW/RogXXXnstxpgyxy+44AKmT59OUVER77zzToV9/vKXvyxJnQIIDQ3lwQcfBGD16tUlx3/88UdOnjxJ586dS1K03Bo0aMCAAQMq7H/06NHccMMNJd9HRUWVfBJ/4MABnn32WcLCwkra77//fowxpFU1ba2aGjRowIgRI8qMAeC8887j8ccfB+Ctt96q8NqrrrqK//mf/ylzbMaMGUDZ1w7gT3/6E9Zapk+fzjXXXFNuDNdff33J9+5Zoj59+vDMM8/QsGHDkrbo6GheeuklAF588UVvnmpQUsqTiIiIiJ8YY7jpppvOes6aNWtYtWoVBw8eJD8/H2st3377LQA7duyo8JohQ4aUO9a1a1eAkmsBWrVqRbt27UhPT2fq1KkkJSXRqVOnSsddUf/u63r37k3z5s3LtDVt2pRzzz23zL39IT09nffff5+9e/dy4sQJrLXk5OQA3r12LVq0KDf+U6dOldSFJCUlVTqW999/H4ARI0YQElL+M3x3TcVnn31WaV/BTgGFiIiIiJ+0atWqzCfVpR07doxbbrmFFStWeLze/eb4TO3atSt3LCoqCnDVF5T28ssvM2rUKJ544gmeeOIJYmJiiI+PZ9SoUQwbNqzC/iuaUWncuLHHNnd7dna2x+fipIKCAsaNG8frr7/u8RxvXjtwvX4//PBDyffZ2dnk5eVx7rnnlgugKrJnzx4Apk+fzvTp0z2el5+fX2lfwU4BhYiIiDiraVM4csT36596yre0pxkz4L77fL+vH/ZGiDjLfhkPPPAAK1asID4+nkceeYTu3bvTrFkzQkNDef/997nuuuvcdaDlVPQJuCeDBg1i586dvPvuuyxbtoxVq1bxyiuv8MorrzBy5EgWLlzoVf/e3LumPPPMM7z++utcfPHFPPnkk/Tq1YvmzZsTFhbG9u3b6datmyOvnTfcK0v179+fzp0718g9goUCChEREXGWMdCsme/X33UXPPlk+c3sziYszHVdde4bYIsXLyY0NJS3336bJk2alGnbvXu3o/dq0qQJo0ePZvTo0QB88skn3HbbbSxatIglS5aUqQWoDRYvXgzA66+/TlxcXJk2p167li1bEhkZyQ8//MDRo0dpVsm/NffMx4gRI7ivOoFuLRD4kFJERESktNatodQ+BVWSmBiwJWPdwsPDAVeuvS+OHDlCkyZNygUTAP/617+qNbbK9O3bl7FjxwKwZcuWGr1XRSp77SprP1I8I1ZR+pJTr11oaCgJCQkAJQXVZzN48GDgp2CnLlNAISIiIsEnJQXi46t2bny86/wAa9myJWFhYezatYvTp097fX3Xrl05cuQIb7zxRpnjs2fPZuXKlY6Mcd++faSmpnLixIkyx/Pz80vu0b59e0fu5Y02bdoAsG3bNp/a3QXo//d//1fm+MKFC3nllVecGiYPPPAAxhiSk5PL/T85deoUS5YsKfn+iiuuYPDgwaxdu5aJEyfy448/lutv06ZNLFu2zLHxBYoCChEREQk+kZGwdCkkJbnSmSoSFuZqX7YMzlKb4C/h4eEMHTqUzMxMevTowa9+9SsmTJjAggULqnS9e5nXUaNGMXDgQEaPHk1cXBz3338/99xzjyNj/OGHH0hMTOS8884jPj6eMWPGMGLECDp06MAnn3xC7969ueWWWxy5lzeGDBlCREQEs2fPZtiwYYwfP54JEyaUBBB9+/alVatWLFy4kISEBH79618zYcKEkt2sp0yZQmhoKFOnTqV3796MHj2aPn36cNtttzn22gHEx8fz5JNPkpOTw6BBg+jTpw+jR49myJAhtG3btiSFzO3VV1/l0ksv5YUXXiAmJoarr76aMWPGcOONN9KhQwd69uypgEJERESkxkRGwpw5sH+/q0h78GDo29f1mJzsOj5nTlAEE27z5s1j7NixZGdn849//IP58+dXeS+GMWPG8N5779G3b1/S09NZunQpbdq0YcWKFQwfPtyR8XXu3Jmnn36ahIQE9u3bx5tvvsmaNWuIiYlh9uzZpKWleVyFqia1adOGt956i759+7JmzRr+9re/MX/+/JJlWyMiInjvvfcYPHgw6enppKamMn/+fLZv3w7AwIEDWbNmDYMGDWL37t28++67hIeHs2jRogp3/q6O+++/n7S0NG6++Wb27dvHwoUL2bJlCxdffDFPP/10mXNbtWrFunXrePbZZ4mNjWXjxo0sXLiQzZs306lTJ/785z9z//33Ozq+QDCeKt7rOmPM1tjY2NitW7cGeigiIiK1SlFRUcknx926dQuKVX5E6jJvfubi4uLIyMjIsNbGeTzJYfoNICIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPlNAISIiIiIiPmsQ6AGIiIiInE1WbhbzNswjbW8aOQU5RIVHkdAxgfGXjie6cXSgh1fCGENMTAx79uwJ9FBE/EoBhYiIiASlvMI8Ji2bRGp6KoVFhWXalu9ezsxVM0nsmUjKsBQiGkQEaJQiooBCREREgk5eYR7DXhtG2t40j+cUFhXy0oaX2Ja9jaVjlhIZFunHEYqIm6M1FMaYSGPMLGPMdmNMvjHmkDHmb8aYtj70NdgY854x5rAxptAYk22Med8Yc7OTYxYREZHgM2nZpLMGE6Wl7U1j8rLJNTwiEfHEsYDCGBMBrAAeAhoDbwH7gURgozGmkxd9TQbeB4YB24FFwNfAtcCbxphkp8YtIiIiwSUzN5PU9FSvrlmQvoCs3KyaGZADlixZwuDBg2nevDkRERF069aNqVOncvTo0XLnWmt57bXX6N+/P9HR0URERNC+fXuuvfZann/++TLnFhQU8MILL9CnTx9atGjBOeecQ8eOHbnxxhv55z//6a+nJ/WckylPM4C+wMfAEGttLoAx5l7gaeBvQEJlnRhjzgP+BBQCg621aaXaBuIKNB40xsy31u52cPwiIiLiAGstx04e8/n65z57rlzNRGUKiwp57vPnuK/ffT7ft2nDphhjfL7ek8cff5xp06bRoEED4uPjadmyJWvXruWJJ55g8eLFrF69mujon4rLp0yZwlNPPUXDhg0ZOHAgLVu2JDMzk82bN7Nz504mTpxYcu6YMWNYuHAhUVFRDBgwgCZNmnDw4EHWrFlDbm4uo0aNcvz5iJzJkYDCGBMO3FX87UR3MAFgrX3GGHMHEG+Mucxa+0Ul3V0BNAT+WzqYKO5rtTHmv8BwoDeggEJERCTIHDt5jOZPNPf7fR9d/SiPrn7U5+uPPHCEZhHNHBwRfP7558yYMYPGjRvzwQcfcMUVVwBw8uRJxo4dy7///W8mTpzIwoULAcjPz+evf/0rUVFRbNq0iQsuuKCkr1OnTvHxxx+XfP/NN9+wcOFCYmJi+OKLL2jRokVJW35+Phs3bnT0uYh44lTK01VAU2CXtbaif70Lix9vqkJfJ6t4z+wqniciIiISEM899xxFRUXcfffdJcEEQMOGDXnuueeIjIxk8eLF7N+/H4Aff/yRkydP0rlz5zLBBECDBg0YMGBAyfeHDx8G4NJLLy0TTABERETQr1+/mnpaImU4FVD0KH7c4KHdffySKvT1GXAUGGSMiS/dUJzydB2wA/jIh3GKiIiI+M1HH7nerowZM6ZcW6tWrRgyZAhFRUWsXbu25Fi7du1IT09n6tSp7N7tORnjoosuolGjRrz33nv8+c9/5tChQzXzJEQq4VQNRYfixwMe2t3HYyrryFp7zBgzHvgHsNIYs674+nbAlcBa4FfW2oKqDMwYs9VDU+eqXC8iIiLeadqwKUceOOLz9U+te4rkj7xff2XGwBnVrqFwmvtNfseOHStsdx8/ePBgybGXX36ZUaNG8WKk8pMAACAASURBVMQTT/DEE08QExNDfHw8o0aNYtiwYSXnNWnShLlz55KUlMSUKVOYMmUKXbt25eqrr2bs2LFcddVVjj8fkYo4NUPRuPjxhIf248WPUVXpzFr7Jq4VnrJxpVP9ovgxB1dR9kHPV4uIiEggGWNoFtHM56+7Lr+LsJAwr+4ZFhLGXX3uqtZ9a6IguzIV3XPQoEHs3LmT1157jbFjx1JUVMQrr7zC9ddfz6233lrm3F/+8pfs3r2buXPnctttt3H06FHmzJlD//79ue8+34MrEW84ug+FU4wx9wEfAKtxpUk1Ln5cAcwC3qxqX9bauIq+gF01MHQRERGpptaNWzOu5zivrknsmUh04+jKT/SzNm3aALB3794K2/fs2QNA27Zlt+xq0qQJo0eP5pVXXmHfvn18/PHHtGvXjkWLFrFkyZIy55533nlMmDCBf/3rX2RmZrJ06VKaNGnCM888w9atnhI1RJzjVEDhXtXpHA/tjYofcyrryBiTADwFpAO3WWu/tNYet9Z+CdxafPwGY8yws3QjIiIitVjK0BTiY+IrPxGIj4knZVhKDY/IN+4i6tdff71c2+HDh/nvf/+LMabS9KS+ffsyduxYALZs2eLxPGMMQ4cO5YYbbgBQQCF+4VRAsa/4sZ2HdvfxisPzssYWPy621haVbrDWnuan2YmBXo1QREREao3IsEiWjllKUq8kj+lPYSFhJPVKYtnty4hoEOHnEVbNxIkTCQkJ4dlnn2X9+vUlxwsKCrj77rvJy8vjlltuoX379gDs27eP1NRUTpwom0Wen5/PypUrAUrO3bhxI2+++SYFBWXLSn/44Qc+/fTTMueK1CSnirI3FT/28tDuPr65Cn25gw9PO+K4j/t/gWsRERHxm8iwSObcNIdZV89i/sb5rNqzipyCHKLCo0jomMD4S8cHZZpTaZdffjl//OMfmT59Ov369SMhIaFkY7v9+/dz4YUXltn9+ocffiAxMZGJEyfSu3dv2rVrx/Hjx1m3bh2HDx+md+/e3HLLLYArjWrkyJE0bdqU3r1707p1a44ePcrq1avJycnhpptu0tKx4hdOBRRrcb3R72yM6WmtTT+j3V1B9E4V+sosfuztob1P8eMer0YoIiIitVJ042imDZjGtAHTAj0Un0ybNo0ePXowe/ZsPv/8c/Ly8ujQoQNTpkxh6tSpNG/+02eknTt35umnn+bDDz8kIyODzz77jEaNGnHBBRcwbdo0kpKSaNiwIeBKg3r00UdZsWIF27Zt46OPPqJ58+ZccskljB8/nttvvz1QT1nqGWOtdaYjYx4FpgPrgCHW2uPFx+8FngbSrLUJpc6/C9fu2outtQ+WOn4zrrSm08AIa+27pdp+zk8pT7HW2m3VGO/W2NjYWOUWioiIeKeoqIht21x/grt160ZISFCu8SJSZ3jzMxcXF0dGRkZG8SJEfuHUDAXAo8C1uPaK2GGM+QjXvhNXAIeBX59xfkugG3D+Gcf/A/wbuA14xxizHvgGuICfZi2mVyeYEBERERERZzj2kYK1Nh+4Gvgjrv0oRuAKKFKBXtZaz1s9lu3H4tp3YjyuZWO7ADcDHYElwDBr7WNOjVtERERERHzn5AwF1to84H+Lvyo7dyYw00ObBf5W/CUiIiIiIkFKSY8iIiIiIuIzBRQiIiIiIuIzBRQiIiIiIuIzBRQiIiLiFWNMyX8XFRUFcCQi9UPpn7PSP3/BQgGFiIiIeMUYQ3h4OADHjx8P8GhE6j73z1l4eHhQBhSOrvIkIiIi9UNUVBTZ2dlkZWUB0KhRI21wJ+KwoqIijh8/XvJzFhUVFeARVUwBhYiIiHitRYsWHD9+nPz8fA4dOhTo4YjUeREREbRo0SLQw6iQAgoRERHxWmhoKB06dCA7O5ucnBwKCgoCPSSROik8PJyoqChatGhBaGhooIdTIQUUIiIi4pPQ0FBatWpFq1atsNbi2pdWRJxijAnKmokzKaAQERGRaqstb3xExHmqnhIREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ8poBAREREREZ81CPQARKRuycrNYt6GeaTtTSOnIIeo8CgSOiYwvOtw3tr2VpWPj790PNGNowP9dERERKQSxlob6DEEhDFma2xsbOzWrVsDPRSRWunMwKFRWCOOnTxG+rfpnLKnqt1/WEgYPVr3oEl4E06cOqFAQ0REpAri4uLIyMjIsNbG+euemqEQkbOq6cDBk8KiQtYfWl/m2PLdy5m5aiajuo+iU7NOrDuwTrMaIiIiAaYZCs1QiACBCxycpFkNERGp7zRDISJ+l1eYx6Rlk0hNT6WwqDDQw6mWs81qJPZMJGVYChENIgI0OhERkbpJAYVIPVN6JuLYyWPsyN7BkfwjgR5WjSosKuSlDS/x/u736dSsk2YvREREHKSAQqSeqEszEb7ac3QPe47uKflesxciIiLVp4BCpA5zz0as3LOSzw5+Rk5BTqCHFHTcsxf/zvg3F7a4kKYNm2rmQkRExAsKKETqIM1GeO9I/hE+O/gZoJkLERERbyigEKkD6mNdRE1T3YWIiEjVKKAQqcU0E1HzVHchIiJydgooRGoRzUQEB/fsxbbsbSwds5TIsMhAD0lERCRgFFCI1AKaiQhOaXvTuPCvF9K+aXulQ4mISL2lgEIkyOUV5jHstWGk7U0L9FCkAgdzDnIw5yCgdCgREamfFFCIBCl3etOL618secNaG4SFhNGzdU+iwqPKFDIP7zqct7e/zao9q8gpyPF4vFFYI348+SMbv93IKXsq0E/Ha0qHEhGR+sZYawM9hoAwxmyNjY2N3bp1a6CHIlJGbUlv8hQ4OJXyk5WbxfyN88sEIFe1v4qdP+zkja1vBPVr49Y2qq3SoURExK/i4uLIyMjIsNbG+eueCigUUEgQCcb0ppoOHHxxZrBRW2Y1wkLClA4lIiI1SgGFHymgkGCU9E4SczfMDegYmkc0r7U7RlcUaOw6sqvMsq/BID4mXulQIiJSIxRQ+JECCgkW7lqJ93e/z0d7P8ISmJ/JuvrpeV5hHpOXTWZB+oKgSpNSOpSIiNSEWh1QGGMigQeBUUAH4AdgGfCQtdbrilJjTEdgKnAd0AbIAXYAb1pr/+zAeBVQSEAFQ61EVHgUfdr04ZpO19T5N7SlZy+CcQ+PuhrQiYiIf9XagMIYEwGsBPoC3wIfAR2By4HDQF9r7W4v+hsGLAQigQ24AokWwMXAcWttFwfGrIBCAibQtRJ68xq8MxdKhxIRkeoIREDh1LKxM3AFEx8DQ6y1uQDGmHuBp4G/AQlV6cgYcxHwJq4ZicHW2nWl2kKAXg6NWSRgJi2b5NdgojbXRdSUyLBI5tw0h1lXzwqqugttliciIrVNtWcojDHhwHdAU6CXtXbjGe2bgEuA3tbaL6rQ3xJgGHCDtXZJtQZ39vtohkICIjM3kw6zO/jlU3HNRPgm2GYv9P9RRESqqrbOUFyFK5jYdWYwUWwhroDiJuCsAYUxpj2umondNRlMiARKVm4WY98cW2NvUjUT4QxPsxf7j+0PyCaD2ixPRESCmRMBRY/ixw0e2t3HL6lCXwlACLDOGNMAuAVXwBIKbAHesNYGTxWlSBXVdAG2PsGuGdGNo5k2YBrTBkwDAl/7onQoEREJRk4EFB2KHw94aHcfj6lCX7HFj7m4Crv7ntGebIy51Vq7sqqDM8Z4ymnqXNU+RKqjJt+Eto1qy+/7/F5vKv0kMiySpWOWBjQd6mDOwZJZkuW7lzNz1UwFkyIiElAhDvTRuPjxhIf248WPUVXoq3nx4wTgImA0cC7QDXi1+L8XG2Pa+jZUEf+rqQLs+Jh4dv7PTqYNmKZgwo/c6VD779lP8qBkBncaTN92fWkbFZhfS+50qKGvDiWvMC8gYxARkfrNqVWenOIOcBoAv7HW/qv4+yPAWGNMN6AP8HtgelU69FSQUjxzEVtRm4hTMnMzSU1PdbRPpTcFB6VDiYiIuDgxQ5Fb/HiOh/ZGxY85XvSVC/y7gvYFxY/xVRuaSGDN3zDfkbQYg2Fgh4EkD0pm/z37mXPTHAUTQcadDpXUK4mwkLCAjOFgzkE+OfAJy3cvZ/qK6bSf3Z7fvPMb8k/lB2Q8IiJSPzgxQ7Gv+LGdh3b38b1V6Mt9zj5b8Xq2e4ofW1VtaCKB5dSn1Xf2upM5N81xpC+pOVodSkRE6iMnAopNxY+eNpxzH99chb7cy84299B+bvFjrod2kaCRlZvF199/Xe1+4mPiSRmW4sCIxF+CMR1qYOpAmkc0J6cgRylRIiLiKCdSntYCx4DOxpieFbTfWvz4ThX6WgdkA62L6yXO5E51qmi/C5GgkFeYR9I7SbSf3Z79P+73uZ+wkDCSeiWx7PZlSm+q5YIhHWr9ofUs371cKVEiIuK4agcU1toC4Lnib583xrhrJjDG3Itr/4m00rtkG2PuMsZ8bYx5/Iy+TgHPAKa4ryalrrkWGAdYQLkfEpTcn0TP3TC3WrUTgzsNVq1EHRNsq0OBVogSERFnmIpLFbzsxJgIYBVwBfAtrj0kYoq/Pwz0tdbuLnX+TOBh4GVr7bgz+goDlgDXAlnAJ0BLXHtShALTrbWPOTDmrbGxsbFbt3rapkLEe0nvJDF3w9xq9REWEsb+e/YrFaWeCHQ6lFvbqLZaIUpEpA6Ii4sjIyMjw9NKpzXBiZQnrLX5wNXAH3HtRzECV0CRCvQqHUxUoa9C4HrgAeB74DrgYiANuMmJYEKkJji1RGxiz0S9katHgiEdCrRClIiI+M6RGYraSDMU4rTk1cnMWDmjWn3Ex8SrZqIey8rNCorVoUqLj4nXClEiIrVIIGYogm1jO5FaqzopK9qsTiD4VocCrRAlIiKVU0Ah4pCcgqrs3Vhe+ybt+fzOz/XmTMpxp0NNXjaZBekLHNkk0RfrD60v8/3y3cuZuWqmgmAREQEcqqEQqe+ycrM4fPywT9de1PIiBRPiUTCuDgVaIUpERH6iGQqRasgrzGPSskmkpqf6/OlxQscEZwcldVIwpkOBUqJERERF2SrKFp858YZOS8RKdeQV5gU8HcoT1QWJiASGirJFapFJyyZV+9NhLREr1eFOh5p19awyq0NFhUdxJP9IudoHf3KnRG3L3qZVokRE6jjNUGiGQnyQmZtJh9kdqvWpsJaIlZoULClRAEm9kphz05xAD0NEpF7QDIVILTF/w3yfgwmlgog/BMsKUQDzNs7j6++/pqCoQDUWIiJ1kGYoNEMhPhjy9yEs373c6+u6NO/Cml+v0Rsp8atg3DBPgbWISM0IxAyFlo0V8YGve060bNRSwYT4nXuFqPfHvs/H4z9mx907iI+JD+iYtOysiEjdoZQnER9EhUf59ToRJwVTOpSWnRURqf0UUIh4ISs3i3kb5rH7yG6frteeExIsgmmFKO3ELSJSu6mGQjUUUgVObGCnPSektgimFaLiY+K17KyIiBdUQyEShNxvruZumFut1BDtOSG1hTslKqlXEmEhYQEdS9reNCYvmxzQMYiIyNlphkIzFFKJpHeSmLthbrX60J4TUluduUJUIFKiQkwI/dv317KzIiJVEIgZCgUUCijkLKq7gZ2WxpS6KBhSovSzJSJSMaU8iQQZXzew69K8C8mDktl/z37m3DRHb3ikTgmGlCgtOysiEjwUUIicha+fwF7Q/AKmDZimlAyps9yrRO2/Zz/Jg5IZ3Gkwfdv1ZWDMQEKM//60qMZCRCTwtGysyFn4uoGdr9eJ1DbuTfOmDZhWcsyJuiNvzNs4j6+//1o1FiIiAaIZCpGz0AZ2It5LGZri1524i2wRq/et5pMDn7B893Kmr5hO+9nt+c07vyH/VL7fxiEiUl8poBDxICs3C18XLdAGdlKfqcZCRKR+0SpPWuVJzlDdTey0gZ3IT4Jh2dnebXrTPKJ5yf2VEiUidZmWjfUjBRRSESeWw0zqlcScm+Y4OCqRukXLzoqI1BwtGysSYJOWTarWm5z4mHhShqU4OCKRukcpUSIidYsCCpFimbmZpKan+nRtWEgYSb2StBu2SBVp2VkRkbpDy8aKFPN1E7vBnQbz95v/rnxsER8Ew7KzC9IXcG7kuXzx7ReqsxAR8YECCpFi1Ul10psOEeekDE1he/Z2v9VYFBYV8qe1fypzbPnu5cxcNVN1FiJSe2Rlwbx5sHcvneACf95aKU8ixbSJnUhwCIYaC1CdhYjUEnl5kJQE7dvDjBlw/DgREOnPIWiGQgTX0paHjx/26VptYifiPHeNxayrZ5VZdjY8NJw1+9ZQZIv8Npa0vWkMTB2opWdFJPjk5cGwYZAWuFXzQMvGatnYeq66e04AJA9KLpP/LSI1y981Fp5o6VkRqRHu1KW0NMjJgagoSEiA8eMh+owPMZKSYG7Z34futWK3Wmv8Ml4UUCigqMecWAtfm9iJ+F8w7GNRWnxMPEvHLCUyzK8ZBiJS1+TlwaRJkJoKhRV8yBkWBomJkJICERGQmQkdOpQ7NxABhVKepN6q7p4TAIk9ExVMiPiZu8Zi8rLJLEhf4PPsolPS9qbxm3d/Q7cW3Ujbm6a0KBHxXlVSlwoL4aWXYNs2WLoU5s+vOPAIAM1QaIaiXsrMzaTD7A7VeiMSHxOvfSdEAiwrN6tMjUVUeBRH8o+w/tD6QA8NUFqUiFRRBalLZzVmDGzZAps2lWtSypMfKaCo35JXJzNj5QyfrtUbBJHgFmwpUaC0KJF6qyr1EB5Sl3ylgMKPFFDUX1m5WVz1t6vYdWSX19d2ad6FNb9eoxQGkSCXV5gXNClRbmMvGau0KJH6wpt6iKefdi336hAFFH6kgKL+cWJFp77t+vLx+I8dHpmI1JSKUqIua3MZT619ilP2VKCHB2jWU6TO8WYp1/h4CAmBlSu9v0+DBnCq/O8xFWWL1BCnUiC054RI7RLdOJppA6aVW9o5+0R2UCw9Cz9toLcte5vSokTqgkmTqr4vRFoaGB/f9/foAZs3B0VhtgIKqRecWNEJIKFjQvUHIyIBlzI0he3Z24OqzkKrRYkEsaruDZGZ6Upz8oav2ULnngvjxnlXzF1DlPKklKc6LSs3i9mfzObJtU9iqd6/de05IVK3BGOdhSdKixIJEG/3hkhOdrQe4qySk+Gee8qlV6mGwo8UUNRtTtRLnCmpVxJzbprjSF8iEjyCfenZ0rRalIgfeVsLsXQpDBkCa9Z4fy9jvJupCAuD/ftdsyN5eTB5MixYAIWFCij8SQFF3ZSVm8WL61/k2U+f5Uj+Ecf61Z4TIvVLMC4966YPN0T8xNu9IVq0gOxs3+7VqhV8913Vz09Kgjln/B7IyoL584l77DHyjx/P22XtOb4NxnsKKBRQ1Ak1MSMBSjMQqc+CNSUqLCSM+/rdxxfffqE6CxFfBGBviEpdc41rxaaqzoYsW+ZKsapAXFwcGRkZGdbauApPqAEKKBRQ1FpZuVnM2zCPlXtW8tnBz8gpyHG0f+05ISJQcUpUQscEMg5n8NqXrwV6eCX0AYhIJQK4N0Sl3PUQpVKXzjo+D8EEKKDwKwUUtYs7eEjbm8axk8c4+ONBvs39liJbVGP3TB6UXG6pSRERt2BNi1KdhUgFvK2HCA2FFSu8v0/v3rBpk3czG6XrIaAkdYlVq86+opQHCij8SAFF7VBTqUyV0YpOIlIVwZoWpToLkTN4Ww8REQH5+d7fp29fuPhi7+5VUT1ENSig8CMFFMEvkJ/+6Y+xiHgj2NKiwkLC2JC0gbe2vaU9LaRuC7Z6iMGD4a23vJsNOUs9hC8UUPiRAorg5U5venH9ixzMOej3+2tFJxFxSiA/GAkxIRWmharWQuqEYK2HSE6GadPKLeV61vE5GEyAAgq/UkARfAKV3uSmP7IiUhOCNS1KtRZSa/mrHsJbZ9ZCQLXrIXyhgMKPFFAEl0B+imcwTO0/lUlXTFIagIjUmGBLiwKld0oQqkoKk7f1EOecAydOeD8WJ/aGCAAFFH6kgCK4JL2TxNwNXvxycPLe+oMqIgEUyA9UtKeFBI2qpjA9+CB07eqfeggH94bwJwUUfqSAInhk5mbSYXaHgKQCqF5CRIJBsKVFKQVU/MqbFKaOHWHPnpoekYuDe0P4kwIKP1JAETwe/OBB/rT2T369p/5Yikgwqigt6rI2l/HU2qc4ZU/5fTyqs5BqqUr6EnifwuQPDu8N4U8KKPxIAUXguYuw522Yh8U//w6bRzRn0hWT+G3v32o6X0RqDaWFSq3izQpMR4/6Z0nXWloP4YtABBQN/HUjkdL8nTOsGQkRqc1ShqawPXt7QOosFqQv4O7L79aeFlI1VUlfKiyEl16Cbdtcn/L7ox7i4ou9q4dISan5MdUhmqHQDEVA+OPTtqjwKPq06cM1na7RHz0RqfUqq7PwtO+EE7SnhQA1swJT27Zw0A97TtXSeghfKOXJjxRQBE5NF2HrD5yI1GWelp+9qetNXPbSZQFb4EK1FnVYMK7A5I1aXA/hCwUUfqSAIjCycrO4/c3b+eCbDxzpL9SE0rpxa9o2aUvThk01BS8i9ZpqLcRxwboCkzdqcT2EL1RDIXWW07tgt41qy+/7/F7Bg4hIKaq1EK9UJYVp0qSqBRPg32CiqsGL6iH8QjMUmqGocU4XYGvvCBERzwK5p4VqLWqJupDCtH07PP54na+H8IVSnvxIAYX/ODUFbzDc2etO/UESEamCYNvTAlRrERSCMYXJ28Ls0ilMdbwewhcKKPxIAYV/OFmA/WD/B3nsmsccGJWISP0VyDoLcNVazLp6FvM2zFNqVCAE4yZyM2fCypVVX9J12bJ6N+vgDQUUfqSAwj+SVyczY+WMavcTFhLG/nv26w+NiEg1+XsfoDOFmBBCCKlwlkSpUdVUWU1EZqZ/NpHzhnsFpiZN6sWSrv6gomypc97f9b4j/ST2TFQwISLigMiwSJaOWRqwPS2KbBFFVNx3YVEhL214iW3Z25Qa5Y2z1UQsX+6aAUhMhPPPD65gAlzjcqcmzZkDs2YphakW0gyFZihqhHtVp3kb5mGp3r8xFWGLiNSMYNzTwk2pUVXkTU1E8+Zw5EjNj8mbFZiUvuQ4pTz5kQKKmuPUdLqmvkVEAifQtRZKjaqiYKuJ0ApMAaeUJ6kTJi2bVO1gYnCnwfz95r/rEygRkQAJ5J4WoNSoEmeri7DWleYUTBITXTMUSl+qVzRDoRkKRzmxqpMKsEVEgkNle1rUZK1FVdXZ1Kiq7BXRowesX1/zY1EKU62ilCc/UkDhvKzcLG5/83Y++OaDavWT1CuJOTfNcWhUIiJSXcFca1EnU6O8qYuoaUphqnUUUPiRAgrnuAuwU9NTq/1HRQXYIiK1S6BrLaqi1m2oF0x1EdpErtZRQOFHCiic4VQBtnbBFhGpnQK9r0VV1ZrZb3/sFVHV1Z6UwlQrqSg7mFW2WUw95UQBNsDAmIG14xe9iIiUUZV9LcJCwjhddNpjkbU/LEhfwKyrZwEEvt7ibO8p5s+v+b0iJk2CQ4eUwiSO0QzF2WYosrLgxRddP/SHDrlWUyjfEZx3HnTvDtdcU68CjC+zvqTnnJ6OFOQlD0pm2oBpDoxKREQCxVOtxfhLx/PQyocCnhrVu01vNmVu8hj01Hi9RVUKrVu1goMHa+b+7nvs3+96r6IUpjpJKU9+VGFA4f7EYOVK+PJLOHy44iDibOpBVO+umZi/cb4jwYRWdRIRqftqS2pUjdVbBEuhdemaCKmTAhFQhPjrRkEtL8/1A9a+PcyYAR9+CN99530wAa5PHF56Cdq0gUcecQUpdUBWbhbJq5O59pVriX4qmrkb5jq2VGBiz0QFEyIidZw7NSqpVxJhIWEVnhMWEkZIgN+apO1NY/Kyyc53PGlS4IOJ+HjXB54iDnNshsIYEwk8CIwCOgA/AMuAh6y1Ps/dGWMuBDYDEcCH1tprHRiua4aia9fYrb/8JTz7bM1tRe/EjEUA6jeycrOYt2EeK/es5MvvvuTw8cNYnJ/N0qpOIiL1T7CnRrlnzsGhegt/FFr36QPp6aqJkNqb8mSMiQBWAn2Bb4GPgI7A5cBhoK+1drePfa8E4gGD0wEFxPptjacLLoAVK1ybw1SmqqlXDtZv+CuAcKu1a4OLiEiNCpbUKK/qLSr74C852ZUBUVPcdRGgmgip1QHFo8B04GNgiLU2t/j4vcDTQJq1NsGHfscD84CXgCRqc0Dh1qpV+QDAidoN8OkTCCf3kKiqwZ0G8/eb/640JxERqVBlO3QHw6pRAPHtB7D0ky5Epr569pmBnTtdHyrWFNVFSCm1MqAwxoQD3wFNgV7W2o1ntG8CLgF6W2u/8KLfaOArYD3wGK4ZkNofUJQdBDRsCPn5zvbbvLkrV/O3v63wE4nSsxGfHfyMnIIcZ+9/FirAFhGRqgr21CiApPUw591KToqKcs0Y1ATtFSFnqK0BxdXACmCXtbZLBe0PAbOAR6y1M73o95/Az4GLrQc6/wAAIABJREFUgXbUxYCipp0xYxGI2Ygz1ZqNhUREJKgFS2pU2GnY/wxEH6+Bztu2dS0So7oI8UJt3diuR/HjBg/t7uOXVLVDY8z1wC+A/7XW7jTGtKvG+Oqv4hWnsnZ/yYtTrubZDS9yJL+Gis+rID4mnpRhWl1CRESqr6ob6vVo3YP1h9bX2DgKQyHlCmhUCGkdISccogogYQ+M31DNQOP3v/9pszvVRUgQc2KG4hngHmC2tfbeCtp7AOnABmvtZVXorxGwFcgHLrHWFhhjEvBxhsIY42kSonMsNKzLMxR5DWDSMEjt6fqFFyghJoQJl05QAbaIiNSIs6VGWSwdZneo0Zl5Y8Ga8sfDTkPiRkhZBhGnvOy09AZ0Il6orTMUjYsfT3hod8fmUVXs71EgBrjaWltQnYHVZ3kNYNjtrk9LAinUhJL+m3S6R3cP7EBERKTOim4czbQB05g2YFqF7eN6jqvReouKgglwfZj3Um/Y1hKWvgqR3gQViYkKJqTWCKqN7YwxvYH/AV6x1q5yok9rbVxFX8AuJ/oPVpOGFQcTAd4Iffyl4xVMiIhIQKUMTSE+Jj4wN7euv8eTh+JKV6oKbUAntYwTAUVu8eM5HtobFT+edXkDY0wDYC5wFLjfgXHVW1+2gvmXFn/j4VMTf1DNhIiIBIOq7tLdu01v529e/Hd4waWQdVUP1xKvYRWPgbAwV7tWbZJaxomUp33Fj54Kp93H91bSTzugJ5AJ/NuYMu+EmxU/XmaMWQXgy74WdZ27ZmL+pVAUwLknbVonIiLBJjIskjk3zWHW1bPOXm/xVFsKjfP7WxSGwvwrI5n20Byypt7NvFfvIS17AzmmgCgbTkKLyxh/+zNEX6BZfal9nAgoNhU/9vLQ7j6+uYr9tS7+qkgzXLtm+1dFO1ID/N//uaYkjwRu5SS3QNdMRIVH0adNH67pdA3jLx2vfSZERCQoRR+Haast09Jw5U5EAQlAV8DCuHTL3Etd/+30LP+K8/PY807ST8u3N/+pbXnRcma+2ksfyEmt5PTGdpdaa9PPaPdpY7sz+kggEPtQVLJBHAB5eTB5smun66LA7NiZ1QhuHA3r21IjvwDPRrMRIiJSK+Tluf6mp6Z63tehRw/y0tfX2Ad0UeFRVdpMNj4mnqVjlhIZFun8IKTOC8QqT9VOjCleiem54m+fL172FQBjzL24gom00sGEMeYuY8zXxpjHq3v/GhES4sphPHQIHn747KssREb+//buPV6qut7/+OuzYbiIiKZAioKmqWfjSVNQTG2ThJeTesw6laY/b7nzVD9Bq99J05NamKejJadzLG9HPan9spuXCryE4CVNyWvwMw3DEGWDoYiwucnn98d3jQzDzN4za9asNbPn/Xw85rHca61Z68sWZuY93+/n+w3L3S9YALvumloTIfRKdB4Lu5wXhQmoe5gwhxFvw+EvwbTfwqLvOtf8KsZ0eCIiImnp7oajj4brrisdJiDsnzuXwRvCjEydc8O0r6VYzO9iKwkTAHNensPUmVPj3UQkAzX3UACY2SBgNnAQ8BrwEGHq14OAZcAEd3+p4PyLgW8AN7v7aRVcfyJp9VDsths88ACMGVP9RfO9FTfeWP4FqxZmMHAgrFmT2hAncxi+CvZZCpP+0sMiPUOHwvjxm4aEaao7ERFpFJ2dIUxUqWsI3LA/zN518wXrVg6Ayw9LvJWbybXlWHTuIgCuf/J65rw8Z4uaDw0vllKy6KFIJFAAmNlg4HzgJGAXYDkwE7jI3V8pOvdiGi1QJLmEfVdXWNVy1ix47jlYtgzi/J5L1W6MHEnXX/7IMf9zBHN5rS5DnIaugfGv9hIgepLk71JERKQWS5bA6NGJftG3ZGsYfW79F40dt9M4nlnyTNlVwDXkWEpp6kDRbMxsXvvgwe3zDj20/kvYFwWMrtVLuf6DoXdh5QAY+A7kNsC6/rBuUH8GDtqa3Hbbs+69I1jXzxnYbyC5thxr3lnDC397gWWrluF1WGCiphU9S+nogBkzwrAwERGReunqCrWMc+bAypWh1zz/3n799XDhhYnfsvNYuO4AUq9dLKZ6CymmQJEiM5vX3t7ePm9ej2XZieh6u4vrn7yeBxY+wHNLn6tbIIhru26Y8hicPTdGb0Rvxo+Hu+/WECgREUleJYXWI0bA4sXJ37qKoceVFmPH1bl/J9cce03dri/NRYEiRWkEiu713UyZOWXT9HANpm0jfO7JBHskyunfH844Q0OgREQkOflC6zlzsmtCf5h67t7cOHRBj8OS/rz8z8xaOKtu7VC9hRRSoEhRPQNF19td/GDuD/iP3/8Hb6zJfo2KUvpthKd/APssS/GmGgIlIiLV6Gko00UXxSq0TlRHB8ycSdeGFWUXyxu59UimPTiNCx9IfthVIdVbSJ4CRYrqESgavUcCeHesZ+dcuOZXGdy/szNMsysiIlJOJUOZ3nmnvus/jR8PTz9d/v5VTD6y5O0ljP7e6Mw/G6jeojUoUKQo6UCx8I2FfOR/PsLCNxcmcr166lgIM2/JaO2IXA4WLVJNhYiIlNYAQ5nefa+CMKnK7Nlb9pBU+T7WeXcn1z2ZcY8Km9db5Gs8NUSqb1GgSFFSgSLfK3HDUzew0bNZKbtS/d+BM5KcxSmuyZPhRz9SqBARkS3FXDMi8TYk3Jvevb6bo289mjkvZxiUCMOfXvjSC1z28GVlR1RoiFRzU6BIkZnNG7zT4PZDLz80dhpvlBeHSoxnJ+6+bhUjF6/IuimB1qoQEWld5Wojjj0WDjigPovDViqqi6jHe1P3+m6mzpzKjU/fWPaD/L7v3Ze5r85N/N6Fdt1214pGVGiIVHNSoEiRmc1jOO18MfoZY6ehO3HW/mdx9rizew0XXW93ccxtxzD3tfr+o09Cx5gOZp48k0Hrvb4recehQm0RkdbRW21EW1t96yJGjYKlSxOpi6hF19tdZYu4HW+Ieos8DZFqPgoUKSoOFMVGDBnBPsP3YdL7Jm32D6UpCq8j/a0/Z3zwjC27LPML7RWOCR03Lhx7/PHaVveOQ4XaIiJ9XyPURkybFuofEqqLqJdGqbcADZFqRgoUKeotUGx2Lsb2W23P1gO2ZsnbS1izYU39G1ij8TuN5+4T747/jUFx6Bg4EAYMCG8IzzwT9iVFhdoiIn1f1rURTfRe02hDquMOkVKPRjYUKFJUTaBoNu8OcarXtwTd3ckPnTr/fLjssmSuJSIi2Wjk2ogm6w1vlHqLanXu38lVR13V42gO9WjUlwJFivpioEj9H2hXFxxzDMxN4MXMDM46S0XaIiLNKOvaiN7UsdC63pqp3gLCZ5FxO43j0Vce7fXcUkXf6tWonQJFivpCoDCM4UOGl6z1SE3SY2JVpC0i0lyyro1oawuPDSXmQ2+BGQUbqd4ijnzRd281qurVqJwCRYoaIVDkA8Fe2+/FoP6DWLthLes2rmNgv4EM6Deg7M8Nl9aTHgLVZN3SIiItLevaiM5OuPTShi+0rpdq6i0qrYVIU77o+7Q7T6voz6CpbHunQJGiLAPFdoO2Y8pBUyqanrapdHXBySfD/ffXdp0mKpwTEWkJjVob0cRDmZJUSb3F6fudzvmHnc+e39+zoYZIAYzbaVxVtSCayrZnChQpyiJQtNHG5/b/XN/urluyBEaPrv3NRatpi4hkL+vaiHLXb4GhTHH0VG+R/3DdiEOkDMOp/POoprLtmQJFitIOFLttuxsPnPoAY7Ydk84Ns5RU97feMEREspN1bUQuB08+CXfd1ZJDmeql2YdI5Wm17/KyCBRtad2oVfW3/nTu38n8L85vjTABIQB0dNR+nfXr4dpr4aijwhubiIikZ8qUbBehO/102GcfuOACuPdeePTRsL3gAoWJGgzODWbGZ2fQuX8nubZcyXNybTk69+/kgVMfKHtO1ioNOnNensPUmVPr2xhRD0U9eyhqXlyumeULta+7LpkVt8ePh7vv1puIiEiSVBvR0pp1iFS1cm05Fp276N0/U1+vu9CQpxTVO1DUfXG5ZnH++XD55clcq39/OOMMDYESEamVaiOkQtUMkTp454N5YvETbPASU/iWUW39RFzTDp/GuRPObYmpaRUoUlSvQNFX/jImJqki7UJDh4Yei0mTNI5WRKRaqo2QKlU6i9T0o6dzzoxzqurRqHaGp7gm7TaJDRs3tMTUtAoUKTKzeXvuvWf7qdeeyqy/zOK5pc+xbNWy2Cm5z04Fm4R6zlHe1gY77gijRsGwYXpDEhHpTSOsG6G1hppSJUOkqunR6BjTwU3H35TKVLbbDNyGt9a+VfH5hVPTNhsFihSZ2bz29vb2efPmvbsv/w+lmoChHokKpP1tmEKGiEhp9eg1roZqI1pCNT0ag/oPasg6jWauu1CgSFGpQFGsMIkv717OijUrANhm0DZsP3j7hv2L1JCSXk27WmYwfHiYMURDpUSkFZQquHavffHRnqg2QgpU0qMBjTuVbbPWXShQpKiSQCF1kNRq2rXSm5uI9FW9FVzXi2ojpAaNuNp3s9ZdKFCkSIEiQ1l3uRfq6IAZM2Bw9i8AIiI1y7LgWrURkoB6TGU7YsgIlq5aWnVbBrQNYN3GdRWfX6ruIouhUgoUKVKgyFjWRYGFRo2CL3xB36CJSPPL6rVVtRGSomoLv/tZP2YtnFX3dhXWXXSv785sqJRWypbWkdRq2klYvBi+/nXYZRf4/OdhzZqsWyQiUl5XF0ybBkccAQcfHLaXXQbPPReGOdVDW5mPC7lcCDEKE5Kialb7nnnyTA7f7fBU2rV+43pueOqGdwPPdU9eV3Zo1vqN67n2yWs56paj6F7fnUr76kk9FOqhyE7WhdrlaBiUiDSirBajU22ENLBKhkgteXsJo783uqq6i7gL7n14zIfZa/u9qhqSVTxUqtZhUhrylCIFigbS1QXHHANz67+wTcU0FlhEGolqI0RqklbdRRz5oVLbDNwmkWFSGvIkrWnkSHjwwcYZAgWh16SrK+tWiIgEU6ZkEyY6OsIQVZEmN/2o6XSMqexzRseYDvYZvk+dW7TJ+o3r+eHcHzb1MCkFCmkMgweHYUadnaF7PWvr18MNN2TdChGRMDNevWojylFthPQxjVp3kXfdk9dVVGQOMOflOUydOXWL/V1vdzHtwWm8vOJl2I7dkm5jTzTkSUOeGk9XV/gwP2sWPP54GK+bhT32gIcf1vhgEUlPFovRTZ4ctqqNkBZRr7qLXFuOPd+zJ/Ner/9nyx5nlPqvcI4vdat7QyIKFAoUjS3rwm0tgCciachyMbpFixQeREqotu6ic/9O/vLmX7jvpfvq2KpN8it5bzGFbgaBQkOepLENHhyKARctCtMkTp4MBx4Y1o4oN41hktavh2uvhaOOCm/4IiJJyxdcX3dd+l+cnH66woRIGdXWXUw/uvLzk/DAwgeYMnNKxUOl6kk9FOqhaF75oVGzZ8OKFWE9iddeq8+0iaCZTkSkPrQYnUjD6l7fzdSZU7nx6RsrmnkpzlCpuIzQAbHF9LYa8pQeBYo+qrD+4rnnYNmyMP44CRoaICK1KFUfccABcMUVsGFD8vcrty6FhnKKVK2Suou8aodKjRo6isUrFyfXWAWK9ChQtIh8wLj66tCDUatp0+CCC2q/joi0jizqI7QYnUhm8itlVzIUqWNMBxPHTOSSBy9JrgEKFOlRoGgxSS0KNXky3HtvMm0Skb4vqwXpNERTJFPVDJV6c82bsWaU2v09u/P8689veVCBIj0KFC0oiRmjdtkFnnhC3+6JyOZKDWeaOBHmz4dbb023LaqNEGkYlQ6VSnRGKQWK9ChQtLCuLjj0UPjzn+M9X+OPRSQvq+leS9Frk0jTqnaY1MyTZ3Ll767kwgcu3PIETRsrkoKRI+G00+I/X1PJighkO93r5MnhMWFC2E6bFiaNuOYahQmRJlTtSt6D+g/izP3PLHtu2tRDoR6K1rRkCYweXfuHAI1TFmldWU33qhnnRPq0mmeU0pCn9ChQSCIfBvTGLtL3pT3da2/0RYaIREoOlVKgSI8ChSQ2+4qmkhXpmxqpPiJPBdciUmSLGaVUQyGSosGDYcaM8G1froYxiLNnJ9YkEWkQWdZHlJLLhdcqhQkRKTI4N5hrjr2GRecuYtrh0xgyYAhsINUiT/VQqIdCIAxpGD8+DF+q1oQJ8OijybdJRNJRakjTG2/A3LnptyWXg698Jdxbi9GJSAxjx45l/vz58919bFr37J/WjUQa2siRsPfe8QLF66+HDyR6sxdpLo04pOn00+Gyy7JuhYhIVVp6yNNf/gJHHBFeu7u6sm6NZK6jI97z/vznsODd5z8Pa9Yk2yYRqY9GG9IE4TVo+vSsWyEiUrWWDhTd3XDfffD1r8OOO8LOO8MllyhctKwzz4xfS6G1KUSay5QptU/IkBTVR4hIk2vpGgpob4fSNRQjRsBee8HAgbBuXXhoKGsLSGIqWU3pKNI4GmnK11NOCUMrZ89WfYSI1E0WNRQKFGUCRc/PhZ12grPOgrPP1vtAn5LEVLJam0Ike41WH6HpXkUkJVkEipYe8hSXOyxeDBdfDO99b/jcOGmSajH6hCSmkl2/Hm64Idl2iUjlGqk+QsOZRKQFqIciRg9Fz9eF7bcPPdn9+sE224Sf1avdhLq64NBDQ9F1tSZPhnvvTb5NIrK5RpnyVdO9ikiD0LSxfYB7mEX09dc333/ffXDhhTB8OOyzT+jR0PtMgxs5EnbYIV6gWLky+faIyCaNNqRJ072KSAvTkKcUucPSpTBr1qaZpTRcqsENHRrveS+8oP+pIvXSSEOaQNO9ikjL05CnhIc81cIs9GBodqkGMm1a6FqKK5cL31xOn67x0yJJSWI2tiTo37eINCDN8pSiRgwUvVHgyMCSJTB6dO3fgnZ0hGLvwYOTaZdIqyiukRgwAB5+GDZuTLcd48fDttuqPkJEGp4CRYqaMVD0Jh84VKORsKS+DdX6FCKVa6QaCU35KiJNRIEiRX0xUBRTj0ZCklibArQ+hTSkUpMkFb4+ZHH8o4d0M/WeoxnwaLYrWXsuh2lIk4g0GQWKFLVCoOhNYY/GQQeFovE//EE9+iV1d8PUqXDjjbV9WzptGlxwQXLtEiHeh/ZDDoEXX4Tbby/9VzqXg913D5OclVpQup7Hr6GTTq7DAYvzC6nSOnJcwVcYx1yGspKVDGU2E/mf/mfysTNGcsEFcMst5X+/0Pv/AxGRtChQpEiBojIaRlWkqyssWnfllbB8efXP1/oU0oOePpRC8qGgEYygi89xPR3MYSgrWcsADuVh+lP/Gol8YLmGTs6m/HBEs/CFS7FcDk45JRy75Zby/w/ynRwrVvQeOhRMRKRWChQpUqCIR70akYMPhsceq/55EybAo48m3x5pKEkHg7Zogu+065DraRDdTGcKp3ETA8iuRmI2HRzFTNZS3yFNO+4Iy5aV761JOpiAwolIq1KgSJGZzdtzz/b2U0+dx6xZ8Nxz4cW+RX8diWqJXo0jjgirFVZrjz3CDDV97hfSeuIMI+qLwSCOQXQzg6OZSHY1EuvIcSOnM4XpdQ8TSeotmJx+Olx+OfzLv5SvZ1c4EenbFChSZGbz2tvb2+fN29RDkR/NMnt2eMEcOBD69YP588OxFv1V1aytLbwJjhoVZk3N5fpAgXgt61No7vqmESc0SO/SqpHIX/9xxvMm225WH3EDZ7KUZnrRqdywYSEo9CbpcDJoUOXBQwFFpH4UKFJUKlD0JB821JtRH4WhY9iwJnhjSWJ9Cq1N0RAUGuoryxoJSG9IU19VaTg57LDQAdvbkC0FFJH6U6BIUbWBopgCRv01/NCpJNan0NoUqSn+wDFkSPig9PTTjVuw3MyyrpFo1iFNfV3WAQWqCx8KKtKMFChSVGugKFY4XGr58vCC6R62f/ubwkYSGq4XI4n1KbQ2RSzl3uSPOw7uvFO9DVlLu0ZiA/14mEMYyLqWGNIkm1QaUDo64Be/gK99rbLw4d7zuopxe1LyFFSknhQoUpR0oOiJejPqpzhkjBuX8sxTSaxPobUpStJQpOZQPKRpJUPZljc4kLmp1Uj0Nu2rCIT3itde6/28ww4L7yMPP1zZuZX0pOSDR28LwKtHRZKgQJGiNANFMQWM9NV9+FRXFxx6aJjwv1otvjaFhiI1p6yHNOWpRkKaQb6H5IQTKuvUrnePSp56VvomBYoUZRkoipWaXWrAgPBNxgsvKHDUQ12GT8Vdm2KXXeCJJ/r8q7GCQ9+haV9Fqrf33vD885WfX48elfxcIFBdAInTswLxAohCS+0UKFLUSIGiNwoc6ag5ZMRdmwKacirZSusYFByaW5ZDmvJUIyGSnDPOgAULquspqaZnJU5gAYWWJClQpKiZAkVvNISqfqqq0bi+hrUp8hpwKln1LLSmof27uX7IFI5fkd2Qpp5qJHI52H33MMqw3DoKvR0/5ZTw3z/6UekPMGZ6LZW+p62tusU1q+1ZqTaw5ANIpXOcZBFaIH4IySK8NHWgMLPBwPnAZ4DRwHJgJnCRuy+u8BrbAv8AHAtMAEYBa4H5wG3A1e6eyDtbXwoUxdSjkb5cDg5vX8Ldz44mV+tf0QSnki31QlYuFKlnoe+p5EP1Zz4TznnkkU1/Hz56SDdT7zmaAY9mN6Qpb93BHXzvyJn89pFBJd+Mi1/vqj0O5c856ST49rfLz7lQSSjJ5WCHHSobuiLSF1QbWDo7w3tSNbOwpxVaBg+OH0KyCC95TRsozGwQ8AAhBLwGPATsChwILAMmuPtLFVznW8DXCV9MPQ28AAwHDgEGAg8DR7r76gTa3GcDRW9KvXmOGxeOPf64ejlqkcQKwOstx9VfW8Skk0ZuMQWqwkBramsL21Jv0uVCQcUfqinxzvXGGzB3bop/whIaaBhgLaHkzDNhm216ngwu6WBS6VSqIo0glwvva9W8V6URWjo74aqr4oeQNMNLXj6IfPvbY1m1ak23+4Ktem9BQty95gfwLUII+B2wdcH+86L9syu8zvnAvwGji/a/H3g5utZlCbV5Xnt7u0tpS5a4T5vmPmmS+4gR7mbu4Z+jHj09BrHaH6Cj5gvNZLKPYEnmfx49kn20tYVHqWO5nPspp7hffLH75MnuEyaE7bRp4d9j/t9kqWOxrF7tftZZ4cZZ/2L69XP/8IcT+oM1rkr+H/Z0zurV7p2d5f+X5XLh+PLlvZ+3446V/a8ZNiz7vx566FHrI5dz79+/+ud89rPVPaezM/w7Puus6p+3erV7R0dl53d0hPPztnw5b3dod/faPy9X+qi5h8LMBgBLgWHA/u7+VNHxZ4APAOPc/Q813OdEwrCnhe6+Ww1Nzl+vZXso4lCvRuUG0c1VTOV0bqxp7LlmsWlOvfUYQO/fdtddEosyJkkrxlelkh6T3s6rpNckvxL1177We+/Kiy/CQw/13nb1nkgzqbYnJJcLIwgOOKC6tZJyOfjUp+DWWyt/Tv5ls/TLeRjp5D4vjXkzgASGPJnZR4BZwAJ336PE8YuAS4FL3P3iGu7TDswD1rn7wLjXKbieAkXCCt+8VqyAxYtD93w1/xj7khF08QTjGc2imq4zmw6OZgZraJxC7VZXyTCjhlI8IHfRovAPtBF0dMDMmZkPa2pVSYSTkSN7X+Mz64BS6RSsIrUaOxbifLSME14WLYKLLio1lKs5A8VU4HvAT939UyWOfwz4FfBLdz+hhvscA9wN/NXdx8S9TsH1FChSUK5AfO1aWLOm74eOeziCI4g5lSxoJeAMNV1oKKW3AblZaqAaCUlOIwaU/Hkf/3hlnXKHHRa2lQSaalQ7W5JIb049NfRqbFl/0pyB4rvAucD33P28Esf3JRRYP+nuB9Rwn/uAjwLfd/dzqnheucSwe3t7+0AFiuz15WlvL2Aa06hxKllgI8aDHMa9HKn59xOWy8F++4UPNKtXN2FogNJTghxyCNxzDzz6aNatg/HjYdttmzSVSZaSCihQeUjJrzxdSSF9pT0pHR3wy19WHmrq2aNS7Tfh0oyaM1BcC5wFTHP3LT45mdkewIvAi+6+Z8x7nA38AHgTGOvur1bxXAWKJtOXhk6NZAl/ZXSi8/irtqJ6faK3oZRG7oHI05AmaTCVhpRKzq0mpORXnK6056UePSodHfC+94X7V6ranpU4Q3fc+86MhP37N8KfRYGi1PUPA+4HcsAn3P2XNTV403U15KmJNHPISGIq2VJatbaip3Bw3HFw110ZFzynpdEKq4tpSJO0kGpCSqXnJ92jUnhuNVOaVtOzEiewxJnSNY3QErcn56Mfhfvvr/55cZRffLM5A0XdhjyZ2T7Ag8B2wDnu/v2aGrv5tRUomlyzzDw1iG5mcDQTSe6DXz6cLGIUP+ALfXIYVJ8YipSkRi6s1pAmkbpJskclr149K3ECy8yZ1T8njdBy8snwk59UP1tTnFme4oaXPfYIC5duKf1AUfO8s8BUwueb28sc/1h0/BdVXnc34NXoud9Ier5ctA5FSyic0/3AA91HjSq/DkA9H4NY7T+k09dSnzn/15LzH9LpA+lO/c9W6yOXcx8/3v3ww/v8MgTxNNJ6EaUeHR3u3d1Z/5ZEJIZq17ep9PxK10wpfOmo5jnVrtkQ9zlx1pNwr/55J59c/Ut8Luf+ta+VO96c61AkPm2sme1IWG17d2C6u0+tqZGl76EeihaV5fCpEXTxI07mCOrTH7qc7XiR97OCYcxmYkP1XKjHoQLFvRBDhsCCBbBwYdYt25KGNIlIL6odClbNc6rtZYnznGpXvM6Xi8V53jnnVL+S9yWXwOjRpf4szTnkqXBhuw+6+9NFx6ta2M7MtgPmAH8P3Aic6bU2svR9FCjkXWmGjHoUapdTzwLucgGhpeoYktLoxdUa0iQiDaqeoQXiBZc4z4sbXjo7+8g6FABm9i3g68DvgCPcfVW0/zzgSmCOu08sOP9LwJeGk+thAAAajElEQVQIa1OcX7B/K0IB9sHA7cBJ7v5OzQ0s3WYFCulRTzUac+du+gL5rbfgqaeqK/aqV6F2ObX0XKhnoc4avbhaszSJiMQKLtU+L26vS59YKRvAzAYBs4GDgNcIw5XGRD8vAya4+0sF518MfAO42d1PK9j/PUJNxjvAT6D0V7iFz6mhzQoUkphKwkfhvmd/380lTxzNASuz+RC5oS3HPTuezndGTWfgsEHqWUhTIxdXF9KQJhGRTFQbXrYMIk0aKADMbDBwPnASsAuwHJgJXOTurxSdezGlA8VNwKm93cvda/4FKVBI5qJXAL/xRiyrYS5Dh4bhLJMmKTnUW6MPaxo1CnbZRUlSRKRJ5YPIZZeNZdWqNd3uC7ZK696JBYpmo0AhDSP/CnDPPWFloqz+Teob6eQ0U3E1aFiTiEgfMnbsWObPnz/f3cemdU8FCgUKaSSlq6vStd128P73w7Bh+qa6Wo3eC1FMIVJEpM/JIlD0T+tGIlKB6dPhhReyLdB9442wMiDAfffBxRfrQ2cpzdQLcfDBcOSRWy4trrAoIiIJUKAQaSSDB8OMGT1P85C29evh2mvhpz9VzwU0Vy+EeiBERCQFGvKkIU/SqPK1FVdfrVmAstJMvRCgwmoREVENRZoUKKRpNPo6BX2x5qKZeiHyVFgtIiKohkJESmnEYVCFStVc7LsvbLNNc6yE12y9EMVaoadIREQamnoo1EMhzaRwtZsVK+DFF8MH+maQy2UbNEoFhxUr4Omnq1vmPEu77RYezRDUREQkExrylCIFCukTtlwes/mUCxrHHQd33rkpABR+gIbNw0FPzznkkBC8br+9uX9H6oUQEZEKKFCkSIFC+pRm7rmoVltb2G7cmG076knF1SIiEpNqKEQknpEj4YILwgP6Rs9FOX05SICKq0VEpOm0Zd0AEamDwYPhmmtg0SKYNg0mTQrfdkvjyuXCSukKEyIi0mTUQyHSlxX2XPTlXotmo+JqERHpQxQoRFpFvtfi0ktbp96i0ai4WkRE+iAFCpFW00r1FllSL4SIiLQIBQqRVqeei2SpF0JERFqMAoWIBOq5qJ56IURERBQoRKQM9VyUp14IERGRdylQiEjPWr3nQr0QIiIiPVKgEJHqlOq5WLkShgyBt96Cp56CDRuybmXt1AshIiJSEQUKEYmnuOcir6uruYJGLgf77Rd6H9QLISIiUjUFChFJVr2DRltb2G7cWPlzcjn4zGdg993hkUfC/RUcREREEqFAISLpqDRo5D/oH3cc3HXXlvvPPDM8r9rnKDSIiIjUhbl71m3IhJnNa29vb583b17WTRERERERScTYsWOZP3/+fHcfm9Y929K6kYiIiIiI9D0KFCIiIiIiEpsChYiIiIiIxKZAISIiIiIisSlQiIiIiIhIbAoUIiIiIiISmwKFiIiIiIjEpkAhIiIiIiKxKVCIiIiIiEhsChQiIiIiIhKbAoWIiIiIiMSmQCEiIiIiIrEpUIiIiIiISGwKFCIiIiIiEpsChYiIiIiIxKZAISIiIiIisSlQiIiIiIhIbAoUIiIiIiISmwKFiIiIiIjEpkAhIiIiIiKxKVCIiIiIiEhsChQiIiIiIhKbAoWIiIiIiMSmQCEiIiIiIrEpUIiIiIiISGwKFCIiIiIiEpsChYiIiIiIxKZAISIiIiIisSlQiIiIiIhIbAoUIiIiIiISmwKFiIiIiIjEpkAhIiIiIiKxKVCIiIiIiEhsChQiIiIiIhKbAoWIiIiIiMSmQCEiIiIiIrEpUIiIiIiISGwKFCIiIiIiEpsChYiIiIiIxKZAISIiIiIisSlQiIiIiIhIbAoUIiIiIiISmwKFiIiIiIjEpkAhIiIiIiKxKVCIiIiIiEhsChQiIiIiIhKbAoWIiIiIiMSmQCEiIiIiIrEpUIiIiIiISGwKFCIiIiIiEpsChYiIiIiIxJZYoDCzwWZ2qZm9YGZrzOxVM/tvMxsV41rbmdl0M3vZzNZG26vMbNuk2isiIiIiIrVLJFCY2SBgFnARsDVwJ7AIOB14yszeV8W1dgAeB84BNgB3ACuBKcDvzew9SbRZRERERERql1QPxYXABOBRYE93/7S7HwR8GRgO/HcV17oK2AP4BbBXdK19gO8DewLfTajNIiIiIiJSo5oDhZkNAL4U/fhFd387f8zdvws8C3SY2QEVXGtH4ERgHfAFd99QcPirwDLgZDMbUWu7RURERESkdkn0UBwCDAMWuPtTJY7/LNoeW8G1jora9JC7dxUecPe1wN1AP+Af4jdXRERERESSkkSg2DfaPlnmeH7/B1K+loiIiIiI1Fn/BK4xOtq+UuZ4fv+YlK8FgJnNK3No7wULFjB27NhKLyUiIiIi0tAWLFgAsEua90wiUGwdbVeXOb4q2g5N+Vq9aVu7du3G+fPnP5/AtaRv2z3aLsi0FdIs9PdFKqW/K1IN/X2RSu0NDE7zhkkEiobm7iW7IPI9F+WOi+Tp74pUQ39fpFL6uyLV0N8XqVQPo3PqJokaivysTluVOT4k2q5M+VoiIiIiIlJnSQSKv0bbncscz+9/OeVriYiIiIhInSURKJ6JtvuXOZ7f/2zK1xIRERERkTpLIlA8AqwAdjez/Uoc/2S0vbuCa80ENgKHFS9eZ2YDCWtZvAP8Jn5zRUREREQkKTUHCndfB/xn9ON/mVm+zgEzO4+wZsQcd/9Dwf4vmdnzZvbtomu9BvwYGABcbWaFRePfAYYDt7j70lrbLSIiIiIitTN3r/0iZoOA2cBBwGvAQ4S1Ig4ClgET3P2lgvMvBr4B3OzupxVdawfgMcL0aAuAucBYYB/gxehay2tutIiIiIiI1CyJIU+4+xrgI8A3CWtIHE8IFDcB+xeGiQqu9TpwIPB9Qk/Fx4FhwH8ABypMiIiIiIg0jkR6KEREREREpDUl0kMhIiIiIiKtSYFCRERERERiU6AQEREREZHYFChERERERCQ2BQoREREREYlNgUJERERERGJrqUBhZrPNzHt4HJV1G6WxmNlwM7vCzP5kZt1mttzMnjSzf8+6bdIYzGxiL68r+ce/Zt1WaQxmNt7MbjezV81svZm9aWYPmdnpZmZZt08ai5n9nZndamavmdlaM1toZv8ZLQQsLcTMDjCzr5nZL8zslfz7SwXPO83MHjezt6PPMb8xsw8l2rZWWofCzGYDHcDPgbdLnHKluz+XaqOkYZnZAcA9wPbAPOCPwDZAO7Czu/fPsHnSIMxsb+BrZQ73A06O/vtwd38gnVZJozKzTwA/IfzdeBL4MzAcOAzoD9zm7p/NroXSSMzscOBuYCvgeWA+sA+wJ/AKcLC7v5JdCyVNZnYH8I/F+9297BcRZnYVMAXoBu4FBgGTAAM+6e53JNK2Fg0Uu7n7wmxbI43MzIYTXri3Ak5097uKjh/o7o9n0jhpGmZ2NPAbYBEwxlvpBVe2YGb9gcXACOCz7n5bwbG/Ax4G3oPCpwBmthXwEjASuNTdvxHtN+A7wFeAe939yOxaKWkys38BhgBPRI+FwMBygcLMPgrcB/yNED5fjPYfDMwGVhM+E79Za9taasiTSBUuAXYAvlocJgAUJqRC+d6JWxUmBNibECb+VBgmANz9/wG3RD+OT7th0pBOIISJPxHekwCIXksuIHyYPMLM9s2kdZI6d/83d/9Xd7/b3ZdU8JTzou238mEius6jwA+BbYEzk2ibAoVIETMbTPgguAq4MePmSJMysyFs6pr+UZZtkYaxtsLz/lbXVkizOCDaPujuGwsPuPt64JHoxy2GwIhEn2UOj378WYlT8vuOTeJ+rToG/Ewz2x7YCLwA3OHuf824TdI4xgFDgYfdvTsatjKZMO7wBeB2d381ywZKUziB0DX9lLvPz7ox0hBeAhYAe5nZSSWGPJ0MvAH8MqP2SWMZEm3fKHM8HzzVQyGl7AUMBJaVqbN5Mtp+IImbtWqguLDo5yvM7Jvu/s1MWiONpj3aLi1TAHWZmZ3p7j9OuV3SXPLDndQ7IQC4+ztmdirwK+BWM/sy8CJhGNRhhLqt09x9eYbNlMaxLNqOKXN8t16OS2sbHW1LFu27+yozexPYzsyGuvvKWm7WakOeHgROAXYnFNvuBXwd2ABcamZTMmybNI7tou1xwFHAFwlv+LsCVwCDgZvNbL9MWicNz8x2JMyi8Q6g4CnvcvdHCJODvATsD3wa+Aihx/y+aL8IhM8sAB8rniLWzEYRes4h9KiLFNs62q7u4ZxV0bbmv0MtFSiiQpZb3P0ld+929xfc/TLg+OiUi6MxZ9La8v8u+gP/6u5Xu/syd3/Z3b8K/BTIAV/NrIXS6E4kTAt6X4WFc9IizOxE4HHCzF8HEd709wRuAr4MzDKzgZk1UBrJvYRhKVsDM8zsQDPbOpqhZwabRplsLHcBkbS0VKAox93vBeYSqt0Pyrg5kr3CNUpKFWXn93Wk0BZpThruJFsws/cDNwOvA8e4++PuvsrdX3T3zxOGQu0PnJFlO6UxRLM5nUBYB2kc8HtgJfA7Qq/5xdGp5WospLXlP8ts1cM5+TqdmoY7gQJFofx0Wjtm2gppBC9H29XuvqzE8YXRdkQ6zZFmEhXXfpDwYp7IgkHSZ3yG0Ls5091LLa56e7T9cHpNkkbm7i8D+wGfAq4CrgHOIdT6dUWnzcumddLg8pMN7VzqYDQT4bbAG7XWT0DrFmWXkh83v6rHs6QVPBVtB5vZQHcvnurxPdG21AcCkVOi7S/cvaexq9J68m/sK8ocz+/frsxxaUHuvoEw1PanhfvN7EPRf85Ou03SFP5EmKp6uJmNcvfFRcf3j7bPJnEz9VDw7qrIh0U/PtnTudL3RVMIP0NYlr7UsKb8vqdKHJMWFq1ge1L0o4Y7SbF8Pc24MsfzC9otrH9TpJmZ2XuBTxKmjv1Fxs2RBuTu3cCs6Md/KnHKJ6Pt3Uncr2UChZl9yMyON7N+Rft3Jcz5PQS4q8xcvdJ6vhNtr4hm7AEgmtnpy9GPP0y9VdLoDiNM4biYTS/kInl3RtsPm9k/Fx4wswnAudGPpRahkhZkZvuY2aCifTsT/i4NBb4cfXAUKeW70fbCqIYLgKiw//PAm8ANSdyolYY87Ukopl1iZk8SfoljCCtRDiKMQTwru+ZJI3H328zsCOBUYL6Z/Y4wXeyHCAvFXOfuP+3pGtKS8sXYtxWvbCvi7k+a2RXAV4CrzeyLhLUndgIOJnzJd627359hM6WxfAX4ePS55TVC7d6hhPehb7r7zVk2TtJlZh8DLirYNSDa/1jBvm+6+68B3P1+M5sOTAGeNrP7oudMJozCON3d30ykbWESgb4vKpT834RZnHYhjFFdBfw/wrjEHyjlS6Fo+MrnCCn+7wAnjDW8Ri/iUiya6vM1wmvLvu6eyLhU6XvM7OPA2YQvtIYRZlh5mvBFhdYtkXeZ2fGEvyv7AtsTZnR6FLjK3Wdn2DTJgJmdRunZJwud7u43lXjelwifZdYBjxGCx+8Sa1urBAoREREREUley9RQiIiIiIhI8hQoREREREQkNgUKERERERGJTYFCRERERERiU6AQEREREZHYFChERERERCQ2BQoREREREYlNgUJERERERGJToBARERERkdgUKEREREREJDYFChERERERiU2BQkSkRZjZQjPzrNuRJDObbWZuZrsmfN3/FV33Y0leNwlmdnzUtk9l3RYREVCgEBGRBpZFCDKzQcC3gCfc/ddp3rtCdwLPAJeZWS7rxoiIKFCIiIhs7p+BXYBvZ92QUtzdgcuB3YHPZdwcEREFChERkSL/DCwHfpV1Q3pwJ7ASODvrhoiIKFCIiBCGuZjZGjNbWOLYHdGY9YdLHJtrZhvNbHjBvsPM7D/N7Fkze8PMus3seTO73My2LXr+CdG1f9JD266MzjmnaP9WZna+mT1lZm9Hj8fM7NQYf/5dojYviH4Py83sV2b2oRLnTozac5OZvcfMfmBmr5nZWjP7o5md0cN9TojauNrMXjezn5rZHmZ2cXTN0wrvAYyJfvaCx8Iy1z4+uvaqqP0/NrOdq/w9dADvB37u7utLHO/p/qdFxy8u2n9TtH+imX3UzB40s5VmttTMrjOzYdF5I8zsGjNbHP0/eNzMJpa6l7t3A3cAHzCzg6r5M4qIJE2BQkQEcPc1wO+BMYUFvmbWBnw4+nG8mW1VcGwY8EFgvrsvK7jcvwNnAt3Ab6PHNsC/AA+b2dYF5/4aWAEcW7S/8P6fAd4B/m/B/hHAo8BlwHuBOcCDwN7ATWb2/Ur/7GZ2MGFM/heB9VGb/ggcCTxoZp8u89RtozYcBzwEPBLd/wYz22IojplNAX4OjCf8ru8DDgAeB3YrOn0JcDOwKvr55oLHz0q05QvR/m7gN8DbhN/bLDMb3OMvYHPHRNvZVTynUh8HZgIWbdcShizdaWY7EH6XRxJ+l08Tfk8zzezvy1wv38aGKxwXkRbj7nrooYceergDXAI4cFrBvg9G+/4YbT9acOzYaN9/Fl3naGBY0b6BwDXR+f9adOz6aP8pJdo0KTo2o2j/r6P9VwEDC/aPBJ6Ijh1V9JyFREPwC/ZtA7wKbAA+W3RsHGHoz0pgeMH+idH1Hfhx0f2Pj/a/XHSt9xE+QK8FPlKwvz/w3wXXO623Nhcdnx09bxVwcMH+rQgBx4Ezqvg78Fj0nPeXOe7AwjLHTouOX1y0/6Zo/zvAxwr2DwWei47NA34E5AqOfzM6dnOZ++0THZ+T9b8dPfTQo7Uf6qEQEdlkdrSdWLAv/9+X9nBsTuFF3H2Gu68o2rcWmEr44P6PRfe9Jdp+tkSb8vtuze8ws/2AfyAEh/Oia+fv0wV0Rj/+c4nrFTsD2BG4yt1vLTzg7nMJH2q3Bk4u8dy3gC8V3f8OQvgaXTSV6xnAAOBH7v5AwfkbgPMIPQq1+J67P1pw3dXAd6MfP1z6KSV9ANgI/LnG9pRymxfMGuXuK4Hroh93Bs7xzYdZXUEIDB1lrvd8tN0v6YaKiFSjf9YNEBFpII8RvkGfWLBvIuEb+p8DL5c4BiWGx5jZKEIPxt6EXoD8FzjrCGP0Cz0IvAJMMrMR7r40usYg4BOEb99/WXD+EdH2DnffWHxvd3/KzN4GDizz5yyUv9Yvyhx/KNqWutYf3P1vJfa/QPj2fEdCDwPAIdH2pyXa+6aZ3QucUEF7y7m3TDuI2tGraMjZYGC5u9djqtpSbXwp2s519zcKD7j7CjNbTpn2u/sGM1sJbGNmA9x9XbLNFRGpjAKFiEjE3bvN7HHgsOjb9b8ChwEPufs7ZjYbODGqo8gRvhkurp/AzM4jTOtZ0RoB7r7RzH4MfBX4NJCvfziGEEZuc/dVBU/ZNdpOM7NpPVx6UAW3z1/rETPr6bwdSux7pcy5K6PtwIJ9+Q/Fi8o856893bwCpdpSqh09GVb0vKQtLrHv7R6O5Y9v38M13yIMndoWWBq/aSIi8SlQiIhsbjYhREwkFCpvx6YeiNnAqcCHCN9kt1E03MnMJgBXEgqtp0TPWZIfFmRmr1L6G+dbCIHiJDYFii2GO0XyvR0PAwuq+LOVkr/Wz9hUAF3K8yX2bdE7kqEk2pIfpjY0gWuV0lMb47Y/H4LejPl8EZGaKVCIiGxuDnARIVBsF+2bXbSdSAgUhfvyPh5tv+7uNxceiGYbem+pm7r7s2b2R2CCmb0PeINQJ7GMLYfK5L+Nv8Pdr+zlz9ObV4C9gMvd/Q81Xqsnr0X32QWYX+L4LnW8d0Xc/W0z6wa2NbO2UsPJIuV6fkr14tSNhVWytwbe0nAnEcmSirJFRDb3O0Kdw8To8RbwJIC7L2RTHcXE6Pw5mz/93RBSagjOPxGmDC0n3xNxEvBJQhHzT6LC5UL3RduPU7skr9WTR6LtJ4oPRNPvHlG8P7IuOietL8CeIbw37tHDOSPya0cU+WB9mlTW3tH26ZTvKyKyGQUKEZECHhYMe5ywoNoRRPUTBafMJhQo7wc8H82qVChfCHxm9A0yAGbWDvxbL7e/jTCrz0mUH+6Eu+fXcDjEzP7LzLYpPsfM9jWzo3q5H4SpbJcC/8fMOqN1Lwqv09/MjjSzfSq4Vk9uJISD/2Vm7866ZGb9CEPEyg0zejXa7lXj/SuVL0If38M5BlxqBUUnZvYR4FPRj0Pq1LZi+UL54lArIpIqBQoRkS3lP6ANYsshTbMJxdZtJY5B+OC8hDDD05/M7Cdmdh/hW+SHCD0cJbn7X6Nz/o4wVegCd3+szOknA08RFnR72cweMLNbo9Wt/xrdr9dA4e5vEqaxXUEIFwvN7DfRtX5LGHI1k56/se+Vuy8A/g+hQPoBM5sVFaK/QOi1yE+dWzx0565o+9to5evrzezyWtrSi/y0rhN7OGcFYc2J+WZ2u5k9AtxP+P8BIUz+tm4t3GRitP11TyeJiNSbAoWIyJZml/nv3o4RTaM6ntDbMICwivQoQl3GiRXc+9Yy/118n6WE4vBzCDUJHyQMk/oAYSrSrxLWMehVFFr+HvgOYYhXB2GBujGEcHUa4QNzTdx9etTGucAEwqrQTwMHAWui04qnof0P4FuE2Y4+QViB/DO1tqWHNs4hCjlmNqDMaW8SFi9cQQiOewDTCb+3nxGCaF2HaEX1OMcDz0Y9ViIimbH6TLUtIiJSmWjY07OEnpmd3H1Jxu2ZQliB/JPu/vOiY/lVwHfNom0F7TiREFq/4O4/yLItIiIKFCIikgoz2x34WzTMKr9vIHAZYbXs+919clbty4sWFPwTsNTdxxcdyzxQRLUbTxFmeGrXDE8ikjVNGysiImn5J+ASM/sDYYG7bYB9CetyvA58KcO2vcvd15jZRcDNZnaMu/8q6zYV+UfC7+3TChMi0gjUQyEiIqkws/GEnogJwHDCl1qLgXuAb7t7uVW0G0Yj9FCIiDQaBQoREREREYlNszyJiIiIiEhsChQiIiIiIhKbAoWIiIiIiMSmQCEiIiIiIrEpUIiIiIiISGwKFCIiIiIiEpsChYiIiIiIxKZAISIiIiIisSlQiIiIiIhIbAoUIiIiIiISmwKFiIiIiIjEpkAhIiIiIiKxKVCIiIiIiEhs/x/jCfQXK+UvGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "wl = []\n",
    "Rs = []\n",
    "Ts = []\n",
    "for i in range(nfreq):\n",
    "    wl = np.append(wl, 1 / flux_freqs[i])\n",
    "    Rs = np.append(Rs, -bend_refl_flux[i] / straight_tran_flux[i])\n",
    "    Ts = np.append(Ts, bend_tran_flux[i] / straight_tran_flux[i])\n",
    "\n",
    "if mp.am_master():\n",
    "    plt.figure(dpi=150)\n",
    "    plt.plot(wl, Rs, \"bo-\", label=\"reflectance\")\n",
    "    plt.plot(wl, Ts, \"ro-\", label=\"transmittance\")\n",
    "    plt.plot(wl, 1 - Rs - Ts, \"go-\", label=\"loss\")\n",
    "    plt.axis([5.0, 10.0, 0, 1])\n",
    "    plt.xlabel(\"wavelength (μm)\")\n",
    "    plt.legend(loc=\"upper right\")\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}