phcr-cavity-Control.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial on MEEP simulation run with phcr-cavity.py\n",
    "The tutorial is based on __[documentation example](https://meep.readthedocs.io/en/latest/Python_Tutorials/Resonant_Modes_and_Transmission_in_a_Waveguide_Cavity/)__. \n",
    "\n",
    "To run the simulation we have to address the python translator by running system command with <b>!</b> symbol. \n",
    "\n",
    "We run simulation code <b>phcr-cavity.py</b>\n",
    "\n",
    "The console output is saved at <b>Output.txt</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! python phcr-cavity_sim.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using MPI version 3.1, 1 processes\n",
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000931978 s\n",
      "Working in 3D dimensions.\n",
      "Computational cell is 16.1 x 6 x 6 with resolution 20\n",
      "     prism, center = (0,0,0.129862)\n",
      "          height 12.2654, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-6.13272,0,0)\n",
      "          (-6.13272,0.112464,0.194793)\n",
      "          (-6.13272,-0.112464,0.194793)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.06248,0,0)\n",
      "          (-0.06248,0.37488,0.649311)\n",
      "          (-0.06248,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0.30976,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (0.24728,0,0)\n",
      "          (0.24728,0.37488,0.649311)\n",
      "          (0.24728,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-0.30976,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.37224,0,0)\n",
      "          (-0.37224,0.37488,0.649311)\n",
      "          (-0.37224,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0.61952,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (0.55704,0,0)\n",
      "          (0.55704,0.37488,0.649311)\n",
      "          (0.55704,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-0.61952,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.682,0,0)\n",
      "          (-0.682,0.37488,0.649311)\n",
      "          (-0.682,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0.92928,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (0.8668,0,0)\n",
      "          (0.8668,0.37488,0.649311)\n",
      "          (0.8668,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-0.92928,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.99176,0,0)\n",
      "          (-0.99176,0.37488,0.649311)\n",
      "          (-0.99176,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (1.23904,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (1.17656,0,0)\n",
      "          (1.17656,0.37488,0.649311)\n",
      "          (1.17656,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-1.23904,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-1.30152,0,0)\n",
      "          (-1.30152,0.37488,0.649311)\n",
      "          (-1.30152,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     ...(+ 32 objects not shown)...\n",
      "time for set_epsilon = 3.53367 s\n",
      "-----------\n",
      "creating output file \"./phcr-cavity_sim-eps-000000.00.h5\"...\n",
      "creating output file \"./phcr-cavity_sim-ey-slice.h5\"...\n",
      "on time step 1 (time=0.025), 10.1679 s/step\n",
      "on time step 16 (time=0.4), 0.336156 s/step\n",
      "on time step 32 (time=0.8), 0.332429 s/step\n",
      "on time step 48 (time=1.2), 0.322165 s/step\n",
      "on time step 64 (time=1.6), 0.332871 s/step\n",
      "on time step 80 (time=2), 0.332094 s/step\n",
      "on time step 96 (time=2.4), 0.333506 s/step\n",
      "^C\n"
     ]
    }
   ],
   "source": [
    "! python phcr-cavity_sim.py -defect false"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 16 python phcr-cavity_sim.py -defect false -r false -pixavg true > Output-norm-ph-cav-flux-pixel-avg.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 16 python phcr-cavity_sim.py -defect true -r false -pixavg true > Output-defect-ph-cav-flux-pixel-avg.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 16 python phcr-cavity_sim.py -defect true -r true -pixavg true > Output-defect-ph-cav-Q-factor-pixel-avg.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 10 python phcr-cavity_sim.py -defect true -r true > Output-defect-ph-cav-Q-factor.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 10 python phcr-cavity_sim.py -defect true > Output-defect-ph-cav-transmission.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 10 python phcr-cavity_sim.py -defect false > Output-norm-ph-cav-transmission.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using MPI version 3.1, 1 processes\n",
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000907898 s\n",
      "Working in 3D dimensions.\n",
      "Computational cell is 16.1 x 6 x 6 with resolution 20\n",
      "     prism, center = (0,0,0.129862)\n",
      "          height 12.2654, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-6.13272,0,0)\n",
      "          (-6.13272,0.112464,0.194793)\n",
      "          (-6.13272,-0.112464,0.194793)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.06248,0,0)\n",
      "          (-0.06248,0.37488,0.649311)\n",
      "          (-0.06248,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0.30976,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (0.24728,0,0)\n",
      "          (0.24728,0.37488,0.649311)\n",
      "          (0.24728,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-0.30976,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.37224,0,0)\n",
      "          (-0.37224,0.37488,0.649311)\n",
      "          (-0.37224,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0.61952,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (0.55704,0,0)\n",
      "          (0.55704,0.37488,0.649311)\n",
      "          (0.55704,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-0.61952,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.682,0,0)\n",
      "          (-0.682,0.37488,0.649311)\n",
      "          (-0.682,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (0.92928,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (0.8668,0,0)\n",
      "          (0.8668,0.37488,0.649311)\n",
      "          (0.8668,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-0.92928,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-0.99176,0,0)\n",
      "          (-0.99176,0.37488,0.649311)\n",
      "          (-0.99176,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (1.23904,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (1.17656,0,0)\n",
      "          (1.17656,0.37488,0.649311)\n",
      "          (1.17656,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     prism, center = (-1.23904,0,0.432874)\n",
      "          height 0.12496, axis (1,0,0), sidewall angle: 0 radians, 3 vertices:\n",
      "          (-1.30152,0,0)\n",
      "          (-1.30152,0.37488,0.649311)\n",
      "          (-1.30152,-0.37488,0.649311)\n",
      "          dielectric constant epsilon diagonal = (4.84,4.84,4.84)\n",
      "     ...(+ 32 objects not shown)...\n",
      "time for set_epsilon = 3.61723 s\n",
      "-----------\n",
      "creating output file \"./phcr-cavity_sim-eps-000000.00.h5\"...\n",
      "creating output file \"./phcr-cavity_sim-ey-slice.h5\"...\n",
      "on time step 1 (time=0.025), 10.206 s/step\n",
      "on time step 16 (time=0.4), 0.346382 s/step\n",
      "on time step 32 (time=0.8), 0.34019 s/step\n",
      "on time step 48 (time=1.2), 0.331899 s/step\n",
      "on time step 64 (time=1.6), 0.32357 s/step\n",
      "on time step 80 (time=2), 0.334761 s/step\n",
      "on time step 96 (time=2.4), 0.321274 s/step\n",
      "on time step 112 (time=2.8), 0.324389 s/step\n",
      "on time step 128 (time=3.2), 0.320441 s/step\n",
      "on time step 144 (time=3.6), 0.331495 s/step\n",
      "on time step 160 (time=4), 0.332482 s/step\n",
      "on time step 176 (time=4.4), 0.332585 s/step\n",
      "on time step 192 (time=4.8), 0.332406 s/step\n",
      "on time step 208 (time=5.2), 0.335781 s/step\n",
      "on time step 224 (time=5.6), 0.332243 s/step\n",
      "on time step 240 (time=6), 0.320944 s/step\n",
      "on time step 256 (time=6.4), 0.332424 s/step\n",
      "on time step 272 (time=6.8), 0.33216 s/step\n",
      "on time step 288 (time=7.2), 0.336336 s/step\n",
      "on time step 304 (time=7.6), 0.324784 s/step\n",
      "on time step 320 (time=8), 0.331819 s/step\n",
      "on time step 336 (time=8.4), 0.320829 s/step\n"
     ]
    }
   ],
   "source": [
    "! python phcr-cavity_sim.py -defect false"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For optional parallel run on several core we use mpi.\n",
    "\n",
    "We may execute  <b>mpirun -np 6 </b>  routine for running the code with 6 cores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 6 python phcr-cavity.py > Output-Parallel-6.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mpirun -np 6 python phcr-cavity_sim.py > Output-Parallel-6.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize the simulation results\n",
    "The results of simulation (electric and magnetic vectors and given times) can be saved in h5 file. To look at them we can use program <b>h5topng</b>\n",
    "\n",
    "First we visualize the structure for simulation, which is stored in <b>phcr-cavity-eps-000000.00.h5</b> by executing the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm create*.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm phcr*.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm norm-phcr*.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "! h5topng -S 6 -y 55:65 create_geometry_v3-eps-000000.00.h5 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "! h5topng -S 6 -z 64:74 phcr-cavity_sim-eps-000000.00.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "! h5topng -S 6 -x 45:55 phcr-cavity_sim-eps-000000.00.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "! h5topng -S 6 -y 55:65 phcr-cavity_sim-eps-000000.00.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is <b>phcr-cavity-eps-000000.00.png</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dynamics of the field during simulation can be visualized as an image series combined into a single video."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to make a series of png frames  from h5 file by a following command:\n",
    "```Bash\n",
    "! h5topng -t 1:100 -S2 -Zc  dkbluered -a yarg -A phcr-cavity-eps-000000.00.h5 phcr-cavity-hz-slice.h5\n",
    "```\n",
    "\n",
    "For normal PC structure, then for a defected PC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h5topng error: invalid slice of HDF5 data\n"
     ]
    }
   ],
   "source": [
    "! h5topng -t 1:500 -z 70 -S2 -Zc  dkbluered -a yarg -A norm-phcr-cavity_sim-eps-000000.00.h5 norm-phcr-cavity_sim-ey-slice.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h5topng error: invalid slice of HDF5 data\n"
     ]
    }
   ],
   "source": [
    "! h5topng -t 1:500 -z 70 -S2 -Zc  dkbluered -a yarg -A phcr-cavity-defect_sim-eps-000000.00.h5 phcr-cavity-defect_sim-ey-slice.h5 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's convert all files into the video file video.mp4 and remove images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ffmpeg version 4.2.3 Copyright (c) 2000-2020 the FFmpeg developers\n",
      "  built with gcc 7.5.0 (crosstool-NG 1.24.0.123_1667d2b)\n",
      "  configuration: --prefix=/home/cloud-user/anaconda3/envs/pmp --cc=/home/conda/feedstock_root/build_artifacts/ffmpeg_1590573566052/_build_env/bin/x86_64-conda_cos6-linux-gnu-cc --disable-doc --disable-openssl --enable-avresample --enable-gnutls --enable-gpl --enable-hardcoded-tables --enable-libfreetype --enable-libopenh264 --enable-libx264 --enable-pic --enable-pthreads --enable-shared --enable-static --enable-version3 --enable-zlib --enable-libmp3lame\n",
      "  libavutil      56. 31.100 / 56. 31.100\n",
      "  libavcodec     58. 54.100 / 58. 54.100\n",
      "  libavformat    58. 29.100 / 58. 29.100\n",
      "  libavdevice    58.  8.100 / 58.  8.100\n",
      "  libavfilter     7. 57.100 /  7. 57.100\n",
      "  libavresample   4.  0.  0 /  4.  0.  0\n",
      "  libswscale      5.  5.100 /  5.  5.100\n",
      "  libswresample   3.  5.100 /  3.  5.100\n",
      "  libpostproc    55.  5.100 / 55.  5.100\n",
      "Input #0, image2, from 'norm-phcr-cavity_sim-ey-slice.t%03d.png':\n",
      "  Duration: 00:00:04.92, start: 0.000000, bitrate: N/A\n",
      "    Stream #0:0: Video: png, rgb24(pc), 644x240, 25 fps, 25 tbr, 25 tbn, 25 tbc\n",
      "Stream mapping:\n",
      "  Stream #0:0 -> #0:0 (png (native) -> h264 (libx264))\n",
      "Press [q] to stop, [?] for help\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0musing cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mprofile High, level 2.1\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0m264 - core 152 - H.264/MPEG-4 AVC codec - Copyleft 2003-2017 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=7 lookahead_threads=1 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=10 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00\n",
      "Output #0, mp4, to 'video-normal-ph-cav.mp4':\n",
      "  Metadata:\n",
      "    encoder         : Lavf58.29.100\n",
      "    Stream #0:0: Video: h264 (libx264) (avc1 / 0x31637661), yuv420p, 644x240, q=-1--1, 10 fps, 10240 tbn, 10 tbc\n",
      "    Metadata:\n",
      "      encoder         : Lavc58.54.100 libx264\n",
      "    Side data:\n",
      "      cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: -1\n",
      "frame=   51 fps=0.0 q=-1.0 Lsize=      90kB time=00:00:04.80 bitrate= 154.0kbits/s dup=0 drop=72 speed=14.2x    \n",
      "video:89kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 1.587057%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mframe I:2     Avg QP:15.54  size:  5240\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mframe P:13    Avg QP:17.35  size:  2411\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mframe B:36    Avg QP:18.09  size:  1347\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mconsecutive B-frames:  3.9%  3.9%  5.9% 86.3%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mmb I  I16..4: 22.8% 60.7% 16.5%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mmb P  I16..4:  2.5% 17.7%  0.9%  P16..4: 22.8%  7.7%  3.9%  0.0%  0.0%    skip:44.4%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mmb B  I16..4:  0.5%  2.2%  0.3%  B16..8: 20.3%  4.7%  0.9%  direct: 4.0%  skip:67.0%  L0:59.2% L1:36.2% BI: 4.6%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0m8x8 transform intra:73.9% inter:93.5%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mcoded y,uvDC,uvAC intra: 42.6% 69.0% 44.0% inter: 10.2% 17.4% 3.2%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mi16 v,h,dc,p: 74%  7% 14%  4%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mi8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 40%  7% 14%  5%  7% 12%  3% 11%  3%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mi4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 56% 10% 10%  2%  6%  8%  3%  5%  2%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mi8c dc,h,v,p: 50%  7% 34%  9%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mWeighted P-Frames: Y:15.4% UV:7.7%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mref P L0: 43.0%  8.5% 17.0% 25.0%  6.4%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mref B L0: 73.0% 20.9%  6.1%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mref B L1: 93.0%  7.0%\n",
      "\u001b[1;36m[libx264 @ 0x556ba123d040] \u001b[0mkb/s:141.66\n"
     ]
    }
   ],
   "source": [
    "! ffmpeg -y  -start_number 01 -i norm-phcr-cavity_sim-ey-slice.t%03d.png -r 10 -pix_fmt yuv420p -c:a libvo_aacenc -vcodec libx264 video-normal-ph-cav.mp4 > video_out-normal-ph-cav.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ffmpeg version 4.2.3 Copyright (c) 2000-2020 the FFmpeg developers\n",
      "  built with gcc 7.5.0 (crosstool-NG 1.24.0.123_1667d2b)\n",
      "  configuration: --prefix=/home/cloud-user/anaconda3/envs/pmp --cc=/home/conda/feedstock_root/build_artifacts/ffmpeg_1590573566052/_build_env/bin/x86_64-conda_cos6-linux-gnu-cc --disable-doc --disable-openssl --enable-avresample --enable-gnutls --enable-gpl --enable-hardcoded-tables --enable-libfreetype --enable-libopenh264 --enable-libx264 --enable-pic --enable-pthreads --enable-shared --enable-static --enable-version3 --enable-zlib --enable-libmp3lame\n",
      "  libavutil      56. 31.100 / 56. 31.100\n",
      "  libavcodec     58. 54.100 / 58. 54.100\n",
      "  libavformat    58. 29.100 / 58. 29.100\n",
      "  libavdevice    58.  8.100 / 58.  8.100\n",
      "  libavfilter     7. 57.100 /  7. 57.100\n",
      "  libavresample   4.  0.  0 /  4.  0.  0\n",
      "  libswscale      5.  5.100 /  5.  5.100\n",
      "  libswresample   3.  5.100 /  3.  5.100\n",
      "  libpostproc    55.  5.100 / 55.  5.100\n",
      "Input #0, image2, from 'phcr-cavity-defect_sim-ey-slice.t%03d.png':\n",
      "  Duration: 00:00:04.92, start: 0.000000, bitrate: N/A\n",
      "    Stream #0:0: Video: png, rgb24(pc), 640x240, 25 fps, 25 tbr, 25 tbn, 25 tbc\n",
      "Stream mapping:\n",
      "  Stream #0:0 -> #0:0 (png (native) -> h264 (libx264))\n",
      "Press [q] to stop, [?] for help\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0musing cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mprofile High, level 2.1\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0m264 - core 152 - H.264/MPEG-4 AVC codec - Copyleft 2003-2017 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=7 lookahead_threads=1 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=10 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00\n",
      "Output #0, mp4, to 'video-defect-ph-cav.mp4':\n",
      "  Metadata:\n",
      "    encoder         : Lavf58.29.100\n",
      "    Stream #0:0: Video: h264 (libx264) (avc1 / 0x31637661), yuv420p, 640x240, q=-1--1, 10 fps, 10240 tbn, 10 tbc\n",
      "    Metadata:\n",
      "      encoder         : Lavc58.54.100 libx264\n",
      "    Side data:\n",
      "      cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: -1\n",
      "frame=   51 fps=0.0 q=-1.0 Lsize=      88kB time=00:00:04.80 bitrate= 150.7kbits/s dup=0 drop=72 speed=13.3x    \n",
      "video:87kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 1.617090%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mframe I:1     Avg QP:13.02  size:  1555\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mframe P:14    Avg QP:17.48  size:  2548\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mframe B:36    Avg QP:18.08  size:  1419\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mconsecutive B-frames:  3.9%  3.9%  5.9% 86.3%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mmb I  I16..4: 42.7% 43.8% 13.5%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mmb P  I16..4:  2.7% 19.3%  1.2%  P16..4: 25.1%  7.8%  3.4%  0.0%  0.0%    skip:40.4%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mmb B  I16..4:  0.5%  1.9%  0.3%  B16..8: 21.2%  5.4%  1.1%  direct: 4.0%  skip:65.7%  L0:62.5% L1:32.9% BI: 4.5%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0m8x8 transform intra:73.3% inter:94.4%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mcoded y,uvDC,uvAC intra: 36.8% 60.5% 35.1% inter: 11.0% 17.8% 3.7%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mi16 v,h,dc,p: 77%  5% 15%  3%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mi8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 40%  5% 18%  6%  8% 10%  3%  8%  2%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mi4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 56%  8% 13%  2%  6%  6%  3%  5%  2%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mi8c dc,h,v,p: 56%  6% 31%  6%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mWeighted P-Frames: Y:28.6% UV:14.3%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mref P L0: 37.2%  9.1% 22.5% 22.1%  9.0%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mref B L0: 72.1% 20.2%  7.7%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mref B L1: 91.9%  8.1%\n",
      "\u001b[1;36m[libx264 @ 0x5646fd473d00] \u001b[0mkb/s:138.53\n"
     ]
    }
   ],
   "source": [
    "! ffmpeg -y  -start_number 01 -i phcr-cavity-defect_sim-ey-slice.t%03d.png -r 10 -pix_fmt yuv420p -c:a libvo_aacenc -vcodec libx264 video-defect-ph-cav.mp4 > video_out-defect-ph-cav.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d9405bb543004973ba1c574fdbd28aa2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Video(value=b'\\x00\\x00\\x00 ftypisom\\x00\\x00\\x02\\x00isomiso2avc1mp41\\x00\\x00\\x00\\x08free\\x00\\x01csmdat\\x00\\x00\\…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import Video, Image\n",
    "video = Video.from_file(\"video-normal-ph-cav.mp4\",play=True)\n",
    "video"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8553190fbbbb4ce582df39e4516e2ad0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Video(value=b'\\x00\\x00\\x00 ftypisom\\x00\\x00\\x02\\x00isomiso2avc1mp41\\x00\\x00\\x00\\x08free\\x00\\x01[\\xa4mdat\\x00\\x…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import Video, Image\n",
    "video = Video.from_file(\"video-defect-ph-cav.mp4\",play=True)\n",
    "video"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9b3/8deHJICyiQhcBRGwFFlEQAQVFdSqaL3uVqtWtPqzam1te9VSe2tprVdtvVqt2+VaxR1rXWrV1l4XFK2KUFFBRQFR4soiiGAgy+f3x3cmGSaTyZlkJpPkvJ+PxzxmOct8ZpI5n/Ndj7k7IiISXx2KHYCIiBSXEoGISMwpEYiIxJwSgYhIzCkRiIjEXGmxA8jVdttt5wMHDix2GCIibcr8+fNXuXvvTMvaXCIYOHAg8+bNK3YYIiJtipm939AyVQ2JiMScEoGISMwpEYiIxFybayPIpLKykvLycioqKoodirRTnTt3pn///pSVlRU7FJG8axeJoLy8nG7dujFw4EDMrNjhSDvj7qxevZry8nIGDRpU7HBE8q5dVA1VVFTQq1cvJQEpCDOjV69eKnFKu9UuEgGgJCAFpf8vac/aTSIQkeZ780147rliRyEtTYkgT0pKShg9enTtbfny5cUOqUlmz57N4YcfnvH1Hj16MGbMGIYNG8avfvWr2mVz585lv/32Y+jQoeyyyy6ceeaZbNy4sdH3Wrt2LTfeeGPeYl++fDn33HNP3vYXRyNGwKRJxY5CWlrBEoGZ3Wpmn5nZwgaWm5ldZ2ZLzOx1MxtbqFhawlZbbcWCBQtqb1GnwaiqqspbDO5OTU1N3vaXbt999+XVV19l3rx53HXXXcyfP59PP/2U448/niuvvJLFixfz1ltvMWXKFNavX9/o/rIlgurq6pzja2oiaMp7ibQnhSwRzASmZFl+KDAkcTsLuKmAsRTFggUL2HPPPRk1ahRHH300n3/+OQCTJ0/m4osvZtKkSVx77bVbbDN9+nS++93vMnnyZAYPHsx1111Xu+zqq69m5MiRjBw5kt///vdAOPgNGzaMc889l7FjxzJnzpzas/KRI0dy8skn8+STTzJx4kSGDBnC3LlzgXAWv/feezNmzBj23ntvFi9eHPlzdenShd13352lS5dyww03MHXqVPbaay8g1KUfd9xx9O3bd4ttFi1axPjx4xk9ejSjRo3i3XffZdq0aSxdupTRo0dz4YUXMnv2bPbff39OOukkdt11V5YvX87IkSNr93HVVVcxffp0AJYsWcI3vvENdtttN8aOHcvSpUuZNm0ac+bMYfTo0VxzzTXMnDmT8847r3b7ww8/nNmzZwPQtWtXLrnkEiZMmMCLL77IXXfdVRvf9773PSUHiZWCdR919+fMbGCWVY4E7vBwrcyXzGwbM9ve3T9uzvv+6EewYEFz9lDf6NGQOO426KuvvmL06NEADBo0iIceeohTTz2VP/zhD0yaNIlLLrmEX/3qV7UH8LVr1/Lss89m3Nfbb7/NM888w/r16xk6dCjnnHMOr7/+Orfddhsvv/wy7s6ECROYNGkSPXv2ZPHixdx2223ceOONLF++nCVLlnD//fczY8YM9thjD+655x6ef/55HnnkEf7rv/6Lhx9+mF122YXnnnuO0tJSnnzySS6++GIeeOCBSN/H6tWreemll/jFL37BrFmzmDp1aqPb3HzzzZx//vmcfPLJbN68merqaq644goWLlzIgsQfbPbs2cydO5eFCxcyaNCgrNVrJ598MtOmTePoo4+moqKCmpoarrjiCq666ioeffRRAGbOnNng9hs2bGDkyJH8+te/5q233uLKK6/khRdeoKysjHPPPZe7776bU089NdL3IdLWFXMcQT9gRcrz8sRr9RKBmZ1FKDUwYMCAFgkuV8mqoaR169axdu1aJiUqXKdOncrxxx9fu/yEE05ocF/f/OY36dSpE506daJPnz58+umnPP/88xx99NF06dIFgGOOOYY5c+ZwxBFHsNNOO7HnnnvWbj9o0CB23XVXAEaMGMGBBx6ImdWeZSfjmzp1Ku+++y5mRmVlZaOfcc6cOYwZM4YOHTowbdo0RowYEfn72WuvvbjssssoLy/nmGOOYciQIRnXGz9+fKN99devX8+HH37I0UcfDYTBXrkqKSnh2GOPBeCpp55i/vz57LHHHkBI6n369Ml5nyJtVTETQab+eJ5pRXefAcwAGDduXMZ1kho7c28tkgf0TDp16lT7uKSkhKqqKkLBKdq+Urfv0KFD7fMOHTrUtkn84he/YP/99+ehhx5i+fLlTJ48udGY991339qz7aQRI0Ywf/58jjzyyKzbnnTSSUyYMIHHHnuMQw45hFtuuYXBgwdn/SylpaVbtHkk+/Fn+y5SNbQ9hORRUlJSu7+pU6dy+eWXR9qvSHtTzF5D5cCOKc/7Ax8VKZa869GjBz179mTOnDkA3HnnnbWlg6bYb7/9ePjhh9m4cSMbNmzgoYceYt99923y/tatW0e/fv2A7FUojTnvvPO4/fbbefnll2tfu+uuu/jkk0+2WG/ZsmUMHjyYH/7whxxxxBG8/vrrdOvWLWujct++ffnss89YvXo1mzZtqk1C3bt3p3///jz88MMAbNq0iY0bN9bb38CBA1mwYAE1NTWsWLGitn0k3YEHHsif//xnPvvsMwDWrFnD++83OGOvSLtTzETwCHBqovfQnsC65rYPtDa33347F154IaNGjWLBggVccsklTd7X2LFjOe200xg/fjwTJkzgzDPPZMyYMU3e30UXXcTPfvYzJk6c2KyG0b59+zJr1iwuuOAChg4dyrBhw5gzZw7du3ffYr377ruPkSNHMnr0aN5++21OPfVUevXqxcSJExk5ciQXXnhhvX2XlZXVNugefvjh7LLLLrXL7rzzTq677jpGjRrF3nvvzSeffMKoUaMoLS1lt91245prrmHixIm11WQXXHABY8dm7pg2fPhwfvOb33DwwQczatQoDjroID7+uF39K4pkZVGL2Tnv2OxeYDKwHfAp8EugDMDdb7YwVPN6Qs+ijcDp7t7oFWfGjRvn6Remeeuttxg2bFhe4xdJF4f/s+QA6gIdFqSIzGy+u4/LtKyQvYa+3chyB75fqPcXEZFoNLJYROop4LhEaYXaTSIoVBWXCMTv/0vj6eKlXSSCzp07s3r16tj9WKVlJK9H0JTxCm1VHmc+kTagXVyYpn///pSXl7Ny5cpihyLtVPIKZXGhRBAv7SIRlJWV6cpRInmkRBAv7aJqSETyS4kgXpQIRKQeJYJ4USIQkXqUCOJFiUBE6lEiiBclAhGpR4kgXpQIRKQeDSiLFyUCEalHJYJ4USIQkXqUCOJFiUBE6lEiiBclAhGpR4kgXpQIRKQeJYJ4USIQkXqUCOJFiUBE6lEiiBclAhGpR4kgXpQIRKQeDSiLFyUCEalHJYJ4USIQkXqUCOJFiUBE6lEiiBclAhGpR4kgXpQIRKQeJYJ4USIQkVodEkcEJYJ4USIQkVpm4V6JIF6UCESklkoE8aREICK1SkrCvQaUxYsSgYjUUokgngqaCMxsipktNrMlZjYtw/IeZvZXM3vNzBaZ2emFjEdEslMiiKeCJQIzKwFuAA4FhgPfNrPhaat9H3jT3XcDJgP/bWYdCxWTiGSXrBpSIoiXQpYIxgNL3H2Zu28GZgFHpq3jQDczM6ArsAbQv6BIkSRLBJWVxY1DWlYhE0E/YEXK8/LEa6muB4YBHwFvAOe7e036jszsLDObZ2bzVq5cWah4RWIv2X20pt6vUNqzQiYCy/Capz0/BFgA7ACMBq43s+71NnKf4e7j3H1c79698x+piAB1CcDTf6nSrhUyEZQDO6Y870848091OvCgB0uA94BdChiTiGSRTAQqEcRLIRPBK8AQMxuUaAA+EXgkbZ0PgAMBzKwvMBRYVsCYRCSLZElAiSBeSgu1Y3evMrPzgCeAEuBWd19kZmcnlt8MXArMNLM3CFVJP3X3VYWKSUSyU9VQPBUsEQC4++PA42mv3Zzy+CPg4ELGICLRqWoonjSyWERqJUsCKhHEixKBiNRSiSCelAhEpJYSQTwpEYhILVUNxZMSgYjUUokgnpQIRKSWxhHEkxKBiABbVgepaihelAhEBNiyFKASQbw0KRGY2Yx8ByIixZVaClAiiJemlgj+J69RiEjRpR78VTUUL01KBO4+P9+BiEhxqWoovhqda8jMnqH+dQRw9wMKEpGIFIUai+MryqRzF6Q87gwciy4nKdLuqEQQX40mggzVQC+Y2bMFikdEikSJIL6iVA1tm/K0A7A78G8Fi0hEikJVQ/EVpWpoPqGNwAhVQu8BZxQyKBFpeSoRxFeUqqFBLRGIiBSXxhHEV1MHlKlqSKSdaco4gooKqKwsTDzScpo6oOyPeY1CRIquKVVDW20Fu+9emHik5TR1QNk38x2IiBRXU6uG3ngj/7FIy4p88Xoz60MYRwCAu39QkIhEpCg0xUR8NVoiMLMjzOxdQm+hZ4HlwN8KHJeItDD1GoqvKFVDlwJ7Au8kehAdCLxQ0KhEpMVpHEF8RUkEle6+GuhgZh3c/RlgdIHjEpEWphJBfEVpI1hrZl2B54C7zewzNNeQSLujRBBfUUoERwIbgR8DfweWAv9eyKBEpOWpaii+oows3pB4WAPcXthwRKRYVCKIL12zWESA5k0xoRJE26ZEICJA7uMIUtffvDn/8UjLUSIQESD3qqGqlC4jX32V/3ik5eScCMzsdjO7ycxGRlh3ipktNrMlZjatgXUmm9kCM1ukC96IFE+uVUOpiWDjxvzHIy2nKSWC64Enge9kW8nMSoAbgEOB4cC3zWx42jrbADcCR7j7COD4JsQjInmQa9WQSgTtR+S5hpLc/RXgFeCBRlYdDyxx92UAZjaL0BX1zZR1TgIeTM5b5O6f5RqPiORHc6qGVCJo26JcqvLrwIXATqnru/sBjWzaD1iR8rwcmJC2zteBMjObDXQDrnX3OzLEcBZwFsCAAQMaC1lEmiDXcQSp1yFQiaBti1IiuB+4GfhfoDqHfVuG19L/vUoJ10A+ENgKeNHMXnL3d7bYyH0GMANg3Lhx6qgmUgAqEcRXlERQ5e43NWHf5cCOKc/7Ax9lWGdVYtDaBjN7DtgNeAcRaVHqNRRfURqL/2pm55rZ9ma2bfIWYbtXgCFmNsjMOgInAo+krfMXYF8zKzWzrQlVR2/l9AlEJC9yrRpSImg/opQIpibuL0x5zYHB2TZy9yozOw94AigBbnX3RWZ2dmL5ze7+lpn9HXidMIXFLe6+MNcPISLN15wSQXUulcbS6kSZa2hQU3fu7o8Dj6e9dnPa898Bv2vqe4hIfjRnHEGV5iNu06L0GioDzgH2S7w0G/gfd69scCMRaXNyHUeQ2mtIiaBti1I1dBNQRhj4BWEg2U3AmYUKSkRaXjIRlJSoaihuoiSCPdx9t5TnT5vZa4UKSESKI1kKKCnJvbFYJYK2LUqvoWoz2zn5xMwGk9t4AhFpA5KlgNJStRHETZQSwYXAM2a2jDBIbCfg9IJGJSItTlVD8RWl19BTZjYEGEpIBG+7+6aCRyYiLUpVQ/HVYCIwswPc/WkzOyZt0c5mhrs/WODYRKQF5Vo1pF5D7Ue2EsEk4GkyX6jeASUCkXZEVUPx1WAicPdfJu7VHiASA6oaiq9Gew2Z2flm1t2CW8zsX2Z2cEsEJyItpzklAiWCti1K99HvuvsXwMFAH0KPoSsKGpWItLjUEoESQbxESQTJ6wocBtzm7q+R+VoDItKGpTYW51o1pDaCti1KIphvZv8gJIInzKwbYaZQEWlHVDUUX1EGlJ0BjAaWufvGxLUI1IAs0s7k2lis7qPtR5QSwV7AYndfa2anAP8JrCtsWCLS0tR9NL6iJIKbgI1mthtwEfA+UO8C8yLStjV1riEzlQjauiiJoMrdHTgSuNbdrwW6FTYsEWlpybP6XBuLO3VSImjrorQRrDeznwGnAPuZWQnh+gQi0o4kSwFlZbmVCDp3ViJo66KUCE4ANgFnuPsnQD90aUmRdqepbQSdO6uNoK2LMvvoJ8DVKc8/QG0EIu1OruMIkr2GVDXU9mWbffR5d9/HzNYTJpmrXQS4u3cveHQi0mJS2wiilgjMoGNHJYK2Ltukc/sk7tUwLBIDTek1VFYW1lfVUNsWpbEYM+sJ7Ji6vrv/q1BBiUjLa8oUE6Wl4aYSQdvWaCIws0uB04Bl1E0t4cABhQtLRFpa8qw+l8bi0tKwvhJB2xalRPAtYGd331zoYESkeJpTIlDVUNsWpfvoQmCbQgciIsWV6ziCykpVDbUXUUoElwOvmtlCwngCANz9iIJFJSItrinjCFQ11D5ESQS3A1cCb6Dpp0XaraZMMZHsNaRE0LZFSQSr3P26gkciIkXVlCkmklVDFRWFjU0KK+qFaS43s73MbGzyFmXnZjbFzBab2RIzm5ZlvT3MrNrMjoscuYjklaqG4itKiWBM4n7PlNca7T6amJzuBuAgoBx4xcwecfc3M6x3JfBE1KBFJP80jiC+osw1tH8T9z0eWOLuywDMbBZhKus309b7AfAAsEcT30dE8iDXKSZSew2p+2jb1mjVkJmdb2bdLbjFzP5lZgdH2Hc/YEXK8/LEa6n77gccDdzcSAxnmdk8M5u3cuXKCG8tIrlqyhQTKhG0D1HaCL7r7l8ABwN9CNcrviLCdpbhtfQC5++Bn7p71vMJd5/h7uPcfVzv3r0jvLWINOSLL+DLL+u/npoIoPHqoWSvIbURtH1R2giSB/TDgNvc/TUzy3SQT1dOmJ8oqT/wUdo644BZid1tBxxmZlXu/nCE/YtIE5xwAmy3Hdx555avp04xASERZPula2Rx+xElEcw3s38Ag4CfmVk3oo0neAUYYmaDgA+BE4GTUldw90HJx2Y2E3hUSUCksD75JPOBuyklAlUNtQ9REsEZwGhgmbtvNLNehOqhrNy9yszOI/QGKgFudfdFZnZ2YnnWdgERKYzKyrqLyqRKTwQ1NXWlg0yqqmDrrVU11B5E6TVUY2afAsPNLNK01SnbPg48nvZaxgTg7qflsm8RaZqqqoYTgRl06FD3vLH9qETQPkSZhvpKwnWL3wSSBUoHnitgXCJSIJWVmQ/c1dUhCSQTQWNVQ+o+2n5EOcM/Chjq7psaXVNEWr1sVUMlJbmXCFQ11PZF6T66DCgrdCAi0jKqqjIfuGtqQhJI9hSKkgg06Vz7EKVEsBFYYGZPseU01D8sWFQiUjDZSgS5VA2p+2j7ESURPJK4iUg70FBjcbKNIJcSQTIRZNqftB1Reg3d3hKBiEjLaKixuKltBGVlqhpq66L0GhpCuErZcKBz8nV3H1zAuESkQLJ1H21Kr6GOHcO21dXZxx1I6xWlsfg24CagCtgfuAO4M+sWItJqNdZ9NNeqoY4dw/PNm/Mbp7ScKIlgK3d/CjB3f9/dp9PItQhEpHWqqQm3fDUWl5UpEbQHURqLK8ysA/BuYsqIDwmzkIpIG5MsCeRzHEEyEWzSSKM2K0qJ4EfA1sAPgd2BU4CphQxKRAojmQDyNY5AVUPtQ9YSQeIykt9y9wuBL4kw2ZyItF7ZSgS5TjGhRNB+NFgiMLPSxAVjdo94/QERaeVSSwTpB/r0NoJsJQL3ukTQqVN4TYmg7cpWIpgLjAVeBf5iZvcDG5IL3f3BAscmInmWWiVUXV035TTUtRFEqRpKvb6xSgRtX5TG4m2B1YSeQk64YpkDSgQibUxqlVByHEBSLr2GkglFvYbah2yJoI+Z/QRYSF0CSGqk9lBEWqPURJDeYJzLOILktioRtA/ZEkEJ0JVoF6EXkTYg9eCf3mCcS/fRTIlA3UfbrmyJ4GN3/3WLRSIiBZetRNCUqiGVCNqHbOMI1FNIpJ3JViIoRtVQnz7wrW/lvp3kV7ZEcGCLRSEiLSK9sThVskSQbEDONqNoctvmdh9duRLuvz/37SS/GkwE7r6mJQMRkcJLPbhnqhoqKak7w892jYF8lAg2bsxtfSmcKFNMiEg7EaVEEKXxNx/dRz/+uO5xY6OYpbCUCERiJEr30ShVPcltU0sQzUkE69bltq3klxKBSIw01n00tUSQ7cBeURHuO3duevfRTz6pe1xentu2kl9KBCIx0ljVUNQz/GT9/tZb56dEoERQXEoEIjHSWGNx1BJBPhJBaong889z21byS4lAJEaylQiSbQRRDuxffRXum5MIUtsF1EZQXEoEIjESZYqJXEsEJSXhlmsiWL8eevcOj7/4IrdtJb+UCERiJMoUE8leQ9kaf1MTAYTkkWsi+OIL6Ns3JBElguJSIhCJkXxVDeUjEaxfD926QffuSgTFVtBEYGZTzGyxmS0xs2kZlp9sZq8nbv80s90KGY9I3BWisRiUCNq6giWCxPWObwAOBYYD3zaz4WmrvQdMcvdRwKXAjELFIyL57z7auXO479gx93EESgStRyFLBOOBJe6+zN03A7OAI1NXcPd/unuy49hLQP8CxiMSe/kaUPbVVyEJJKesbmoi6N5diaA1KGQi6AesSHlennitIWcAf8u0wMzOMrN5ZjZv5cqVeQxRJF6iTDFRVhaeN1YiSFYLAXTpkvskcioRtB6FTASRr2xmZvsTEsFPMy139xnuPs7dx/VO9jcTkZxFKREkp6LOJRF07Qpffhk9DnclgtYkysXrm6oc2DHleX/go/SVzGwUcAtwqLuvLmA8IrHXWPfRkpLwuFOnxruPNicRbNwY3k+JoHUoZIngFWCImQ0ys47AicAjqSuY2QDgQeA77v5OAWMREaJNQw2N9wJqbiJYvz7cKxG0DgUrEbh7lZmdBzwBlAC3uvsiMzs7sfxm4BKgF3CjhevjVbn7uELFJBJ3US5VCS2fCDZsCO+fLJFIyypk1RDu/jjweNprN6c8PhM4s5AxiEidqFVDhU4EyRJAMhFASA7bbBN9H5I/GlksEiNVVbDVVuFxc6qGvvqqbj+QeyJIzja67bZ1iUATzxWPEoFIjFRW1h3AG+o+Co0ngvXrw8E/qWvXkByqq6PFsTrRLaRXr7pEoHaC4lEiEImRqqowTqC0NHuJoFOn7Ingiy+gR4+658mksGFDtDiUCFoXJQKRGKmsDEmgoUSQ2kaQrfvoF1/UHcChLhFErR5KJoLUqiElguJRIhCJkcrKUCIoK2t40jnIXjW0eXO4ZnFzE0HXruF9lAiKT4lAJEaqqhouEURtI0h2/WxOIlizJlQLpe5HiaB4lAhEYiQfJYLkAbu5JYJkIki2NSgRFI8SgUiMJEsEZWXZ2wi6dGn4oJ7s5pmvRNClC5gpERSTEoFIjCRLBI31Gtp+e/j448z7yFQiSB7UV62KFsfq1aGhGMJ7duumRFBMSgQiMZLsPpqpaii1jWCHHcKgr6++qr+PTIng3/4t3H/ySbQ4PvmkbpvkvpQIikeJQCRGGuo+6h4SQ2oigMylgobaCLbeGj79tPEY1q8PVUjJ94AwtcSaNbl9FskfJQKRGGmosXjp0nD2v/PO4XnyIP1RvYnjMycCCGf4UUoEyeSy/fa5byuFoUQgEiMNNRY//XS4P/DAcJ8tEWRqLAbo2zdaiSCZCFJLBDvskPm9pGUoEYjESEONxa+8AtttB1//eniePFv/8MP6+1i6NDQOp84+CiERRDmrTx7wU0sEO+wQEoRnvIahFJoSgUiMNNRYvHYt9OkTunEC9OwZ1vnss/r7eOMN2HXXunWTcq0aSi0RbL99SEyrdY3ColAiEImRhhqL163bsqrHDHr3hpUrt9zeHRYuhJEj6+97xx3DgbyxiedWrAgzoKZOWpetKkoKT4lA2rTKynDWqiqFaBoqEaRPIgchEaSXCD74IPT4yZQIvva1cL90afYY3nwThg3bskShRFBcSgTS5tTUwGOPwZQpYVRq374weDD85S/Fjqz1a6hEkD6tNISqovQSweLF4X7YsPr7TiaCd9/NHkOyainVjjuG+/ffz76tFIYSgbQJNTWwaBFcfTUMHQqHHx4OKOefD9dcE85mjzoKLr+82JG2bqndR7NVDUHmEsE774T7ZKNyqmTX0yVLGn7/1atDG0F6iaJfP+jcufEkIoVR0GsWi0S1ejXMmgWPPBIO+MmDlFl4vH593Wt77w2XXgrHHhsOaADnngunnw4XXxzOZM84ozifo7XLtWoovUTwzjth8FjfvvX33aNH+O7feqvh9//Xv8L9qFFbvt6hQ0gk2ZJIYz74AK6/Hn784y17JEnjlAikqBYuDGf5d98dZrscNgy+8Y3QmJis9y8tDVVAw4fDPvvUnXmm6tgRbrstJJSzzgoHsSOOaNnP0hZkqhqqrg71/ulVQ717hwT84INwzDHhtbffDqWB9B5DSfvsE8YkuGde55lnwnvvtVf9ZUOG1JU4clVVBZMnw3vvwcaNISFIdEoE0uLc4amn4Kqr4IknwkH/jDPge9+D3XZr+n47doQ//zkMijrhhLDv/fbLX9ztQaYSQabrC0Bdwj3xxJCkX301/N3OP7/h/U+ZEhLHm2/CiBH1lz/1FOyxR5hkLt2QIfC3v9UNesvFAw+EJNCjB9x+O/zud3XXZpbGqY1AWsxHH8Fvfxvqhw86CBYsgN/8JnQnvPHG5iWBpK5dQ0PyoEHwzW/CCy9E2666Gh5+GL7zHZg6Ff75z+bH0hqllgiWLoXnnqubMiK9RHD88XDSSWGbRYtg+vQwJ9AllzS8/8MOC/cPPVR/2QcfwNy5deuk22OPcHnM+fNz/ljcdRcMGAD33BNKN889l/s+4kyJQAqqpiacIU6ZEnqG/PSn4WBy662hh8jPf143hXG+bLcdPPlkqCc+8MBQZdRQ91L3cNDadVc4+mj4xz/gr3+FiRNDe0N765aaWiIAmDSp4bmDSkrg1FPD45EjQ/vNeeeFv19D+vWDffcNVX3V1Vsuu+WWcH/SSZm3nTw53D/zTOSPA4RxC08+GToLTJ4MnTqF0qDkwN3b1G333Xd3af1qatwfe8x9zBh3cB8wwP0//9P9nXdaLoZVq9wnTw7vP2mS+0MPua9dW7fsT39y33PPsHzoUPf77nOvrHRfv979zDWL3wYAAA77SURBVDPD69Ont1y8LaFDB/ef/9z9rLPC5wP3qVPD/RNP1F9/9eq69fbf3/3LLxt/j1mzwvq//nXda6+95l5W5n788dm3HTXKfZ99cvpIfsst4f2efTY8P+QQ9513Dv+DUgeY5w0cV4t+YM/1pkTQ+j37rPvEieG/a9Ag9zvucK+qKk4sVVXuN9zgvsMOdQe00tK6x/37hwNJZeWW21VXu59+eljnyiuLE3u+VVfXJbfjjqv7DsB9/Phw0M9k9epwUI16YK2pcT/55Lokc8YZ7l26uPft675yZfZtL7ssbPfee9E/09Ch4YQjGd8f/xj2MXdutH3EhRKBFNznn7vfeaf7AQeE/6rtt3e/8Ub3TZuKHVmwebP7//2f++WXu0+b5n7VVe6zZ9dPAKmqqty//e3wea6+uuViLZRNm8Jnuewy9/32q0sCW23lvnFjft+roiKUqjp1Cu9x0EHRDszLl7uXlLiff36093noobD/e++te+3zz0Pp4yc/aVrs7ZUSgeTVRx+5z5zp/tOfuh91lPuQIXUHlQED3H/7W/cNG4odZX5UVrofe2z4bBdcEBJKW/Xll3UlnF12qfubjRlTuPesrnZfsSK3apozzwwH8qefzr7e55+HEufXvlY/oR9xhHu/fsUribZG2RKBGoslsldeCf3J+/WD004L/f8XLw4NiZddBnPmhC58F15Yf4ritqq0FO69NwxYu+oqGD06dFHdvDn6PlatgtdeC10qKyoKF2tjkt1F02cVTV6DoBA6dID+/Rsed5DJb38bupIedVToDeQZGuwXLgxjFsrLQ2eA9O6mp58eptC++urmxR8X5pm+5VZs3LhxPm/evGKHERs1NfDss+HH+fe/hx4j554b+ukPH557f++27C9/gf/4j7r5+A84AMaMCYPgttkmHPTWrw/X+n3zzTAFxoIF4WCV1KkTHHlk6KJ68MEt+/19/HGY3O2668JnuPZauOOO8Lfs2LHl4ohixYrQ/feNN2CnncJ33atXSKSvvhq6BffoEXqkHXBA/e3d4bjj4NFH4fnnQ9fUuDOz+e4+LuOyQiYCM5sCXAuUALe4+xVpyy2x/DBgI3Cau/8r2z6VCPKjoiIc4J99NvTv7tIljCTt0ycc0NatCwezZ58NZ1a9esEFF8D3v595MFBcVFeHLqZ33QUvvhhKQJmUloY5kUaNgrFjw7iGTZvCNvfeG0ZAb789nHIK/Pu/hwNV586Fjf3SS8MYgHnzQgKrqGjdJTd3mDkzjO/45z/DiGGzkHiPOgq++93sU0msWhW++5UrwyC4s84Kf4eopZOamjDlxWOPhbEN1dVh/Mu3vhXGq7Q1RUkEZlYCvAMcBJQDrwDfdvc3U9Y5DPgBIRFMAK519wnZ9hvHRPDll2FQT5cuoVif+o/sHs5CP/ssnCmtWBHOlHr2DGepHTuG/uDV1eGH8cEH4UD29NPhGrVlZaGqZ+PGcHBK7fs9YADsvnsYWHTEEeH9ZUtr14az6y++CAeObt1Cf/zBgxs+y968ORxcZs4M99XVYd2ddw63Pn3CWIhttw3JIXlpydSpISoqQmKpqAi3yspQ9VNSUrd+MpaKinAG/ac/hQPoAw+06FdUVJ9+Gk5g7r47/Fb69oUJE0LJaJtt6n4bNTV1t88/DwPoFi2qu7ZC//7hvrw8JIFDDw3jJXbaKfxOdtwx/K06dKh/y6VarJCKlQj2Aqa7+yGJ5z8DcPfLU9b5H2C2u9+beL4YmOzuHze036YmgieeCJNRNSbq15HLesl10x9D+JGuWxeqDJL/SCUlYfnmzWH52rV1+yspCeskJf+JczF4cCh2H3poGICTHIpfUwNr1oT37tpVQ/Rbwpo1oW3lxRdDe8uyZSFhr1qVWztEMlHU1ISEkD6YC2D//cNguTgm9HfeCdNbvPhiKBGtXBl+V+51v7nkgbtr11DtOWJEuB1yCAwcGNZ96aUwGPKxx+qutBaFWeYkkZossi1PrnP22XDRRU37DrIlgkLWUPYDVqQ8Lyec9Te2Tj9gi6/YzM4CzgIYMGBAk4Lp3j3zxTQyiZrBc1kvuW76444dw5lJRUX44ace2Dt1CsuTU/Ru2BDO3FOTUElJ2L5XrzA6dvDgcHb6+echwSTPFDt0CGeZffuGs6FMsSfXkZaz7bahzeDII7d83T38rTdtCn+/5N+xsjL8T3TuHG7J/5H0v2d1dfgf+PTTUGLs0SP8b3Tq1HKfrTX5+tfD7Zxzmr4PszBZ3l57hb/PypVhdPz774eSQmXlliWL1Jt7w8tyWaeJh79GFTIRZDpMpp9HR1kHd58BzIBQImhKMMk/YBz07BmKrNJ2mYUz96aevZeUhKSuxF4YZqEKr0+f9tEQXcjuo+XAjinP+wPpF6KLso6IiBRQIRPBK8AQMxtkZh2BE4FH0tZ5BDjVgj2BddnaB0REJP8KVjXk7lVmdh7wBKH76K3uvsjMzk4svxl4nNBjaAmh++jphYpHREQyK+hwFnd/nHCwT33t5pTHDny/kDGIiEh2mmJCRCTmlAhERGJOiUBEJOaUCEREYq7NzT5qZiuB94sdR5rtgFXFDqIBiq3pWnN8iq1p4hzbTu7eO9OCNpcIWiMzm9fQHB7FptiarjXHp9iaRrFlpqohEZGYUyIQEYk5JYL8mFHsALJQbE3XmuNTbE2j2DJQG4GISMypRCAiEnNKBCIiMadE0Agzm2Jmi81siZlNy7C8h5n91cxeM7NFZnZ64vUdzewZM3sr8fr5rSW2lOUlZvaqmT3ammIzs23M7M9m9nbi+8vrJYWaGduPE68tNLN7zSyvl5yPEFtPM3vIzF43s7lmNjLqtsWKrZX8Fhr83hLLi/lbyPY3LehvoZa769bAjTB99lJgMNAReA0YnrbOxcCVice9gTWJdbcHxiZe7wa8k75tsWJLWf4T4B7g0dbyvSWe3w6cmXjcEdimNcRGuIzqe8BWiWV/Ak5r4dh+B/wy8XgX4Kmo2xYxttbwW8gYWyv5LTQYWyF/C6k3lQiyGw8scfdl7r4ZmAWkXV0WB7qZmQFdCQeNKnf/2N3/BeDu64G3CAeSoscGYGb9gW8Ct+QxpmbHZmbdgf2APwK4+2Z3X9saYkssKwW2MrNSYGvye0W9KLENB54CcPe3gYFm1jfitkWJrZX8Fhr63lrDbyFjbC3wW6ilRJBdP2BFyvNy6v8DXw8MIxwQ3gDOd/ea1BXMbCAwBni5FcX2e+AioIb8a05sg4GVwG2JovotZtbEK/fmNzZ3/xC4CvgA+JhwRb1/tHBsrwHHAJjZeGAnwiVeo2xbrNhqFfG3kC22Yv8WGoqt0L+FWkoE2VmG19L72x4CLAB2AEYD1ycyediBWVfgAeBH7v5Fa4jNzA4HPnP3+XmMJy+xEc64xwI3ufsYYAOQz/ru5nxvPQlnc4MSy7qY2SktHNsVQE8zWwD8AHiVUFqJsm2xYgs7KO5vIWNsreS30ND3VujfQq2CXqGsHSgHdkx53p/6VQGnA1d4qMRbYmbvEer55ppZGeEf/253f7AVxTYROMLMDgM6A93N7C53z9dBrTmxfQCUu3vyjPHP5Pefvzmx7QS85+4rAczsQWBv4K6Wii1xAE12SDBCm8V7hGqqxj5XsWKj2L+FLLGdSJF/C438TQv5W9giCN0abugpBZYRzgCTDT0j0ta5CZieeNwX+JAwi6ABdwC/b22xpa0zmfw3kDUrNmAOMDTxeDrwu9YQGzABWET4gRqhIe8HLRzbNtQ1qv8/4I6o2xYxttbwW8gYWyv5LTQYWyF/C1vEUIidtqcbcBihl8NS4OeJ184Gzk483gH4B6EueSFwSuL1fQhFwNcJ1QwLgMNaQ2xp+8j7P39zYyNUx8xLfHcPAz1bUWy/At5OvH4n0KmFY9sLeDcRw4Op302mbVtDbK3kt9Dg99YKfgvZ/qYF/S0kb5piQkQk5tRYLCISc0oEIiIxp0QgIhJzSgQiIjGnRCAiEnNKBNJmmNk1ZvajlOdPmNktKc//28x+ksf3m2lmx+Vrfyn7vTjl8UAzW5jv9xDJhRKBtCX/JIzkxcw6EAZ5jUhZvjfwQhHiytXFja8i0nKUCKQteYFEIiAkgIXA+sR87p0IE8W9amaXmNkrFq4ZMMOCYWY2N7mjxJn464nHu5vZs2Y2P1HK2D79jRtax8xmm9mViXnk3zGzfROvb21mf0rMMX+fmb1sZuPM7ArC7KULzOzuxO5LzOx/LczV/w8z2yrbl2Bm083sgpTnCxOfZ2Bi3vpbEq/dbWbfMLMXzOzdxIRmIvUoEUib4e4fESYKG0BICC8SZrHcCxgHvO5hqt/r3X0Pdx8JbAUc7u5vAR3NbHBidycAf0rMgfMH4Dh33x24Fbgs9X0jrFPq7uOBHwG/TLx2LvC5u48CLgV2T3yGacBX7j7a3U9OrDsEuMHdRwBrgWOb8TV9DbgWGEWYH+kkwsjeC1BJRBqgSeekrUmWCvYGriZM6bs3sI5QdQSwv5ldRJgTaFvC/EB/JVxI5luE2R5PSNyGAiOB/wvzfVFCmGI6VWPrJCdRmw8MTDzeh3BAxt0XJksfDXjP3Rdk2EdTvOfubwCY2SLCRU7czN5o5n6lHVMikLYm2U6wK6FqaAXwH8AXwK0WLh15IzDO3VeY2XTCrJIA9wH3J2YNdXd/18x2BRa5e7ZLAFoj62xK3FdT95vKNP1wQzalPK4mlGKyqWLL0nzq5TJT91WT8rwG/d6lAaoakrbmBeBwYI27V7v7GsLsjXsRqoqSB8VVifnva3v9uPtSwoH2F4SkALAY6G2Ja8GaWZmZpTZAR10n3fOE0gdmNpyQuJIqE9VNWZnZeWZ2XoZFywnz1GNmYwkzW4o0mRKBtDVvEHoLvZT22jp3X+XhUn7/m3jtYeCVtO3vA04hVBORaFM4DrjSzF4jzIy5d+oGUdbJ4EZC8ngd+Clh9sh1iWUzgNdTGosbsguwOsPrDwDbJi5kcg5hZkuRJtPsoyIFYGYlQJm7V5jZzoRr0n49kVSi7uNR4JhcthFpCtUZihTG1sAziSogA87J9YDu7ocXJDKRNCoRiIjEnNoIRERiTolARCTmlAhERGJOiUBEJOaUCEREYu7/Awvbf1Oi/WIqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9b3/8deHkEBA9q1V1KB1Q5YIQaq1Ii5orVXb4tJbW7S1lFZbl+p16a1LrY9qwfvTagulSqvW64IL11ZbrSuKG4uAoKIIKhG9bIZ9C3x+f3xnyJBMJmeSmUwy834+HvOYOet8Mpk5n/NdzveYuyMiIoWrTa4DEBGR3FIiEBEpcEoEIiIFTolARKTAKRGIiBS4trkOIF09e/b0srKyXIchItKqzJ49e5W790q2rNUlgrKyMmbNmpXrMEREWhUz+6i+ZaoaEhEpcEoEIiIFTolARKTAtbo2gmS2b99OZWUlW7ZsyXUokkfat29P3759KS4uznUoIlmVF4mgsrKSTp06UVZWhpnlOhzJA+7O6tWrqayspF+/frkORySr8qJqaMuWLfTo0UNJQDLGzOjRo4dKmVIQ8iIRAEoCknH6TkmhyJtEICIZ8M47MH16rqOQZqZEkCFFRUWUl5fvenz44YcZ2W9ZWRmrVq1Kuc7UqVM55JBDGDlyZNr7/+tf/8ry5cuTLjv33HPp168f5eXlDBkyhFdffXXXsgkTJnDwwQczYMAABg8ezD333BPp/V544QVeeeWVtOOsz7Rp03j77bcztr+C178/jBiR6yikmSkRZEhpaSlz587d9Yg6DEZ1dXWT3/uuu+7ij3/8I88//3za26ZKBADjx49n7ty53HTTTfz4xz8GYNKkSfz73//mjTfeYMGCBUyfPp2oNzhKlQga81k0JhFk4jMXySvu3qoeQ4cO9drefvvtOvOaW8eOHevMe/PNN3348OE+cOBAP/30033NmjXu7j5ixAi/6qqr/Oijj/YJEybsts2qVav8hBNO8PLych87dqzvs88+vnLlSnd3v/fee33YsGE+ePBgHzt2rFdXV/v111/vHTt29AMPPNAvu+wyr66u9ssuu8wrKip84MCBPmnSpF37vvnmm33AgAE+aNAgv+KKK3zq1Km7th08eLBv2rRpt1jGjBnjU6dOdXf3zZs3e2lpqbu777333r548eIGP5PbbrvNDznkEB84cKCfddZZvnTpUu/Tp4/vueeePnjwYJ8+fbqPGTPGL7nkEj/mmGP80ksv9WuvvdbHjx+/ax+HHnqoL1261N3d7777bh84cKAPGjTIzznnHJ8xY4Z369bNy8rKfPDgwb548WIfMWKEz5w5093dV65c6fvuu6+7u//lL3/x0aNH+ymnnOIjR470DRs2+HnnnecVFRVeXl7u06ZNS/o3tITvVrOC8JC8A8zyeo6redF9dDcXXwxz52Z2n+XlcOutKVfZvHkz5eXlAPTr14/HHnuM73//+9x+++2MGDGCa665huuvv55bY/upqqrixRdfrLOf66+/nqOOOoprrrmGJ554gsmTJwPwzjvv8OCDDzJjxgyKi4v56U9/yn333cc111zDc889x4QJE6ioqGDy5Ml06dKFmTNnsnXrVr7yla8watQo3n33XaZNm8brr79Ohw4dWLNmDd27d+eOO+7YtW0qf//73xk4cCDr169n/fr17L///g1+bDfddBNLly6lXbt2VFVV0bVrV8aNG8cee+zBZZddBoTSzHvvvcczzzxDUVER1113XdJ9LVy4kBtvvJEZM2bQs2fPXfGfeuqpnHLKKYwePbrBeF599VXmz59P9+7dufrqqzn22GOZMmUKVVVVHH744Rx//PF07Nixwf2I5Jv8SwQ5Eq8ailu7di1VVVWMiNW3jhkzhjPOOGPX8rPOOivpfqZPn86jjz4KwNe//nW6desGwLPPPsvs2bMZNmwYEBJP796962z/9NNPM3/+fB5++OFdcbz//vs888wznHfeeXTo0AGA7t27R/q7Lr/8cn7zm9/Qq1cv7rrrLtw9cm+aQYMG8d3vfpfTTz+d008/vd71zjjjDIqKilLu67nnnmP06NH07NkzrfgTnXDCCbu2e/rpp3n88ceZMGECELogf/zxxxxyyCFp71ektcu/RNDAmXtLkerMM9mB1t0ZM2YMv/3tb1Pu1925/fbbOfHEE3eb/69//atR3SHHjx9f52y7Y8eOLFmyhP322y/ltk888QTTp0/n8ccf54YbbmDhwoVJ10v8LNq2bcvOnTt3Tcf78UdNQInb174GIPF93J1HHnmEgw46qMF9iuQ7NRZnSZcuXejWrRsvvfQSAPfee++u0kEqRx99NPfddx8A//znP/n8888BOO6443j44YdZsWIFAGvWrOGjj+qOKnviiScyceJEtm/fDsB7773Hxo0bGTVqFFOmTGHTpk27tgfo1KkT69evT+tvu+qqq7jgggtYt24dAOvWrdtVhRW3c+dOli1bxsiRI/nd735HVVUVGzZsaPD9ysrKmDNnDgBz5sxh6dKlu/7+hx56iNWrV6eMv6ysjNmzZwPsKhUlc+KJJ3L77bfvauR+88030/oMRPJJ1hKBmU0xsxVmtqCe5WZmvzezxWY238yGZCuWXLn77ru5/PLLGTRoEHPnzuWaa65pcJtrr72W6dOnM2TIEJ5++mn22WcfAPr3789vfvMbRo0axaBBgzjhhBP49NNP62x//vnn079/f4YMGcKAAQP48Y9/THV1NSeddBKnnnoqFRUVlJeX76oSOffccxk3bhzl5eVs3rw50t/1k5/8hJEjRzJs2DAGDBjAiBEjdlU5xe3YsYNzzjmHgQMHcthhh3HJJZfQtWtXvvGNb/DYY49RXl6+K0km+va3v82aNWsoLy9n4sSJHHjggQAceuih/PKXv2TEiBEMHjyYSy+9FICzzz6b8ePHc9hhh/HBBx9w2WWXMXHiRI488siU3W5/9atfsX37dgYNGsSAAQP41a9+FelvF8lHFj8jyviOzY4GNgD3uPuAJMtPBn4GnAwMB25z9+EN7beiosJr35jmnXfeUd2uZEXBfbfi1W9ZOi5I7pjZbHdP2iskayUCd58OrEmxymmEJOHu/hrQ1cy+mK14REQkuVy2EewFLEuYrozNq8PMxprZLDObtXLlymYJTkSkUOQyESTrApK0POruk929wt0revVKeu/lyFe2ikSl75QUilwmgkpg74TpvkD9Yx2k0L59e1avXq0frmSMx+5H0L59+1yHIpJ1ubyO4HHgQjN7gNBYvNbd63aDiaBv375UVlaiaiPJpPgdygqSe03DseS9rCUCM7sfOAboaWaVwLVAMYC7TwKeJPQYWgxsAs5r7HsVFxfrLlIimbRjB7TNv+tNJbms/afd/TsNLHfggmy9v4g0gRJBQdGVxSJS144duY5AmpESgYjUpURQUJQIRKQuJYKCokQgInUpERQUJQIRqUuJoKAoEYhIXUoEBUWJQETqUiIoKEoEIlKXEkFBUSIQkbqUCAqKEoGI1KVEUFCUCESkrp07cx2BNCMlAhGpSyWCgqJEICJ1KREUFCUCEalLiaCgKBGISF1KBAVFiUBE6lIiKChKBCJSlxJBQVEiEJG6lAgKihKBiNSlRFBQlAhEpC4lgoKiRCAidSkRFBQlAhGpS4mgoCgRiEhdSgQFRYlAROpSIigoSgQiUpcSQUFRIhCRupQICooSgYjUpURQUJQIRKQuJYKCktVEYGYnmdkiM1tsZlcmWd7FzP5uZvPMbKGZnZfNeEQkIiWCgpK1RGBmRcAfgK8B/YHvmFn/WqtdALzt7oOBY4BbzKwkWzGJSERKBAUlmyWCw4HF7r7E3bcBDwCn1VrHgU5mZsAewBqgOosxiUgUSgQFJZuJYC9gWcJ0ZWxeojuAQ4DlwFvARe6uu2aL5JpuXl9QspkILMk8rzV9IjAX2BMoB+4ws851dmQ21sxmmdmslStXZj5SEQnaxA4JKhEUlGwmgkpg74TpvoQz/0TnAY96sBhYChxce0fuPtndK9y9olevXlkLWKTgKREUpGwmgpnAAWbWL9YAfDbweK11PgaOAzCzPsBBwJIsxiQiqSgRFKS22dqxu1eb2YXAU0ARMMXdF5rZuNjyScANwF/N7C1CVdIV7r4qWzGJSAMsVqOrRFBQspYIANz9SeDJWvMmJbxeDozKZgwikgaVCAqSriwWkRpKBAVJiUBEaqhqqCApEYhIDZUICpISgYjUUCIoSI1KBGY2OdOBiEgL4L77sxSExpYI/pTRKESkZYgnAA0xUVAalQjcfXamAxGRFiCeAJQICkqD1xGY2fPUHSMIdz82KxGJSO6oRFCQolxQdlnC6/bAt9FQ0SL5SYmgIDWYCJJUA80wsxezFI+I5FI8AaixuKBEqRrqnjDZBhgKfCFrEYlI7qhEUJCiVA3NJrQRGKFKaCnww2wGJSI5okRQkKJUDfVrjkBEpAVQIihIjb2gTFVDIvlIbQQFqbEXlN2V0ShEpGVQiaAgNfaCsq9nOhARaQHSSQQbN8KiRdmNR5pF5ERgZr3NbJ/4I5tBiUiOpJMITj8dDj5Y1Uh5oMFEYGanmtn7hN5CLwIfAv/Mclwi0twSD+hREsEzz4Tnal1f2tpFKRHcAHwZeC/Wg+g4YEZWoxKR5peYCNI5y9+2LfOxSLOKkgi2u/tqoI2ZtXH354HyLMclIs0t3RJBnBJBqxflgrIqM9sDmA7cZ2Yr0FhDIvkn8eCvRFBQopQITgM2AZcA/wI+AL6RzaBEJAdUIihYUa4s3hh7uRO4O7vhiEjOKBEULN2zWEQCNRYXLCUCEQnURlCwlAhEJFDVUMFKOxGY2d1mNtHMBmQjIBHJESWCgtWYEsEdwDPA9zIci4jkUuLBX20EBSXKdQS7cfeZwEzgkcyHIyI5oxJBwYoy1tCBZvZnM3vazJ6LP6Ls3MxOMrNFZrbYzK6sZ51jzGyumS3UvZBFckiJoGBFKRFMBSYBfwZ2RN2xmRUBfwBOACqBmWb2uLu/nbBOV+CPwEnu/rGZ9U4neBHJICWCghUlEVS7+8RG7PtwYLG7LwEwswcIVym/nbDOfwCPuvvHAO6+ohHvIyKZoO6jBStKY/HfzeynZvZFM+sef0TYbi9gWcJ0ZWxeogOBbmb2gpnNNrPvJ9uRmY01s1lmNmvlypUR3lpE0qYLygpWlBLBmNjz5QnzHNivge0sybza3662wFDC0NalwKtm9pq7v7fbRu6TgckAFRUVuguGSDakUzWUuK4SQasXZayhfo3cdyWwd8J0X2B5knVWxcYz2mhm04HBwHuISPNKp2oo8eCvRNDqRek1VGxmPzezh2OPC82sOMK+ZwIHmFk/MysBzgYer7XO/wJfNbO2ZtYBGA68k+4fISIZkE6JYMuWmtdKBK1elKqhiUAxoXcPhAvJJgLnp9rI3avN7ELgKaAImOLuC81sXGz5JHd/x8z+BcwnjG56p7svaNyfIiJNkk4bwdatNa+VCFq9KIlgmLsPTph+zszmRdm5uz8JPFlr3qRa0+OB8VH2JyJZpBJBwYrSa2iHme0fnzCz/UjjegIRaSXSaSNQiSCvRCkRXA48b2ZLCD2B9gXOy2pUItL80ikRKBHklSi9hp41swOAgwiJ4F1339rAZiLS2qSTCKoTbluuRNDq1ZsIzOxYd3/OzL5Va9H+Zoa7P5rl2ESkOaUz+uj27TWvlQhavVQlghHAcyS/Ub0DSgQi+UQlgoJVbyJw92tjz2oPECkE6SQClQjySpQLyi4ys84W3Glmc8xsVHMEJyLNqLElgq1qMmztonQf/YG7rwNGAb0JPYZuympUItL80uk+mlgiSEwK0ipFSQTxweNOBv7i7vNIPqCciLRm6VxZnJgIEl9LqxQlEcw2s6cJieApM+tEGA5CRPJJY6qGiopUIsgDUS4o+yFQDixx902xexGoAVkk3zSmsbi0VCWCPBClRHAEsMjdq8zsHOC/gLXZDUtEml06bQTxUkBpqUoEeSBKIpgIbDKzwcB/Ah8B92Q1KhFpfo1pI1CJIC9ESQTV7u6E+w3f5u63AZ2yG5aINLvGtBGoRJAXorQRrDezq4BzgKPNrIhwfwIRySfxg39RkdoICkyUEsFZwFbgh+7+GeEG9Lp/gEi+iZcIoiQClQjySpTRRz8D/jth+mPURiCSf9JJBIklgrXqO9LapRp99GV3P8rM1hMGmdu1CHB375z16ESk+SQmgnQai1UiaPVSDTp3VOxZDcMihSCdNoLEqiG1EbR6URqLMbNuwN6J67v7nGwFJSI50NiqIZUIWr0GE4GZ3QCcCyyhZmgJB47NXlgi0uwa01jcrp1KBHkgSongTGB/d9eg4yL5LLFqKEobQXFxeKhE0OpF6T66AOia7UBEJMfSLRHEE4FKBK1elBLBb4E3zWwB4XoCANz91KxFJSLNLzERNHTXse3boW3b8FCJoNWLkgjuBm4G3kLDT4vkr3giaNsWtmxJva5KBHklSiJY5e6/z3okIpJb6Q4xoRJB3oiSCGab2W+Bx9m9akjdR0XySbxE0KZNeo3FO3eGR5soTY7SEkVJBIfFnr+cME/dR0XyTWLVUJTG4niJID5dUpLd+CRroow1NLI5AhGRHEu3aiheIgAlglauwbKcmV1kZp0tuNPM5pjZqCg7N7OTzGyRmS02sytTrDfMzHaY2eh0gheRDEq3+2hiiUANxq1alEq9H7j7OmAU0Jtwv+KbGtoodt+CPwBfA/oD3zGz/vWsdzPwVBpxi0imJbYRNKZEIK1WlERgseeTgb+4+7yEeakcDix29yWxq5IfINzlrLafAY8AKyLsU0SyJbGNoKHG4nj3UZUI8kKURDDbzJ4mJIKnzKwT0a4n2AtYljBdGZu3i5ntBXwTmJRqR2Y21sxmmdmslStXRnhrEUlbY7qPqkSQF6Ikgh8CVwLD3H0TUEKoHmpIslJD7dOMW4Er3H1Hqh25+2R3r3D3il69ekV4axFJW2OGmFCJIC9E6TW008z+D+hvZpGGrY6pJAxdHdcXWF5rnQrgATMD6AmcbGbV7j4tjfcRkUxIdxhqlQjyRpRhqG8m3Lf4bSB+5u7A9AY2nQkcYGb9gE+As4H/SFzB3fslvM9fgX8oCYjkSDp3KKuuhvbtVSLIE1HO8E8HDnL3rQ2umcDdq83sQkJvoCJgirsvNLNxseUp2wVEpJmpjaBgRUkES4BiEoaXiMrdnwSerDUvaQJw93PT3b+IZFDtEoE7WD0dBOPdR1UiyAtREsEmYK6ZPcvuYw39PGtRiUjzS+w+Gp9uKBGoRJAXoiSCx2MPEcln8eqg+OBxqQaS27o13KZSJYK8EKXX0N3NEYiI5NiOWF+Q+Fl+qgbjeCJQiSAvROk1dADhLmX9gfbx+e6+XxbjEpHmFk8E8cHjUjUYb90a1lOJIC9EuaDsL8BEoBoYCdwD3JvNoEQkB+Jn9fGz/FSJYNs2lQjySJREUOruzwLm7h+5+3XoXgQi+SedRKA2grwSpbF4i5m1Ad6PXRfwCWEUUhHJJ1GrhtzrthE0dLN7adGilAguBjoAPweGAucAY7IZlIjkQO0SQX2NxdXVYVm7duEBITFIq5WyRBC7V8CZ7n45sIFog82JSGsUtUQQP+grEeSNeksEZtY2NiroULP6rioRkbwRLxGkkwjat999nrRKqUoEbwBDgDeB/zWzqcDG+EJ3fzTLsYlIc4raWBw/6JeUqESQJ6I0FncHVhN6CjnhPgMOKBGI5JOoF5TFG4ZVNZQ3UiWC3mZ2KbCAmgQQ18AYtSLS6jSmakiJIC+kSgRFwB5Eu9OYiLR26VYNtWsXRiotKlIiaOVSJYJP3f3XzRaJiORW7aqhKIkg/qxE0Kqluo5APYVECkljSgTxZyWCVi1VIjiu2aIQkdyrrg7DThcVhen6GouVCPJOvYnA3dc0ZyAikmM7doSxgxLvR5BMYvdRUCLIA1GGmBCRQlBdHUoDDSWCxO6j8WclglZNiUBEgniJID6QgNoICoYSgYgE1dW7Vw2pjaBgKBGISBCvGmpoaGklgryjRCAiQbxqqEOHML15c/L1lAjyjhKBiATxEkE8EWzalHw9JYK8o0QgIkG8jSBqIoh3H23fXomglVMiEJGgdtVQfYlg06bQoKwSQd5QIhCRIGrV0IYNsMceNd1MlQhaPSUCEQmilgjWr4dOnWqmm5oIPv4Ypk5t/PbSZFlNBGZ2kpktMrPFZnZlkuXfNbP5sccrZjY4m/GISApR2wjiJYK4piaCo46CM8+s//0k67KWCGI3vv8D8DWgP/AdM+tfa7WlwAh3HwTcAEzOVjwi0oCoVUOZLhEsWxaelyxp/D6kSbJZIjgcWOzuS9x9G/AAcFriCu7+irt/Hpt8DeibxXhEJJV41VBxcUgIUUsEpaXhmoP6rkRu6D3jPvgg/e0lI7KZCPYCliVMV8bm1eeHwD+TLTCzsWY2y8xmrVy5MoMhisgu8aohs1AqiFoi6NQpJIGNG9N/z8RSwOLF6W8vGZHNRBD5FpdmNpKQCK5IttzdJ7t7hbtX9OrVK4Mhisgu8aohSJ0IapcI4klh/fr037Oysua1EkHOZDMRVAJ7J0z3BZbXXsnMBgF3Aqe5++osxiMiqcSrhiC9EkHnzjXz0/V//xeeS0p2TwrSrLKZCGYCB5hZPzMrAc4GHk9cwcz2AR4Fvufu72UxFhFpSFNLBOvWpf+e8UQwYACsWpX+9pIRqW5e3yTuXm1mFwJPAUXAFHdfaGbjYssnAdcAPYA/Wrg4pdrdK7IVk4ikEG8jgPoTwc6doS2gdhsBNK5EsGJFeM8DDoDZs9PfXjIia4kAwN2fBJ6sNW9SwuvzgfOzGYOIRBSlaijeIJxYImhq1VCvXtC7t0oEOaQri0UkqF01lKwXUPxgn6xE0NiqoT59oGdPqKqC7dvT34c0mRKBiASJJYJOnZKf4ccP9pmsGoonAoDV6i+SC0oEIhIkthF07pz8DD9+HU/v3jXz4lVDjS0R9O4dqodA1UM5okQgIkFi1VCXLrB2bd11VqwIz4mJoLQ0DEudbonAffeqIVAiyBElAhEJEquGunQJjcW16+yTJQKz+quSUtmwAbZsCftSIsgpJQIRCRJLBPX1BIongviBO66+qqRU4tcQJJYINIRMTigRiEiQ2EbQpUt4rl09tGIF9OhRs15c9+7pN/QmJoIePcJrlQhyQolARILaVUOQPBEkVgvF9e6d/tl8PBH07h2GmOjSRYkgR5QIRCRIVjVUu7onVSKIVxtFFV+/T5/w3LOnEkGOKBGISBClamjZMthzz7rb9uqVfiKIlwjiXUeVCHJGiUBEgtrXEcDuiWD79nB/4f32q7tt796hF1A6t5usrAzbFReH6Z491VicI0oEIhIGk9u0CTp2DNPJSgTLloV2hP33r7t9vLoonQP50qXQr1/NdK9eKhHkiBKBiISzeagpCXTrFtoLPv20Zp343cSSlQji1TvpJIIPP4SyspppVQ3ljBKBiNQ0CscTQUkJHHwwzJtXs078nsL1VQ3B7okjlR07QjVTYiLo0yfc+7iqKq3QpemUCESk5sKxeCIAOOwwmDu3Zvqll8KZ/15Jbj3+pS+F5/ffj/Z+n34a2hwSE0H89UcfRY1aMkSJQESSjypaXh4adD/7LJzB/+tf8LWvhXGFauvZM1wU9u670d4vnjASSxf77huelQianRKBiNStGgI48cTw/Oij8MYb4crhk0+ufx+HHALvvBPt/eIljcGDa+YpEeRMVu9QJiKtRLKqoQED4NBDYeLE0GbQpk1Nckjm4IPhscfCqKLh1rP1e/NN+OIXay4mg1DtVFoaGpGlWalEICLJq4YAfvQjWLAA5syBkSOha9f693HYYaHUEOWM/s03Q9VTIrPQTrB0aVqhS9MpEYhI8qohgB/8AI4/HiZMgIcfTr2P4cPD8+uvp15v1aqQXI44ou6y/v3DMmlWSgQikvxexPHpf/8bfvGL1KUBgEGDoH17ePXV1Os9/3x4PuGEussGD4bFi5PfL1myRm0EkntbtoSqhzfeCH3Vt2wJB5Rhw0IvlfjFSg3Ztg2eeiqcuT7zTDi4HXssTJkShkmW+q1bFz7zkpLG76O4GI47Lnz+t9xSM4BdbY8+GpJKRUXdZYMGhTaGt96CL3+58bFIWpQIpPl9/nnok/7SS/DyyzB7ds2dsLp0CcMcbNgAd9wR6o2HD4dTTgmPQYPqNkQuXw533hkaNT/7LFwVO3JkuDDpiSdg6FCYOjX5gae1mTYNvvlNeOghOOOMzO137dq6pYHGGDMmfOb//Gf4f9W2fDk88ghccEHdexoADBkSnl97TYmgGZm75zqGtFRUVPisWbNyHYakY9UqeOEFmD4dXnwxnO25h7PPww+Hr3wl1BcPHw5f+ELYZufO0MXwH/8Ij5kzw/wePUKjZI8eoQQwZ05N4+TJJ4cDzAkn1Axk9sYb4YD52Wfw+9/D2LEN92hpyf7rv+DGG0Oie+65zO33q18Nn/mMGU3bz/btcOCBoQT2+ut1D/bjxsFdd8GiRcmvUIbQDXWffULprqncw8VrffrUX0IpEGY2292Tnw25e6t6DB061KUV+OQT9zvucB850r1NG3dw79DB/fjj3W+4wf3FF903b46+v08/db/rLvfzz3evqHA/6CD3gw92P/NM91tucV+0qP5tV61y/9rXQgznnOO+bl3T/75cGTcu/B3l5Znb5/bt7qWl7hddlJn9PfBAiHHcOPdt22rm33JLmP/zn6fe/tJL3UtK3D//vGlx7Njh/r3vhfccPNh97dqm7a+VA2Z5PcfVnB/Y030oEbRAO3e6L1ni/sgj7pdf7j5sWPhqQThY//KX7q+84r51a+5i3LEjJCAz9733dn/wwTAvlZ073efMcb/2WvcjjnDv2zc8T5nS8LbZMnp0+Fz79s3cPufODfu8777M7fMXvwj7PPRQ97Fj3b/85TB9xhm7J4dkZs8O6/7+902L4cYbw35OOcW9bVv3c89t2v5aOSUCybx588JBf8QI9y5dag78xcXuRx3l/utfuy9cmOso63rllXBwAvf998PdFBgAAA3MSURBVHe/+mr3J55wf/999+XL3efPd582zf3KK8NyCMnj8MPdv//9cGYJ7sOHu8+a1fzxjxgR3r9du5CoMmHy5LDP99/PzP7iHnjA/cgjw/dj+HD3CRNC6SOKI45w33PPxp/Fv/hiKImefXb4nK64IvyNL73UuP3lASUCyYyNG93vvLPmjL+kJPzAx41z/9Of3F9/3X3TplxH2bDqavf773c/9tiaaqvaj6Ii91Gj3P/8Z/cVK2q23bnT/e673fv0CdtedFHTqzDS0b9/TYzr12dmnz/6kXu3bplLLJnw2mvh8/3mNxsuQdT29tshiXzpSzWJZN0697KykNwz9bm1MkoEUr8dO8KBsT7btrlPn+5+ySXhYBEv7t92W6h7b+02bHB/4QX3e+5xnzgxnMXOnOleVZV6u6oq9wsuCKWFLl1C9dHSpU2LZdOmcFa+enX9B+XevUOpC5r+fu7uCxaEapNRo5q+r0y79daa0tfLLzecqLZsCSckXbqEz2nevN2Xv/BC+H+ddVb0kkkeSZUIstpryMxOAm4DioA73f2mWssttvxkYBNwrrvPSbVP9RrKgBUrQjfExx4LvUQ2bAhdLuMjSJaUhNEm16wJF/ds2xZ6f3zzm3DhhaGHSWvueZNJc+fCr38dPkuAgw4K3R4PPDB8lj16hB40PXpAu3bhLmDr14e7fX34Yc1j6dLwiP8eS0rCVbYnnRSupRg8OHTvbNs2dKGdNw/+9jf47ncbF7d76KJ51lnhZjL33gujR2fgA8mw+++HSy8Nvb722y/0MjvggNDNuG3b8HdUVoZeSC+9FLrBHnEEPPBA6HlU2803w5VXhv/RDTeE3lcF0psoVa+hrCUCMysC3gNOACqBmcB33P3thHVOBn5GSATDgdvcfXiq/RZMIti6NRyMO3RIvV51dThgv/VW6KL5+uvhILLPPqErZteuoVvg5s3hBzN/PrzySpi3//4walS4YGv16tDNc/XqmgN/587hRzd8eLhQqKErSwvZkiXhQqkXXgjXRXz2WbTtvvCFML5OWVlIIv36hesfli8P/8uXXw7fg6Ki8P/4/PMw7MOUKWH7UaPCY489woBt8RvE7NwZnrduDc9VVeFCvffeC0lo/vwwFHRRUYj71FMz+GFk2IYN8D//E65PWLAgJM743wfhQrh+/eCoo+DMM8N3NdWJyv33w8UXhxOiHj3CmEf9+4cupj161CSZ4uLw+WzcGD6/qqqQaOKvq6pCYu/TJ/ze9t47PPr2Dcm/W7eafbUAuUoERwDXufuJsemrANz9twnr/Al4wd3vj00vAo5x93pvc9ToRPDUU3DJJelvV590Prf4uom10PVNV1eHH3v8Evt27cIPvF27sHzHjvAjiD9v3FizfZs24WwRQt/6zz/fPY4uXcLB5qST4NvfhoEDdWafLZs3hwS9Zk1IrmvWhANxx47hoN23bzh4lJam3k9VVUgGr79ec1b8ox+Fg938+eF/nnhQbEhJSfgeDB0Kp50G3/pWTfJoLbZtC59ldXX47nfvnv73ePPmkFiefDKcRC1aVDPMRirFxeGEKP4oLQ3/l2XLwj6TKSkJ8dX3aNMm+rILL4Srrkrvb41JlQiymar2ApYlTFcSzvobWmcvYLdEYGZjgbEA+yQr7kXRuXMYVjeT0vnyxddN/Kcmmy4qCl/s7t3D6/gBZOvW8KVo0ybMj7/u3DlU6XzpS+HCrMSrQ7dvD2cwRUU1yUQH/uZRWhru5JXsbl7p6Nq15qrqRPPmhZMBs5As1q8PJwVr1oQzULOQIEpKwv8/ftDq1q1pw0i0BCUlTf8bSktDVVhiddjWrSFpr1sXkkz80aHD7gf+ZL8h97DtsmXwySfhJCz+2Lw5WXeEmsfOndGXH3hg0/7uemQzESQ74tQ+jY6yDu4+GZgMoUTQqGiOOCL5aIf5rLg4JAnJT/G67fiJgzRNu3aw557hkS6z8Fvr2TNc+d7KZHP00Upg74TpvsDyRqwjIiJZlM1EMBM4wMz6mVkJcDbweK11Hge+b8GXgbWp2gdERCTzslY15O7VZnYh8BSh++gUd19oZuNiyycBTxJ6DC0mdB89L1vxiIhIclnt1+TuTxIO9onzJiW8duCCbMYgIiKp6Q5lIiIFTolARKTAKRGIiBQ4JQIRkQLX6m5VaWYrgY9yHUctPYFVuQ6iHoqt8VpyfIqtcQo5tn3dvVeyBa0uEbREZjarvjE8ck2xNV5Ljk+xNY5iS05VQyIiBU6JQESkwCkRZMbkXAeQgmJrvJYcn2JrHMWWhNoIREQKnEoEIiIFTolARKTAKRE0wMxOMrNFZrbYzK5MsryLmf3dzOaZ2UIzOy82f28ze97M3onNv6ilxJawvMjM3jSzf7Sk2Mysq5k9bGbvxj6/jN5RqImxXRKbt8DM7jez9s0cWzcze8zM5pvZG2Y2IOq2uYqthfwW6v3cYstz+VtI9T/N6m9hF3fXo54HYfjsD4D9gBJgHtC/1jpXAzfHXvcC1sTW/SIwJDa/E/Be7W1zFVvC8kuB/wH+0VI+t9j03cD5sdclQNeWEBvhNqpLgdLYsoeAc5s5tvHAtbHXBwPPRt02h7G1hN9C0thayG+h3tiy+VtIfKhEkNrhwGJ3X+Lu24AHgNNqreNAJzMzYA/CQaPa3T919zkA7r4eeIdwIMl5bABm1hf4OnBnBmNqcmxm1hk4GrgLwN23uXtVS4gttqwtUGpmbYEOZPaOelFi6w88C+Du7wJlZtYn4rY5ia2F/Bbq+9xawm8haWzN8FvYRYkgtb2AZQnTldT9At8BHEI4ILwFXOTuOxNXMLMy4DDg9RYU263AfwI7ybymxLYfsBL4S6yofqeZdWwJsbn7J8AE4GPgU8Id9Z5u5tjmAd8CMLPDgX0Jt3iNsm2uYtslh7+FVLHl+rdQX2zZ/i3sokSQmiWZV7u/7YnAXGBPoBy4I5bJww7M9gAeAS5293UtITYzOwVY4e6zMxhPRmIjnHEPASa6+2HARiCT9d1N+dy6Ec7m+sWWdTSzc5o5tpuAbmY2F/gZ8CahtBJl21zFFnaQ299C0thayG+hvs8t27+FXbJ6h7I8UAnsnTDdl7pVAecBN3moxFtsZksJ9XxvmFkx4Yt/n7s/2oJi+wpwqpmdDLQHOpvZ39w9Uwe1psT2MVDp7vEzxofJ7Je/KbHtCyx195UAZvYocCTwt+aKLXYAjXdIMEKbxVJCNVVDf1euYiPXv4UUsZ1Njn8LDfxPs/lb2C0IPepv6GkLLCGcAcYbeg6ttc5E4LrY6z7AJ4RRBA24B7i1pcVWa51jyHwDWZNiA14CDoq9vg4Y3xJiA4YDCwk/UCM05P2smWPrSk2j+o+Ae6Jum8PYWsJvIWlsLeS3UG9s2fwt7BZDNnaaTw/gZEIvhw+AX8bmjQPGxV7vCTxNqEteAJwTm38UoQg4n1DNMBc4uSXEVmsfGf/yNzU2QnXMrNhnNw3o1oJiux54Nzb/XqBdM8d2BPB+LIZHEz+bZNu2hNhayG+h3s+tBfwWUv1Ps/pbiD80xISISIFTY7GISIFTIhARKXBKBCIiBU6JQESkwCkRiIgUOCUCaTXM7P+Z2cUJ00+Z2Z0J07eY2aUZfL+/mtnoTO0vYb9XJ7wuM7MFmX4PkXQoEUhr8grhSl7MrA3hIq9DE5YfCczIQVzpurrhVUSajxKBtCYziCUCQgJYAKyPjefejjBQ3Jtmdo2ZzbRwz4DJFhxiZm/EdxQ7E58fez3UzF40s9mxUsYXa79xfeuY2QtmdnNsHPn3zOyrsfkdzOyh2BjzD5rZ62ZWYWY3EUYvnWtm98V2X2Rmf7YwVv/TZlaa6kMws+vM7LKE6QWxv6csNm79nbF595nZ8WY2w8zejw1oJlKHEoG0Gu6+nDBQ2D6EhPAqYRTLI4AKYL6HoX7vcPdh7j4AKAVOcfd3gBIz2y+2u7OAh2Jj4NwOjHb3ocAU4MbE942wTlt3Pxy4GLg2Nu+nwOfuPgi4ARga+xuuBDa7e7m7fze27gHAH9z9UKAK+HYTPqYvAbcBgwjjI/0H4crey1BJROqhQeektYmXCo4E/pswpO+RwFpC1RHASDP7T8KYQN0J4wP9nXAjmTMJoz2eFXscBAwA/h3G+6KIMMR0oobWiQ+iNhsoi70+inBAxt0XxEsf9Vjq7nOT7KMxlrr7WwBmtpBwkxM3s7eauF/JY0oE0trE2wkGEqqGlgG/ANYBUyzcOvKPQIW7LzOz6wijSgI8CEyNjRrq7v6+mQ0EFrp7qlsAWgPrbI0976DmN5Vs+OH6bE14vYNQikmlmt1L84m3y0zc186E6Z3o9y71UNWQtDYzgFOANe6+w93XEEZvPIJQVRQ/KK6KjX+/q9ePu39AOND+ipAUABYBvSx2L1gzKzazxAboqOvU9jKh9IGZ9SckrrjtseqmlMzsQjO7MMmiDwnj1GNmQwgjW4o0mhKBtDZvEXoLvVZr3lp3X+XhVn5/js2bBsystf2DwDmEaiJibQqjgZvNbB5hZMwjEzeIsk4SfyQkj/nAFYTRI9fGlk0G5ic0FtfnYGB1kvmPAN1jNzL5CWFkS5FG0+ijIllgZkVAsbtvMbP9CfekPTCWVKLu4x/At9LZRqQxVGcokh0dgOdjVUAG/CTdA7q7n5KVyERqUYlARKTAqY1ARKTAKRGIiBQ4JQIRkQKnRCAiUuCUCERECtz/B8jJHNZJAIgNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xb9bn48c8jecgrXnGmk9jZezphFxLKKKP0toUWCjSMckOh0LIv97bQQn/AbSij7N0WWqAtLWkvBUopK6wMnL3jeGY5dryn9P39oSNXVmRbsnUsK37er5dfkY6Ojh8rOuc53y3GGJRSSg1ejmgHoJRSKro0ESil1CCniUAppQY5TQRKKTXIaSJQSqlBLi7aAYRr6NChJi8vL9phKKVUTFmzZk2lMSYn2Gsxlwjy8vJYvXp1tMNQSqmYIiLFXb2mVUNKKTXIaSJQSqlBThOBUkoNcjHXRqCUsk9bWxtlZWU0NzdHOxTVSy6Xi9zcXOLj40N+jyYCpVSHsrIy0tLSyMvLQ0SiHY4KkzGGQ4cOUVZWRn5+fsjvs61qSESeE5EDIrKxi9dFRB4WkZ0isl5E5tsVi1IqNM3NzWRnZ2sSiFEiQnZ2dtglOjvbCF4Azuzm9a8Ak6yfq4DHbYxFKRUiTQKxrTf/f7YlAmPMB0BVN7ucB/zGeH0KZIjISLviUepoUNnayh8OHIh2GOooE81eQ6OBUr/nZda2I4jIVSKyWkRWHzx4sF+CU2ogunjLFi7YvJnio7gxV0S48cYbO54vX76cO++8s19jOOWUUzoGrp511lkcPny4T8d77733OOecc4K+duGFFzJ79mweeOABli5dyh//+Mc+/a7eiGYiCFZ+CbpKjjHmKWNMgTGmICcn6AhppQaFyrY2APa3tkY5EvskJiby2muvUVlZ2av3t7e3RzSeN954g4yMjIge02ffvn18/PHHrF+/nh/96Ee2/I5QRDMRlAFj/J7nAhVRikWpmJBldQmsjvDFbiCJi4vjqquu4oEHHjjiteLiYk499VRmz57NqaeeSklJCQBLly7lhhtuYPHixdx6660sXbqUq6++msWLFzN+/Hjef/99Lr/8cqZNm8bSpUs7jnf11VdTUFDAjBkzuOOOO4LGk5eXR2VlJU888QRz585l7ty55Ofns3jxYgDefvttjjvuOObPn8/5559PfX09AG+++SZTp07lxBNP5LXXXgt67NNPP50DBw4wd+5cPvzww6C/F2D16tWccsopAFx33XX87Gc/A+Ctt97iS1/6Eh6PJ8RPN7hodh9dAVwrIi8DxwA1xpi9UYxHqQEvK857ylZZJQM7/XDHDgqti1qkzE1N5cFJk3rc75prrmH27NnccsstnbZfe+21XHrppXz3u9/lueee47rrruMvf/kLANu3b+edd97B6XSydOlSqqureffdd1mxYgXnnnsuK1eu5JlnnmHhwoUUFhYyd+5cfv7zn5OVlYXb7ebUU09l/fr1zJ49O2hMy5YtY9myZbS1tbFkyRJuuOEGKisrufvuu3nnnXdISUnhvvvu45e//CW33HIL3/ve93j33XeZOHEi3/rWt4Iec8WKFZxzzjkUFhYC8Oyzz/b42dx7770sXLiQk046ieuuu4433ngDh6Nv9/R2dh/9PfAJMEVEykTkChFZJiLLrF3eAHYDO4Gnge/bFYtSR4tMq0RQdRSXCACGDBnCpZdeysMPP9xp+yeffMJFF10EwCWXXMJHH33U8dr555+P0+nseH7uueciIsyaNYvhw4cza9YsHA4HM2bMYM+ePQC8+uqrzJ8/n3nz5rFp0yY2b97cY2zXX389S5Ys4dxzz+XTTz9l8+bNnHDCCcydO5df//rXFBcXs3XrVvLz85k0aRIiwsUXXxyBT8UrOTmZp59+mtNOO41rr72WCRMm9PmYtpUIjDEX9vC6Aa6x6/crdTTqzxJBKHfudvrhD3/I/Pnzueyyy7rcx7+rZEpKSqfXEhMTAXA4HB2Pfc/b29spKipi+fLlrFq1iszMTJYuXdpj//sXXniB4uJiHnnkEcA7gOu0007j97//faf9CgsL+9wNNy4urqPKJzCuDRs2kJ2dTUVFZGrTda4hpWJImnXHezS3EfhkZWVxwQUXdKouOf7443n55ZcBeOmllzjxxBN7ffza2lpSUlJIT09n//79/P3vf+92/zVr1rB8+XJefPHFjqqYY489lpUrV7Jz504AGhsb2b59O1OnTqWoqIhdu3YBHJEoQpGXl8eaNWsA+NOf/tSxvbi4mPvvv58vvviCv//973z22WdhHzuQJgKlYlB/lAgGghtvvLFT76GHH36Y559/ntmzZ/Pb3/6Whx56qNfHnjNnDvPmzWPGjBlcfvnlnHDCCd3u/8gjj1BVVcXixYuZO3cuV155JTk5ObzwwgsdXUCPPfZYtm7disvl4qmnnuLss8/mxBNPZNy4cWHHd8cdd3D99ddz0kkndVR5GWO44oorWL58OaNGjeLZZ5/lyiuv7PPcUOKtoYkdBQUFRhemUYPVPcXF3F5UxLnZ2ayYNSvix9+yZQvTpk2L+HFV/wr2/ygia4wxBcH21xKBUjHEd9vW0sfugkr500SgVAyJrfK7ihWaCJSKIb6qXDunhYu16mLVWW/+/zQRKBVD7L5Eu1wuDh06pMkgRvnWI3C5XGG9TxemUSqG2H15zs3NpaysDJ3cMXb5VigLhyYCpWKI3YkgPj4+rJWt1NFBq4aUiiFaZaPsoIlAqRhiAv5VKhI0ESgVQzQRKDtoIlAqhmgiUHbQRKBUDPG1EWhbgYokTQRKxSBNAyqSNBEoFUO0akjZQROBUjFEE4GygyYCpWJIRyLQNgIVQZoIlIohHY3FUY5DHV00ESgVQ7RqSNlBE4FSMUQTgbKDJgKlYoi2ESg7aCJQKoZoAlB20ESgVAzRqiFlB00ESsUQTQTKDpoIlIohmgiUHTQRKBVDdNI5ZQdNBErFEC0RKDtoIlAqhmgiUHbQRKBUDNFEoOygiUCpGKJtBMoOtiYCETlTRLaJyE4RuS3I6+ki8lcRWScim0TkMjvjUSrWaYlA2cG2RCAiTuBR4CvAdOBCEZkesNs1wGZjzBzgFOB+EUmwKyalYp0mAmUHO0sEi4CdxpjdxphW4GXgvIB9DJAmIgKkAlVAu40xKRXTNBEoO9iZCEYDpX7Py6xt/h4BpgEVwAbgemOMJ/BAInKViKwWkdUHDx60K16lBjxtI1B2sDMRSJBtgd/eM4BCYBQwF3hERIYc8SZjnjLGFBhjCnJyciIfqVIxQksEyg52JoIyYIzf81y8d/7+LgNeM147gSJgqo0xKRXTNBEoO9iZCFYBk0Qk32oA/jawImCfEuBUABEZDkwBdtsYk1IxTROBskOcXQc2xrSLyLXAW4ATeM4Ys0lEllmvPwHcBbwgIhvwViXdaoyptCsmpWKdthEoO9iWCACMMW8AbwRse8LvcQVwup0xKHU00RKBsoOOLFYqhmgCUHbQRKBUDNKEoCJJE4FSMUSrhpQdNBEoFUO0sVjZQROBUjFESwTKDpoIlIohmgiUHTQRKBVDNBEoO2giUCqGaBuBsoMmAqViiJYIlB00ESgVQzQRKDtoIlAqhmgiUHboVSIQkaciHYhSqmfaRqDs0NsSwZMRjUIpFRItESg79CoRGGPWRDoQpVTPNBEoO/Q4DbWI/Isg3ztjzBJbIlJKdUkTgbJDKOsR3OT32AV8A2i3JxylVHe0jUDZocdEEKQaaKWIvG9TPEqpbmiJQNkhlKqhLL+nDmABMMK2iJRSXdJEoOwQStXQGrzfO8FbJVQEXGFnUEqp4DQBKDuEUjWU3x+BKKV61tFGEOU41NGltwPKtGpIqSjoqBrSxmIVQb0dUPZsRKNQSoVE2wiUHXo7oOzsSAeilOqZJgJlh1AaiwEQkWF4xxEAYIwpsSUipVSXtI1A2aHHEoGIfFVEduDtLfQ+sAf4u81xKaWC0DYCZYdQqobuAo4Ftls9iE4FVtoalVIqKK0aUnYIJRG0GWMOAQ4RcRhj/gXMtTkupVQQmgiUHUJpIzgsIqnAB8BLInIAnWtIqajSRKAiKZQSwXlAI/Aj4E1gF3CunUEppYLTSeeUHUIZWdxgPfQAv7Y3HKVUd7RqSNnB1jWLReRMEdkmIjtF5LYu9jlFRApFZJPOaqpU9zQRKDuEPI4gXCLiBB4FTgPKgFUissIYs9lvnwzgMeBMY0yJNVZBKdUFTQTKDnaWCBYBO40xu40xrcDLeNsb/F0EvOYbnGaMOWBjPErFPG0jUHYIOxGIyK9F5HERmdnDrqOBUr/nZdY2f5OBTBF5T0TWiMilXfzOq0RktYisPnjwYLghK3XU0BKBskNvSgSPAO8Al/SwnwTZFvj9jcO70M3ZwBnAj0Vk8hFvMuYpY0yBMaYgJyenFyErdXTQRKDsEHYbgTFmFbAK+FMPu5YBY/ye5wIVQfaptHomNYjIB8AcYHu4cSk1GGgCUHYIZanKycDNwDj//Y0xS3p46ypgkojkA+XAt/G2Cfh7HXhEROKABOAY4IGQo1dqkNFJ55QdQikR/AF4AngacId6YGNMu4hcC7wFOIHnjDGbRGSZ9foTxpgtIvImsB7vOIVnjDEbw/0jlBosdNI5ZYdQEkG7Mebx3hzcGPMG8EbAticCnv8C+EVvjq/UYKNtBMoOoTQW/1VEvi8iI0Uky/dje2RKqSNoIlB2CKVE8F3r35v9thlgfOTDUUp1R9sIlB1CmWsovz8CUUr1TNsIlB1C6TUUD1wNfMna9B7wpDGmzca4lFJBaNWQskMoVUOPA/F45wQC70Cyx4Er7QpKKRWcJgJlh1ASwUJjzBy/5++KyDq7AlJKdU3bCJQdQuk15BaRCb4nIjKeMMYTKKUiRxOAskMoJYKbgX+JyG688weNAy6zNSqlVFD+icAYg0iwKb2UCk8ovYb+KSKTgCl4E8FWY0yL7ZEppY5gAh5rGlCR0GUiEJElxph3ReTrAS9NEBGMMa/ZHJtSKoB/t1GtJlKR0l2J4GTgXYIvVG8ATQRK9bPAqiG0akhFQJeJwBhzh/WvtgcoNUAEVg0pFQk99hoSketFZIh4PSMia0Xk9P4ITinVmSYCZYdQuo9eboypBU4HhuHtMXSvrVEppYLSNgJlh1ASga8S8izgeWPMOrSzglJRoRd/ZYdQEsEaEXkbbyJ4S0TS8C4io5SKIp14TkVKKAPKrgDmAruNMY3WWgTagKxUFGgbgbJDKCWC44BtxpjDInIx8D9Ajb1hKaWC0USg7BBKIngcaBSROcAtQDHwG1ujUkoFpY3Fyg6hJIJ24/32nQc8ZIx5CEizNyylVDBHDChTKgJCaSOoE5H/Ai4GviQiTrzrEyil+plWDSk7hFIi+BbQAlxhjNkHjAZ+YWtUSqmgNBEoO4Qy++g+4Jd+z0vQNgKlokLbCJQdupt99CNjzIkiUkfn75wAxhgzxPbolFKdaBuBskN3k86daP2rDcNKDRBaNaTsEEpjMSKSCYzx398Ys9auoJRSwWkiUHboMRGIyF3AUmA3/55awgBL7AtLKRWMthEoO4RSIrgAmGCMabU7GKVU97SNQNkhlO6jG4EMuwNRSvVMq4aUHUIpEdwDfCEiG/GOJwDAGPNV26JSSgWliUDZIZRE8GvgPmADOv20UlGlbQTKDqFUDVUaYx42xvzLGPO+7yeUg4vImSKyTUR2isht3ey3UETcIvLNkCNXahDSNgJlh1BKBGtE5B5gBZ2rhrrtPmrNSfQocBpQBqwSkRXGmM1B9rsPeCvM2JUadPTSr+wQSiKYZ/17rN+2ULqPLgJ2GmN2A4jIy3hnMN0csN8PgD8BC0OIRalBTdsIlB1CmWtocS+PPRoo9XteBhzjv4OIjAb+A29S6TIRiMhVwFUAY8eO7WU4SsU+bSNQduixjUBErheRIeL1jIisFZHTQzh2sAXuA7+7DwK3GmPc3R3IGPOUMabAGFOQk5MTwq9W6uikJQJlh1Aaiy83xtQCpwPD8K5XfG8I7yvDOy2FTy5QEbBPAfCyiOwBvgk8JiJfC+HYSg1K2lis7BBKG4Hvzv4s4HljzDoRCXa3H2gVMElE8oFy4NvARf47GGPyO36JyAvA34wxfwklcKUGI4P37s2DlghU5ITaa+htIB/4LxFJI4TxBMaYdhG5Fm9vICfwnDFmk4gss15/og9xKzUoGWNwiOAxRhOBiphQEsEVwFxgtzGmUUSy8VYP9cgY8wbwRsC2oAnAGLM0lGMqNZj5SgS+x0pFQii9hjwish+YLiIhTVutlLKHARwiYIy2EaiICWUa6vvwrlu8GfD17jHABzbGpZQKwvDvRjtNAypSQrnD/xowxRjT0uOeSinbadWQirRQuo/uBuLtDkQp1TNfYzFoIlCRE0qJoBEoFJF/0nmuoetsi0opFVSnxmJtI1AREkoiWGH9KKWirKOxGC0RqMgJpdfQr/sjEKVUzwzeQTm+x0pFQii9hibhXaVsOuDybTfGjLcxLqVUEMYYHA5v5ZAmAhUpoTQWPw88DrQDi4HfAL+1MyilVHDaRqDsEEoiSDLG/BMQY0yxMeZOel6LQPXRu9XVPFpeHu0w1ACjbQTKDqE0FjeLiAPYYc0dVI53FlJlo1PXrQPgmtGjoxyJGkj8B5QpFSmhlAh+CCQD1wELgIuB79oZlFIqOB1HoOzQbYnAWk/4AmPMzUA9IU42p5Syh046p+zQZYlAROKslcMWhLj+gFLKZp3aCLSxWEVIdyWCz4H5wBfA6yLyB6DB96Ix5jWbY1OAx68qQCktESg7hNJYnAUcwttTyNdWZQBNBP2gxeMhyenseUc1KGgbgbJDd4lgmIjcAGzkyM4K+h20UbPb3fFYE4HypyUCZYfuEoETSCV4bzX9DtrocHt7x+MWrQdWfrSNQNmhu0Sw1xjzs36LRHWo8SsRNHt6XB5aDSJaIlB26G4cgbZQRklTQNWQUj7aRqDs0F0iOLXfolCd+FcHaSJQ/rREoOzQZSIwxlT1ZyDq31r9Lv6aCJQ/bSPonVcOHGB7Y2O0wxiwQpliQvUz/4u/NhYrf1oiCF+j2823N29m0Zo10Q5lwNJEMAD5X/y1sVgF0jaC8Kyvrwc6d8JQnWkiGIBatGpIBeGrCtISQXi+sBIBQIMmg6A0EQxAmghUML4Lv7YRhKfQLxGsqauLYiQDlyaCAUgTgQqmIxEEPFfdK2tpYYLLhRN4s0r7wASjiWAA6tR9VO/6lMVXAvAN8NFvRmgOtrUxMSmJE9PT+b9Dh6IdzoCkiWAA0hKB6o7OCh+eg21tDI2P59ghQ9jS2GhLlVpjjLc9aCIYgPzHEWivIeWjVUO9c7C1lZyEBIYnJNBmTKe5vCKhpLmZlA8/5JmKiogetz/ZmghE5EwR2SYiO0XktiCvf0dE1ls/H4vIHDvjiRU6slgFo43F4Wtyu2nweMiJj2d4QgIA+1tbI/o7fAPVnty7N6LH7U+2JQJrmctHga8A04ELRWR6wG5FwMnGmNnAXcBTdsUTS7RqSAWjJYLwVba1AXROBNa2SPEdb3VdHXfv2RPRY/cXO0sEi4CdxpjdxphW4GXgPP8djDEfG2OqraefArk2xhMzWjweEkUQtLFY/VvHOAIdUBayg36JYFh8PBD5EkFpc3PH4x9rIjjCaKDU73mZta0rVwB/D/aCiFwlIqtFZPXBgwcjGOLA1OLxkOhwkOhwaIlAddASQfg6EoHVRgBwIMKJoKylpdNzTwzevNmZCEJe0EZEFuNNBLcGe90Y85QxpsAYU5CTkxPBEAemFmO8iUBEE8EgZYzBHXBB0TaC8B2yEsHQ+Hiy4+NxEPmqodKWFkYkJHBMWhrw7+QTS+xMBGXAGL/nucARzeoiMht4BjjPGKOdfOlcItBeQ4PPnqYmjl27lvM3beq0XUsE4auyLsqZcXE4RciJj6ci4A6+r0pbWpiXmsrt48YB3l5EscbORLAKmCQi+SKSAHwbWOG/g4iMBV4DLjHGbLcxlgHnzHXrGPXxx0Ffa/F4SBDBpVVDg06j203BmjV8XlfHnysrWec3PYK2EYSv2uoqmhnnXYxxXloaH9fWRuz4de3tbGxoYFZKCmMTEwEoiXCi6Q+2JQJjTDtwLfAWsAV41RizSUSWicgya7efANnAYyJSKCKr7YpnoHmrupq9XdRVtvqqhhwObSweZNbW1XGovZ1fjB9PggjP+3VJ1BJB+Kra20l1Ool3eD+1L2dmsrWxkbII3bW/d/gwbcZwRlYWY10uAPZoiaAzY8wbxpjJxpgJxpifW9ueMMY8YT2+0hiTaYyZa/0U2BnPQFQTZHCLNhYPXp9bk6JdMmIEZ2dn88rBgx1tBdpGEL7qtraO0gDAKRkZAHwSoVLBBzU1JIpwQno6mXFx5LlcvBWD8xnpyOIo29XUdMS22vZ20pxObSwehD6qqWFcYiLDExK4cNgw9rW28v7hw4CWCHqjur2dLL9EkG/dtQf29OmtipYWRicmkuhwICJcMnw471RXR7wdwm6aCKKsKEgxssbtJj0uTksEg0yT283bVVWclZ0NwDnZ2aQ6nfzB6jKtbQThq2pvJ9MaPwDetoIkh4PyCF2o97e2dnRLBTg3OxsDfBbBdoj+oIkgSnwffLCFMmra20l3OrXX0CCzsqaGBo+Hc61EkOR0sigtrWMOfS0RhC+wakhEyE1MjFiJYH9bW8dANYAZKSkIsL6hISLH7y+aCKLEd3kPdsdf095Oelyct9eQ1gMPGpusOWsWWP3RAWanprKpoQG3MUe0Ebx64EB/hxhzAquGAEZHMhEElAiSnU4mJSV1LI8ZKzQRRFlgIjDGdCQCrRoaXLY0NJAVF0eO3x3m7JQUGj0edjc1dTQaO63Xntu3LwpRxpbAqiEgYiUCtzFUtrV1SgQAs1JStESgwhN4x9/gduMGMuLitLF4kNnS2MjU5ORO6w0cO2QIAH+urOSANThqmN+Fp12/H11qdrtp9ng6VQ2BNxFUtLbS1sfPrrKtDQNHJILZqansamqKqfWRNRFEWeCFvsb68miJYHAxxrC5sZFpycmdtk9LSWFJRgZPVlSwzxp3MtLvwqNtSF0LHEzmsygtjTZj+Kimpk/H901eNyygxDE7JQUDbIqhUoEmgihwd7PegG9cga+xWNsIBofi5mYq29o6tQ/4nJ6Vxe7mZrZYF5ZR1ghWgCZNBF2qss6lrIAL9WmZmSSKsKKPy1b6EvOIICUCIOLtBFsaGvpciumKJoIo6G69gY5EYJUI9I5vcFhl9QxaGCQRzEpJAeAf1d4Z20f5XXiikQher6xk5uef88vS0gG9RGO13zxD/lLj4jg+PZ2P+1gi8M0p5BtR7JPncpHicPBhH4/vr83jYdHatdywa1fEjulPE0EU+F/cA+/4/ROBzjU0eKyqqyNBpONu0p9/IkhxOEhzOjtei0Yi+P3+/WxqbOTGXbt4vbKy339/qLqqGgJv9c2mhoY+TRm9p7kZJzA6oETgEGHpiBH8Zv/+PicbnzV1ddS73Zycnh6R4wXSRBAFzd2VCKw7rIy4OJIdDpo8npic31yFp7C+npkpKSQ4jjwlcxMTGZWQQJPHw8jExE6NyU1RuCPf3tTEXCth7Yvw3P6R1FXVEMDMlBQaPJ4+zQtUbI0qjgvyf/bT/HwgclNZvGuNLvdNkRFpmgiioLtEcNivRJBi3fkN5OK36jtjDIX19R0X10AiwvdHe9d0WpKR0Wmhj/6uOvzPbdv4or6ek9PTcTKw597vqmoIYJb1WW/oQ4NucXMz4wKqhXyy4+PJjIsLOoVMb3xaW8v05GSGBpQ+IkUTQRR0WyLwayxOtRJBg1YPHdX2trZysK2ty0QAcPOYMfzfrFk8Nnlyp0TQn1VD7R4PT1mzoZ6Qnk5OQsLATgR+N1WBZqWkkCDCyj5U3XSXCAAmJCWxM0KJoKylpWOeJDtoIoiCntoInECK09lRIoil/sgqfL41B+Z0kwgSHA7Oys7Gaa1l7dOfiWCNFee948fzjZwccuLjORjBqqHS5uaIfter29vJsBakCZTsdHJiejpv9nKm0HaPh/KWFsb59eAKNDEpic9qaztWSeuLipaWTr3FIk0TQRQ09VAiGBIXh4iQYtU91vfh5Gj1eFhfX8+KykreOHSIygFcpxsLPMbwdlUVqyM4qVhhCInAX5vfzUN/JYLd1qppAF8fOhSHtdpXpEoEK2tqGPfpp0z87LOOVcX6qipgnqFAZ2RlsaGhodPi86Eqb23FDd2XCFwuat1uvrFxY9jH99fu8XCgra3T+JFI00QQBXV+axAEayz2FWVT+1AicBvDI2VlTPjsM+asXs15Gzdy9oYNDP/4Y67dvl17I/WCMYYLN2/mjPXrWbh2LT8pKorIcQvr68l3uYJWYQTj/33or8bi+0tLAbhu9GgmJiUBkBMf3zHaua9+UlSEwdv4/Kvy8ogcs6q9vdtE8PWhQwF4pRdzNvkambtLBN8dMQKA92tq+rR2xAFrBLOWCI4yvp5BQ+PjgzYWZ1hfXl/VULglAo8xXLVtGz/YuZN8l4tnp0zhn3Pm8M85c1g2ahSPVlRw+rp1mgzC9NL+/bx68CD/NXYsS0eM4K7iYl6zpojui8/r6pgfYmkAOn8femosXllTg7z3Xq8bLRvdbnJWruSxigquHDmShyZN6ui1NCwhISJVQ8YYPq+r4wejR3NedjYPlZVRG2TBpnDtbW3t9i56YnIyC9PSeLUX/4fFViLI6yYRTEpO5uGJEwH6lDB9axtoieAo4/uSD4uPD9pGkG4lgHDaCH67bx8nrl3LMxUV3FtSwnP79vGTceP4YN48Lh85kiWZmSzJzOTRyZP5zdSpfFBTw7U7dkT4Lzt6GWO4r7SUWSkp3J2fz1OTJzMvNZXrduzoU712WXMze5qbOTGM/uE5fheEJ/fu7TRSPdBzVuPuu9ZgtHC8VVVFyocfUmldxG4aM6ZzHPHx1LjdtPbxhmJfayv1bjeTk5L4n3HjqG5v5+fFxX06JkC51b2zO1/JymJNXR2Hw7xQ+xLB2B6OP9WaMmSrNbNsb+wNMrVIpGkiiAJfz6BhCVAhrlAAABsESURBVAlH3JXvampijHWXEWqvofX19VyxbRsra2v53vbt/HdREefn5HBnXl7Q/S8ZMYLbx47lmb17ua+kRJc8DME71dVsbGjghtxcHCLEOxw8OmkS5a2t3NOHi5ZvvpuTwugffsyQIbw7Zw4Aq+vq+HUIs5DeW1LC5jC7Sj5YVtbxePWCBUwJmAfJN0tqZR+rh3ZYpZXJyckUDBnC5SNGcH9pacfI3d5o8Xg42NbWYyJYkpmJB2/1TTiKm5sZHh+Py29wXzC+RHBKYSH7ejnjqa80p1VDR5latxsBsuPiOiWCA62tVLS2Ms+qJgilasgYw7Lt28mMi2P/8cfz0MSJPDhxIi9Om9Zp4FGgn+Xn882cHG7bvZuZq1Yxd9Uqxn3yCQtWr+a7W7awKsZWWLLbw+XlDI+P58Lhwzu2HZeezneGDWN5aSlFvax6+aimhlSnkznW6OFQneyXOEK5EO9ubua0detCPv671dW841eKCFZ15SuZHOhj9dB26255ktX28JO8PAxwd3Fxr29SfAPdRvVwF33skCFkx8Xx46IimsMo2e1ubibfirc7uYmJHefzW70olbV5PDxcXs781NQe/5a+0EQQBTXt7QxxOo+YQsLXe6QjEVi9hrqrevhndTWf1NZyZ14ewxISuC43l+tzc4OOUPXnFOGV6dN5aOJExrpcjHO5OCUjg5z4eFYcOsSitWv51qZN7I2xtVeDqetjfXN5SwtvHDrEFSNHkhjwud43YQJOEW7q5RwwH9bUcNyQIUFHp3bH4Zfka7r4+94/fLjTmgW1YVzoTl23jnbrIvzfY8cGvanwzbrZ155DO5qaSBDpmLNnnMvF1aNG8fTevTxZUdGrY/qWouypRJDocPDklClsaGgI60K9vbGRySEkAhFh9YIFZMTF9WrMQmF9PUXNzdw8Zky3N3Z9FVo3BdUrxhieqKjgvKFDOxXraq2eQRlxcRzyO4l9RWTfVMQ9lQiMMfy0uJjRCQlcPnJk2PE5RLguN5frcnM7ba9rb+f+0lLuKy3lraoqlk+YwOUjR3a6+IC3EesvlZVUtLQwJC6OBWlpnJSe3mMS6g++BvMVhw5xsK2Neamp3JWfz9nWMpDheGHfPjwQ9DMenZjI7ePG8T9FRbxVVcUZWVkhH/dAaysbGho4Pycn7Jj8dTVNwimFhZ2ex4dwIWl2uzt1SZ2enMzd48cH3TcnQolge1MTE5KSOvX3/9WkSaypr+fu4mK+kp3dbe+cYEJNBOBdZzjF4eCtqirOs3oSdae+vZ3y1lYmB1SVdcUhwglDhvCP6mrq2ttJC7F3GHg/Gwi9a3FvRf+MPYrtbGri+zt2MPqTTzpt940VGONycbi9nXorGexuaiLJ4ehY6CLe4SDV6eyyX/WbVVV8VFPD7ePGHXGn2hdpcXHcmZ/P+oIC5qam8r3t25m/ejW/27+fz2treaaigsWFheR+8gnX7NjBPSUl3Lp7N19et45RVvfUDw8fjtocScYYfrZnD8/u28fijAzuzMuj2ePhnA0b+OmePWFVN3iM4dm9e1mSkcGELu4Ab8zNZUpSEsu2bw+r4fiPBw9igK+FcPHpTqjz5fSUCG7auZPcTz4ha+XKjm03BDQQ+/NVDXWVCEL9nIPdXYsId+Xlcai9ndPWrQt7+uUSKxGMCSERJDgcLMnM5JUDBzoG93XHN1o4lBKBz/dHj6a0uZll27eH/B6AbY2NOIDxYfyu3tBEYPng8GEei1D/ZZ/dfidomd9jX88g35e01PrSFjU3k+dydSoCjkhICDqxV7PbzS27dzPe5eLKXpQGQjEpOZl3587ld9OmUed2850tWzhm7Vq+t3075S0t3JGXx5aFC2k/+WSqTziB12fO5MuZmTy7bx9fKixk2MqVfG3DBm7auZN7i4t5du9eXq+s5NOaGva2tNjWSP1YRQU/LS7mO8OG8fL06dyRl8faBQv47vDh3LlnD7fs3h3y7363upqi5uZuP2OX08nTU6awp7mZO8IYW/DqgQNMT05mZi/v9jYtXMj5OTl8XFvLhZs391jH3V0iaPF4uL+srFMJde2CBVzRzd+dGRfnnW/I7/u5o7GR72/fzvhPPyXhgw8Y/+mnnPTFFzxaXh60q6nbGHY1NTEpyN31l7OyeHX6dHY0NXHFtm1hrca2o7GRrLi4I5ap7Mp948cTJ8ItIVTxbfNr3A7VWdnZ3DRmDL8/cIBtYfQg2t7YSJ7LFdEbvWAGXdVQbXs7Py4q4mBbGz8YPZpFQ4Z4p3e1itGr6uq4ICeHlTU13uJfUhLHp6fjsb6wHrzduCYmJVHndrO/tRW3MQxLSCBOhJr2drY2NlLa0sKL+/d3/N4Pamq4yCreFjc3Mz8trVMimJaSwu6mpiPmExmRkMD+IHdc95SUsLGhgb/OnGlrVYxDhAuHD+ebOTkU1tezt7WVPJeLWSkpnRJWRnw8Xx06lK8OHUpdezt/O3SIf1RXs7Kmhrerq4OOgHU5HOS5XORbP8MTEhgaH09OfDzjXC4mJyWREeKJ7PN2VRXX79jBudnZ/MavwdzldPLc1KmkOJ0sLy2l1ePhgYkTj6juCnRfaSk58fH8Rw937SdlZHDVyJE8UFbG6VlZnN5DFVFFSwsf1NR02bMrFNNTUnh00iQ+rKnh5QMHePnAAd6YNYslmZkcZ40C9tddmvinX/34V7OzuXj48G7nPgLvd2Oo36Cy/zt0iPM3baLF42FJZiZnZ2dzoLWVlVZX5XtLSrhv/Hgu8mtwL21upsWYLu+uz8nO5o5x4/hpcTGpTieP+o1j6M6OpqaOxudQTEtJ4XsjR3JPSQklzc1HrDHg74PDh0lxOI5YTa4nN4wZw8Pl5dxdXMxvp00L6T1bGxvDSji9NWgSgTGGnxcX8+M9ezq2/T7IiMIX9u3jhQgvCp7kcPCdLVvIiovjhPR0djU3c9nIkR1ftpLmZmra29nY0MA5AXXYw+Pj2RxwB7Guvp7lpaWcn5PDOX2sVghVvMPBQmv93J6kxcVx4fDhnXrYNLrdHGxro7KtjX2trRQ3N1PU3ExRUxNFzc18WlvbMUmYjwPv5GbnDR3KednZTOzhhPi0poavb9zIjJQUXpw27YiLvEOERyZNIsHh4MGyMqrb23l2yhTiu0ik71k9Z345YUKP3QQB7p8wgU9qa7lg0yY+nj+f6d30BPplaSkCXDhsWI/H7U5OQgLbFi0i/aOPADhrw4Yu9z3U1oYx5ogL6YrKSs7zmwbhP4YO5fwQ48pzudjS2MjfKiv52saNzE1N5bWZMztdSJvdbv5cWcnNu3bxnS1b2NLYyJ15eThF+NxakCfYOgzgrSK6Mz+fRo+HX5SW4nI4WD5hQo8JfEdTU6eeVaG43Erk523cyOfz5wf9Xhhj+HtVFadmZoZ9AzYsIYH/HDWKB8vKSHc6eWTy5G73//DwYdY1NPCtPn5HQjFoEsFz+/Z1SgJduWT4cMa7XKysreVweztjExM5OSODJo+H1XV1uI1hREICc1NTaXC72djQQJwIeS4XU5OTqXW78RhDbmIiL+3fz77WVnITE3l23z5+VV7eMTZgTkoKoxISELwlgnerq3HDEY2NIxISOuYih383gqbHxfGgNWoxFiQ7nYxzOrtt9GvzeDjU1sb+tjaKm5tZXVfH65WV3LRrFzft2sX05GRvUhg6lIK0tI7GxTaPh2f27uWGXbsYnZDAm7NnM6SLBjkR4ZcTJni7DO7Zw97WVn43bVqnQVrgrSr54c6djE5I4OpRo0L6G1Pj4vjbrFkcs3YtZ61fz6fz5zMiSB11SXMzj5SXc+mIEUGrRMLV1d86LD6eS4YPJ8/l4sX9+/msro46t/uI/f3rxU8YMoRvhtF4vSQzk3tKSviopoYFqam8O3fuEcd3OZ1cOHw4FwwbxhVbt3J3cTHlLS08M2VKx911TyOr7x0/nka3mwfKykh0OLg7Pz/oZHLgnXajrKWlYyqMUI1PSuLFadP4+qZN/Kq8PGj7yE7rxuXmbtpOunNPfj5FTU08WlHB90eP7vZm4RelpYxISOD6gM4cdhg0ieA7w4bhNoZZKSneEZNOJ2dnZzM9JYUnysu5escOHpw4MaIf+mnWRb3Z7eZPlZU0WHdG8SIcl55OvMPByIQESltaaDcGJ7Ao4K57eEIC1e3ttHg8JDocPFJezud1dbwwdaqtA0yiId7hYERiIiMSE5mTmspXhw7lZ9aJs+LQIf5SWcn/lpRwT0kJKQ4Hc1JTcTkcrG9ooLKtjcUZGbw6fXqPc7aLCP+Tl0duYiLLtm9n7urV/GLCBL41bBhOEera27l061bWWVVvoZQGfMa6XPx15kxOLizkuC++4M8zZjDXb/nJZreb727dioE+VQsFen3mTL62cSMG7x39/44f36kElREXx2dbt1La0sKMgAu1f2Pv23PmkBzG33tGVhb3lJQA8McZM7pMSuDtsvz81KnkuVz8tLiY/a2t/OvwYU6yzoXuOET41aRJ1Lvd3FtSwv7WVp6dMiVoNdGqujoMhD02A7wN9+dkZ3P77t2clJ5+RCnYN1tpOL3D/LmcTpZPmMDrhw4xY9UqPp43j+OCjCqvaW/nraoqrhk9Oqz/j94aNInA5XRylXVnF/jBLxs9mjOzsrqtF+zr7/7a0KH8o6qKspYWzsjKItuq+x7rclHa0sKhtjamJCcf0Sjku6vZ3tjI0Ph4/ruoiK9kZXGpX7XL0S4/KYnrrfERVW1tvFlVxSe1tayvr6fR4+HsrCy+mZPD2dnZYfW1XjpyJHNSU7li2za+s2ULP9q5k8nJyRTW19PgdvPQxIm9qnorGDKEf82dy9c2bmTBmjV8b+RIrhg5klq3m/9XXMx7hw/z26lTw+4S2Z2vDh3K4RNP5NPa2qDtE74bjJU1Nczwu0CWt7R0mnQt3IvOl9LTeX3mTOanppIbwt8jItyRl0eS08ltu3cjwK1jx4b0u8RKJLmJify8pISqtjZemTHjiHPmH9XVOOjdal4iwrNTpnDM2rWctm4dqxcs6Eioxhj+XFnJ5KSkPvXimZiczO+mTeOiLVu4dOtW/jB9eqebBYAnKypoNYYL+ti1OFSDJhH0JM/m7lkTk5J4weo1cYk1KyF4u7etq6+nzRgWBVm4fKZ10m5oaODtqipaPR4enjjR1sElA1lWfDwXDR/eqcGxL+alpbF6wQL+UlnJHw8eZF9rK+fn5LBs1KgjSmfhWDRkCOsLCvhZcTGPV1TwpDXnT4rDwXNTpnCx33cgUobExXXZSD0pKYlh8fG8f/gwV40aRXlLC68eOMCPi4po8HiYl5rKQ72oahQRvhpmshQRbh07lm8PG8b+1tawPmcR4a78fLLj47lh1y5OX7eOV6ZP76iCM8aworKSRUOGhN3RwGdYQgLvzplDwZo1nFJYyAtTp/LlrCx+u38//zp8mP/tYlxFOC4cPpxmj4crt21j8bp1vDRtGqdnZhLncPDXykp+UlTE14YO5Zg+fAfDoYmgn8zzqwP1722Q73LxB2v2w2us5Qj9TUlOJk6Ex8rLWVlby21jx/bYaKrC4xDh6zk5fD3Cd19DExJ4eNIkbh4zhrX19cSJcHJ6OqlhDCiKFBHhvKFDeX7fPm6rr+eaHTv40Brp+oPRo1k2alS39dV2GGeNaA+XiPCjMWMYnpDAldu2MX/NGp6bMoUzsrJ4u7qa9Q0NPN1DQ2xP8pOSeGv2bJZu3cpp69eT73Kxp7mZk9PT+VGEqo8vGzmSUzIyOHXdOs7esIF0azGqitZWCtLSeHry5H674RM7JxwTkTOBhwAn8Iwx5t6A18V6/SygEVhqjDmy35ufgoICs3r1apsitk91W1vHQJ2NCxd2FM/fqaritPXrASgsKAg6gvC8DRtYcegQIxMS2L5oUVQuJCr2HWprY+THH3NDbi6/O3CgY/yK5+STY7aEub6+nm9s2sTOpiamJydT1NzMqIQENi5cGFbbTlcOt7XxWEUFa+rqyHO5+EleXsjrRoSqxePh9cpKnqioIEGEiUlJ3DN+fFgjkEMhImuMMQXBXrPtiiIiTuBR4DSgDFglIiuMMZv9dvsKMMn6OQZ43Pr3qJMZH8+NubkMiYtjut8d/ckZGSxKS2NSUhKzurgje2zyZKaUlXHhsGGaBFSvZcfHU5CWxoc1NTRag8+uGjkyZpMAeLudbly4kKcrKlhx6BAL09K4Kz8/IkkAvONjbh83LiLH6kqiw8EFw4ZxQT90E+2KbSUCETkOuNMYc4b1/L8AjDH3+O3zJPCeMeb31vNtwCnGmL1dHTdWSwRKDQQ379rFcmu1sUcmTQpaHamOTlEpEQCjgVK/52UcebcfbJ/RQKdEICJXAVcBjA2xh4FS6kj/OXIklW1tZMTF8e0o3oGqgcXORBCsvBlY/AhlH4wxTwFPgbdE0PfQlBqcJiYn8/zUqdEOQw0wds5kVAb4D7/LBQInFw9lH6WUUjayMxGsAiaJSL6IJADfBlYE7LMCuFS8jgVqumsfUEopFXm2VQ0ZY9pF5FrgLbzdR58zxmwSkWXW608Ab+DtOroTb/fRy+yKRymlVHC29kU0xryB92Lvv+0Jv8cGuMbOGJRSSnVPF6ZRSqlBThOBUkoNcpoIlFJqkNNEoJRSg5ytk87ZQUQOAsXRjiPAUKAy2kF0QWPrvYEcn8bWO4M5tnHGmKBT7MZcIhiIRGR1V3N4RJvG1nsDOT6NrXc0tuC0akgppQY5TQRKKTXIaSKIjKeiHUA3NLbeG8jxaWy9o7EFoW0ESik1yGmJQCmlBjlNBEopNchpIuiBiJwpIttEZKeI3Bbk9XQR+auIrBORTSJymbV9jIj8S0S2WNuvHyix+b3uFJEvRORvAyk2EckQkT+KyFbr8ztuAMX2I2vbRhH5vYi4+jm2TBH5s4isF5HPRWRmqO+NVmwD5Fzo8nOzXo/mudDd/6mt50IHY4z+dPGDd/rsXcB4IAFYB0wP2Od24D7rcQ5QZe07EphvbU8Dtge+N1qx+b1+A/A74G8D5XOznv8auNJ6nABkDITY8C6jWgQkWa+9Cizt59h+AdxhPZ4K/DPU90YxtoFwLgSNbYCcC13GZue54P+jJYLuLQJ2GmN2G2NagZeB8wL2MUCaiAiQivei0W6M2WuMWQtgjKkDtuC9kEQ9NgARyQXOBp6JYEx9jk1EhgBfAp4FMMa0GmMOD4TYrNfigCQRiQOSieyKeqHENh34J4AxZiuQJyLDQ3xvVGIbIOdCV5/bQDgXgsbWD+dCB00E3RsNlPo9L+PIL/AjwDS8F4QNwPXGGI//DiKSB8wDPhtAsT0I3AJ4iLy+xDYeOAg8bxXVnxGRlIEQmzGmHFgOlAB78a6o93Y/x7YO+DqAiCwCxuFd4jWU90Yrtg5RPBe6iy3a50JXsdl9LnTQRNA9CbItsL/tGUAhMAqYCzxiZXLvAURSgT8BPzTG1A6E2ETkHOCAMWZNBOOJSGx477jnA48bY+YBDUAk67v78rll4r2by7deSxGRi/s5tnuBTBEpBH4AfIG3tBLKe6MVm/cA0T0XgsY2QM6Frj43u8+FDrauUHYUKAPG+D3P5ciqgMuAe423Em+niBThref7XETi8X7xXzLGvDaAYjsB+KqInAW4gCEi8qIxJlIXtb7EVgKUGWN8d4x/JLJf/r7ENg4oMsYcBBCR14DjgRf7KzbrAurrkCB42yyK8FZT9fR3RSs2on0udBPbt4nyudDD/6md50KnIPSn64aeOGA33jtAX0PPjIB9HgfutB4PB8rxziIowG+ABwdabAH7nELkG8j6FBvwITDFenwn8IuBEBtwDLAJ7wkqeBvyftDPsWXw70b17wG/CfW9UYxtIJwLQWMbIOdCl7HZeS50isGOgx5NP8BZeHs57AL+29q2DFhmPR4FvI23LnkjcLG1/US8RcD1eKsZCoGzBkJsAceI+Je/r7HhrY5ZbX12fwEyB1BsPwW2Wtt/CyT2c2zHATusGF7z/2yCvXcgxDZAzoUuP7cBcC50939q67ng+9EpJpRSapDTxmKllBrkNBEopdQgp4lAKaUGOU0ESik1yGkiUEqpQU4TgYoZIvKAiPzQ7/lbIvKM3/P7ReSGCP6+F0Tkm5E6nt9xb/d7nCciGyP9O5QKhyYCFUs+xjuSFxFx4B3kNcPv9eOBlVGIK1y397yLUv1HE4GKJSuxEgHeBLARqLPmc0/EO1HcFyLyExFZJd41A54Sr2ki8rnvQNad+Hrr8QIReV9E1liljJGBv7irfUTkPRG5z5pHfruInGRtTxaRV6055l8Rkc9EpEBE7sU7e2mhiLxkHd4pIk+Ld67+t0UkqbsPQUTuFJGb/J5vtP6ePGve+mesbS+JyJdFZKWI7LAmNFPqCJoIVMwwxlTgnShsLN6E8AneWSyPAwqA9cY71e8jxpiFxpiZQBJwjjFmC5AgIuOtw30LeNWaA+dXwDeNMQuA54Cf+//eEPaJM8YsAn4I3GFt+z5QbYyZDdwFLLD+htuAJmPMXGPMd6x9JwGPGmNmAIeBb/ThY5oIPATMxjs/0kV4R/behJZEVBd00jkVa3ylguOBX+Kd0vd4oAZv1RHAYhG5Be+cQFl45wf6K96FZC7AO9vjt6yfKcBM4B/e+b5w4p1i2l9P+/gmUVsD5FmPT8R7QcYYs9FX+uhCkTGmMMgxeqPIGLMBQEQ24V3kxIjIhj4eVx3FNBGoWONrJ5iFt2qoFLgRqAWeE+/SkY8BBcaYUhG5E++skgCvAH+wZg01xpgdIjIL2GSM6W4JQOlhnxbrXzf/PqeCTT/clRa/x268pZjutNO5NO+/XKb/sTx+zz3o+a66oFVDKtasBM4BqowxbmNMFd7ZG4/DW1XkuyhWWvPfd/T6Mcbswnuh/THepACwDcgRay1YEYkXEf8G6FD3CfQR3tIHIjIdb+LyabOqm7olIteKyLVBXtqDd556RGQ+3pktleo1TQQq1mzA21vo04BtNcaYSuNdyu9pa9tfgFUB738FuBhvNRFWm8I3gftEZB3emTGP939DKPsE8Rje5LEeuBXv7JE11mtPAev9Gou7MhU4FGT7n4AsayGTq/HObKlUr+nso0rZQEScQLwxpllEJuBdk3aylVRCPcbfgK+H8x6lekPrDJWyRzLwL6sKSICrw72gG2POsSUypQJoiUAppQY5bSNQSqlBThOBUkoNcpoIlFJqkNNEoJRSg5wmAqWUGuT+P+C5Bgzs5T1qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "186\n",
      "499\n",
      "499\n",
      "Quality factor of the 1D PC resonator (from the left sideband)= 1386.3636363635583\n"
     ]
    }
   ],
   "source": [
    "from subprocess import call\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "simfileName_n ='Output-norm-ph-cav-flux-pixel-avg.txt'\n",
    "fileName_n = 'flux-norm-ph-cav-pixel-avg.dat'\n",
    "\n",
    "simfileName_d ='Output-defect-ph-cav-flux-pixel-avg.txt'\n",
    "fileName_d = 'flux-defect-ph-cav-pixel-avg.dat'\n",
    "\n",
    "ex_str_n =f'grep flux1: {simfileName_n}  | cut -d , -f2,3- | sed 1d > {fileName_n}'\n",
    "ex_str_d =f'grep flux1: {simfileName_d}  | cut -d , -f2,3- | sed 1d > {fileName_d}'\n",
    "\n",
    "call(ex_str_n, shell=\"True\")\n",
    "call(ex_str_d, shell=\"True\")\n",
    "\n",
    "f_n = np.genfromtxt(fileName_n, delimiter=\",\")\n",
    "f_d = np.genfromtxt(fileName_d, delimiter=\",\")\n",
    "\n",
    "#print(f)\n",
    "\n",
    "#g = np.genfromtxt(\"flux-nh.dat\", delimiter=\",\")\n",
    "\n",
    "nzd_f = f_d[:,1]/f_n[:,1] # Normalized flux\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(1/f_n[:,0],f_n[:,1], 'b-', label='For normal PC structure')\n",
    "plt.xlabel(\"Wavelength, um\")\n",
    "plt.ylabel(\"Transmission, a.u.\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(1/f_d[:,0],f_d[:,1], 'r-', label='For defect PC structure')\n",
    "plt.xlabel(\"Wavelength, um\")\n",
    "plt.ylabel(\"Transmission, a.u.\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(1/f_d[:,0],nzd_f/max(nzd_f), 'c-', label='Normalized flux')\n",
    "plt.xlabel(\"Wavelength, um\")\n",
    "plt.ylabel(\"Transmission, a.u.\")\n",
    "#plt.xlim(left=0.88, right=0.94)\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# Quality factor of the 1D PC resonator\n",
    "\n",
    "\n",
    "# Search of the flux's global maximum and its cell's number\n",
    "for i in range(len(nzd_f)):\n",
    "    if nzd_f[i] == max(nzd_f):\n",
    "        i = i\n",
    "        break\n",
    "    else:\n",
    "         i = i + 1\n",
    "print(i) #Cell's number of the flux's global maximum\n",
    "\n",
    "# Search of the half of the maximum flux and its cell's number (left sideband)\n",
    "for j in range(len(nzd_f)):\n",
    "    if nzd_f[j] >= 1*max(nzd_f)/2 and nzd_f[j] <= max(nzd_f)/2:\n",
    "        j = j\n",
    "        break\n",
    "    else:\n",
    "         j = j + 1\n",
    "print(j) # Cell's number of the half flux's global maximum (left sideband)\n",
    "\n",
    "# Search of the half of the maximum flux and its cell's number (right sideband)\n",
    "for k in range(len(nzd_f)):\n",
    "    if nzd_f[k] >= 1*max(nzd_f)/2 and nzd_f[k] <= max(nzd_f)/2 and k > i:\n",
    "        k = k\n",
    "        break\n",
    "    else:\n",
    "         k = k + 1\n",
    "print(k) # Cell's number of the half flux's global maximum (right sideband)\n",
    "\n",
    "w_r = f_d[i,0] # Resonance frequency\n",
    "w_s_l = f_d[i-1,0] # HWHM's left frequency\n",
    "#w_s_r = f_d[k,0] # HWHM's right frequency\n",
    "delta_w_ls = 2*(w_r - w_s_l) # FWHM (full width at half maximum) (from lef sideband) \n",
    "#delta_w_rs = -2*(w_r - w_s_r) # FWHM (full width at half maximum) (from right sideband)\n",
    "Q_ls = w_r/delta_w_ls # Quality factor of the 1D PC resonator (from left sideband)\n",
    "#Q_rs = w_r/delta_w_rs # Quality factor of the 1D PC resonator (from right sideband)\n",
    "\n",
    "print('Quality factor of the 1D PC resonator (from the left sideband)=', Q_ls) \n",
    "#print('Quality factor of the 1D PC resonator (from the right sideband)=', Q_rs)\n",
    "\n",
    "#print(w_r, w_s_l, delta_w_ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8725552532540507"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/1.1460592280782966"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9120598801138553"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/1.096419239354292"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9511747426384649"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/1.0513315326541353"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8725552532453101"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/1.146059228089777"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}