simulation-overview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "33592e1c-344d-46ea-a22b-996b4c08d820",
   "metadata": {},
   "source": [
    "# About This Notebook\n",
    "\n",
    "The purpose of this notebook is to break up the Smith Purcell simulation in meep step by step.  There may have been modifications to the main .py file for the simulator since this notebook was written, but the core simulation routine should be the same.  \n",
    "\n",
    "This is hopefully helpful in familiarizing a new meep user with how such a simulation works, and inspire them on how to modify the code for their own purposes if needed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29e2177e-097a-47a2-ae7c-8e836e649cc7",
   "metadata": {},
   "source": [
    "## Headers, Functions and Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "100e4e44-fb9d-4547-857f-5cb1ec235d39",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Headers\n",
    "import meep as mp\n",
    "import numpy as np\n",
    "import h5py\n",
    "from mpi4py import MPI\n",
    "import meep.materials as meepmat\n",
    "\n",
    "#Constants\n",
    "# -- Units noted.  \n",
    "# -- Fundamental unit length of 1-um assumed for meep calculations (see meep manual)\n",
    "tcon = 3.33564 #time conversion to meep units (fs/meep time)\n",
    "q =  1.6e-19 # electron charge in coulomb\n",
    "m =  9.10938356e-31 # electron mass in kg \n",
    "\n",
    "#Functions\n",
    "\n",
    "#Functions for calculating index of refraction -- meep style\n",
    "def calc_sig_d(n, k, fcen):\n",
    "    eps = (n + 1j*k)**2\n",
    "    eps_r = np.real(eps)\n",
    "    return 2*np.pi*fcen*np.imag(eps)/eps_r\n",
    "\n",
    "def calc_eps_r(n, k):\n",
    "    eps = (n + 1j*k)**2\n",
    "    return np.real(eps)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04345446-6f7d-4107-bb13-a2a7adc6ce8b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Settings for our Calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b76e41cb-b058-4464-8be8-5fea6f4f819e",
   "metadata": {},
   "source": [
    "### Output File Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "58297573-a705-4e1d-9fbb-64dc39fbff60",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_prefix = 'test_sim'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d736f00f-3fff-4b99-bdee-afe001667f6c",
   "metadata": {},
   "source": [
    "### Grating Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "b569da95-67c9-410c-8ad6-2f0ceab5ddfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Grating Structure Properties\n",
    "grating_N = 10 #No. of grating elements\n",
    "grating_period = 0.1 #grating period in um\n",
    "grating_height = 0.02 #grating height in um\n",
    "grating_fill_factor = 0.5 #fill-factor from 0 to 1\n",
    "\n",
    "#If you want another grating under the underlayer...\n",
    "lower_grating_height = 0.02\n",
    "\n",
    "#Grating Material Properties\n",
    "grating_material_eps = 13.5\n",
    "want_grating_material_by_name = True\n",
    "grating_material_name = 'Au'\n",
    "\n",
    "#Grating fill material properties\n",
    "# -- This is the material between each \"tooth\" of the grating\n",
    "grating_fill_material_eps = 1.0\n",
    "want_grating_fill_material_by_name = False\n",
    "grating_fill_material_name = 'Au'\n",
    "\n",
    "#This underlayer is a layer of the grating type of whatever thickness. This is more for traditional\n",
    "#metalic gratings on-top of a planar film.  \n",
    "# -- Can set to 0 thickness if you do not want it\n",
    "underlayer_thick = 0.03\n",
    "want_underlayer_material_by_name = True\n",
    "underlayer_material_eps = 1.0\n",
    "underlayer_material_name = 'Au'\n",
    "\n",
    "\n",
    "#Substrate Material (semi-infinite film underneath grating):\n",
    "substrate_material_eps = 1.55**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58dbe6aa-be56-42fc-8d79-dd7f3e0064cd",
   "metadata": {},
   "source": [
    "### Electron Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "e10dabd6-8913-498f-9fe5-958306594c61",
   "metadata": {},
   "outputs": [],
   "source": [
    "e_height = 0.01 #Height of electron above grating surface (um)\n",
    "W_e = 4e3 # Kinetic energy of electron (eV)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81c6b95b-c834-4e41-a51a-dd978f17716a",
   "metadata": {},
   "source": [
    "### FDTD-Specific Settings\n",
    "\n",
    "These are important simulation parameters for running the FDTD simulation.\n",
    "\n",
    "There are six basic components to every simulation.  I find that if I just note to take care of these three components one by one, I create simulations very quickly without any major issues.\n",
    "\n",
    "1. The cell.  This is the world the simulation runs in.  I like to define the cell as being composed of an \"active region\" and a \"boundary region\".  The active region is defined by sizes `sx`, `sy`, and `sz`.  For a two-dimensional simulation you can simply make `sz = 0`.  The full cell is then defined by `sX`, `sY`, `sZ`. The full cell is typically slightly enlarged as you often place PMLs around the boundary of the active region. Note these are personal convention -- meep sees everything as being inside one cell.\n",
    "2. Geometry. This is the definition of the optical properties of each voxel inside the simulation.  There are convenience shapes (blocks, spheres, cones, etc.) that can be used to help you define this.  There is also a way to simply pass your own functions to define optical properties at each voxel.   \n",
    "3. Sources.  Inside the active region you typically place sources and monitors (see point 3).  The sources are essentially all comprised of some oscillating current either at a point, over a line, or over a sheet.  By manipulating these sources with phase-profiles you can create plane-waves, dipoles, etc.  \n",
    "4. Monitors.  The monitors do exactly as they say -- they monitor the fields over a certain region in space.  These can be at a point, over a line, over a sheet, or over a volume.  You can place as many as you like, and use them to monitor fields of any type and in any specified polarization.  You choose this at the output.  \n",
    "5. Fluxes. This is a compliment to monitors.  You can technically monitor the fields and use these to calculate power spectral densities, but meep has these flux components to do this for you and in a more sophisticated and accurate way.  You specify a sheet over which to monitor the flux.  This is essentially the power spectral density through a specified region calculated internally using the pointing vector.  \n",
    "6. The run sequence.  This is the final step.  The simulator is calculating a lot of data on each iterative step of the FDTD process.  The run sequence allows you to tell it what you want to do at each step, how you want the simulation to be run, and what data to actually pull out and save.  This is extremely useful as it gives you the opportunity to tell the code what is worth keeping and what is not so that you don't save things you do not need.  You can also use this to do hack the simulator at each simulation step and do extra calculations beyond those typically performed in FDTD analysis.  \n",
    "\n",
    "My strategy is to typically draw out my simulation by hand on paper.  I define my geometry using the basic provided shapes (or function if needed).  I then place my sources, monitors, and fluxes where they should go depending on what I want to calculate and where.  Finally, I think of my boundary conditions and where I would like to place PMLs, etc.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "9ce91aef-88aa-4f1e-8890-a31749135e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Simulation Parameters\n",
    "# -- These parameters are important to running the simulation\n",
    "\n",
    "#Is this the reference simulation?\n",
    "reference = False\n",
    "\n",
    "#PML Thicknesses\n",
    "dpml = 0.5;\n",
    "\n",
    "#Monitor settings\n",
    "monitor_width = 5 #Width in no. of gratings\n",
    "monitor_height = 1.0 #height as a ratio of cell size\n",
    "\n",
    "# Sampling time for monitors in fs\n",
    "t_sample = 0.5 \n",
    "\n",
    "#EPS Averaging in Meep Calculation\n",
    "avgeps = True\n",
    "\n",
    "#Source Terms\n",
    "src_offset = 0.001 # offset of source to material edge\n",
    "\n",
    "#Perfectly Matched Layers\n",
    "dpml = 1.0 #PML thickness around simulation\n",
    "\n",
    "#Cell height (without PMLs)\n",
    "sy = 5 # size of y-dimension\n",
    "\n",
    "#Spatial Resolution\n",
    "res = 200 # Resolution: points per meep-unit\n",
    "        # e.g. if using 1 meep unit = 1 um,\n",
    "        # 1/res um spacing\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "558935fb-4037-4af2-a76b-b4865cd2c47d",
   "metadata": {},
   "source": [
    "### Frequency Output Range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "309c7dff-17e6-4efe-9aa6-878d6b1e75ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Frequency range for calculations\n",
    "fcen  = 2.0 #Central frequency (inverse wavelength in um)\n",
    "nfreq  = 1000 #Number of frequency samples\n",
    "df  = 1.9*fcen #Bandwidth of simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a6e0603-8914-4f02-a1c2-a0b7d2fb8732",
   "metadata": {},
   "source": [
    "## Calculated Parameters\n",
    "\n",
    "These are remaining parameters that are calculated based on user inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce3865ac-0839-4509-8e3f-049c12efdd51",
   "metadata": {},
   "source": [
    "### Calculate Remaining Cell Sizes\n",
    "\n",
    "We want to calculate the active cell width `sx` as it depends on the number of grating periods we include in our simulation and it is all well-defined.  We will also then calculate the total cell sizes `SX`, `SY`.  The z-sizes are all 0 for our 2D simulation here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "4f2057a3-0c75-4a49-bf4b-e206c06475ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Determine overall simulation cell size\n",
    "sx  = grating_period*grating_N*2.0 #Width of the cell in um\n",
    "\n",
    "#Full simulation cell sizes accounting for PMLs\n",
    "sX = 2.0*dpml+ sx\n",
    "sY = 2.0*dpml + sy\n",
    "sZ = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da41bda9-78e9-4f8a-b22f-84ef300c67a3",
   "metadata": {},
   "source": [
    "### Electron Velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "id": "a38768b2-0e0b-4295-aac4-88d53d0667d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_e = 1e6*1e-15*tcon*np.sqrt(2.0*q*W_e/m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4014e243-352b-4aed-84c6-ca6948361a32",
   "metadata": {},
   "source": [
    "## Setup the FDTD Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f35c8b5-76e6-4bee-95d2-adbeee9c622c",
   "metadata": {},
   "source": [
    "### Define Materials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "id": "44ce0530-1097-442f-960e-f5fd0d685f14",
   "metadata": {},
   "outputs": [],
   "source": [
    "if(want_grating_material_by_name):\n",
    "    grating_material = eval('meepmat.' + grating_material_name)\n",
    "else:\n",
    "    grating_material = mp.Medium(epsilon=grating_material_eps)\n",
    "    \n",
    "if(want_grating_fill_material_by_name):\n",
    "    grating_fill_material = eval('meepmat.' + grating_fill_material_name)\n",
    "else:\n",
    "    grating_fill_material = mp.Medium(epsilon=grating_fill_material_eps)\n",
    "    \n",
    "if(want_underlayer_material_by_name):\n",
    "    underlayer_material = eval('meepmat.' + underlayer_material_name)\n",
    "else:\n",
    "    underlayer_material = mp.Medium(epsilon=underlayer_material_eps)\n",
    "\n",
    "\n",
    "substrate_material = mp.Medium(epsilon=substrate_material_eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65d1b160-7b3e-4ed7-8f7f-dda7952fa6c7",
   "metadata": {},
   "source": [
    "### Define the Geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "0a8189f2-71f1-4005-a0e7-58dd52de93e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -- Cell --\n",
    "cell_size = mp.Vector3(sX, sY, sZ)\n",
    "\n",
    "# -- Geometry --\n",
    "geometry = []\n",
    "\n",
    "#Substrate Block\n",
    "geometry.append(mp.Block(center=mp.Vector3(0.0,\n",
    "                                           float(-0.25*sY),\n",
    "                                           0.0),\n",
    "                         size=mp.Vector3(mp.inf,\n",
    "                                         float(0.5*sY),\n",
    "                                         mp.inf),\n",
    "                         material=substrate_material))\n",
    "\n",
    "#Grating Fill Block\n",
    "geometry.append(mp.Block(center=mp.Vector3(0.0,\n",
    "                                           float(0.5*grating_height),\n",
    "                                           0.0),\n",
    "                         size=mp.Vector3(mp.inf,\n",
    "                                         float(grating_height),\n",
    "                                         mp.inf),\n",
    "                         material=grating_fill_material))\n",
    "\n",
    "#Underlayer Block\n",
    "geometry.append(mp.Block(center=mp.Vector3(0.0,\n",
    "                                           float(-0.5*underlayer_thick),\n",
    "                                           0.0),\n",
    "                         size=mp.Vector3(mp.inf,\n",
    "                                         float(underlayer_thick),\n",
    "                                         mp.inf),\n",
    "                         material=underlayer_material))\n",
    "\n",
    "\n",
    "if(not(reference)):\n",
    "\n",
    "    #Loop through and create all desired objects...\n",
    "    for i in range(grating_N):\n",
    "\n",
    "        #Upper Grating\n",
    "        geometry.append(mp.Block(center=mp.Vector3(float(grating_N*grating_period*-0.5 + grating_period*0.25 + i*grating_period),\n",
    "                                                   float(grating_height*0.5),\n",
    "                                                   0.0),\n",
    "                                 size=mp.Vector3(float(grating_period*grating_fill_factor),\n",
    "                                                 float(grating_height),\n",
    "                                                 mp.inf),\n",
    "                                 material=grating_material))\n",
    "        \n",
    "        #Lower Grating\n",
    "        geometry.append(mp.Block(center=mp.Vector3(float(grating_N*grating_period*-0.5 + grating_period*0.25 + i*grating_period),\n",
    "                                                   float(lower_grating_height*-0.5 - underlayer_thick),\n",
    "                                                   0.0),\n",
    "                                 size=mp.Vector3(float(grating_period*grating_fill_factor),\n",
    "                                                 float(lower_grating_height),\n",
    "                                                 mp.inf),\n",
    "                                 material=grating_material))       \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4491c54-1ccb-4c44-96c7-40cae01582d7",
   "metadata": {},
   "source": [
    "### Boundary Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "id": "74e396f6-f669-4da9-a13b-3b314708cfef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -- Periodic Bdr. Cond -- \n",
    "k_point = mp.Vector3(0, 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4befe65-20f9-462a-a838-d9265ce558f7",
   "metadata": {},
   "source": [
    "### Symmetries\n",
    "\n",
    "This simulation does not have symmetry.  However, you should note that if your simulation has symmetry, then you can speed up the simulation significantly by including it. \n",
    "\n",
    "This symmetry is NOT symmetry of the objects, but of the fields.  If it is inversion symmetric, then you can include a negative pre-factor to account for it.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32967455-db9c-4fd2-b0a8-65fda6f1fb2d",
   "metadata": {},
   "source": [
    "### Define the PML Layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "id": "5f180fd9-3b60-4b1b-aef5-ca8e18f3231b",
   "metadata": {},
   "outputs": [],
   "source": [
    " # -- PML Layers --\n",
    "pml_layers = [mp.PML(thickness=np.double(dpml), direction=mp.X),\n",
    "              mp.PML(thickness=np.double(dpml), direction=mp.Y),\n",
    "              mp.PML(thickness=np.double(dpml), direction=mp.Z)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b6da666-c06a-49c7-9ece-a258493b883e",
   "metadata": {},
   "source": [
    "### Define the Field Monitors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "d3666fa2-b49f-4305-9033-86ce2c812b62",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a \"volume\" monitor (technically 2D as z = 0) for fields:\n",
    "monitor_xy = mp.Volume(center=mp.Vector3(0, 0),\n",
    "                                   size=mp.Vector3(monitor_width*grating_period, monitor_height*sy, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a841f574-ad4a-452f-9659-27c9ea38e7b5",
   "metadata": {},
   "source": [
    "### Sources\n",
    "\n",
    "Here will will have a delta source that is a point current source with zero size.  This reflects a moving charge in space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "id": "a81eb0b8-f4a7-42a9-b4af-df39251bb5bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Y-Polarized Gaussian source next to the dpml at -0.5 sz\n",
    "sources = [mp.Source(mp.ContinuousSource(frequency=1e-10),\n",
    "                     component=mp.Ex,\n",
    "                     center=mp.Vector3(sx*-0.5, grating_height+e_height),\n",
    "                     size=mp.Vector3(0,0,0))]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a96e9f52-bdb5-42c0-a39f-62e6cd99dde2",
   "metadata": {},
   "source": [
    "### Define our Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "d50d3443-1910-4972-a94a-81989d1088b1",
   "metadata": {},
   "outputs": [],
   "source": [
    " # -- Load Settings Into Simulation -- \n",
    "sim = mp.Simulation(resolution=res,\n",
    "                    cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    sources=sources,\n",
    "                    k_point=k_point,\n",
    "                    geometry=geometry,\n",
    "                    dimensions=2,\n",
    "                    filename_prefix=save_prefix)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e934088d-3231-4813-9789-e9b17a52fc5d",
   "metadata": {},
   "source": [
    "### Define and Add the Flux Regions\n",
    "\n",
    "We define the flux regions and add them to our simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "2ed7c832-6465-4cd9-9a58-1d42820534df",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Air Side\n",
    "air_fr = mp.FluxRegion(center=mp.Vector3(0, sy*0.05, 0), \n",
    "                       size=mp.Vector3(sx, 0, 0),\n",
    "                       direction=mp.Y)\n",
    "air_flux = sim.add_flux(fcen, df, nfreq, air_fr)\n",
    "\n",
    "#Substrate Side\n",
    "substrate_fr = mp.FluxRegion(center=mp.Vector3(0, -1*sy*0.05, 0), \n",
    "                       size=mp.Vector3(sx, 0, 0),\n",
    "                       direction=mp.Y)\n",
    "substrate_flux = sim.add_flux(fcen, df, nfreq, substrate_fr)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2f12285-f102-4e73-b9f5-7beafab8fa4d",
   "metadata": {},
   "source": [
    "## Run the Simulation\n",
    "\n",
    "Now we have it all setup -- we just have to run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "id": "390c5e60-f214-4122-be7a-ce5bb7ceeea3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.00882316 s\n",
      "Working in 2D dimensions.\n",
      "Computational cell is 4 x 7 x 0 with resolution 200\n",
      "     block, center = (0,-1.75,0)\n",
      "          size (1e+20,3.5,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (2.4025,2.4025,2.4025)\n",
      "     block, center = (0,0.01,0)\n",
      "          size (1e+20,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (0,-0.015,0)\n",
      "          size (1e+20,0.03,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.475,0.01,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.475,-0.04,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.375,0.01,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.375,-0.04,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.275,0.01,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.275,-0.04,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     block, center = (-0.175,0.01,0)\n",
      "          size (0.05,0.02,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "     ...(+ 13 objects not shown)...\n",
      "time for set_epsilon = 4.47161 s\n",
      "lorentzian susceptibility: frequency=10.7433, gamma=1.78571\n",
      "lorentzian susceptibility: frequency=3.47141, gamma=2.01155\n",
      "lorentzian susceptibility: frequency=2.39466, gamma=0.701702\n",
      "lorentzian susceptibility: frequency=0.66944, gamma=0.278261\n",
      "lorentzian susceptibility: frequency=0.33472, gamma=0.19438\n",
      "drude susceptibility: frequency=1e-10, gamma=0.0427474\n",
      "-----------\n",
      "creating output file \"./test_sim-eps-000000000.h5\"...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dkeathley/anaconda3/envs/pymeep/lib/python3.9/site-packages/meep/simulation.py:1343: RuntimeWarning: Note: your sources include frequencies outside the range of validity of the material models. This is fine as long as you eventually only look at outputs (fluxes, resonant modes, etc.) at valid frequencies.\n",
      "  warnings.warn(warn_src, RuntimeWarning)\n",
      "/home/dkeathley/anaconda3/envs/pymeep/lib/python3.9/site-packages/meep/simulation.py:1347: RuntimeWarning: DFT frequency 0.10000000000000009 is out of material's range of 0.16131113692089302-4.0327458966810505\n",
      "  warnings.warn(warn_dft_fmt.format(dftf, min_freq, max_freq), RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "creating output file \"./test_sim-Ex_xy.h5\"...\n",
      "creating output file \"./test_sim-Ey_xy.h5\"...\n",
      "Meep progress: 0.1075/15.995224211855438 = 0.7% done in 4.1s, 599.8s to go\n",
      "on time step 46 (time=0.115), 0.0900326 s/step\n",
      "Meep progress: 0.2475/15.995224211855438 = 1.5% done in 8.1s, 516.3s to go\n",
      "on time step 104 (time=0.26), 0.0694881 s/step\n",
      "Meep progress: 0.3925/15.995224211855438 = 2.5% done in 12.1s, 482.9s to go\n",
      "on time step 162 (time=0.405), 0.0695776 s/step\n",
      "Meep progress: 0.5375/15.995224211855438 = 3.4% done in 16.2s, 465.3s to go\n",
      "on time step 220 (time=0.55), 0.0695529 s/step\n",
      "Meep progress: 0.675/15.995224211855438 = 4.2% done in 20.2s, 458.5s to go\n",
      "on time step 275 (time=0.6875), 0.0729669 s/step\n",
      "Meep progress: 0.8175/15.995224211855438 = 5.1% done in 24.2s, 450.1s to go\n",
      "on time step 331 (time=0.8275), 0.0737854 s/step\n",
      "Meep progress: 0.9525/15.995224211855438 = 6.0% done in 28.3s, 447.1s to go\n",
      "on time step 387 (time=0.9675), 0.0721615 s/step\n",
      "Meep progress: 1.0975/15.995224211855438 = 6.9% done in 32.3s, 438.7s to go\n",
      "on time step 445 (time=1.1125), 0.0691151 s/step\n",
      "Meep progress: 1.245/15.995224211855438 = 7.8% done in 36.4s, 430.8s to go\n",
      "on time step 504 (time=1.26), 0.0689113 s/step\n",
      "Meep progress: 1.3900000000000001/15.995224211855438 = 8.7% done in 40.4s, 424.2s to go\n",
      "on time step 563 (time=1.4075), 0.0686487 s/step\n",
      "Meep progress: 1.5375/15.995224211855438 = 9.6% done in 44.5s, 418.0s to go\n",
      "on time step 621 (time=1.5525), 0.0703532 s/step\n",
      "Meep progress: 1.6775/15.995224211855438 = 10.5% done in 48.5s, 413.9s to go\n",
      "on time step 678 (time=1.695), 0.0709821 s/step\n",
      "Meep progress: 1.825/15.995224211855438 = 11.4% done in 52.5s, 407.9s to go\n",
      "on time step 737 (time=1.8425), 0.0686307 s/step\n",
      "Meep progress: 1.9725000000000001/15.995224211855438 = 12.3% done in 56.6s, 402.4s to go\n",
      "on time step 796 (time=1.99), 0.0688966 s/step\n",
      "Meep progress: 2.1175/15.995224211855438 = 13.2% done in 60.6s, 397.2s to go\n",
      "on time step 854 (time=2.135), 0.069062 s/step\n",
      "Meep progress: 2.255/15.995224211855438 = 14.1% done in 64.7s, 394.0s to go\n",
      "on time step 907 (time=2.2675), 0.0757006 s/step\n",
      "Meep progress: 2.3925/15.995224211855438 = 15.0% done in 68.7s, 390.6s to go\n",
      "on time step 962 (time=2.405), 0.0731499 s/step\n",
      "Meep progress: 2.5175/15.995224211855438 = 15.7% done in 72.7s, 389.3s to go\n",
      "on time step 1013 (time=2.5325), 0.0785889 s/step\n",
      "Meep progress: 2.64/15.995224211855438 = 16.5% done in 76.7s, 388.1s to go\n",
      "on time step 1063 (time=2.6575), 0.0813446 s/step\n",
      "Meep progress: 2.7575/15.995224211855438 = 17.2% done in 80.8s, 387.8s to go\n",
      "on time step 1108 (time=2.77), 0.089251 s/step\n",
      "Meep progress: 2.8875/15.995224211855438 = 18.1% done in 84.8s, 384.9s to go\n",
      "on time step 1162 (time=2.905), 0.0752005 s/step\n",
      "Meep progress: 3.02/15.995224211855438 = 18.9% done in 88.8s, 381.7s to go\n",
      "on time step 1215 (time=3.0375), 0.0767113 s/step\n",
      "Meep progress: 3.1550000000000002/15.995224211855438 = 19.7% done in 92.9s, 378.0s to go\n",
      "on time step 1269 (time=3.1725), 0.074589 s/step\n",
      "Meep progress: 3.295/15.995224211855438 = 20.6% done in 96.9s, 373.5s to go\n",
      "on time step 1325 (time=3.3125), 0.0723764 s/step\n",
      "Meep progress: 3.4375/15.995224211855438 = 21.5% done in 100.9s, 368.6s to go\n",
      "on time step 1383 (time=3.4575), 0.0700939 s/step\n",
      "Meep progress: 3.58/15.995224211855438 = 22.4% done in 104.9s, 363.9s to go\n",
      "on time step 1440 (time=3.6), 0.0714301 s/step\n",
      "Meep progress: 3.7225/15.995224211855438 = 23.3% done in 109.0s, 359.3s to go\n",
      "on time step 1498 (time=3.745), 0.0697877 s/step\n",
      "Meep progress: 3.8625000000000003/15.995224211855438 = 24.1% done in 113.0s, 355.1s to go\n",
      "on time step 1553 (time=3.8825), 0.0730724 s/step\n",
      "Meep progress: 3.9975/15.995224211855438 = 25.0% done in 117.1s, 351.4s to go\n",
      "on time step 1607 (time=4.0175), 0.0742494 s/step\n",
      "Meep progress: 4.11/15.995224211855438 = 25.7% done in 121.1s, 350.1s to go\n",
      "on time step 1651 (time=4.1275), 0.0919766 s/step\n",
      "Meep progress: 4.23/15.995224211855438 = 26.4% done in 125.1s, 348.0s to go\n",
      "on time step 1700 (time=4.25), 0.082572 s/step\n",
      "Meep progress: 4.355/15.995224211855438 = 27.2% done in 129.1s, 345.2s to go\n",
      "on time step 1751 (time=4.3775), 0.0790105 s/step\n",
      "Meep progress: 4.4750000000000005/15.995224211855438 = 28.0% done in 133.2s, 342.8s to go\n",
      "on time step 1797 (time=4.4925), 0.087554 s/step\n",
      "Meep progress: 4.61/15.995224211855438 = 28.8% done in 137.2s, 338.8s to go\n",
      "on time step 1853 (time=4.6325), 0.0723907 s/step\n",
      "Meep progress: 4.7475000000000005/15.995224211855438 = 29.7% done in 141.2s, 334.6s to go\n",
      "on time step 1908 (time=4.77), 0.0727974 s/step\n",
      "Meep progress: 4.8925/15.995224211855438 = 30.6% done in 145.3s, 329.7s to go\n",
      "on time step 1966 (time=4.915), 0.0697177 s/step\n",
      "Meep progress: 5.0375000000000005/15.995224211855438 = 31.5% done in 149.3s, 324.8s to go\n",
      "on time step 2024 (time=5.06), 0.0696059 s/step\n",
      "Meep progress: 5.1850000000000005/15.995224211855438 = 32.4% done in 153.4s, 319.8s to go\n",
      "on time step 2083 (time=5.2075), 0.068659 s/step\n",
      "Meep progress: 5.33/15.995224211855438 = 33.3% done in 157.4s, 315.0s to go\n",
      "on time step 2141 (time=5.3525), 0.0696223 s/step\n",
      "Meep progress: 5.4775/15.995224211855438 = 34.2% done in 161.5s, 310.1s to go\n",
      "on time step 2199 (time=5.4975), 0.0695822 s/step\n",
      "Meep progress: 5.62/15.995224211855438 = 35.1% done in 165.5s, 305.5s to go\n",
      "on time step 2257 (time=5.6425), 0.0695874 s/step\n",
      "Meep progress: 5.765/15.995224211855438 = 36.0% done in 169.5s, 300.8s to go\n",
      "on time step 2315 (time=5.7875), 0.069226 s/step\n",
      "Meep progress: 5.91/15.995224211855438 = 36.9% done in 173.6s, 296.2s to go\n",
      "on time step 2371 (time=5.9275), 0.0717046 s/step\n",
      "Meep progress: 6.055/15.995224211855438 = 37.9% done in 177.6s, 291.6s to go\n",
      "on time step 2430 (time=6.075), 0.0684785 s/step\n",
      "Meep progress: 6.202500000000001/15.995224211855438 = 38.8% done in 181.7s, 286.9s to go\n",
      "on time step 2489 (time=6.2225), 0.0688777 s/step\n",
      "Meep progress: 6.3500000000000005/15.995224211855438 = 39.7% done in 185.7s, 282.1s to go\n",
      "on time step 2548 (time=6.37), 0.0684183 s/step\n",
      "Meep progress: 6.495/15.995224211855438 = 40.6% done in 189.8s, 277.6s to go\n",
      "on time step 2606 (time=6.515), 0.0694275 s/step\n",
      "Meep progress: 6.6425/15.995224211855438 = 41.5% done in 193.8s, 272.9s to go\n",
      "on time step 2663 (time=6.6575), 0.0708676 s/step\n",
      "Meep progress: 6.7825/15.995224211855438 = 42.4% done in 197.8s, 268.7s to go\n",
      "on time step 2721 (time=6.8025), 0.0696011 s/step\n",
      "Meep progress: 6.9275/15.995224211855438 = 43.3% done in 201.9s, 264.3s to go\n",
      "on time step 2779 (time=6.9475), 0.0699265 s/step\n",
      "Meep progress: 7.07/15.995224211855438 = 44.2% done in 205.9s, 259.9s to go\n",
      "on time step 2836 (time=7.09), 0.0706054 s/step\n",
      "Meep progress: 7.215/15.995224211855438 = 45.1% done in 209.9s, 255.5s to go\n",
      "on time step 2895 (time=7.2375), 0.0688381 s/step\n",
      "Meep progress: 7.3625/15.995224211855438 = 46.0% done in 214.0s, 250.9s to go\n",
      "on time step 2954 (time=7.385), 0.0685991 s/step\n",
      "Meep progress: 7.5075/15.995224211855438 = 46.9% done in 218.0s, 246.5s to go\n",
      "on time step 3012 (time=7.53), 0.0693077 s/step\n",
      "Meep progress: 7.6525/15.995224211855438 = 47.8% done in 222.0s, 242.0s to go\n",
      "on time step 3070 (time=7.675), 0.0690058 s/step\n",
      "Meep progress: 7.8/15.995224211855438 = 48.8% done in 226.1s, 237.5s to go\n",
      "on time step 3128 (time=7.82), 0.0697583 s/step\n",
      "Meep progress: 7.9475/15.995224211855438 = 49.7% done in 230.1s, 233.0s to go\n",
      "on time step 3187 (time=7.9675), 0.0687527 s/step\n",
      "Meep progress: 8.09/15.995224211855438 = 50.6% done in 234.1s, 228.8s to go\n",
      "on time step 3245 (time=8.1125), 0.0700344 s/step\n",
      "Meep progress: 8.235/15.995224211855438 = 51.5% done in 238.1s, 224.4s to go\n",
      "on time step 3303 (time=8.2575), 0.0693148 s/step\n",
      "Meep progress: 8.38/15.995224211855438 = 52.4% done in 242.1s, 220.0s to go\n",
      "on time step 3362 (time=8.405), 0.0689097 s/step\n",
      "Meep progress: 8.5275/15.995224211855438 = 53.3% done in 246.2s, 215.6s to go\n",
      "on time step 3420 (time=8.55), 0.0694345 s/step\n",
      "Meep progress: 8.6725/15.995224211855438 = 54.2% done in 250.2s, 211.3s to go\n",
      "on time step 3480 (time=8.7), 0.0683235 s/step\n",
      "Meep progress: 8.82/15.995224211855438 = 55.1% done in 254.2s, 206.8s to go\n",
      "on time step 3539 (time=8.8475), 0.0678418 s/step\n",
      "Meep progress: 8.9675/15.995224211855438 = 56.1% done in 258.3s, 202.4s to go\n",
      "on time step 3598 (time=8.995), 0.0686682 s/step\n",
      "Meep progress: 9.115/15.995224211855438 = 57.0% done in 262.4s, 198.0s to go\n",
      "on time step 3656 (time=9.14), 0.069414 s/step\n",
      "Meep progress: 9.26/15.995224211855438 = 57.9% done in 266.4s, 193.7s to go\n",
      "on time step 3715 (time=9.2875), 0.0691863 s/step\n",
      "Meep progress: 9.405/15.995224211855438 = 58.8% done in 270.4s, 189.5s to go\n",
      "on time step 3774 (time=9.435), 0.0688953 s/step\n",
      "Meep progress: 9.55/15.995224211855438 = 59.7% done in 274.4s, 185.2s to go\n",
      "on time step 3833 (time=9.5825), 0.0686471 s/step\n",
      "Meep progress: 9.695/15.995224211855438 = 60.6% done in 278.4s, 180.9s to go\n",
      "on time step 3891 (time=9.7275), 0.0694123 s/step\n",
      "Meep progress: 9.84/15.995224211855438 = 61.5% done in 282.4s, 176.7s to go\n",
      "on time step 3947 (time=9.8675), 0.0727485 s/step\n",
      "Meep progress: 9.955/15.995224211855438 = 62.2% done in 286.4s, 173.8s to go\n",
      "on time step 3993 (time=9.9825), 0.0873192 s/step\n",
      "Meep progress: 10.0525/15.995224211855438 = 62.8% done in 290.6s, 171.8s to go\n",
      "on time step 4029 (time=10.0725), 0.11258 s/step\n",
      "Meep progress: 10.175/15.995224211855438 = 63.6% done in 294.6s, 168.5s to go\n",
      "on time step 4082 (time=10.205), 0.076998 s/step\n",
      "Meep progress: 10.307500000000001/15.995224211855438 = 64.4% done in 298.6s, 164.8s to go\n",
      "on time step 4134 (time=10.335), 0.0777748 s/step\n",
      "Meep progress: 10.4375/15.995224211855438 = 65.3% done in 302.7s, 161.2s to go\n",
      "on time step 4188 (time=10.47), 0.0750952 s/step\n",
      "Meep progress: 10.5725/15.995224211855438 = 66.1% done in 306.7s, 157.3s to go\n",
      "on time step 4242 (time=10.605), 0.0744089 s/step\n",
      "Meep progress: 10.6875/15.995224211855438 = 66.8% done in 310.8s, 154.4s to go\n",
      "on time step 4281 (time=10.7025), 0.105605 s/step\n",
      "Meep progress: 10.77/15.995224211855438 = 67.3% done in 314.9s, 152.8s to go\n",
      "on time step 4317 (time=10.7925), 0.111276 s/step\n",
      "Meep progress: 10.8925/15.995224211855438 = 68.1% done in 319.0s, 149.4s to go\n",
      "on time step 4369 (time=10.9225), 0.0787911 s/step\n",
      "Meep progress: 11.025/15.995224211855438 = 68.9% done in 323.0s, 145.6s to go\n",
      "on time step 4423 (time=11.0575), 0.0745123 s/step\n",
      "Meep progress: 11.165000000000001/15.995224211855438 = 69.8% done in 327.0s, 141.5s to go\n",
      "on time step 4479 (time=11.1975), 0.0718143 s/step\n",
      "Meep progress: 11.307500000000001/15.995224211855438 = 70.7% done in 331.1s, 137.3s to go\n",
      "on time step 4536 (time=11.34), 0.0711064 s/step\n",
      "Meep progress: 11.4525/15.995224211855438 = 71.6% done in 335.1s, 132.9s to go\n",
      "on time step 4594 (time=11.485), 0.0697675 s/step\n",
      "Meep progress: 11.595/15.995224211855438 = 72.5% done in 339.1s, 128.7s to go\n",
      "on time step 4652 (time=11.63), 0.0700051 s/step\n",
      "Meep progress: 11.735/15.995224211855438 = 73.4% done in 343.1s, 124.6s to go\n",
      "on time step 4708 (time=11.77), 0.0716763 s/step\n",
      "Meep progress: 11.875/15.995224211855438 = 74.2% done in 347.2s, 120.5s to go\n",
      "on time step 4764 (time=11.91), 0.0723048 s/step\n",
      "Meep progress: 12.0225/15.995224211855438 = 75.2% done in 351.2s, 116.1s to go\n",
      "on time step 4823 (time=12.0575), 0.0686204 s/step\n",
      "Meep progress: 12.17/15.995224211855438 = 76.1% done in 355.3s, 111.7s to go\n",
      "on time step 4881 (time=12.2025), 0.0689912 s/step\n",
      "Meep progress: 12.3175/15.995224211855438 = 77.0% done in 359.4s, 107.3s to go\n",
      "on time step 4939 (time=12.3475), 0.0689882 s/step\n",
      "Meep progress: 12.4575/15.995224211855438 = 77.9% done in 363.4s, 103.2s to go\n",
      "on time step 4995 (time=12.4875), 0.0719165 s/step\n",
      "Meep progress: 12.602500000000001/15.995224211855438 = 78.8% done in 367.4s, 98.9s to go\n",
      "on time step 5053 (time=12.6325), 0.069433 s/step\n",
      "Meep progress: 12.75/15.995224211855438 = 79.7% done in 371.5s, 94.6s to go\n",
      "on time step 5111 (time=12.7775), 0.0697191 s/step\n",
      "Meep progress: 12.895/15.995224211855438 = 80.6% done in 375.6s, 90.3s to go\n",
      "on time step 5168 (time=12.92), 0.0710693 s/step\n",
      "Meep progress: 13.035/15.995224211855438 = 81.5% done in 379.6s, 86.2s to go\n",
      "on time step 5224 (time=13.06), 0.0718595 s/step\n",
      "Meep progress: 13.18/15.995224211855438 = 82.4% done in 383.7s, 82.0s to go\n",
      "on time step 5282 (time=13.205), 0.0700622 s/step\n",
      "Meep progress: 13.3125/15.995224211855438 = 83.2% done in 387.8s, 78.1s to go\n",
      "on time step 5334 (time=13.335), 0.0769529 s/step\n",
      "Meep progress: 13.4525/15.995224211855438 = 84.1% done in 391.8s, 74.1s to go\n",
      "on time step 5391 (time=13.4775), 0.0705171 s/step\n",
      "Meep progress: 13.5975/15.995224211855438 = 85.0% done in 395.8s, 69.8s to go\n",
      "on time step 5449 (time=13.6225), 0.0698352 s/step\n",
      "Meep progress: 13.7425/15.995224211855438 = 85.9% done in 399.8s, 65.5s to go\n",
      "on time step 5507 (time=13.7675), 0.069379 s/step\n",
      "Meep progress: 13.887500000000001/15.995224211855438 = 86.8% done in 403.9s, 61.3s to go\n",
      "on time step 5565 (time=13.9125), 0.0696919 s/step\n",
      "Meep progress: 14.0325/15.995224211855438 = 87.7% done in 407.9s, 57.1s to go\n",
      "on time step 5623 (time=14.0575), 0.0692 s/step\n",
      "Meep progress: 14.18/15.995224211855438 = 88.7% done in 412.0s, 52.7s to go\n",
      "on time step 5682 (time=14.205), 0.068815 s/step\n",
      "Meep progress: 14.3275/15.995224211855438 = 89.6% done in 416.0s, 48.4s to go\n",
      "on time step 5741 (time=14.3525), 0.0689857 s/step\n",
      "Meep progress: 14.47/15.995224211855438 = 90.5% done in 420.1s, 44.3s to go\n",
      "on time step 5798 (time=14.495), 0.0706006 s/step\n",
      "Meep progress: 14.61/15.995224211855438 = 91.3% done in 424.1s, 40.2s to go\n",
      "on time step 5854 (time=14.635), 0.0718805 s/step\n",
      "Meep progress: 14.7575/15.995224211855438 = 92.3% done in 428.1s, 35.9s to go\n",
      "on time step 5912 (time=14.78), 0.0691972 s/step\n",
      "Meep progress: 14.89/15.995224211855438 = 93.1% done in 432.1s, 32.1s to go\n",
      "on time step 5966 (time=14.915), 0.0749979 s/step\n",
      "Meep progress: 15.035/15.995224211855438 = 94.0% done in 436.2s, 27.9s to go\n",
      "on time step 6024 (time=15.06), 0.0694761 s/step\n",
      "Meep progress: 15.1775/15.995224211855438 = 94.9% done in 440.2s, 23.7s to go\n",
      "on time step 6081 (time=15.2025), 0.0707235 s/step\n",
      "Meep progress: 15.3225/15.995224211855438 = 95.8% done in 444.2s, 19.5s to go\n",
      "on time step 6140 (time=15.35), 0.0688503 s/step\n",
      "Meep progress: 15.46/15.995224211855438 = 96.7% done in 448.2s, 15.5s to go\n",
      "on time step 6195 (time=15.4875), 0.0737883 s/step\n",
      "Meep progress: 15.605/15.995224211855438 = 97.6% done in 452.2s, 11.3s to go\n",
      "on time step 6254 (time=15.635), 0.0681949 s/step\n",
      "Meep progress: 15.7525/15.995224211855438 = 98.5% done in 456.3s, 7.0s to go\n",
      "on time step 6313 (time=15.7825), 0.0685689 s/step\n",
      "Meep progress: 15.9/15.995224211855438 = 99.4% done in 460.3s, 2.8s to go\n",
      "on time step 6372 (time=15.93), 0.0687192 s/step\n",
      "run 0 finished at t = 15.9975 (6399 timesteps)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "#Step function for moving the source with time given the electron velocity:\n",
    "def move_source(sim):\n",
    "    sim.change_sources([mp.Source(mp.ContinuousSource(frequency=1e-10),\n",
    "                                  component=mp.Ex,\n",
    "                                  center=mp.Vector3(-0.5*sx + v_e*sim.meep_time(), \n",
    "                                                    grating_height + e_height,\n",
    "                                                    0),\n",
    "                                  size=mp.Vector3(0,0,0))])\n",
    "\n",
    "# -- Run Simulation -- \n",
    "sim.run(move_source,\n",
    "        mp.at_beginning(mp.output_epsilon),\n",
    "        mp.in_volume(monitor_xy,\n",
    "                     mp.to_appended(\"Ex_xy\",\n",
    "                                    mp.at_every(t_sample/tcon, mp.output_efield_x))),\n",
    "        mp.in_volume(monitor_xy,\n",
    "                     mp.to_appended(\"Ey_xy\",\n",
    "                                    mp.at_every(t_sample/tcon, mp.output_efield_y))),\n",
    "        until = sx/v_e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ef86404-4464-4fb0-99e3-dda14ca6d3c6",
   "metadata": {},
   "source": [
    "## Analyze the Outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cd2abc6-ec93-4036-8af3-da7b78ef61b2",
   "metadata": {},
   "source": [
    "### Plot Power Out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "b33f2cdc-0398-4e4f-94e4-148c9e57496b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4025000000000003"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "-566.2872510454871"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAHJCAYAAADXbg7bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABdwUlEQVR4nO3dd5xcdb3/8ddn2vaWZNMbCaGEDgESioCAoiLYRRGx/OSq1+5VsXvvVS/qtV/L5dpAuQoiCFe6FOmBhBJKQnov23uZ9v39cWY2u5vd7MzulJ3M+/l4rDtzZuacz47L7Dvfas45RERERGRy8+W7ABEREREZm0KbiIiISAFQaBMREREpAAptIiIiIgVAoU1ERESkACi0iYiIiBSAQL4LGI9p06a5hQsX5rsMERERkTGtXr26yTlXP9HzFGRoW7hwIatWrcp3GSIiIiJjMrNtmTiPukdFRERECoBCm4iIiEgBUGgTERERKQAKbSIiIiIFQKFNREREpAAotImIiIgUAIU2ERERkQKg0CYiIiJSABTaRERERAqAQpuIiIhIAVBoExERESkACm0iIiIiBUChTURERKQAKLSJiIiIFACFNhEREZECoNAmImPqCUfzXYKISNFTaBORUcXiji/f+gLHfv0efvbgxnyXIyJS1BTaRGRUtz+/ixtWbueIGVV8755XWLOzLd8liYgULYU2ERmRc45fPbKFw6dXcuM/raAs6Od/V27Pd1kiIkVLoU1ERrRubycv7e7gyjMWUlMW5JITZnP787sJR+P5Lk1EpCgptInIiB5e3wjABUdPB+C8o+rpCcd4cXd7PssSESlaCm0iMqKHNzRy5IwqZtWUAXDKgikArNraks+yRESKlkKbiBwgHI3z9NZWzl4ybeBYfVUJC6eWs2prax4rExEpXgptInKAdXs7CEfjnLygbsjxkxfU8cz2tvwUJSJS5BTaROQAa3Z649aOm1Mz5PhRM6to6uqnrSecj7JERIqaQpuIHOCFne3UlQeZW1c25PiS6VUAbGzoykdZIiJFTaFNRA6wZlc7x82txcyGHD98eiUAGxTaRERyTqFNRIaIxuJsauji6FlVBzw2p7aM0qBPLW0iInmg0CYiQ2xv6SEciw90hQ7m8xmL6ysV2kRE8kChTUSGSHZ9Lkl0hQ63qL6SLU3duSxJRERQaBORYZKtaItHCW1z68rY095LPO5yWZaISNFTaBORITbs62RObRmVJYERH59TW0Yk5mjo7M9xZSIixU2hTUSG2NzUzaL6ilEfn5NYBmRna0+uShIRERTaRGSYbc09LJhaPurj8xKhbVdbb65KEhERFNpEZJC2njDtvREWTBm9pW12bbKlTaFNRCSXFNpEZMC2Zq/Lc/5BWtrKQwGmVIQU2kREckyhTUQGbGvxQtvBukfBm4yg7lERkdxSaBORAdubvfXX5k85eGibUV3Kvva+XJQkIiIJCm0iMmBHSy/1VSWUh0Ze7iNpRnUJDZ0KbSIiuZSz0GZmR5rZc4O+OszsU2Y2xczuM7MNie91uapJRIba3d47MNHgYKZXldLaEyEcjeegKhERgRyGNufcK865E51zJwKnAD3ArcDVwP3OuSXA/Yn7IpIHe9r7mFVdOubzpleXANDYpQV2RURyJV/do+cDm5xz24BLgesSx68D3pSnmkSKmnOOPW29zKpNIbRVeaGtoUNdpCIiuZKv0HYZ8MfE7RnOuT0Aie/T81STSFHr6IvSHY4xuya17lFAW1mJiORQzkObmYWAS4A/p/m6q8xslZmtamxszE5xIkVsT7u3hMfMmtS7RxXaRERyJx8tba8DnnHO7Uvc32dmswAS3xtGepFz7lrn3DLn3LL6+voclSpSPPYklvCYnUL36NSKEGbQqO5REZGcyUdoexf7u0YBbgeuTNy+Ergt5xWJCHvavAA2K4Xu0YDfx9SKErW0iYjkUE5Dm5mVAxcCtww6fA1woZltSDx2TS5rEhHPnvZefLZ/ksFY6qtKaNLsURGRnDn4CpoZ5pzrAaYOO9aMN5tURPJod1sf06tKCfhT+7fc1IoQzd3hLFclIiJJ2hFBRADY25Hach9JUypCtCi0iYjkjEKbiADemLZUlvtImlIRoqVLoU1EJFcU2kQE5xy723uZlcJyH0lTK0J09kfpj8ayWJmIiCQptIkI7b0R+iLxlNZoS5pSGQKgtTuSrbJERGQQhTYRYXdbco221LtHp1Z4oa25WzNIRURyQaFNRNjbkfpuCElTKrylQTQZQUQkNxTaRISGDq+1bEZ1erNHQaFNRCRXFNpEZGCR3GSXZyoGukc1g1REJCcU2kSExs5+qksDlAb9Kb+mpiyI32dqaRMRyRGFNhGhqSvMtBS3r0ry+Yy68qB2RRARyRGFNhGhsbOf+sr0QhtAbXmIth6FNhGRXFBoExGauvrTbmkDqC0L0tajddpERHJBoU1EJtDSFqStV6FNRCQXFNpEilxfJEZnf5T6cbS01ZSFaFf3qIhITii0iRS5xk5vuQ+1tImITG4KbSJFLrlG27Sq1NdoS6otC9ITjmnTeBGRHFBoEyly+1vaUt8NIam2PAh4G86LiEh2KbSJFLmmxI4G42ppK/de064ZpCIiWafQJlLkki1tUyvGN6YN0Lg2EZEcUGgTKXJNXf3UlgcJBdL/OKgt81ratFabiEj2KbSJFLnGzn6mjWPmKAxqadOyHyIiWafQJlLkmrrGt7AuQI0mIoiI5IxCm0iRaxznFlYAVSUB/D5T96iISA4otIkUuaZxbmEFYGZUlQZo61X3qIhItim0iRSxvkiM7nCMqZXpL/eRVF0apLMvmsGqRERkJAptIkWsudtrIZtaMf7QVlUaUGgTEckBhTaRItaaCG11Ew5tGtMmIpJtCm0iRSwzLW3qHhURyQWFNpEilomWturSIB1a8kNEJOsU2kSKWEsitE0p15g2EZHJTqFNpIi1dIfxGdSUBcd9jurSAF3hKPG4y2BlIiIynEKbSBFr6QlTVx7C57Nxn6OqNIhz0BVWa5uISDYptIkUsZauMFMmMJ4NoLosAKAuUhGRLFNoEyliLT3hCU1CAK+lDdBkBBGRLFNoEylird3hCS33Ad5EBFBLm4hItim0iRSxlu7MtbRpgV0RkexSaBMpUvG4o7UnPKHlPsCbPQpqaRMRyTaFNpEi1d4bIe6Y8ESEgTFtamkTEckqhTaRItXSk1hYV2PaREQKgkKbSJEa2A1hgqGtNOgn5PeppU1EJMsU2kSKVKZCG3hrtamlTUQku3Ia2sys1sxuNrN1ZrbWzFaY2RQzu8/MNiS+1+WyJpFi1ZrB0FZVGlRoExHJsly3tP0YuNs5dxRwArAWuBq43zm3BLg/cV9Esqw5EdrqJjh7FLxxbVpcV0Qku3IW2sysGngV8GsA51zYOdcGXApcl3jadcCbclWTSDFr7Q5TFvRTFvJP+FxVpQGt0yYikmW5bGlbBDQCvzWzZ83sV2ZWAcxwzu0BSHyfPtKLzewqM1tlZqsaGxtzV7XIIaqle+L7jiZVq3tURCTrchnaAsDJwC+ccycB3aTRFeqcu9Y5t8w5t6y+vj5bNYoUjdaeMHUVwYycy2tpU2gTEcmmXIa2ncBO59zKxP2b8ULcPjObBZD43pDDmkSKVmtPJCPj2SA5EUHdoyIi2ZSz0Oac2wvsMLMjE4fOB14GbgeuTBy7ErgtVzWJFLP23gg1ZZlraesOx4jG4hk5n4iIHCiQ4+t9HLjBzELAZuD9eMHxJjP7ILAdeHuOaxIpSq094Yy1tFUntrLq6o9Sm6FziojIUDkNbc6554BlIzx0fi7rECl28bijvTdCXXnmWtrA28pKoU1EJDu0I4JIEeroi+Ac1GRwTFvyvCIikh0KbSJFqK3HC1eZammr1qbxIiJZp9AmUoRae7zdEGozFdoSExq0K4KISPYotIkUobZEuMrU+LMqtbSJiGSdQptIEWpLtrRlbMkP7zxaq01EJHsU2kSK0P4xbZlpaassUUubiEi2KbSJFKHWnghm+8eiTVQo4CPk99EdjmXkfCIiciCFNpEi1N4Tpro0iN9nGTtnRYmfnrBa2kREskWhTaQItfZEMjZzNKk8FKCrX6FNRCRbFNpEilBbbyTjOxdUlgTo6Vf3qIhItii0iRShtp5wxmaOJpWX+OlW96iISNYotIkUIW+z+MyGtsqSAN3qHhURyRqFNpEi1NaT+e7R8pCfbnWPiohkjUKbSJGJxuJ09kUzPhGhoiSg7lERkSxSaBMpMu29mV1YN6kipO5REZFsUmgTKTKtPcl9R7PR0qbuURGRbFFoEyky7b2JfUcz3tLmJxyNE4nFM3peERHxKLSJFJnW7kRLW4aX/KhI7D+qtdpERLJDoU2kyLRla0xbiR+ALk1GEBHJCoU2kSLT1uN1j9ZkYUwbQI8mI4iIZIVCm0iRaeuJ4PcZ1aWBjJ63IuSdT/uPiohkh0KbSJFp7QlTUxbEzDJ63oGWNs0gFRHJCoU2kSLjbRaf2a5R8HZEALW0iYhki0KbSJHJxmbx4O09CtCjiQgiIlmh0CZSZNp6IhmfOQpQnpw9qiU/RESyQqFNpMi09UQyPnMUBrW0qXtURCQrFNpEikxbTzgrLW1lQT9maP9REZEsUWgTKSLhaJzucCwrY9rMzNs0XrNHRUSyQqFNpIi0Jfcdrch8Sxt4M0g1EUFEJDsU2kSKSFtPdvYdTaooCWgigohIlii0iRSRZGjLxpg28PYf1UQEEZHsUGgTKSKtiX1Hs7G4LkB5KKDFdUVEskShTaSItCe7R7MU2ipLAtrGSkQkSxTaRIrIwESELHWPlof8WvJDRCRLFNpEikhrT4SAz6hI7BOaaZUlAbo1e1REJCsU2kSKSFtPhNryEGaWlfOXhwJ0a/aoiEhWKLSJFJH23nDWxrMBVJb46Q5Hcc5l7RoiIsVKoU2kiLR2R7K2RhtAeUkA56A3otY2EZFMU2gTKSJtvZGsTUIAb3FdQF2kIiJZoNAmUkTaerLbPZqc4KAZpCIimRfI5cXMbCvQCcSAqHNumZlNAW4EFgJbgXc451pzWZdIsWjriVCXzdCWbGnTDFIRkYzLR0vbec65E51zyxL3rwbud84tAe5P3BeRDOuLxOiNxLLbPRpS96iISLZMhu7RS4HrErevA96Uv1JEDl3tvd5uCDVZnIhQUZLoHlVLm4hIxuU6tDngXjNbbWZXJY7NcM7tAUh8n57jmkSKQrY3i4fBExEU2kREMi2nY9qAM51zu81sOnCfma1L9YWJkHcVwPz587NVn8ghqy3Lm8XD/tDWo+5REZGMy2lLm3Nud+J7A3ArcBqwz8xmASS+N4zy2mudc8ucc8vq6+tzVbLIIaO1Jwfdo4nZo11qaRMRybichTYzqzCzquRt4DXAi8DtwJWJp10J3JarmkSKSXtis/i6iux1j5YnJiL0aEybiEjG5bJ7dAZwa2LPwwDwv865u83saeAmM/sgsB14ew5rEikayZa2bO6IEAr4CPl9dKl7VEQk43IW2pxzm4ETRjjeDJyfqzpEilVbT4Sg3yhPdGFmS0WJXy1tIiJZMBmW/BCRHPA2iw+RaO3OmvJQQGPaRESyQKFNpEhke7P4pMqSgGaPiohkgUKbSJFo6w1ndY22pPISvxbXFRHJAoU2kSLR1hOhJotrtCVVlgS0uK6ISBYotIkUibae3HSPlof82ntURCQL0g5tifXWsjv9TEQyrq03nNU12pIqSgLqHhURyYIxQ5uZ+czs3WZ2h5k1AOuAPWb2kpl9z8yWZL9MEZmIvkiMvkg8q7shJFWE1D0qIpINqbS0PQgsBr4IzHTOzXPOTQfOBp4ErjGz92SxRhGZoORm8dncdzTJa2lT96iISKalsrjuBc65yPCDzrkW4C/AX8ws+38JRGTc2pJbWOVg9mhFyE84GicSixP0a9isiEimjPmJOlJgG89zRCR/Wruzv4VVUkVJYv9RTUYQEcmolLexMrNlwJeBBYnXGeCcc8dnqTYRyZDkZvG1uWhpK/HmKXWFozlZYkREpFiks/foDcDngBeAeHbKEZFsyPWYNoAeTUYQEcmodEJbo3Pu9qxVIiJZ05rL0BbyPla0/6iISGalE9q+bma/Au4H+pMHnXO3ZLwqEcmott4woYCPsmD2l1gcaGnTDFIRkYxKJ7S9HzgKCLK/e9QBCm0ik1x7YjcEM8v6tcpDiTFtamkTEcmodELbCc6547JWiYhkTWtPOCddo+DtPQrQo10RREQyKp1FlJ40s6VZq0REsqatJ5KTmaMA5cnZo1ryQ0Qko9JpaTsLeJ+ZbcYb06YlP0QKRFtPhAVTy3NyrUrNHhURyYp0QttFWatCRLKqrTfMieW1OblWWdCPGdp/VEQkw9IJbVeOcvzfMlGIiGSP1z2amzFtZuZtGq/ZoyIiGZVOaOsedLsUuBhYm9lyRCTTesMx+qPxnO5OUB7yq6VNRCTDUg5tzrnvD75vZv8JaLFdkUkul5vFJ1WWqKVNRCTT0pk9Olw5sChThYhIdgxsYZWDzeKTykvU0iYikmnpbBj/At5iugB+oB6NZxOZ9Fp7vJa2XHaPVoQCCm0iIhk2ZmgzsxXAk3hj2JKiwD7nnD6VRSa59oGWttx1j1aUBGjo7MvZ9UREikEq3aNXAquB7wDnAf3OuV0KbCKFoSXR0jalIrehrUeL64qIZNSYLW3OuQ8DmNlRwOuA35lZDfAgcDfwmHNOn84ik1Rrd2IiQkUuu0f92ntURCTDUp6I4Jxb55z7oXPuIuDVwKPA24GV2SpORCaupTtCZUmAkoA/Z9esKAnQo9mjIiIZlcqYtk8BjwHPJrtEnXO9wJ2JLxGZxFq6+3PaygZeS1t3OIpzDjPL6bVFRA5VqcwenQv8GDjKzNYAj+OFuCeccy3ZLE5EJq6lJ8KUHK7RBl5Lm3PQG4lRHkpnDW8RERlNKmPa/gXAzELAMuAM4APA/5hZm3NuaXZLFJGJaO0OM60yt6GtPLFpfFd/VKFNRCRD0llctwyoBmoSX7vReDaRSa+lO0xdDmeOAlSWeOPnNINURCRzUhnTdi1wDNCJF9IeB37gnGvNcm0ikgEt3eGcd48mW9c0g1REJHNSaWmbD5QAe4FdwE6gLYs1iUiG9IZj9EZieWhp80KbZpCKiGROKmPaLjJv+tcxeOPZPgsca2YteJMRvp7lGkVknJJbWE3NcWgrD3ndo9rKSkQkc1IaIeycc8CLZtYGtCe+LgZOAxTaRCaploGFdfPT0tYdVmgTEcmUVMa0fQKvhe1MIEJiuQ/gN8ALWa1ORCYkGdpyuYUV7J89qpY2EZHMSaWlbSFwM/Bp59ye7JYjIpmU7B6ty/FEhMpQMrRpTJuISKakMqbtM7koREQyL9nSlvMxbSUa0yYikmnprNMmIgWmpTuMz6C6LLfbWAX9PkIBH92aPSoikjFphTYzqx/pdprn8JvZs2b2t8T9KWZ2n5ltSHyvG895ReRALd1hastD+H253/+zIuSnRxMRREQyJt2Wtj+bWdDMyoA/j/OanwTWDrp/NXC/c24JcH/ivohkQGtPmLry3LayJVWUBLS4rohIBqUb2r4M3AHcDXwh3YuZ2VzgDcCvBh2+FLgucfs64E3pnldERtbSHWZqRUlerl0RCmgbKxGRDEo3tF0M3AX8AnjXOK73I+DzQHzQsRnJWamJ79PHcV4RGYG372i+Wtr8WqdNRCSD0g1ty51zP3TO/Qk4Kp0XmtnFQINzbnWa10y+/iozW2VmqxobG8dzCpGi09IdyfkabUkVJQHNHhURyaCUdkQY5Dtm9j3A4bWapeNM4BIzez1QClSb2R+AfWY2yzm3x8xmAQ0jvdg5dy1wLcCyZctcmtcWKTqxuKOluz9v3aPlIT8NHf15ubaIyKEorZY259zdzrnPOec+75y7O83XftE5N9c5txC4DHjAOfce4HbgysTTrgRuS+e8IjKy1p4wcQfTKvPY0qbuURGRjJkM67RdA1xoZhuACxP3RWSCmrq8Vq5pVfmbiKDuURGRzEm3ezQjnHMPAQ8lbjcD5+ejDpFDWVOntxvCtMo8hbaSgBbXFRHJoDFb2sxszFU5U3mOiOTWQEtbvkJbyE84GicSi4/9ZBERGVMq3aMPmtnHzWz+4INmFjKzV5vZdewfkyYik0QytNXnsaUN0FptIiIZkkr36EXAB4A/mtlhQBve7E8/cC/wQ+fcc9kqUETGp7Grn5DfR3VZXkZBUJHYNL4rHKUmT7syiIgcSsb8NHfO9QE/B35uZkFgGtDrnGvLcm0iMgFNnWGmVobI1+iF/S1tmowgIpIJaf0T3DkXAfZkqRYRyaDm7v68jWcDb/YooP1HRUQyZDIs+SEiWdDU1Z+3NdpgUEubZpCKiGSEQpvIIaqpM5zXlrbyUGJMm1raREQyIuXQZp73mNnXEvfnm9lp2StNRMbLOed1j+ZpYV2AyoGWNoU2EZFMSKel7efACuBdifudwM8yXpGITFh7b4RIzOW3pS05e1RLfoiIZEQ6ExFOd86dbGbPAjjnWs0sfwNmRGRU+xfWzd9/opWaPSoiklHptLRFzMwPOAAzqwe01LnIJNSY2MIqXwvrApQF/Zih/UdFRDIkndD2E+BWYLqZfQt4FPh2VqoSkQnJ92bxAGbmbRqv2aMiIhmRcveoc+4GM1uNt7m7AW9yzq3NWmUiMm753nc0qTzkV0ubiEiGpLu47jpgXZZqEZEMaejsJ+Azasvyu31UZYla2kREMmXM0GZmnznY4865H2SuHBHJhH3tfcyoLsXny88WVknlJWppExHJlFRa2qoS348ETgVuT9x/I/BwNooSkYnZ29HHjOr8do2Ct5WVQpuISGaksmH8vwKY2b3Ayc65zsT9bwB/zmp1IjIuezv6OGpm1dhPzLKKkgANnX35LkNE5JCQzuzR+UB40P0wsDCj1YhIRiS7R/OtoiRAjxbXFRHJiHQmIvweeMrMbk3cfxNwfcYrEpEJ6eyL0B2OMXMyhLaQX3uPiohkSDpLfnzLzO4CzsZbYPf9zrlns1aZiIzLvg6vO3JmzSQIbSUBejR7VEQkI1IObcmN4gd5o5m90Tn3bxmuSUQmYG+7t0bbpOgeDfnpDkdxzmGW35msIiKFLp0xbd2DvmLA69CYNpFJZ2+ypW0yhLaSAM5Bb0StbSIiE5VO9+j3B983s/9k//IfIjJJTKbu0fLEpvFd/VHKQ2mt5S0iIsOk09I2XDmwKFOFiEhm7G3vo6YsSGnQn+9SqCzxatAMUhGRiUtnTNsLeBMQAPxAPfDv2ShKRMZvb0ffpOgaBQZa1zSDVERk4tLpr7h40O0osM85p09ikUlmX0cfMyZB1yh4e48CmkEqIpIB6XSPftQ5ty3xtcs5FzWz72StMhEZl73tfcycBFtYAZSHvO5RbWUlIjJx6YS2C0c49rpMFSIiExeJxWnq6p803aOVJeoeFRHJlDG7R83sI8BHgUVmtmbQQ1XAY9kqTETSt6etj7iDuXXl+S4F2D97tCes0CYiMlGpjGn7X+Au4D+Aqwcd73TOtWSlKhEZlx2tPQDMnVKW50o8+1vaNKZNRGSixgxtzrl2oB14V/bLEZGJ2NHihbZ5k6SlrSIxpq2rTy1tIiITlUr36KPOubPMrJP9S34AGOCcc9VZq05E0rKjtQe/z5g1SWaPBvw+SoM+utU9KiIyYam0tJ2V+F6V/XJEZCJ2tPQyu7aUgH8i62ZnVmVJQBMRREQyYPJ8sovIhO1o7Zk0XaNJFSUBLfkhIpIBqXSPJrtFbYSH1T0qMonsaOnl/KOm57uMISpCCm0iIpmQSveoukVFCkBvOEZTVz/zJsnM0SR1j4qIZEY6e4+W4q3XdhZey9sjwC+dc31Zqk1E0rAzsdzHvCmTrXvUT1NXON9liIgUvHTGtF0PHAP8FPgvYCnw+2wUJSLpG1ijbZKNaassDap7VEQkA9LZMP5I59wJg+4/aGbPZ7ogERmfHS29AJOwe9RPp0KbiMiEpdPS9qyZLU/eMbPT0TZWIpPGtuYeSoM+6isnx2bxSZqIICKSGemEttOBx81sq5ltBZ4AzjGzF4btSToiMys1s6fM7Hkze8nM/jVxfIqZ3WdmGxLf68b1k4gUuY2NXSyur8RspIne+VNREqAnHCMed2M/WURERpVO9+hFE7xWP/Bq51yXmQWBR83sLuAtwP3OuWvM7Gq8/U2/MMFriRSdjfs6Oe2wKfku4wDJ/Ue7w1GqSoN5rkZEpHClHNqcc9smciHnnAO6EneDiS8HXAqcmzh+HfAQCm0iaenqj7K7vY8lMybfCj0VydDWH1NoExGZgHSW/FgGfBlYkHhdcu/R49M4hx9YDRwO/Mw5t9LMZjjn9uCdbI+ZTa6VQUUKwKYG799Di+sr81zJgSpKEpvGa1ybiMiEpNM9egPwOeAFID6eiznnYsCJZlYL3Gpmx6b6WjO7CrgKYP78+eO5vMgha2MitC2ZMflC20D3qEKbiMiEpBPaGp1zt2fios65NjN7CG+c3D4zm5VoZZsFNIzymmuBawGWLVumEc0ig2xo6CLoNxZMsoV1QaFNRCRT0pk9+nUz+5WZvcvM3pL8SvXFZlafaGHDzMqAC4B1wO3AlYmnXQnclkZNIoLX0nbYtAoC/nT+k86N5Jg2rdUmIjIx6bS0vR84Cm8CQbJ71AG3pPj6WcB1iXFtPuAm59zfzOwJ4CYz+yCwHXh7GjWJCLCxoZOls6vzXcaI1NImIpIZ6YS2E5xzx433Qs65NcBJIxxvBs4f73lFil1fJMb2lh4uOWF2vksZUYVCm4hIRqTTl/KkmS3NWiUiMi4v7e4g7uCYOTX5LmVEyZa2rv5YnisRESls6bS0nQVcaWZb8BbKNSA+bD9SEcmxNTvbADhhbm1e6xhNadCHz9TSJiIyUePdEcGAucCXMluOiKTrhZ3tTK8qYWZNab5LGZGZUVES0DptIiITlNaOCGZ2IvBu4B3AFuAvWapLRFL0/M42jp87ObtGkypLtGm8iMhEjRnazOwI4DLgXUAzcCNgzrnzslybiIyhsy/C5qZuLj1xTr5LOajKkgDdYYU2EZGJSKWlbR3wCPBG59xGADP7dFarEpGUvLS7A+fguEne0lZREqCzT6FNRGQiUpk9+lZgL/Cgmf2PmZ2PN6ZNRPIsOQnh+Ek6czRJ3aMiIhM3Zmhzzt3qnHsn3sK6DwGfBmaY2S/M7DVZrk9EDuKpLS0cNq2CqZUl+S7loCpK/HRryQ8RkQlJeZ0251y3c+4G59zFeDNHnwOuzlZhInJwsbjjqS0tnH7YlHyXMibNHhURmbhxbVTonGtxzv23c+7VmS5IRFKzdk8HHX1Rli+amu9SxqSJCCIiEzf5dpcWkZSs3NICwOmLCqOlTWPaREQmRqFNpEA9ubmZBVPLmVVTlu9SxlRZEiASc/RHNa5NRGS8FNpEClC8gMazwf79RzUZQURk/BTaRArQy3s6aO+NsGLx5B/PBl73KECX1moTERk3hTaRAvT4piYAViyaludKUlNZ4gfQDFIRkQlQaBMpQI9vamZRfcWk3SR+uGRLm2aQioiMn0KbSIGJxOI8taWFMxcXRisbDOoeVUubiMi4KbSJFJg1O9voCcc4o0DGs8HgiQgKbSIi46XQJlJgHt/YDFAQi+omVWoigojIhCm0iRSYxzc1s3RWNXUVoXyXkrKqUi+0dSq0iYiMm0KbSAHpi8RYvb21oLpGASpCAcygsy+S71JERAqWQptIAXlmWyvhaJwzDi+s0ObzGZUlATrU0iYiMm4KbSIF5PFNzfh9xqkLC2MnhMGqS4PqHhURmQCFNpEC8vimJo6fW0NVaTDfpaStqjSg7lERkQlQaBMpEP3RGC/u6ijIVjZIhja1tImIjJdCm0iBWLenk3AszonzavNdyrhUlQbp7FdLm4jIeCm0iRSI53e2AXBCwYY2tbSJiEyEQptIgXhuRxv1VSXMLpD9RodTaBMRmRiFNpEC8dyONk6YW4uZ5buUcakqDdLRG8E5l+9SREQKkkKbSAFo742wubGbE+fV5LuUcasqDRCNO/oi8XyXIiJSkBTaRArAi7vaATh+bm1+C5mA5DIlWvZDRGR8FNpECsDaPR0ALJ1dnedKxq86sf+odkUQERkfhTaRAvDK3k6mVYaYVlmS71LGrVotbSIiE6LQJlIAXtnXyVEzC7eVDbwxbYBmkIqIjJNCm8gkF4s7XtnbyZEzq/JdyoTsH9Om0CYiMh4KbSKT3LbmbvqjcY4q+NCWbGlT96iIyHgotIlMcuv2dgKoe1REpMgptIlMcuv2duIzWDKjMt+lTEhFKICZWtpERMZLoU1kklu/t5OFUysoDfrzXcqE+HxGZUlAS36IiIyTQpvIJLepsYtF9YXdypZUXRpU96iIyDjlLLSZ2Twze9DM1prZS2b2ycTxKWZ2n5ltSHyvy1VNIpNdNBZnW3MPi6dX5LuUjPA2jVf3qIjIeOSypS0KfNY5dzSwHPhnM1sKXA3c75xbAtyfuC8iwM7WXsKxOIsPkZY2L7SppU1EZDxyFtqcc3ucc88kbncCa4E5wKXAdYmnXQe8KVc1iUx2m5u6AFhcf6i0tAXp7FdLm4jIeORlTJuZLQROAlYCM5xze8ALdsD0fNQkMhltaugGYNE0tbSJiBS7nIc2M6sE/gJ8yjnXkcbrrjKzVWa2qrGxMXsFikwim5u6mFoRoq4ilO9SMqKqNEBHr1raRETGI6ehzcyCeIHtBufcLYnD+8xsVuLxWUDDSK91zl3rnFvmnFtWX1+fm4JF8mxTQzeLDpGuUUh0j/ZFcc7luxQRkYKTy9mjBvwaWOuc+8Ggh24HrkzcvhK4LVc1iUx2mxq7DplJCOC1tEXjjr5IPN+liIgUnEAOr3UmcAXwgpk9lzj2JeAa4CYz+yCwHXh7DmsSmbTaesI0d4cPsdCW3DQ+QlmosBcLFhHJtZyFNufco4CN8vD5uapDpFBsakxMQjiEukerE/uPdvRFmV7YW6mKiOScdkQQmaQ2NyaX+zh0WtqqB7W0iYhIehTaRCapTY3dBP3G3LqyfJeSMVWJljYt+yEikj6FNpFJalNjFwumVhDwHzr/me4f06bQJiKSrkPnr4HIIWZzYxeHF2rXaMsWeP5PsONpGLS8x/6WNnWPioikK5ezR0UkRZHERvEXHTsz36Wkxzl46Br4x3eARFhb8lp4y7VQVkt1mdfS1qHQJiKSNrW0Zcne9j52tvbkuwwpUNtbeojGXeFNQnj4P+Ef18AJl8GHH4PXfAs2PQA3vx/iMSpCfgI+o61HoU1EJF1qacuCpq5+3vPrlWxs6GLzt1+PzzfaSiciI9vUUIAzR7c/CQ9+E45/J7zpF2AGM4+Fkkr4v0/Ckz/Hzvg4NWVB2rWVlYhI2tTSlmF3v7iXZd/8OxsTf3S3NnfnuSIpRAW3RlssCn/7NFTPhTf8wAtsSae8D5a8Bh76DnQ1UFMepE2hTUQkbQptGbSvo4/P3fz8kGMv7GrPUzVSyDY1djG9qmRgtuWk9+LN0PAyvPZbXsvacK/9NkS64fGfUlMW1KbxIiLjoNCWQZ++8bmBpQymVYbwGazZqdAm6SuoPUfjMfjHd2HGcXD0JSM/Z9oSWPomWP07ZpSENaZNRGQcFNoyaN3eTgAe+fx5rPzSBbz6qBncvHon/dFYniuTQuKcY1NDF4unF0jX6IZ7oWUTnP0Z8B3kI+WMj0N/Bxf2P6AxbSIi46DQliHRWJz23ggff/XhzJtSjt9nvH3ZXNp7I7yg1jZJQ1NXmI6+aOG0tK38JVTNhqPfePDnzTkZZp3Ais67aesJ56Y2EZFDiEJbhuxp7yMWd0O2HDpiRhXgLd8gkqpNhbTnaNt22PyQN9nAn8L4uxPfw+ze9cwNbyIWd2M/X0REBii0ZUBfJMaHrl8FwIKp+7u0ahMLiWr8jqRjILRNL4DQtuYm7/sJ70zt+ce9jbgFuNj3hHZFEBFJk0JbBjywrmFgPNvgJRqqy4KYoeUNJC2bGropDfqYVV2a71IOzjlYcyPMXwF1C1N7TfkUmqaewgW+1RrXJiKSJoW2DFi/zwtspx82hfrKkoHjfp9RXRrU+B1Jy8bGLhZNq5z8izLveR6a1nuL6aahde6FHOHbRc/eDVkqTETk0KTQlgFbm7qZU1vGjf+0ArOhf2hry4PqHpW0rNvTwVGzqvJdxtjW3Aj+EBzzprRe1n/4awEIbbwrC0WJiBy6FNoyYEtT96gr19eUafV3SV1jZz8Nnf0snVWd71IOLh6Hl271djooq0vrpWX1h/FyfAG12+/LUnEiIocmhbYJ2t3Wy/M72zls2sihrSzopy+sddokNWv3dACwdPYkD227n4XOPaMvpnsQNeVB7oufQl3zs9DdnIXiREQOTQptE+Cc44xrHgBg/pTyEZ9TGvTTp8V1JUUvJ0PbZG9pW/c3MD8c8Zq0X1pbFuK+2Mn4iMOGe7JQnIjIoUmhbQLueGHPwO1XHzV9xOeUBn30RRTaJDUv7+5gTm0ZteWhfJdycOvugIVnpd01ChAK+NgaWkJnYCpseiALxYmIHJoU2ibgq399EYDVX7mARaMshFoa9NMXieeyLClgL+/p4OjJ3srWtAGaXoGjLh73KeoqQqwvPwm2POwtHSIiImNSaBun7v4orT0Rrli+gKmDlvkYrjTgV0ubpKQnHGVzYxdLJ/vM0XV3eN+PfN24T1FXHuK5wPHQtc9bNkRERMak0DZOu9t6ATj1sCkHfZ66RyVVz2xrI+7gpAXpdznm1Ct3wszjoXbeuE9RVx7iSXeMd2fLwxkqTETk0KbQNk67EqFtTu3BV633JiKoe1TG9uTmZvw+Y9lkDm3dzbDjqQm1sgFMqQixrn8K1M739i4VEZExKbSN0+62PgBm1ZQd9HklQT/haJy4NseWMazc0syxc2qoKk1h4/V82Xgf4OCI107oNHXlIVq7I3DYq2DroxBXa7SIyFgU2sZhb3sfX7r1BQCmV40+ng287lGAfrW2yUH0hmM8v6Od5WN0t+fd+nugYjrMOmlCp6krD9LVHyUy/2zoa4O9L2SmPhGRQ5hC2zj8/KGNA7cD/oO/haUBP4DGtclBPbu9lXAszvJFU/NdyuhiEdh4v7c2m29iHx11Fd6SJm0zlnsHNK5NRGRMCm3jEE60mlWE/GM+tzSYCG1aYFcO4u9rGwgFfGNObMmr7U9CfzsccdGETzUlEdqafVNg6hKvi1RERA5KoW0c9rR749l+94HTxnxusntUa7XJaOJxx10v7uFVS+qpLAnku5zRrb/b2yB+0bkTPlVtuTdur6U7DPOXw46V3n6mIiIyKoW2cdjV1strj5nBqQvHbhUZaGlT96iM4pGNTexp7+OSE2fnu5SDW3+PtwtCycTXkUu2tLV2R7zQ1tem9dpERMag0JamvkiMjQ1dzKkdea/R4fa3tCm0ych+9chmplaEeO0xM/JdyuiaN0Hzhox0jQJMSWzT1dIThvkrvIPbn8jIuUVEDlUKbWn60i3eLLc5dQdf6iNp/0QEdf3IgR5e38gjG5r4yLmLKQmMPUYyb9YnNnZfkv4G8SNJ7q3a2h2GKYugfJrXRSoiIqNSaEtTMqy98YRZKT2/RBMRZBSRWJxv37mWeVPKuGLFgnyXc3Dr74ZpR8KUwzJyulDAR3VpgOaufjDzukjV0iYiclAKbWn67GuOZOs1b2B61cF3QkgaWKdN3aMyzC8e2sS6vZ185Q1LJ3crW18HbHt8wgvqDjetsoSmrrB3Z/5yaN0KnXszeg0RkUOJQluW7Z+IoO5R2e/l3R385P4NXHLCbF57zMx8l3Nwmx+EeCRj49mSvNDW790ZGNf2ZEavISJyKFFoyzLNHpXhwtE4//Ln56ktD/GvlxyT73LGtv4eKK2Beadn9LTTqkL7Q9vM4yFQqnFtIiIHodCWZaUBzR6VoX724EZe3tPBt9987MDOAJNWPO6FtsMvBH9m15Ab0j0aCMGcZRrXJiJyEAptWbZ/RwR1jwpsbOjiZw9u5C0nzeE1k71bFGD3M9DTlPGuUYCpFSW090YGdhhh/umwZw2EuzN+LRGRQ4FCW5ape1QGu+audZQF/Xz5DUfnu5TUrL8bzAeHn5/xU0+rSqzV1p2cjLACXAx2rsr4tUREDgU5C21m9hszazCzFwcdm2Jm95nZhsT3ulzVkyt+nxH0myYiCCs3N/P3tfv48LmLmVpZku9yUrP+bpi3HMozvyfqtMR7MDCube6pgGlcm4jIKHLZ0vY7YHgfy9XA/c65JcD9ifuHnNKAXy1tRc45x7fvWsesmlI+eFZm1jrLuvZdsPeFjC/1kTSt0mtpa0yGtrJamL5U49pEREaRs9DmnHsYaBl2+FLgusTt64A35aqeXCoJ+unX4rpF7W9r9vD8jjY+c+ERA13mk96Ge73vWRjPBvtb2pqTkxEgsXn80xDXfy8iIsPle0zbDOfcHoDE9+l5ricrSoM+dY8Wsf5ojO/es46jZlbxlpPn5ruc1K2/B2rnQ/2RWTn9Ad2j4IW2cCfseykr1xQRKWT5Dm0pM7OrzGyVma1qbGzMdzlpKQ2qe7SY/eHJ7exo6eWLrz8av8/yXU5qIr2w+SGvlc2yU3NFSYCyoJ/GzmGhDbTIrojICPId2vaZ2SyAxPeG0Z7onLvWObfMObesvr4+ZwVmgtfSlpnQ5pzz9muUSWdPey+/emQzv350C2t2thGJxXnwlQa+e/c6XnVEPeccUUC/t1segWhv1sazJc2oLmFfR9/+AzXzoHoO7FBoExEZLrOrZabvduBK4JrE99vyW052eBMRMtM9+pP7N/KTBzbwl4+cwYnzajNyTpm4Jzc388HfPU13eH849xnEHRw5o4ofvuOEPFY3DuvvhmAFLDgrq5eZUV06NLSZeTsvqKVNROQAOQttZvZH4FxgmpntBL6OF9ZuMrMPAtuBt+eqnlwqDfrpDkcnfB7nHNc9sZVY3HH9E1s5cd6JEy9OJqytJ8zH/vcZZtaU8usrT6Us5GfllhZe2dvBnNpy3nzSHMpCBTL5AMA5bzzb4vMgWJrVS82sKeXZ7W1DD85fAS/dAm07oHZeVq8vIlJIchbanHPvGuWhzK/aOcmUBn00d0+8pW1DQ9fAQqR/f3kfkVicoD/fPdxyzV3raO2JcN0HTmPhtAoALjlhNpwwO8+VjdO+l6BjJ5z7haxfakZ1KXs7+nDOYcmxc/MTe5xuf1KhTURkEP3Fz4GSoJ/+DIxpe3i9NwHjaxcvpaMvyuObmid8TpmYbc3d3LRqB1euWMgxs2sO/uRwN7x0KzzyA3j2D9C2PTdFpmv93d73Ja/J+qVmVJcSjsZp64nsPzj9GAhVab02EZFh8j2mrShkYnHdrU3dXP/ENhbXV/Du0+dz7cOb+eJf1nDnJ8+mtnySbzp+CLv24c0EfD4+fM6igz/x2T/AvV+F3sFLFRosvRRe+y2omURLgay7A2afBFXZ3xt1ZrXX/bq3o4+6isTvsT8A807VzggiIsOopS0HSoO+CW8Y/7XbX2JXWy9fuOgoSoN+rn3vKezr7OeX/9icoSolXa3dYW5evZO3nDyH6dWjjP1yDu79Ctz2zzD9aHjfHfDFnfDPT8NZn/YWsP3Z6fDyJJmD07TB2yT+2Lfl5HIzqr212vYOnowA3tZZ+16C3rac1CEiUggU2nJgouu0RWNxnt7SwhXLF/CaY7zWj+Pn1nLG4qn8fe2+TJUpabp59U76o3Hed+bC0Z/06A/g8Z/Cqf8P3ns7LDwLSqqg/gi44Ovw0Se9rZtuei888C2I53kR5uf/5G0Qf1yuQpsXdhuGh7b5ywEHO5/OSR0iIoVAoS0Hkuu0OefG9fo97X30RmIcPatqyPFzj5zOxoYudrT0ZKJMSUM87vjfp7azbEEdR82sHvlJG++H+/8Njns7vO57XrffcHUL4H1/g5OugIe/CzdeDn0d2S1+NPE4rLkJFp2bk65R2B/a9rYPW3tw7jIwv5b+EBEZRKEtB0oDfuIOIrHxhbbdbb0AzK4tG3L8vCO9xVofWl9YO0QcCh7f1MyWpm4uXz5/5Cf0tMCtH4b6o+CNPwHfQf5TC5TAJT/1gt36e+BXF0DD2uwUfjDbn4D27XD8ZTm7ZCjgY2pFiL0dvcMeqIBZxyu0iYgMotCWA8kNwvvGuWn87vaRQ9th0ypYMLWcW57ZqQ3pc+wPT26jrjzI646dNfIT/v516GmGt/4KQuVjn9AMTr8K3nubN1nh2vNg9e+8MXG5svp33qzNoy/O3TWBOXVl7GrrO/CB+Stg1yqIhg98TESkCCm05UBp0HubxzuubXfiD9rsmqGhzcy46lWLeHZ7Gx+6fvW4u18lPfs6+rhv7T7esWzeQCAfYvtKeOZ6WPFRmHlceic/7Gz48GMw7zT4v0/Cb18PWx/Lfnjr3OctR3LS5V4rVw7NrStjZ+sIXfzzV0C0z5sYISIiCm25UJL4w94/zq2sGjr6qC4NjLiq/uWnL+Bzrz2Sh9c3sm5v54TqlNT86akdxOKOd502QtdoPA53ftbbP/Ocq8d3gaoZcMWt8MYfQ9N6+N3r4b+Wwe2fgOf+CK1bMx/iVv8W4hE49UOZPW8K5taVs6u198B/dCw8CzBvH1QREVFoy4WB7tFxtrTt6+gfGLA9kree7K3x9eRmLbabbdFYnD89vZ2zl0wb2P1giBdvhr0vwIX/BiWV47+Qzw+nvA8+9YIX3qYshpf/Cn/9MPz4BPjB0XDLVfDy7RCLjHW2g4uGYdVv4PALYdrhEzvXOMytK6M/Gqexa9hkhPIpXkvlln/kvCYRkclIi+vmQGkg2T06zpa2zj6mJ9azGsmM6hKmVZbw4q48zTosIg+sa2BPex9ff+MxBz4YDcOD3/KCxjFvycwFQ+VeeDvlfV4rXuM62P44bHsCNtwHa26EmvnellMnXu6NjUvXM9dB1z5Y/pHM1JymuXVet//O1l6mVw37x8lhr4Kn/gcivRAsG+HVIiLFQy1tOZBsaRvvZIGGzv4D/5gNYmYcO6eal3a3j+v8kro/rNzOzOpSLjh6+oEPPnu913X56q8dfLboePl8MGOpt+bb234N/7IB3nWj15162z/D9ZdC5970ztnfBf/4Diw4Cxa/OvM1p2BenTdRY8Slaw47B2L9sOOpHFclIjL5KLTlwP7u0fRb2pxzNHT0H7SlDeDY2TVsaOia8HZZMrrtzT08vL6Ry06bR8A/7D+dcA/843ve4PklF+amIH8AjrwIPnAvXPwj2LnKm3W6a3Xq53jiZ9DdCBd8Y3ytdBkwZ1BL2wEWrPDWa9vycI6rEhGZfBTacmAis0fbeyOEY/GDtrQBHDe3hljc8dJudZFmyw1PbcPvMy47dYQJCKt/C1174fyv5T78+Hyw7P3wwXvBF/BmnL5w89iva9oIj/4Qjr7E2+szT8pDAaZWhEYObSVVMOcUhTYRERTacmIi67Q1dHqDs6dXHbyl7fi5NQA8v6Mt7WvI2PoiMf68aicXHD2dmTXDAnS4Bx79kdeVt+CMvNQHwMxj4aoHvZDzlw/CA98cfVuscA/85QMQCMHrvpvbOkcwb0o525q7R37wsLO91sN+zY4WkeKm0JYDpYHxd4/uS+zJOFZom1ldyuL6Cv6wchtbm0b54yfjducLe2jpDnPF8oUHPrj6t9DdAOeOc4mPTKqYBlf8NbEt1vfgj5dBx+6hz+nvgpuugD1r4M3XQvUoCwTn0KL6CjY3jhbaXgUuBtsez21RIiKTjEJbDkyke7Shw2tpO9iSH+BNRvjqxUtp7OznA797mlhcC+1m0u+f3Mai+grOPHzq0AcivfDYj/PfyjZYILR/W6wtD8N/nQb3fsVbHuTJX8AvzvD2Rb3kJ96YuElgcX0lezv66OqPHvjgvOUQLIcN9+a+MBGRSUShLQdKJrBO20D36BgTEcDbQP6rb1jK5qZuNjV2pX0tGdmLu9p5dnsb7zl9ATZ8vNqq33rLZUyGVrbBkttiffRxWHwePPFzr3Xt7quhtAbedwec/N58Vzlgcb235t2WkVrbgqWw6DxvX1bt+iEiRUzrtOVAsqWtPzq+7tHKkgDlodT+rzp5QS3gjW07YkZV2teTA13/xFbKgn7eesrcoQ9EeuGxH3ndd5OllW24KYvgnb/3xoO1bIayKVAzN28zRUezuN5biHhzUxfHJcZnDnHkRfDKHbDvJW/snohIEVJLWw6E/D7MxtfS1tg59nIfgy2aVkllSYA1O7VmWya090S47bndvOmkOdSUBYc+uPp3XivbeLeryqWSKph1AtTOm3SBDWD+1HJ8BpsaRmkhXvIa7/v6u3NXlIjIJKPQlgNmRmnAP87u0b4xJyEM5vN5C+2u2aXQlrJn/wDfPwq+NRvu+gKE93fR/Xn1Dvqjca5YvmDoayK93nIZC8+GhWfmuOBDT0nAz/wp5WwabTJC1UyYfZLXRSoiUqQU2nKkNOgb1+zRsXZDGMnxc2tZu7uD8Di6Y4vO6uu83QTqDoOj3wgrf+nNuIxFiMTiXP/ENpYtqGPp7Ophr/vd5BzLVsAOn17FK/sOsqzHERfBzqehqzF3RYmITCIKbTlSGky/pc05x76O9FrawFuzLRyL88perWt1UC2b4a7Pe9s3vfc2eMt/w6U/92Zc3vNl/vjUdra39PDhcxYPfV2k11uXbeHZsPCsvJR+KDp2TjWbGrvoCY8wgxS80IaDjffltC4RkclCoS1HSoN++tJs+eroi9IXiR+4mOsYTllQB8DKLc1pva7o3P1F8AW9oBYIecdOuhyW/zM89d9svucXrFg0lfOH7zO6+nfe7gfnfCHnJR/Kjp1dg3Owds8ou3rMOgGqZsPav+W2MBGRSUKhLUfKgn56R2tBGMXAwrpjrNE23KyaMo6YUckP7lvPt+9cqzXbRrL7WW9Q+1mfPGBx2YblX2J14CS+zP/w/eU9Q5f5aNsOD37bW5ftsLNzXPSh7dg53qzRF3eNEtrMYOmlXktbn8ZsikjxUWjLkcqSwMgLhx5EMrTNTDO0AXzt4mOYVlnCtQ9v5p6X9qb9+kPeoz/01is77Z+GHH5hZzuX/mIlH+3/GJGq+cy+833ezgEAvW3w5/eBi3sL00pGzaguYVpliBcPNonmuLdBLKzWNhEpSgptOVJZOp7QltwNIb0xbQBnLZnGg/9yLtMqQwptw3Xs8f7on/xeKN0/weCRDY287ZeP4zPj1x++kPIP3AbBCvjVBd5khV+d7wW4t1wLdQvzV/8hysw4ZnbNwZermXMK1C6AF/6cu8JERCYJhbYcqSwJ0NU3zu7RNGePJvl9xpmHT+OxjU3Ex9tFun0lPHM9dDeN7/WT0bO/9/ayPOX9A4d2tfXy8T8+y2HTKrjtY2d6XXV1C+CfHvZmla79G4Qq4D03w1FvyGPxh7bTDpvCK/s6ae7qH/kJZnDCu2DzQ15XtYhIEVFoy5HxtbT1UV0aoCzkH/d1zzp8Gk1dYdbuHWWc0ME8+B/wm9fA7R+Hn50ODevGXUe6nHO8vLtj9D/e4xWLehMJFr8apu6fFfrtO9fSH4nzy/ecwrTKQS2blfXwtl/D1du8ALfo3MzWI0OsWOzt7frk5pbRn3TS5d73Z2/IQUUiIpOHQluOVJUE6BxHS1u6M0eHO+fIegAeWNsw6nN2t/Vyzvce5EPXr2JPey83r95J+6qb4B/XwImXwwfvA58fbrwcIn0TqidV37n7FV7/k0c493sP8dLuod1l4Wh8/JMrNtwDHbtg2QcHDq3Z2cYda/bwobMPY+G0iomULRN0/JwaKksCPLH5IC27tfO9/VSf/b0XwkfR1R+loaNvXItai4hMRtp7NEdqy0P0R+N090epKBnjbY+G4aVbqG3uYkb1cRO67vSqUk6cV8utz+3iqnMWURI4sNXufx7ZzLbmHrY193Dfy/uop43Xln0eN/sU7I0/Bn8Q3vxL+P2b4eHvwflfnVBNSfG445t3rGX19la++oajWbZwCuAt+fDfD2/idcfO5LkdbXzwd6u45aNnUBb08+93vMxfn91FVWmQq161iCvPWEgs5vj72n1saerm5AW1nHfk9CEzPjc2dLJySwtnH17P/Kd/7S0bccRFA4//9IGN1JQF+dCrFmXk55LxC/h9nHbYFB7b2IxzbujM3cFO/X/wp3fD2tvh2LcMHH5lbyc3rNzG/Wsb2NXWO3B83pQyzl5SzzlH1HPekdMJBfTvVREpPAptOZKcTLCvo49Fic2xR+Qc3Px+WPc3vo2PX9T8GDh9Qtf+5/MO50PXr+Kq61fz6yuXEfAP/YN138v7uHDpDN59+nweWtfAZVt/Railn+dP/Q4n+hP7bS5+tTeW6LEfwbFvhRlLJ1QTeGHxN49toSzo58N/WM0D/3Iu1aVBvn/vK1SXBrnmLcezq62Xd/z3E1z800fxGbT3Rrj89AXsbO3he/e8wo/+vp64Y0jL2+uPm8k333QcVaUB/vOeV/ifRzYTd3B0SRN32f1w7pfA7/3qr9/XyX0v7+OT5y+hqjQ4WqmSQ68+ajpf+euLvLyng2Nmj7B5PMARr4Mpi+CJ/4Jj3kx7b5Rv3vEyNz+zk5Dfx7lH1nP58vlUlQZp7wnz3I52bnt2F/+7cjvTKkt49+nzufz0+cwYx8xsEZF8UWjLkcOne0Ht+/eu53tvP56Azzfyv/afuR7W/Q131mdpfuTXvKHtD8B7J3TtC5fO4N8uPYav3fYStz23m7eeMnfgse3NPexs7eWqVy3ivCOnc178KXj2AX7o3knztlJOPGnQiV7zLW9tszs+C++/M62Nx1dva+HXj26hvrKEz1x4JE3d/Xz/vvW89pgZfOy8JVzys0f50X0beMvJc/j72gY+e+ER1JQHqSkP8ucPr+Cbd7xMRSjAJ85fMrCe1/M72rjzxT0EfT5ec8wMjppZza8f3cIP7nuFR9Y3UVUaYHd7H+86bT6XnTqPddd/imjER8eR72RKoq5fPLSJ8pCf952xcELvsWTOG46bxTdu935XRw1tPh+s+Ge447NsWnkHlz9QRmNXPx86exEfPXcxteWhA14SicV5bGMTv39iGz99YAM/f3AjFx07k/edsZBTFtSN3qonIjJJmHOFt/DqsmXL3KpVq/JdRlricccnfn4r9Xsf4MnSc4hXzODuT509bOHWHfDzFTDnJBrffBO//c4n+XzwRvjwYzDz2Ald3znH6378CAB3fXL/dW9YuY0v3/oif//MORxeHfMmHJRP4+NVP+CJre08+cXzh7bMrb4O/u8T8KZfwInvHvV6N6zcxg/uXc+s2lIW11fyf8/vpqo0SFd/lFk1pVSEAuxp7+Xvnz2H6VWlfOWvL3DDyu3Mrimjsy/Co1e/mupxtnxt2NfJzx/aRHd/lLcvm8eFS2dApI/Y94/i7z1L+PPib3PtFcvY1tLDBT/4B+8/YyFfuXjiLYeSOR/83dO8tLuDR75wHkH/KF2Z0X56vn88G3oq+Gjp9/jFFadw/NzalM6/rbmb3z+xjZtW7aCjL8rSWdW874yFXHLibEqD45/4IyIyEjNb7ZxbNtHzaGBHjvjiEX7Q/zW+Hvg9Pw5/nW37mtja3LP/Cc7h/vZpnIvDJf/Fvs4wf4idT9RfDo//dMLXNzOuWLGAdXs7+dUjW3jvb57itud28cj6JmbXlLK4vgLu+7q3CfolP+H1J8yjqSvMyi3DZvGddAXMPQ3u/Qp07gNgc2MXZ17zAG/62WO0dId5cF0DX771ReZNKcfv83H/2gbedspcHrv61dz0TysI+X3sbO3hu287YWA5ky+/finLFtTR0h3mu287ftyBDWDJjCp++M4Tufa9y7zABrD2dvx9rbhTPsDf1zbwldte5NM3PkdZ0M9VGss26Vy+fD57O/r486qdIz4ejzu+d/8WvtHxRk6wTdz9mpaUAxvAgqkVfOXipTz5pfP59puPIxZ3fP4va1j+H/fzH3etZUdLz9gnERHJMbW05cqLf4GbP8BN9lre4e7hXyNX8NvY6/jRO0/kjhf28NbA41y0/qv8suwqPvgv3+HuF/fy8T8+y9On/J36tdfDJ5+HmrljX+cgOvsinHnNA3QMm8X6zmXz+M4p7XDdxbDiY/Dab9EXiXHqt/7OuUdO56fvOmnoiRrWwbXnwoIVcPlfeO/vVvHEJm+23/wp5TR1hZldW8atHz1jxFYL5xxx560jN1w87vCNcByAjt2w7g5vfa6KelhyIUw/euwf3Dmv3v4O3MdW8Y3/W8t1T2wj6Dd+ctlJvO64WWOeQnLLOcdbf/E4u9p6uffT51BTtj/Ed/dH+cxNz3HPS/u4/NRZ/HvDJ/D1NME/PzVkseR0r7dySwvXP7GVe17aR9w5zj9qBu87YyFnHj5VXaciMiGZamlTaMuVP7wVGl/hlcseY9Zt76B3zzpe1f8j+glRSj8PlnyWfa6Wt4T/DZ/PTzQxsH7tZ5dS9otlcMr74A3fn3AZT29t4Z4X93L58gVc+ZunaGjv4u8XNjD38a97a5J9+DEIlQPwH3et5dqHN/PtNx/HZafOG/jDFYnF8T/zO3x3fJq9p3+J5f84li9cdBRHzqzkE398jtryINd94DQWH2zCRVI87o1POpj2Xd7yI8/e4C2K6w95WxmBt/DtG37o1T6atf8HN74HLv0ZnPQewJtlWFMWnPCSKpI9z2xv5Z3//QSnLKjjZ+8+mbryEI9ubOIbt7/E1uZuvnrxUt53xkJs1zPebhUnXg5v+tmEr7u7rZf/XbmdPz61nebuMIvrK7jyjIW85eS5VI4181tEZAQKbYUU2jr3wg+OhrM+4y2XseVhuO6N3FT/CT6/Yzkf8d/OF4J/4h39X+Upt7/l6PXHzeTnl58C//cpeO4G+MjjMG1JZmrq2E3fyt8QfO73+Lv3wuyT4Z1/gJo5A0/pCUf50PWreGxjM/96yTFcecZCXtrdzpW/eYrKkJ975vya4IY7+Wjsc3zni5+jpjxIXyRG0O8bsRVtwM7V8NS1sP4uCHfDYa+CY94MR10M5VP2P6+nxdsj9KlrIR6DUz8Ip37IWxS3qwFW/xYe+QGU1XkL4C4868BrxSLwy7MhHoWPPjkwa1QKw23P7eLTNz5HwO+jLOinvTfCvCllfOctx3PG4dP2P/H+f4dH/nPMsZbp6I/GuGPNHq57fCvP72ynsiTAW0+ew3vPWJjaP0jSEI87OvuidIejBHxGZWmA8pB+V0UOFQpthRTaHv+pNwbsY6u80OUc/PZ1uJYtvNv9O7/t/RQby0/g7uN+xAnzalmzs41Pnr8EM/PCT+de+Plyb1HR99/lbac0HtEwbH8cnv+Tt3djPAaHn++tebXktSO2eMXijg9dv4pHNjRy3ftP4+pbXqC1O0x3OMpx04P8e+vnOTKwl5J3/wEOv8B7Ubgbdj0DDS9DX4e3MXtptXf8hZu9Gkqq4ehLvOOv3AmtW8H8MPdUb1/Prr2w7XEvdB3/TjjvS962UsPtfcHbxL1lM5xzNZz92aHB7IFvwcPfhXfeAEdfPL73TfJqY0Mnf3pqBz2RGKcurON1x846sNs9FoXfvwm2PwHv+pPXdZ5Bz+1o47rHt3LHmj2EY3GWL5rCW0+ey+uOm5V261tXf5Tnd7Tx3I421u3tZMO+TjY3dROOxoc8r7o0wAnzajlj8TTOWDyVY+fUHPwfQyIyaSm0FUpocw5+eRYESuFD9+8/vn0l/OY1uEApFgvDR56A6UeNfp5X7vIWE51/Bpzzea91qXIGVE7fv/RGPOZ1G8Yi0NfuTSro2OV1L+5YCZsegP4ObxP0k6+A0z8MUw4b80do7Q5z8U8fZVdbLz6DG/9pBZsaurj6lheotzYenf1TSprXwmHnQH8n7F3jtWyNpGY+LP+Id/2Sqv3v0Z7nvYVStzzsTXAonwLzV3ibuo+1Jlx/F9zxGVhzo7eh+IX/BtOXwtO/hge/CSe8G978izF/Tilwfe3wuzd4Yy7f+OP9211lUGNnPzc+vZ2bV+9ka3MPZUE/ZyyeyorFU1m+aCqHT68cEijD0Tg7Wnt4bnsbq7e38sy2Vtbv6yS5rODcujKOmFHF4dMrmVFdSkXIGxrR2RdlZ2sPq7a28sq+TgCqSgOcfthUzlg8lTMPn8YRMyoLYqxdNBansy9KZ1+UcCyGL/GPUZ8ZJQEfUypCB6wdKXKoOaRCm5ldBPwY8AO/cs5dc7Dn5zW07VwFjeu8INSx09tIPR7zxlq5+P4vzAtVkR5Y9zfvj8gp7xt6ruf/5HV7nvReOP7tY197zU1w+ycgun+ld8rqvIDS1QCtW0YPS5Uz4IjXeouSLjon7da6jQ1d/Meda7lg6Qzeddp8AB7f2ERJ0McpM0Pw2I+9FrPSGph3OsxfDrNO9Orra/fCoj8E1XPGHsM2Hs7BS7fAnZ+Dnub9x495M7z5WggcuG6XHIJ6W+Gm93rh/+hL4IJvDNljNlOcczyzvY2/PruLxzY2sbmpe+Cx2vIgpQE//dEY7b2RgYBWVRLgxPm1nDy/jpMX1HHivNohEyxG09jZzxObm3liUxOPb2pmW2LW+dSKECsWT+WMxdNYtrCOxfWVWW+Jc87RE47R0h2mqauflu4wzV1hmrvDNCfuN3WHaenup7krTHtvhJ7wwbcRM4OpFSXMqC5hRnUpC6aWc+SMKo6YWcURM6ryPo6wNxyjobOPhs5+9nX00dDRT2NXPz39UfoiccKxOGYMhNCasiC15UFqyoLUlIWoLQ8yo7qU2bWlI+5II8XhkAltZuYH1gMXAjuBp4F3OedeHu01eQ1tN7wdNtzr3a6o98KQLwDmG/rl4l4XX38XnHCZtzBtJsJK515oWu8FoY7dsO9FaHzFa3GbergXmnwBr/uxcgZUz/ZmnZbVpbUYbsEKd3sLAHfs9sbpLTwz3xVJrsUi3pCEf3wHov3emMklF8KMY73/ZsuneP999ndBuMtrHQ537b8f7vJeF+1PtFyHvf+mAiVei7k/5P13VlE/8LUvVsXKPVG2NvfQ2NlPXyRGadBPXUWIeXVlHD+3lsOnZyZU7Wzt4fFNzTyxqZnHNjbR0NkPQEXIz7Fzalg6u5rZNWXMrCllZmJNxKDfCPh9+Ax6IzF6wjF6+mP0hKPe7fD+293hKL3hGN39MXojUbr6Y7T1JMNZP32R+Ih1lQX9TK0MMbUixJSKEFMqSqgrD1JVGqS6LEBVaZCSgI+4c0Rjjphz9EfjNHb209jZx76Ofva297G1uXtI0JtbV8bSWdUsnV3NMbOTP19pRlsZ+yIxtjZ3s7Wpm81N3vetTT1sbuqmqav/gOcHfEZFSYDSoLdIunPevxt7I15QH21v5PqqEubUljGnrsz7nvyqK2NWTSk1ZcGCaD2V9B1KoW0F8A3n3GsT978I4Jz7j9Fek9fQ1rzJCz/Vc7wPcRGZnLoa4Kn/8Vpgmzem/3pfAPwl3t678RjE+vfPWh7x+cFEiJuWCIdTvYBYPtX7R1P5FCib4rVyB0ogUAbBUu91JP7q4xKt9W5/q3085rWgu8T3eGzgmItH2dPWzcaGHjY29bK+qZetLX10R4wYPqL4ieHbf9v5iSZuRwkQwZ+47Qe8sBAK+KgI+SkPBSgP+SkPeeFzSkWIaZUlTKnwgpkX0BL3K0MZmzgRjzt2tvbyyr5O1u/r5OU9Hazd08GWpm6Sf65qy4MsnVXNkTOrmFtXzuyaUmbXljG7toypFaEDlg2KxR0dvRGau/vZ2tTD1uZutjR5X1ubutnd3jfk+dMqS1g0rYKF08pZMLWC6VVeK+D06hKmV5VSWxYcdWki5xzdYS/otvdGaOuJsKe9j91tvexq7WVX2/6v4eMYS4M+ZtV4AW5mTSkzqr1rDW+5qy4LEvQbIb+PoN9HwG8Efb6BmpxzxOJeMI7HIRqPE49DLHE8nnx80G3vOwO3487hEkszJa8RSnwP+Hz7bw+79qHEOUc07v0jIxyLE0l+RYfdj8WJxNzA7XDUDXksHHO8d8XCQya0vQ24yDn3/xL3rwBOd859bLTXFNSYNhHJv64Gr0W6p9n78gWgpBJCia/k7ZKqRKgqBd8IXVnxuBfe+tqhuzHx1TTo9qD7Pc3Q0wr97bn/ecfBmR/8QcwX9Cbz+IJeYPX5vd6DEY3yh3rU1qIRjqf43Ljz/lD2R70Wuv5ojHA0PhDk3KDnmyVO67zj8WF/5xw2EEZCAR8hf/K2n2DAh/+AmjLwcw477IBozBGJOyJRRyS+/49/NPGHPhpzOA6sfSQOI/EjDzueen2jn3tkQ97zxP/48N58M++YmQ18J3HbZ4njiRcNfE+cw9h/w5KPJG4nr+m8/3MHanOJ24PvJ2t3Q57niDvD4a0X6hLh1LtvA2E1uZ5oqj//WMfP+eZjGQltk2FO+Ug/9QHvkZldBVwFMH/+/GzXJCKHksrp3tdE+XzgK4NgGVTNTO010bA31q63xVvGJtIL0T7vK9IL8Qgk/tB54Sjx3cybUe1LfgUSQzGGHTM/4BKtcENb4oYcczGv6zh5LBbxrh3z7ls8sv/x5GPxqPf4SH+2Rv0H/yjHR3x+6s/1AaU4Bq+s6IBINE5vJEpvOE5fJEZ/NO61FMUdhhuYhR/0G6GAj/Kgn4oSP6EDJj+kU/dBpPi+GBB0jpFHNbrE/3qtX5FYnGii1SYaixON7w8cXsjwgkbcuSEBx2cMCUJeeBp8zIYeG/J4IgQmzw9DAk0y1DiSxwY/N1lfosbE+5J8TSxxPhfff+7BAQznSLxqSPhK3klm5UTMw2du//1BP4tv4HlD34NgMjQmf17zJsZ475kNBMohAXPE24n3MHHbO4cbcu7ksUyZDKFtJzBv0P25wO7hT3LOXQtcC15LW25KExGZoEAIqmZ4X5JRBoQSXzV5riUbDO+P9GT4Qy0T9OHMdB9PhnnWTwNLzOwwMwsBlwG357kmERERkUkl7wHeORc1s48B9+At+fEb59xLeS5LREREZFLJe2gDcM7dCdyZ7zpEREREJqvJ0D0qIiIiImNQaBMREREpAAptIiIiIgVAoU1ERESkACi0iYiIiBQAhTYRERGRAqDQJiIiIlIAFNpERERECoBCm4iIiEgBUGgTERERKQAKbSIiIiIFQKFNREREpAAotImIiIgUAIU2ERERkQJgzrl815A2M2sEtuXh0tOApjxct1Dp/Uqd3qvU6b1Kj96v1Om9So/er9Qd6ZyrmuhJApmoJNecc/X5uK6ZrXLOLcvHtQuR3q/U6b1Knd6r9Oj9Sp3eq/To/Uqdma3KxHnUPSoiIiJSABTaRERERAqAQlt6rs13AQVG71fq9F6lTu9VevR+pU7vVXr0fqUuI+9VQU5EEBERESk2amkTERERKQAKbQlmdpGZvWJmG83s6hEev9zM1iS+HjezEwY9ttXMXjCz5zI1Q2QyS+G9OtfM2hPvx3Nm9rVUX3uoSeG9+tyg9+lFM4uZ2ZTEY8X2e/UbM2swsxdHedzM7CeJ93KNmZ086LGi+r2ClN4vfWYlpPBe6TMrIYX3Sp9Zg5jZPDN70MzWmtlLZvbJEZ6Tuc8u51zRfwF+YBOwCAgBzwNLhz3nDKAucft1wMpBj20FpuX755hE79W5wN/G89pD6Svdnxd4I/BAMf5eJX7eVwEnAy+O8vjrgbsAA5Yn/xsstt+rNN4vfWal/l7pMyvF92rYc4v6MyvxM88CTk7crgLWj/A3MWOfXWpp85wGbHTObXbOhYE/AZcOfoJz7nHnXGvi7pPA3BzXOFmM+V5l6bWFKN2f913AH3NS2STknHsYaDnIUy4FrneeJ4FaM5tF8f1eAWO/X/rM2i+F363RFN3vVprvVVF/ZgE45/Y4555J3O4E1gJzhj0tY59dCm2eOcCOQfd3cuCbPtgH8VJzkgPuNbPVZnZVFuqbTFJ9r1aY2fNmdpeZHZPmaw8VKf+8ZlYOXAT8ZdDhYvq9SsVo72ex/V6NRzF/ZqVKn1lp0GfWgcxsIXASsHLYQxn77CrIHRGywEY4NuK0WjM7D+8D8KxBh890zu02s+nAfWa2LvGvlUNRKu/VM8AC51yXmb0e+CuwJMXXHkrS+XnfCDzmnBv8L9xi+r1KxWjvZ7H9XqVFn1kp0WdW+vSZNYiZVeIF2E855zqGPzzCS8b12aWWNs9OYN6g+3OB3cOfZGbHA78CLnXONSePO+d2J743ALfiNXkeqsZ8r5xzHc65rsTtO4GgmU1L5bWHmHR+3ssY1s1QZL9XqRjt/Sy236uU6TMrNfrMGhd9ZiWYWRAvsN3gnLtlhKdk7LNLoc3zNLDEzA4zsxDeL+Ptg59gZvOBW4ArnHPrBx2vMLOq5G3gNcCIs24OEam8VzPNzBK3T8P7PWtO5bWHmJR+XjOrAc4Bbht0rNh+r1JxO/DexEys5UC7c24Pxfd7lRJ9ZqVOn1np0WfWfonfm18Da51zPxjlaRn77FL3KOCci5rZx4B78GZz/MY595KZfTjx+C+BrwFTgZ8n/tuOOm+j3BnArYljAeB/nXN35+HHyIkU36u3AR8xsyjQC1zmvKkyI742Lz9IDqT4XgG8GbjXOdc96OVF9XsFYGZ/xJvFN83MdgJfB4Iw8F7diTcLayPQA7w/8VhR/V4lpfB+6TMrIYX3Sp9ZCSm8V6DPrMHOBK4AXjCz5xLHvgTMh8x/dmlHBBEREZECoO5RERERkQKg0CYiIiJSABTaRERERAqAQpuIiIhIAVBoExERESkACm0iIiIiBUChTURERKQAKLSJSFaZ2Q/N7FOD7t9jZr8adP/7ZvaZDF+zK8PnqzWzjw66v9DMUlrt3czKzOwfZubPQB0hM3vYzLQwukgRUmgTkWx7HDgDwMx8wDTgmEGPnwE8loe60lELfHSsJ43iA8AtzrnYRItwzoWB+4F3TvRcIlJ4FNpEJNseIxHa8MLai0CnmdWZWQlwNPCsmf3VzFab2UtmdlXyxWb2nWGtXN8ws8+a2XvM7Ckze87M/nuklqzRnpNoKVtrZv+TuN69ZlaWeOyrZrbOzO4zsz+a2b8A1wCLE+f5XuL0/pFeP4LLGbpH40NmdmTi9lQzezFRzzoz+1Xi/g1mdoGZPWZmGxL7YSb9NXFOESkyCm0iklXOud14ezjOxwtvTwArgRXAMmBNogXpA865UxLHPmFmUxOn+BNDW5beAaxKHDvTOXciEGNYkDGzo8d4zhLgZ865Y4A24K1mtgx4K3AS8JZELQBXA5uccyc65z432uuH/+yJTaAXOee2Djp8OLAhcft44IVBx3+cOHYU8G7gLOBf8PYyTHoROHX4tUTk0KdxESKSC8nWtjOAHwBzErfb8bpPwQtqb07cnocXipqdc8+a2XQzmw3UA63AccApwNOJDarLgIZh1zx/jOdscc49l7i9GliI13V7m3OuF8DM/u8gP9NIrx9uGl6gI3G+BcAu51w8ceh4YM2g872QeN5LwP3OOWdmLww+t3MuZmZhM6tyznUepD4ROcQotIlILiTHtR2H11K0A/gs0AH8xszOBS4AVjjneszsIaB00OtvBt4GzMRreTPgOufcFw9yzbGe0z/odgwv1FkaP9NIrx+ul6E/x4nsD2nghcobRzhffND9OAd+VpcAfWnUKiKHAHWPikguPAZcDLQ452LOuRa8wf0r8LpLa4DWRGA7Clg+7PV/Ai7DC2434w3Gf5uZTQcwsymJVqzBUnnOcI8CbzSzUjOrBN6QON4JVKX7QzvnWvHGviWD2wkkQpyZLQEuZX/3aEoS3caNzrlIuvWISGFTaBORXHgBr6vwyWHH2p1zTcDdQMDM1gD/Pux5OOdewgtNu5xze5xzLwNfAe5NvOY+YNaw14z5nOGcc08DtwPPA7fgjZ1rd841A48lJgl872DnGMG9eGPTwGtp85nZ88DXgLXAlWme7zzgzjRfIyKHAHPO5bsGEZFJw8wqnXNdZlYOPAxc5Zx7ZgLnOwn4jHPuCjPbCJw0kbFoZnYL8EXn3CvjPYeIFCaNaRMRGepaM1uK14153UQCG0BiIsWDZlYDxCcY2ELAXxXYRIqTWtpERERECoDGtImIiIgUAIU2ERERkQKg0CYiIiJSABTaRERERAqAQpuIiIhIAVBoExERESkACm0iIiIiBUChTURERKQA/H80V4r+dlh6eAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x540 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "air_flux_out = mp.get_fluxes(air_flux)\n",
    "substrate_flux_out = mp.get_fluxes(substrate_flux)\n",
    "flux_freqs = mp.get_flux_freqs(air_flux)\n",
    "\n",
    "\n",
    "lam = 1.0/np.array(flux_freqs)\n",
    "d_lam = -1*np.diff(lam)\n",
    "d_lam = np.append(d_lam, d_lam[-1])\n",
    "\n",
    "#Plot the wavelength spectrum:\n",
    "fig = plt.figure()\n",
    "plt.plot(lam, air_flux_out/d_lam)\n",
    "plt.plot(lam, -1.0*np.array(substrate_flux_out)/d_lam)\n",
    "plt.xlim(0.15, 2.0)\n",
    "plt.xlabel('Wavelength ($\\mu$m)')\n",
    "plt.ylabel('Amplitude ($\\propto$W/um)')\n",
    "\n",
    "fig.set_size_inches(10, 7.5)\n",
    "\n",
    "display(substrate_material_eps)\n",
    "display((1/0.1 + 1/0.9 - v_e/0.1)/(1/0.9 - v_e/0.1)/v_e)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "e501ac15-1061-4803-b465-5e0ff4284f29",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Wrong number or type of arguments for overloaded function 'dft_flux_save_hdf5'.\n  Possible C/C++ prototypes are:\n    meep::dft_flux::save_hdf5(meep::h5file *,char const *)\n    meep::dft_flux::save_hdf5(meep::h5file *)\n    meep::dft_flux::save_hdf5(meep::fields &,char const *,char const *,char const *)\n    meep::dft_flux::save_hdf5(meep::fields &,char const *,char const *)\n    meep::dft_flux::save_hdf5(meep::fields &,char const *)\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_144171/3808627645.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mair_flux\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msave_hdf5\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'test'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/envs/pymeep/lib/python3.9/site-packages/meep/__init__.py\u001b[0m in \u001b[0;36msave_hdf5\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m   2496\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2497\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0msave_hdf5\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2498\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_meep\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdft_flux_save_hdf5\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2499\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2500\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mload_hdf5\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Wrong number or type of arguments for overloaded function 'dft_flux_save_hdf5'.\n  Possible C/C++ prototypes are:\n    meep::dft_flux::save_hdf5(meep::h5file *,char const *)\n    meep::dft_flux::save_hdf5(meep::h5file *)\n    meep::dft_flux::save_hdf5(meep::fields &,char const *,char const *,char const *)\n    meep::dft_flux::save_hdf5(meep::fields &,char const *,char const *)\n    meep::dft_flux::save_hdf5(meep::fields &,char const *)\n"
     ]
    }
   ],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}