TEMIP_linear_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SimPEG import *\n",
    "from SimPEG import EM\n",
    "from scipy.constants import mu_0\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "from simpegEMIP.TDEM import geteref, Problem3DIP_Linear, SurveyLinear\n",
    "from simpegEMIP.TDEM import Survey, Rx\n",
    "from simpegEMIP.TDEM import Problem3DEM_e, Problem3D_e\n",
    "import matplotlib.pyplot as plt\n",
    "from pymatsolver import PardisoSolver\n",
    "\n",
    "eta, tau, c = 0.1, 0.01, 0.5\n",
    "cs, ncx, ncz, npad = 10., 25, 20, 18\n",
    "hx = [(cs,ncx), (cs,npad,1.3)]\n",
    "hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]\n",
    "mesh = Mesh.CylMesh([hx,1,hz], '00C')    \n",
    "sigmaInf = np.ones(mesh.nC) * 0.001\n",
    "actinds = mesh.gridCC[:,2]<0.\n",
    "layerind = (np.logical_and(mesh.gridCC[:,2]<0, mesh.gridCC[:,2]>-50.)) & (mesh.gridCC[:,0]<100.)\n",
    "sigmaInf[~actinds] = 1e-8\n",
    "sigmaInf[layerind] = 0.1\n",
    "eta = np.zeros(mesh.nC)\n",
    "eta[layerind] = 0.5\n",
    "tau = np.ones(mesh.nC) * 1.\n",
    "c = np.ones(mesh.nC) * 0.5\n",
    "\n",
    "def get_em_data(sigma, eta=None, tau=None, c=None, data_type='em'):\n",
    "    rxloc = np.array([[0., 0., 30.]])\n",
    "    srcloc = np.array([[0., 0., 30.]])\n",
    "    dt = 1.47e-3\n",
    "    tpeak = 2.73e-3\n",
    "    t0 = tpeak + dt\n",
    "    rx_vtem = Rx.Point_dbdt(rxloc, np.logspace(np.log10(2e-5), np.log10(0.009), 51)+t0, orientation='z')\n",
    "    src_vtem = EM.TDEM.Src.CircularLoop([rx_vtem], waveform=EM.TDEM.Src.VTEMWaveform(offTime=t0, peakTime=tpeak, a=3.), loc=srcloc)\n",
    "    survey = Survey([src_vtem])\n",
    "    if data_type == 'em':\n",
    "        prb = Problem3DEM_e(mesh, sigma=sigma)\n",
    "    elif data_type == 'emip':\n",
    "        prb = Problem3D_e(mesh, sigmaInf=sigma, eta=eta, tau=tau, c=c)\n",
    "    prb.timeSteps = [(tpeak/10, 10), ((t0-tpeak)/10, 10), (1e-06, 5), (2.5e-06, 5), (5e-06, 5), (1e-05, 10), (2e-05, 10), (4e-05, 10), (8e-05, 10), (1.6e-04, 10), (3.2e-04, 20)]\n",
    "    prb.Solver = PardisoSolver\n",
    "    prb.pair(survey)\n",
    "    e = prb.fields(sigmaInf)\n",
    "    data = survey.dpred(sigmaInf, f=e)\n",
    "    # waveform\n",
    "    cur = []\n",
    "    for t in prb.times:\n",
    "        cur.append(src_vtem.waveform.eval(t))\n",
    "    cur = np.hstack(cur)\n",
    "    return e, data, cur\n",
    "\n",
    "def get_ip_data(sigma, eref, eta, tau, c):\n",
    "    rxloc = np.array([[0., 0., 30.]])\n",
    "    srcloc = np.array([[0., 0., 30.]])    \n",
    "    rx_ip = Rx.Point_dbdt(rxloc, np.logspace(np.log10(2e-5), np.log10(0.009), 51), 'z')\n",
    "    src_ip = EM.TDEM.Src.CircularLoop([rx_ip], loc=srcloc)\n",
    "    dt = 1.47e-3\n",
    "    tpeak = 2.73e-3\n",
    "    t0 = tpeak + dt\n",
    "    survey_ip = SurveyLinear([src_ip])\n",
    "    t1, t2, t3 = dt, t0-0.001365, t0\n",
    "    prb_ip = Problem3DIP_Linear(\n",
    "        mesh, \n",
    "        sigmaInf=sigmaInf, \n",
    "        eta=eta, \n",
    "        tau=tau, \n",
    "        c=c, \n",
    "        actinds = actinds,\n",
    "        tlags = [0., t1, t2, t3]\n",
    "    )\n",
    "    prb_ip.Solver = PardisoSolver\n",
    "    prb_ip.pair(survey_ip)\n",
    "    prb_ip.set_eref(eref)\n",
    "    ip_approx = survey_ip.dpred([])\n",
    "    return ip_approx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-200.0, 200.0)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAG2CAYAAAC5/Bg9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAq2UlEQVR4nO3df3DTdYL/8VeENpTaxkKhabRAz8XhpIhSWCju2halwPJDThdh2fHoHHJ4ULRbOLHoauGE+gt0Bk5071hAlh2cuRH3B9xtiwLaAdyC4ALuCRyFlqWxiiUBxKSUz/eP3X6+hhYovEtDwvMx85khn8878f32k0yf80maOizLsgQAAICrdlO4JwAAABDpCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMRE1SlpaUaNGiQEhIS1L17d40fP16ff/55yBjLslRSUiKPx6O4uDjl5ORo//79IWMCgYBmzZql5ORkxcfHa9y4cTp27Fh7LgUAAESZiAmqrVu3aubMmdqxY4fKy8t17tw55eXl6cyZM/aYl19+WUuWLNGyZctUWVkpt9ut4cOH69SpU/aYwsJCrV+/XuvWrVNFRYVOnz6tMWPGqLGxMRzLAgAAUcARqX8c+csvv1T37t21detW3XfffbIsSx6PR4WFhZo7d66kv16NSklJ0UsvvaTp06fL5/OpW7duWrNmjSZOnChJOn78uNLS0rRx40aNGDEinEsCAAARqmO4J3C1fD6fJKlLly6SpKqqKnm9XuXl5dljnE6nsrOztW3bNk2fPl27du1SQ0NDyBiPx6OMjAxt27btokEVCAQUCATs2+fPn9fXX3+trl27yuFwXIvlAQCANmZZlk6dOiWPx6ObbmrbN+kiMqgsy1JRUZF+8IMfKCMjQ5Lk9XolSSkpKSFjU1JSdPToUXtMbGyskpKSmo1pun9LSktLNX/+/LZcAgAACJOamhrddtttbfqYERlUBQUF+tOf/qSKiopmxy68YmRZ1mWvIl1uTHFxsYqKiuzbPp9PPXr00A/0I3VUzBXOHgAAhMM5NahCG5WQkNDmjx1xQTVr1iz99re/1YcffhhSl263W9Jfr0Klpqba++vq6uyrVm63W8FgUPX19SFXqerq6jR06NCL/jedTqecTmez/R0Vo44OggoAgIjwt0+NX4uP60TMb/lZlqWCggK9++67+uCDD5Senh5yPD09XW63W+Xl5fa+YDCorVu32rGUmZmpmJiYkDG1tbXat2/fJYMKAADgUiLmCtXMmTP161//Wr/5zW+UkJBgf+bJ5XIpLi5ODodDhYWFWrRokXr37q3evXtr0aJF6ty5syZPnmyPnTp1qmbPnq2uXbuqS5cumjNnjvr166cHHnggnMsDAAARLGKCavny5ZKknJyckP0rV65Ufn6+JOmpp57S2bNnNWPGDNXX12vw4MEqKysLea/0tddeU8eOHfXII4/o7Nmzuv/++7Vq1Sp16NChvZYCAACiTMR+D1U4+f1+uVwu5ehBPkMFAECEOGc1aIt+I5/Pp8TExDZ97Ij5DBUAAMD1iqACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABiKqKD68MMPNXbsWHk8HjkcDr333nshx/Pz8+VwOEK2IUOGhIwJBAKaNWuWkpOTFR8fr3HjxunYsWPtuAoAABBtIiqozpw5o/79+2vZsmUXHTNy5EjV1tba28aNG0OOFxYWav369Vq3bp0qKip0+vRpjRkzRo2Njdd6+gAAIEp1DPcErsSoUaM0atSoS45xOp1yu90tHvP5fFqxYoXWrFmjBx54QJL0q1/9Smlpadq0aZNGjBjR5nMGAADRL6KuULXGli1b1L17d91xxx2aNm2a6urq7GO7du1SQ0OD8vLy7H0ej0cZGRnatm3bRR8zEAjI7/eHbAAAAE2iKqhGjRqltWvX6oMPPtDixYtVWVmpYcOGKRAISJK8Xq9iY2OVlJQUcr+UlBR5vd6LPm5paalcLpe9paWlXdN1AACAyBJRb/ldzsSJE+1/Z2RkaODAgerZs6c2bNighx566KL3syxLDofjoseLi4tVVFRk3/b7/UQVAACwRdUVqgulpqaqZ8+eOnjwoCTJ7XYrGAyqvr4+ZFxdXZ1SUlIu+jhOp1OJiYkhGwAAQJOoDqoTJ06opqZGqampkqTMzEzFxMSovLzcHlNbW6t9+/Zp6NCh4ZomAACIcBH1lt/p06d16NAh+3ZVVZX27NmjLl26qEuXLiopKdHDDz+s1NRUHTlyRPPmzVNycrL+4R/+QZLkcrk0depUzZ49W127dlWXLl00Z84c9evXz/6tPwAAgCsVUUG1c+dO5ebm2rebPtc0ZcoULV++XHv37tXbb7+tkydPKjU1Vbm5uXrnnXeUkJBg3+e1115Tx44d9cgjj+js2bO6//77tWrVKnXo0KHd1wMAAKKDw7IsK9yTiDR+v18ul0s5elAdHTHhng4AAGiFc1aDtug38vl8bf556Kj+DBUAAEB7IKgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYiKqg+/PBDjR07Vh6PRw6HQ++9917IccuyVFJSIo/Ho7i4OOXk5Gj//v0hYwKBgGbNmqXk5GTFx8dr3LhxOnbsWDuuAgAARJuICqozZ86of//+WrZsWYvHX375ZS1ZskTLli1TZWWl3G63hg8frlOnTtljCgsLtX79eq1bt04VFRU6ffq0xowZo8bGxvZaBgAAiDIOy7KscE/iajgcDq1fv17jx4+X9NerUx6PR4WFhZo7d66kv16NSklJ0UsvvaTp06fL5/OpW7duWrNmjSZOnChJOn78uNLS0rRx40aNGDGiVf9tv98vl8ulHD2ojo6Ya7I+AADQts5ZDdqi38jn8ykxMbFNHzuirlBdSlVVlbxer/Ly8ux9TqdT2dnZ2rZtmyRp165damhoCBnj8XiUkZFhjwEAALhSHcM9gbbi9XolSSkpKSH7U1JSdPToUXtMbGyskpKSmo1pun9LAoGAAoGAfdvv97fVtAEAQBSImitUTRwOR8hty7Ka7bvQ5caUlpbK5XLZW1paWpvMFQAARIeoCSq32y1Jza401dXV2Vet3G63gsGg6uvrLzqmJcXFxfL5fPZWU1PTxrMHAACRLGqCKj09XW63W+Xl5fa+YDCorVu3aujQoZKkzMxMxcTEhIypra3Vvn377DEtcTqdSkxMDNkAAACaRNRnqE6fPq1Dhw7Zt6uqqrRnzx516dJFPXr0UGFhoRYtWqTevXurd+/eWrRokTp37qzJkydLklwul6ZOnarZs2era9eu6tKli+bMmaN+/frpgQceCNeyAABAhIuooNq5c6dyc3Pt20VFRZKkKVOmaNWqVXrqqad09uxZzZgxQ/X19Ro8eLDKysqUkJBg3+e1115Tx44d9cgjj+js2bO6//77tWrVKnXo0KHd1wMAAKJDxH4PVTjxPVQAAEQevocKAADgOkZQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGoiqoSkpK5HA4Qja3220ftyxLJSUl8ng8iouLU05Ojvbv3x/GGQMAgGgQVUElSX379lVtba297d271z728ssva8mSJVq2bJkqKyvldrs1fPhwnTp1KowzBgAAkS7qgqpjx45yu9321q1bN0l/vTr1+uuv65lnntFDDz2kjIwMrV69Wt98841+/etfh3nWAAAgkkVdUB08eFAej0fp6emaNGmSDh8+LEmqqqqS1+tVXl6ePdbpdCo7O1vbtm275GMGAgH5/f6QDQAAoElUBdXgwYP19ttv6w9/+IP+4z/+Q16vV0OHDtWJEyfk9XolSSkpKSH3SUlJsY9dTGlpqVwul72lpaVdszUAAIDIE1VBNWrUKD388MPq16+fHnjgAW3YsEGStHr1anuMw+EIuY9lWc32Xai4uFg+n8/eampq2n7yAAAgYkVVUF0oPj5e/fr108GDB+3f9rvwalRdXV2zq1YXcjqdSkxMDNkAAACaRHVQBQIB/fnPf1ZqaqrS09PldrtVXl5uHw8Gg9q6dauGDh0axlkCAIBI1zHcE2hLc+bM0dixY9WjRw/V1dXphRdekN/v15QpU+RwOFRYWKhFixapd+/e6t27txYtWqTOnTtr8uTJ4Z46AACIYFEVVMeOHdNPfvITffXVV+rWrZuGDBmiHTt2qGfPnpKkp556SmfPntWMGTNUX1+vwYMHq6ysTAkJCWGeOQAAiGQOy7KscE8i0vj9frlcLuXoQXV0xIR7OgAAoBXOWQ3aot/I5/O1+eeho/ozVAAAAO2BoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGLphg+qNN95Qenq6OnXqpMzMTH300UfhnhIAAIhQN2RQvfPOOyosLNQzzzyj3bt364c//KFGjRql6urqcE8NAABEoKsKqmHDhmn+/PnN9tfX12vYsGHGk7rWlixZoqlTp+qxxx7T3//93+v1119XWlqali9fHu6pAQCACHRVQbVlyxYtW7ZM48eP15kzZ+z9wWBQW7dubbPJXQvBYFC7du1SXl5eyP68vDxt27atxfsEAgH5/f6QDQAAoEnHq73jpk2bNH36dA0ZMkS/+93v1KtXrzac1rXz1VdfqbGxUSkpKSH7U1JS5PV6W7xPaWlpi1fk1h/Yq8SEDtdkngDCY+yBkeGeAoBrxHEmKP3o2jz2VX+GKjU1VVu3btVdd92lQYMGacuWLW04rWvP4XCE3LYsq9m+JsXFxfL5fPZWU1PTHlMEAAAR4qquUDWFh9Pp1Nq1a/XCCy9o5MiRmjt3bptO7lpITk5Whw4dml2Nqqura3bVqonT6ZTT6WyP6QEAgAh0VVeoLMsKuf3ss89q7dq1Wrx4cZtM6lqKjY1VZmamysvLQ/aXl5dr6NChYZoVAACIZFd1haqqqkrdunUL2ffwww+rT58+2rlzZ5tM7FoqKirSo48+qoEDByorK0u/+MUvVF1drccffzzcUwMAABHoqoKqZ8+eLe7v27ev+vbtazSh9jBx4kSdOHFCCxYsUG1trTIyMrRx48aLrgsAAOBSrvq3/CLdjBkzNGPGjHBPAwAARIEb8pvSAQAA2hJBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwFFVB1atXLzkcjpDt6aefDhlTXV2tsWPHKj4+XsnJyXriiScUDAbDNGMAABANOoZ7Am1twYIFmjZtmn375ptvtv/d2Nio0aNHq1u3bqqoqNCJEyc0ZcoUWZalpUuXhmO6AAAgCkRdUCUkJMjtdrd4rKysTJ999plqamrk8XgkSYsXL1Z+fr4WLlyoxMTE9pwqAACIElH1lp8kvfTSS+ratavuvvtuLVy4MOTtvO3btysjI8OOKUkaMWKEAoGAdu3addHHDAQC8vv9IRsAAECTqLpC9eSTT2rAgAFKSkrSH//4RxUXF6uqqkr/+Z//KUnyer1KSUkJuU9SUpJiY2Pl9Xov+rilpaWaP3/+NZ07AACIXNd9UJWUlFw2ZiorKzVw4ED97Gc/s/fdddddSkpK0o9//GP7qpUkORyOZve3LKvF/U2Ki4tVVFRk3/b7/UpLS9Mjh4YrJj72SpcEAACizHUfVAUFBZo0adIlx/Tq1avF/UOGDJEkHTp0SF27dpXb7dbHH38cMqa+vl4NDQ3Nrlx9l9PplNPpvLKJAwCAG8Z1H1TJyclKTk6+qvvu3r1bkpSamipJysrK0sKFC1VbW2vvKysrk9PpVGZmZttMGAAA3HCu+6Bqre3bt2vHjh3Kzc2Vy+VSZWWlfvazn2ncuHHq0aOHJCkvL0933nmnHn30Ub3yyiv6+uuvNWfOHE2bNo3f8AMAAFctaoLK6XTqnXfe0fz58xUIBNSzZ09NmzZNTz31lD2mQ4cO2rBhg2bMmKF7771XcXFxmjx5sl599dUwzhwAAES6qAmqAQMGaMeOHZcd16NHD/3+979vhxkBAIAbRdR9DxUAAEB7I6gAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYiJqgWLlyooUOHqnPnzrrllltaHFNdXa2xY8cqPj5eycnJeuKJJxQMBkPG7N27V9nZ2YqLi9Ott96qBQsWyLKsdlgBAACIVh3DPYHWCgaDmjBhgrKysrRixYpmxxsbGzV69Gh169ZNFRUVOnHihKZMmSLLsrR06VJJkt/v1/Dhw5Wbm6vKykodOHBA+fn5io+P1+zZs9t7SQAAIEpETFDNnz9fkrRq1aoWj5eVlemzzz5TTU2NPB6PJGnx4sXKz8/XwoULlZiYqLVr1+rbb7/VqlWr5HQ6lZGRoQMHDmjJkiUqKiqSw+For+UAAIAoEjFv+V3O9u3blZGRYceUJI0YMUKBQEC7du2yx2RnZ8vpdIaMOX78uI4cOXLRxw4EAvL7/SEbAABAk6gJKq/Xq5SUlJB9SUlJio2NldfrveiYpttNY1pSWloql8tlb2lpaW08ewAAEMnCGlQlJSVyOByX3Hbu3Nnqx2vpLTvLskL2Xzim6QPpl3q7r7i4WD6fz95qampaPScAABD9wvoZqoKCAk2aNOmSY3r16tWqx3K73fr4449D9tXX16uhocG+CuV2u5tdiaqrq5OkZleuvsvpdIa8TQgAAPBdYQ2q5ORkJScnt8ljZWVlaeHChaqtrVVqaqqkv35Q3el0KjMz0x4zb948BYNBxcbG2mM8Hk+rww0AAOBCEfMZqurqau3Zs0fV1dVqbGzUnj17tGfPHp0+fVqSlJeXpzvvvFOPPvqodu/erffff19z5szRtGnTlJiYKEmaPHmynE6n8vPztW/fPq1fv16LFi3iN/wAAICRiPnahOeee06rV6+2b99zzz2SpM2bNysnJ0cdOnTQhg0bNGPGDN17772Ki4vT5MmT9eqrr9r3cblcKi8v18yZMzVw4EAlJSWpqKhIRUVF7b4eAAAQPRwWXxN+xfx+v1wul4ZvnK6Y+NhwTwcAALRCw5mgyn/0lnw+n/3uVVuJmLf8AAAArlcEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAUMQE1cKFCzV06FB17txZt9xyS4tjHA5Hs+3NN98MGbN3715lZ2crLi5Ot956qxYsWCDLstphBQAAIFp1DPcEWisYDGrChAnKysrSihUrLjpu5cqVGjlypH3b5XLZ//b7/Ro+fLhyc3NVWVmpAwcOKD8/X/Hx8Zo9e/Y1nT8AAIheERNU8+fPlyStWrXqkuNuueUWud3uFo+tXbtW3377rVatWiWn06mMjAwdOHBAS5YsUVFRkRwOR1tPGwAA3AAi5i2/1iooKFBycrIGDRqkN998U+fPn7ePbd++XdnZ2XI6nfa+ESNG6Pjx4zpy5EgYZgsAAKJBxFyhao1/+7d/0/3336+4uDi9//77mj17tr766is9++yzkiSv16tevXqF3CclJcU+lp6e3uLjBgIBBQIB+7bf7782CwAAABEprFeoSkpKWvwg+Xe3nTt3tvrxnn32WWVlZenuu+/W7NmztWDBAr3yyishYy58W6/pA+mXeruvtLRULpfL3tLS0q5glQAAINqF9QpVQUGBJk2adMkxF15RuhJDhgyR3+/XF198oZSUFLndbnm93pAxdXV1kv7/laqWFBcXq6ioyL7t9/uJKgAAYAtrUCUnJys5OfmaPf7u3bvVqVMn+2sWsrKyNG/ePAWDQcXGxkqSysrK5PF4LhluTqcz5HNXAAAA3xUxn6Gqrq7W119/rerqajU2NmrPnj2SpO9973u6+eab9bvf/U5er1dZWVmKi4vT5s2b9cwzz+if//mf7RiaPHmy5s+fr/z8fM2bN08HDx7UokWL9Nxzz/EbfgAA4KpFTFA999xzWr16tX37nnvukSRt3rxZOTk5iomJ0RtvvKGioiKdP39ef/d3f6cFCxZo5syZ9n1cLpfKy8s1c+ZMDRw4UElJSSoqKgp5Ow8AAOBKOSy+JvyK+f1+uVwuDd84XTHxseGeDgAAaIWGM0GV/+gt+Xw+JSYmtuljR933UAEAALQ3ggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMERQAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMAQQQUAAGCIoAIAADBEUAEAABgiqAAAAAwRVAAAAIYIKgAAAEMEFQAAgCGCCgAAwBBBBQAAYIigAgAAMBQRQXXkyBFNnTpV6enpiouL0+23367nn39ewWAwZFx1dbXGjh2r+Ph4JScn64knnmg2Zu/evcrOzlZcXJxuvfVWLViwQJZltedyAABAlOkY7gm0xv/+7//q/Pnzeuutt/S9731P+/bt07Rp03TmzBm9+uqrkqTGxkaNHj1a3bp1U0VFhU6cOKEpU6bIsiwtXbpUkuT3+zV8+HDl5uaqsrJSBw4cUH5+vuLj4zV79uxwLhEAAEQwhxWhl2deeeUVLV++XIcPH5Yk/fd//7fGjBmjmpoaeTweSdK6deuUn5+vuro6JSYmavny5SouLtYXX3whp9MpSXrxxRe1dOlSHTt2TA6Ho1X/bb/fL5fLpeEbpysmPvbaLBAAALSphjNBlf/oLfl8PiUmJrbpY0fEFaqW+Hw+denSxb69fft2ZWRk2DElSSNGjFAgENCuXbuUm5ur7du3Kzs7246ppjHFxcU6cuSI0tPTW/xvBQIBBQKBkP+2JJ37JtjieAAAcP1p+rl9La4lRWRQ/d///Z+WLl2qxYsX2/u8Xq9SUlJCxiUlJSk2NlZer9ce06tXr5AxTffxer0XDarS0lLNnz+/2f7NP15psgwAABAGJ06ckMvlatPHDGtQlZSUtBgq31VZWamBAwfat48fP66RI0dqwoQJeuyxx0LGtvSWnWVZIfsvHNNUqZd6u6+4uFhFRUX27ZMnT6pnz56qrq5u8xNyPfP7/UpLS1NNTU2bXyq9nrFu1n0jYN2s+0bg8/nUo0ePkHe42kpYg6qgoECTJk265JjvXlE6fvy4cnNzlZWVpV/84hch49xutz7++OOQffX19WpoaLCvQrndbvtqVZO6ujpJanZ167ucTmfI24RNXC7XDfVEbJKYmMi6byCs+8bCum8sN+q6b7qp7b/kIKxBlZycrOTk5FaN/ctf/qLc3FxlZmZq5cqVzf5nZGVlaeHChaqtrVVqaqokqaysTE6nU5mZmfaYefPmKRgMKjY21h7j8XiavRUIAADQWhHxPVTHjx9XTk6O0tLS9Oqrr+rLL7+U1+sNudqUl5enO++8U48++qh2796t999/X3PmzNG0adPs+p48ebKcTqfy8/O1b98+rV+/XosWLVJRUVGrf8MPAADgQhHxofSysjIdOnRIhw4d0m233RZyrOkzUB06dNCGDRs0Y8YM3XvvvYqLi9PkyZPt76mS/voWXXl5uWbOnKmBAwcqKSlJRUVFIZ+Pag2n06nnn3++xbcBoxnrZt03AtbNum8ErLvt1x2x30MFAABwvYiIt/wAAACuZwQVAACAIYIKAADAEEEFAABgiKC6Qm+88YbS09PVqVMnZWZm6qOPPgr3lNpUaWmpBg0apISEBHXv3l3jx4/X559/HjImPz9fDocjZBsyZEiYZtw2SkpKmq3J7Xbbxy3LUklJiTwej+Li4pSTk6P9+/eHccZto1evXs3W7XA4NHPmTEnRc64//PBDjR07Vh6PRw6HQ++9917I8dac30AgoFmzZik5OVnx8fEaN26cjh071o6ruHKXWndDQ4Pmzp2rfv36KT4+Xh6PR//4j/+o48ePhzxGTk5Os+fA5b6QOdwud75b87yOtvMtqcXXusPh0CuvvGKPibTz3ZqfWe31+iaorsA777yjwsJCPfPMM9q9e7d++MMfatSoUaqurg731NrM1q1bNXPmTO3YsUPl5eU6d+6c8vLydObMmZBxI0eOVG1trb1t3LgxTDNuO3379g1Z0969e+1jL7/8spYsWaJly5apsrJSbrdbw4cP16lTp8I4Y3OVlZUhay4vL5ckTZgwwR4TDef6zJkz6t+/v5YtW9bi8dac38LCQq1fv17r1q1TRUWFTp8+rTFjxqixsbG9lnHFLrXub775Rp988ol+/vOf65NPPtG7776rAwcOaNy4cc3GTps2LeQ58NZbb7XH9K/a5c63dPnndbSdb0kh662trdUvf/lLORwOPfzwwyHjIul8t+ZnVru9vi202ve//33r8ccfD9nXp08f6+mnnw7TjK69uro6S5K1detWe9+UKVOsBx98MHyTugaef/55q3///i0eO3/+vOV2u60XX3zR3vftt99aLpfLevPNN9tphu3jySeftG6//Xbr/PnzlmVF57mWZK1fv96+3Zrze/LkSSsmJsZat26dPeYvf/mLddNNN1n/8z//025zN3Hhulvyxz/+0ZJkHT161N6XnZ1tPfnkk9d2ctdQS+u+3PP6RjnfDz74oDVs2LCQfZF+vi/8mdWer2+uULVSMBjUrl27lJeXF7I/Ly9P27ZtC9Osrj2fzydJzf6Q5JYtW9S9e3fdcccdmjZtmv03ESPZwYMH5fF4lJ6erkmTJunw4cOSpKqqKnm93pBz73Q6lZ2dHVXnPhgM6le/+pX+6Z/+KeQvB0Tjuf6u1pzfXbt2qaGhIWSMx+NRRkZGVD0HfD6fHA6HbrnllpD9a9euVXJysvr27as5c+ZE/JVZ6dLP6xvhfH/xxRfasGGDpk6d2uxYJJ/vC39mtefrOyK+Kf168NVXX6mxsbHZH1FOSUlp9geXo4VlWSoqKtIPfvADZWRk2PtHjRqlCRMmqGfPnqqqqtLPf/5zDRs2TLt27YrYb90dPHiw3n77bd1xxx364osv9MILL2jo0KHav3+/fX5bOvdHjx4Nx3Sviffee08nT55Ufn6+vS8az/WFWnN+vV6vYmNjlZSU1GxMtLz+v/32Wz399NOaPHlyyB/L/elPf6r09HS53W7t27dPxcXF+vTTT+23hyPR5Z7XN8L5Xr16tRISEvTQQw+F7I/k893Sz6z2fH0TVFfowr/5Z1lW1P4dwIKCAv3pT39SRUVFyP6JEyfa/87IyNDAgQPVs2dPbdiwodmLM1KMGjXK/ne/fv2UlZWl22+/XatXr7Y/rBrt537FihUaNWqUPB6PvS8az/XFXM35jZbnQENDgyZNmqTz58/rjTfeCDk2bdo0+98ZGRnq3bu3Bg4cqE8++UQDBgxo76m2iat9XkfL+ZakX/7yl/rpT3+qTp06heyP5PN9sZ9ZUvu8vnnLr5WSk5PVoUOHZrVaV1fXrHyjwaxZs/Tb3/5Wmzdvbvb3Ey+Umpqqnj176uDBg+00u2svPj5e/fr108GDB+3f9ovmc3/06FFt2rRJjz322CXHReO5bs35dbvdCgaDqq+vv+iYSNXQ0KBHHnlEVVVVKi8vD7k61ZIBAwYoJiYmqp4DFz6vo/l8S9JHH32kzz///LKvdylyzvfFfma15+uboGql2NhYZWZmNrvsWV5erqFDh4ZpVm3PsiwVFBTo3Xff1QcffKD09PTL3ufEiROqqalRampqO8ywfQQCAf35z39Wamqqffn7u+c+GAxq69atUXPuV65cqe7du2v06NGXHBeN57o15zczM1MxMTEhY2pra7Vv376Ifg40xdTBgwe1adMmde3a9bL32b9/vxoaGqLqOXDh8zpaz3eTFStWKDMzU/3797/s2Ov9fF/uZ1a7vr5NPk1/o1m3bp0VExNjrVixwvrss8+swsJCKz4+3jpy5Ei4p9Zm/uVf/sVyuVzWli1brNraWnv75ptvLMuyrFOnTlmzZ8+2tm3bZlVVVVmbN2+2srKyrFtvvdXy+/1hnv3Vmz17trVlyxbr8OHD1o4dO6wxY8ZYCQkJ9rl98cUXLZfLZb377rvW3r17rZ/85CdWampqRK+5SWNjo9WjRw9r7ty5Ifuj6VyfOnXK2r17t7V7925LkrVkyRJr9+7d9m+zteb8Pv7449Ztt91mbdq0yfrkk0+sYcOGWf3797fOnTsXrmVd1qXW3dDQYI0bN8667bbbrD179oS83gOBgGVZlnXo0CFr/vz5VmVlpVVVVWVt2LDB6tOnj3XPPfdE7Lpb+7yOtvPdxOfzWZ07d7aWL1/e7P6ReL4v9zPLstrv9U1QXaF///d/t3r27GnFxsZaAwYMCPk6gWggqcVt5cqVlmVZ1jfffGPl5eVZ3bp1s2JiYqwePXpYU6ZMsaqrq8M7cUMTJ060UlNTrZiYGMvj8VgPPfSQtX//fvv4+fPnreeff95yu92W0+m07rvvPmvv3r1hnHHb+cMf/mBJsj7//POQ/dF0rjdv3tzi83rKlCmWZbXu/J49e9YqKCiwunTpYsXFxVljxoy57v9fXGrdVVVVF329b9682bIsy6qurrbuu+8+q0uXLlZsbKx1++23W0888YR14sSJ8C7sMi617tY+r6PtfDd56623rLi4OOvkyZPN7h+J5/tyP7Msq/1e346/TQgAAABXic9QAQAAGCKoAAAADBFUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAbnhffvml3G63Fi1aZO/7+OOPFRsbq7KysjDODECk4G/5AYCkjRs3avz48dq2bZv69Omje+65R6NHj9brr78e7qkBiAAEFQD8zcyZM7Vp0yYNGjRIn376qSorK9WpU6dwTwtABCCoAOBvzp49q4yMDNXU1Gjnzp266667wj0lABGCz1ABwN8cPnxYx48f1/nz53X06NFwTwdABOEKFQBICgaD+v73v6+7775bffr00ZIlS7R3716lpKSEe2oAIgBBBQCS/vVf/1X/9V//pU8//VQ333yzcnNzlZCQoN///vfhnhqACMBbfgBueFu2bNHrr7+uNWvWKDExUTfddJPWrFmjiooKLV++PNzTAxABuEIFAABgiCtUAAAAhggqAAAAQwQVAACAIYIKAADAEEEFAABgiKACAAAwRFABAAAYIqgAAAAMEVQAAACGCCoAAABDBBUAAIAhggoAAMDQ/wPGHBv0OBI+iQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh.plotImage(np.log10(sigmaInf))\n",
    "plt.xlim(0,200)\n",
    "plt.ylim(-200,200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\81805\\anaconda3\\envs\\simpegEMIP\\lib\\site-packages\\SimPEG\\Props.py:400: UserWarning: Cannot add model as there are no active mappings, choose from: ['cMap', 'etaMap', 'rhoMap', 'sigmaInfMap', 'sigmaMap', 'tauMap']\n",
      "  \"', '\".join(self._all_map_names)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 2.68 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "e_emip, data_emip, cur = get_em_data(sigmaInf, eta, tau, c, data_type='emip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\81805\\anaconda3\\envs\\simpegEMIP\\lib\\site-packages\\SimPEG\\Props.py:400: UserWarning: Cannot add model as there are no active mappings, choose from: ['rhoMap', 'sigmaMap']\n",
      "  \"', '\".join(self._all_map_names)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">> Compute Sensitivity matrix\n",
      " 1 / 1Wall time: 555 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "e_em, data_em, cur = get_em_data(sigmaInf)\n",
    "eref = geteref(e_em[:,0,:], mesh, option=None, tInd=20) \n",
    "ip = get_ip_data(sigmaInf, eref, eta, tau, c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = ip + data_em"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = np.logspace(np.log10(2e-5), np.log10(0.009), 51)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x22d2e8ae2c8>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGnCAYAAABLpnZwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbn0lEQVR4nO3deVzU1f7H8dcwLKYiuYWaS1qZC4qKS1qatmCU9tPU1mva1bqmZWZqLjngUJGlliV6Kyuz1TZtM43K1LLUSMqiXQ3SzCUFwQQZzu+PCa4IIuAw32F4Px+PedR8v+f7/X6mO9f5eM7nnGMzxhhERERE/ESA1QGIiIiIeJKSGxEREfErSm5ERETEryi5EREREb+i5EZERET8ipIbERER8StKbkRERMSvKLkRERERv6LkRkRERPyKkhsRERHxK0puRERExK9UqeRmzpw5tG/fnoiICF544QWrwxEREREfFGh1AGW1detWXnrpJZKTkwG45JJLGDBgAKeffrq1gYmIiIhPqTLJzffff0+vXr2oUaMGAJ06dWLVqlVcd911Zbo+Pz+fXbt2ERoais1mq8xQRURExEOMMRw6dIgmTZoQEFDGASfjIWvXrjUDBgwwjRs3NoBZvnx5sTaJiYnmrLPOMiEhIaZLly5m3bp1Zb7/d999Z9q1a2cOHDhgDhw4YNq0aWPmzJlT5uvT09MNoJdeeumll156VcFXenp6mX/zPdZzk52dTWRkJDfffDNDhgwpdn7ZsmVMmDCBhQsXcsEFF/DEE08QExNDamoqzZs3ByAqKoqcnJxi137wwQe0a9eO8ePHc/HFFxMWFka3bt0IDCx7+KGhoQCkp6dTp06dCn5KERER8abMzEyaNWtW+DteFjZjjPF0IDabjeXLlzNo0KDCYz169KBLly4sWrSo8Fjbtm0ZNGgQCQkJ5X7G6NGjGTx4MFdeeWWJ53NycookSgX/cTIyMpTciIiIVBGZmZmEhYWV6/fbK7OlcnNzSU5OJjo6usjx6OhoNmzYUOb77NmzB4Aff/yRTZs20b9//xO2TUhIICwsrPDVrFmzigUvIiIiVYpXCor37duHy+UiPDy8yPHw8HB2795d5vsMGjSIgwcPUqtWLZ599tlSh6WmTZvGxIkTC98X9NyIiIiIf/PqbKnjZykZY8o1c6k8vTwhISGEhISQmJhIYmIiLperzNeKiIhI1eWVYakGDRpgt9uL9dLs2bOnWG+Op40bN47U1FQ2b95cqc8RERER3+CV5CY4OJioqCiSkpKKHE9KSqJXr17eCEFERESqCY8NS2VlZfHLL78Uvt++fTspKSnUq1eP5s2bM3HiRIYPH07Xrl3p2bMnTz75JGlpaYwZM8ZTIZRIw1IiIiLVi8emgn/yySf069ev2PERI0awZMkSABYuXMhDDz3EH3/8QUREBI888gh9+vTxxONPqiJTyURERMRaFfn9rpR1bnyRkhsREZGqx2fXubFSYmIi7dq1o1u3blaHIiIiIl6gnhsRERHxWeq58aK4uDji4+NLPBcfH09cXJx3AxIRERGgGiQ3lTUsZbfbcTgcxRKc+Ph4HA4HdrtdCZCIiIgF/D65qaxF/GbOnInT6cThcHDttY+Rnn4Up9Od2DidTmbOnFmmBEhEREQ8y6vbL/ibmTNncuSIjQceuJ1XXw0Abue00wbw5pv5/PjjFnr0uIE77qiFwzGpsH1BYlOQAMXFxWG325k5c2ax+8fHx+NyudTDIyIiUg5+33NT2fr3vwn4FXABdfn7786kpETx4oudGT/+bB5/fCJwCIfjSgICnsfh+IvBgx9jzJjpQNmGt0RERKTs/L7nprJXKF679jnAQVBQKEePtiAy8gbs9gh++60WBw82wuU6G6gFdMGYLsBwli+H5cshOPgPmjS5itatz8HhWMbevfOYP38i992n3h0REZGK8vuem8rcOPPYIabc3Eyczmv4+uvpDBqUwr59F3P0aFt27szg6qunAYOw2RzAG7h7eiA3tzE7dkTy00/XAyt4/PGJBAQcwOG4kBYtXic7+2rWrz8IBKl3R0REpKxMNZGRkWEAk5GR4ZH7OZ1OAxin01nq8RO9v/XWKcbh+NBccskK06TJahMY+I2BHAOmhNffpkaNbww8abp2fdq8+eYf5t577y9y39jY2GKxHBtTbGysRz63iIiIN1Xk99vvh6Uqi8vlKhw2OlbBe5fLVax4+Njz7uO1+fBD93un00ls7H3Y7R1xuTpQr94lHD7cmiNH2gB1OHKkA9CBL7+Eq68GmEJw8NW88UY2O3d+TVZWXV58cUqRZ0DR3iUREZHqQCsUV6Ky1socnwQd+/7228ezatVPJCXtJTk5n2++CQCigPASnpgLbKVJk10MHdqCrKyPeeaZu3E641S7IyIiVVKFfr8rrR/JRyxYsMC0bdvWtG7d2qPDUp5S1uGtY48FBQUbONN07XqfiYx8x9Su/amBvScY0jpkwsO/M5ddttn077/QQHiZniUiIuILKjIspZ4bi3mid2fmzJlkZx9m9eofePfd3Tz77Fag2z+v0BKeuoP69X9h4MAGGLOW556bqN4dERHxSeq5KYWnC4q9qSK9O8HBwQYCzIAB08zAgW+bpk3fNwEBqQZcxXp3bLYM06zZN+baa5PNddfNNxCk3h0REfEJFfn99vup4P6gtOJlp9NZuIbPsb05OTk5OJ1xvPtuAt26pZCefjl5eW1ITt7GlCkfYLPdC6wEMjGmDunpHVi2rAuvvDIeOIjD0Ys2bV7m5Zd3Ehd3X7HCaBEREV+l2VJVQGnDQAXJxslmZhW879LlHN5772WMuZ/g4GByc11cfvk95OR0JyWlDgcOdAAaAJfw449www0AYzjjjK4EBrYiMzOfefOcGroSERGfpZ4bP1Hx3p1YVq16gH79vuGvv/qRkVGDp57aQEzM28DLwH6gAXv2XM706a0JCzvKI49chsPxOxMnzinyLC0qKCIiPqESh8l8gq/PlvKmis3MOs1AH1O//rPGZvu1WL1OjRpfm6uuSjZXX/24gdO0qKCIiHiUam5KUJnbL1Q1Fendyc09jNN5Kfv338y99z7Hk09+Ro8ebxEcnAzkc+RIR95+uwtvvnk7cJD586/hssuS2bw5HIdjnraMEBER76vEZMunVOXZUt5U1t6d/Px8s3r118Zmu8XASwbSS1hjx2Xga9Os2QoTH/+jmTx5traMEBGRctH2C3LKyrKtBIDNZmPjxrcw5imCg58jNzeXoUOnkJ9/ARs3BrFr1zkYcy7QkfT0jrgvn0JIyAA++eQIZ5zxG3l5IcyaNb3I/UFbRoiIyCmqxGTLp6jnxrNOtkHo0aNHzapVW8yIEW8bmG8gpYSenRxTq9ZmA/ea669/zBw9qvV0RESkKPXciFeUddp5//6d2LTpHcDxz7Tz04mIuJO//urMrl3tgeZkZ3cFuvLyy/DyyxlAW3r1epp//etGoOwrOIuIiBTw+4Ji8byKTzu/nW+/ncGYMV9y+PAZLF36OZdfvoK6dT/CPeU8DBjKhg3/plWrEMLCtvHii+1xOD4gLu6+Is9SYbKIiJxQJfYk+RQNS3lX+aedBxi7vacBhwkI2FTCNhEHTP36H5rJk78zkyY9pMJkEZFqQlPBxWeUv3cnjry8DTidgeTnd2fkyKkMG/YWDRuuBvYBp7N//yU8/HA75sy5m5o1t5CaehWff55JQIAdh8OhaeciIuJWicmWT9Aifr6rrL07e/bsN/HxHxiYZWBzCQsJ7jRnnfWmgX4mNva+Uu8tIiJVS0V6bmzGGGNpduUlFdoyXSpVeYqFC3ph3IXJDTj77Ans3t2d7OweQI1jrjyIzfY+xrzGvfd2JT5+ermfJSIivqNCv9+Vlmr5GNXcVF2lTTtPSfnF3HzzW6ZRo5UG/izSo2OzZZm2bVPM7Nm/mhkz7i9zDZCIiPiOivx+K7kRn1aewuR77plu4HwDDxvYflyik20aNPjIwFAzY8YDpd5bRER8hwqKxe+UpzB59uwHcDqvwJhJzJjxNNCN+vWfwWbbgTE12bfvYuA17r//Tmy213E4vsThuK/EoSoREam6VHMjVV5Jiwoee3zGjHupVasPS5Zk8/PPkRjTsrCN3X6A88//jWnTmrFp0+MEBqouR0TEl1Tk91s9N1Llnax3JzDQzrRpl/Hjj4OYPv0ZoAswF9iFy1WXzz7rxIAB9UlIGIHDcZQ77nikyH00pVxEpGpRz41UG8f38Nxzz3QeeuhLatb8D4cPXw7UKmxbq9YWbr31NGrUWEFCwrQSkycREal8Ffn91t5SUi2UNHQ1e/YD1K4dj8MxlNtum8yuXT1YvfoMjhzpRXZ2Zx55BOAOGjfuROfO3TAGZs3SlHIREV/nk8NSgwcPpm7dugwdOrTYuXfffZfzzjuPc889l8WLF1sQnVRFJxu6OuOMmqxYMYTDhy9k+fJkbLZJwHdALf7443IGDqxPWFga7757Dg7Ho1oNWUTEh/nksNSaNWvIysriueee4/XXXy88npeXR7t27VizZg116tShS5cubNy4kXr16p30nhqWkrIqSFSCgoI5ejSKmjUncPjwAKDmPy1ygDcYNiyTV14Zw/33l1zQLCIip85vCor79etHaGhoseObNm2iffv2nHnmmYSGhnLFFVewevVqCyIUf3Xs8FVubg5OZwyHD1/L6NEO+vZ9Fbv9ayAEuIHXXhuD3f4zDsdepk9PUGIjIuIjyp3crFu3joEDB9KkSRNsNhsrVqwo1mbhwoW0bNmSGjVqEBUVxfr16z0RK7t27eLMM88sfN+0aVN27tzpkXuLlFSXUzBstXjxXC6++Eeyss4jPn4ljRu/BWQC5wKP8eCD44iO3sqPP+YQFxdXbNjq2GeoJkdEpHKVu6A4OzubyMhIbr75ZoYMGVLs/LJly5gwYQILFy7kggsu4IknniAmJobU1FSaN28OQFRUFDk5OcWu/eCDD2jSpMkJn13SCJrNZivvRxApUWl1OQXna9Sowb33XoEx8TgcN2CzjcCY8eTntyEpqQNt2rho0KAf+/ZNxxhwOIqvu+N0Or36uUREqp1TWRIZMMuXLy9yrHv37mbMmDFFjrVp08ZMnTq1XPdes2aNGTJkSJFjn332mRk0aFDh+/Hjx5sXX3yxxOuPHDliMjIyCl/p6enafkE84vhtG6ZNm2EgxgQGfnTcjuWbzKWXPmNycrTVg4hIRVm+/UJubi7JyclER0cXOR4dHc2GDRtO+f7du3fn22+/ZefOnRw6dIiVK1fSv3//EtsmJCQQFhZW+GrWrNkpP1+kpKGrBx64D6ezJ3l5l9C//yQaNnwbOAJ048MPbyYk5Hccjt04HPerLkdExAs8us7Nvn37cLlchIeHFzkeHh7O7t27y3yf/v3789VXX5GdnU3Tpk1Zvnw53bp1IzAwkLlz59KvXz/y8/OZMmUK9evXL/Ee06ZNY+LEiYXvMzMzleDIKSvL0NWqVVfx/vtfMmXKNr79tg/QFEhk9uw9/P33Nu67rxXBwd6PXUSkuqiURfyOr4MxxpSrNqa0GVBXXXUVV1111UnvERISQkhICImJiSQmJhZusChyKkorBj424YmJ6cqXX77Pt98Ox2YbjTHTyMlpysMPn8Hjj/9J166ruOSS34mLm1HsPloMUETk1Hh0WKpBgwbY7fZivTR79uwp1pvjLePGjSM1NZXNmzdb8nypnv43fHUv+fmJ3HHHY8AdwB8cORLOp5+OYNas67j44ufIyyt+nRYDFBGpOI8mN8HBwURFRZGUlFTkeFJSEr169fLko0R8Vkl1OY899hBO5xlAK5o0mQP8CZzNmjUjqFVrBwkJ6cyapcUARUQ8odzDUllZWfzyyy+F77dv305KSgr16tWjefPmTJw4keHDh9O1a1d69uzJk08+SVpaGmPGjPFo4GWlYSnxtpPX5WTxf/+Xwb///REpKZeRm3sW06cDXMzo0WHMnDne+0GLiPiRcm+/8Mknn9CvX79ix0eMGMGSJUsA9yJ+Dz30EH/88QcRERE88sgj9OnTxyMBV5S2XxBftHnz9/To8QLGTAVCgTyuvPJXXnqpNXXqaA0nEZGK/H775N5SlUHJjfiigiGswMCzyMt7GHBvFhsSspc5c/IYN64xWqdSRKozv9lbypMSExNp164d3bp1szoUkSKOrc05enQ7DsdWIAbYRk5OQ+64ozENG37B998XX81b2ziIiJyY3yc3mi0lvqikouNZs2bhdPYC2lO37gIgl/37z6ddu3xGjdpGwY4lmlElIlI6DUuJWCAuLg673V7irKj4+Hjy8vIIC+vO5Mm1yM/vC0CdOru47LL3eeON0ZpRJSLVhmpuSnDsbKmffvpJyY1UKQcOHKRr17ls2zYOaARA/fqb+OijSCIjQ6wNTkTEC5TclEI9N1KVBQU1JC9vKjAeCAKOcuWVv7F06dnUq6eKYxHxXyooFvFD7mGqfQQHTwfaY7evBIJ4771zaNjwIBdd9GqRVY6PvU5FxyJSHSm5EfFhxxYe5+Tk4HQOx+W6kg4d7sZm+478/LqsW3cNoaHbeO21jGLXqehYRKqjStk405dohWKpqkqaUVXwT4fDwZ132li37jy2bBnMkSOtuOYaaN9+G5de+jbz52sbBxGpvlRzI+KjTjajqmDn8Hff/Yxrr/2Ow4dvxl2Pk8UFF7zP2rXDUMeNiFR1KiguhZIb8Wcul4uQkEhcrkTgIgDOPDOdt98Op0uXYGuDExE5BSooFqmmHnjgAVyu7wgK6g/cCmSwc2czoqJsjBqVXrgAoIhIdaDkRqSKO7Y2Jzf3CE5nM6Addvs7QBDPPNOMunW389FHf5d4rWZUiYi/8fvkRntLiT87UdGx0zkGl+sqzjprMvAnf//dkksvDWHQoHQOHSp6rWZUiYi/Uc2NSBVWlqLjiIg+3HDDLo4e/RcAtWrtZ/Dg5bzwwi2aUSUiPk8FxaVQciPV2aFDh+jefTo//HAX0ApwERX1Pl98MYBAv18QQkSqMhUUi0iJQkND+f77xwkM7AI8DdhJTh5AmzY7SU+vFn+/EZFqRMmNSDXh3sYhg+DgscANQCa//nomZ5+dxYsvHrI6PBERj1FyI1INHL+Nw6xZbYAuwJccPRrKv/4VSkTERyVOGdeMKhGpavw+udFsKanuSppR5X4/AuhFSMgCAL777hIaNvyZH37IL3atZlSJSFWigmIRP3eyGVWHDx9mw4bTWbfu30BDbLZsHnvMxYED84slRSIi3qbZUqVQciNyYsYY5s59mSlTGmNMv3+OPoHDsY9Zs2ZYGpuIVG+aLSUiFWKz2Zg06Qa++qoB4ADygf/w0kuj2b3b4uBERMpJyY2IFHrnnRVAPDbbVcBBfvklnHPPPcjnn+dZHJmISNkpuRERoGjhcV7e23TvfgfwPVlZp3PhhfksWpRtdYgiImWi5EZEis2oCggIYOPG57n66oeAt8nPD2bs2FqMGLGfPHXiiIiPU3IjIrhcrhJnRb3xxrPcdtsHhITMBmDp0vpERe1j//6i12stHBHxJUpuRIS4uLgTTvdeuHAB6en/pmHD/wBZfPNNA8499yBff+2eaKm1cETE1/h9cqNF/EROXcOGDfn998dp33408CsHDpxOVFQO11//gtbCERGfo3VuRKRc+ve/ng8+GAn0B44yaNAKli8fZnFUIuKvtM6NiFS61atfJjBwMPA8EMSKFUN4+OHDVoclIlJIyY2IlIt7d/G/CQwcDTwGBDBlSk2mTs2ievQDi4ivU3IjImV27JTxo0dzGDv2JyAOgNmza3PLLYfIzy/1FiIilU7JjYiUSUm7iycmLuDOOw8C4wF4+ulQhg49pLVwRMRSSm5EpExOtBbOo48+ysSJQdSqdRuQx/LloURHZ3LkSNHrtRaOiHiLkhsRKZPS1sKZO3cuv/wSS1jYv4EjrFlTh969M8nMdJ/XWjgi4k1KbkTEIxo1asS2bY/SoMFwIJMvv6xD9+6HmDZtjtbCERGv8snkZvDgwdStW5ehQ4eW65yIWKtevXps2/YMjRvfCOzlxx9DefDBi5k2bbYSGxHxGp9MbsaPH8/SpUvLfU5ErBcaGsqvv76KzdYX2AN04Z13xnLokLVxiUj14ZPJTb9+/QgNDS33ORHxDXPmzMGYVCAa+Itvv61Nv37ZHNZafyLiBeVObtatW8fAgQNp0qQJNpuNFStWFGuzcOFCWrZsSY0aNYiKimL9+vWeiFVEqoBjp4wfOvQpZ5wxHMggObkW0dGHi82iEhHxtHInN9nZ2URGRrJgwYISzy9btowJEyYwY8YMtmzZQu/evYmJiSEtLa2wTVRUFBEREcVeu3btqvgnERHLHb8WTu3atfnhhxeoW/cGIIvPPqvJgAFHyM21OlIR8WeB5b0gJiaGmJiYE56fN28eo0aNYvTo0YB7DYzVq1ezaNEiEhISAEhOTq5guGWXk5NDTk5O4fvMgjmpIlJpSloLp27duqSmPk3r1jdw6NAyPvroNIYMyWH58hACy/0nkIjIyXm05iY3N5fk5GSio6OLHI+OjmbDhg2efNRJJSQkEBYWVvhq1qyZV58vUh2daC2cRo0asXXr4zRocAuQw7vvhnDjjbm4XN6PUUT8n0eTm3379uFyuQgPDy9yPDw8nN27d5f5Pv3792fYsGGsXLmSpk2bsnnz5jKdO9a0adPIyMgofKWnp1fsQ4mIR7Ro0YL16+8lOPhfwFFefTWYUaPyiuxFpVWMRcQTKqVT2GazFXlvjCl2rDSrV6+u0LljhYSEEBISQmJiIomJibj0V0QRy7Vp04ZRoxqyaNENwCs891wgp53mYuFCO/fd9796HRGRU+HRnpsGDRpgt9uL9dLs2bOnWG+Ot4wbN47U1NQT9vCIiHctXLiQUaPCgJFAPv/9r51LL/1IqxiLiMd4NLkJDg4mKiqKpKSkIseTkpLo1auXJx8lIlXY4sWLGT7cBowD4OOP+zF8+ItKbETEI8qd3GRlZZGSkkJKSgoA27dvJyUlpXCq98SJE1m8eDHPPPMM33//PXfddRdpaWmMGTPGo4GXVWJiIu3ataNbt26WPF9ESrZ06VICAxcDTwABvPnmUH791eqoRMQf2IwxpjwXfPLJJ/Tr16/Y8REjRrBkyRLA3e380EMP8ccffxAREcEjjzxCnz59PBJwRWVmZhIWFkZGRgZ16tSxNBYR+d+aODZbDYz5GOjJeeflkpwcTK1aVkcnIr6iIr/f5U5uqiolNyK+49jF/u655x5atOjJ7t3vAY0YNszFsmV2yjEHQUT8WEV+v31ybylP0rCUiG85fhXj4OBgvvjiTYKDbwSO8tprdubNszpKEanK1HMjIl4VFxeH3W4vVjz83nvvMWDA+8ACbLZ8PvwwgIsvtiZGEfEdGpYqhZIbEd83bdp0HnywNTCSunXz2LIlkBYtrI5KRKykYakSaFhKpOqIj3dy4YUvAckcOBDIoEH5/P231VGJSFWjnhsR8Sl//PEHHTsOZN++94GGjBgBzz6LCoxFqin13IhIlde4cWOuvDICuA5w8dxzsHDh/85r/ykRORklNyLic84++2zgY2AKABMmGNau/d9MK7vdbml8IuLbKmXjTF+ijTNFqp6ZM2dijCE2NhboSl7e9cTEHObvv5/T/lMiclKquRERnzVt2jQefHA+8AnQnfr19/HTTw2oV8/iwETEa1RzIyJ+JSEhgaCgPOAq4Df272/A0KGQm2t1ZCLiy5TciIjPio+P5+jRowQG7gcGAJmsWQO33QbVo89ZRCpCyY2I+KRjt2k4evQo3bvXBK4FXDzzDDz8sNURioiv8vvkRov4iVQ9x+8/BbB+/XqaN08F7gTgnnvgzTctDFJEfJYKikXE55xo/6ndu3fTunVrDh26H7iD004zrFtno2tXa+IUkcqnvaVKoeRGxD9s3LiR3r37cfTo68AVNGoEmzZBs2ZWRyYilUGzpUTE7/Xo0YNFix7HvYLxN+zeDQMHwqFDVkcmIr5CyY2IVDmjRo3ittv+BQzAZvuTr7+G668HrdUpIqDkRkSqqPr169O8uQ1jBmKzHeG998DhcJ/T/lMi1ZvfJzeaLSXin4KDg0lLSyM09AeMGQlAQoJh5MgXtP+USDWngmIRqbKO3UjT5UoE/gP8wZQpLzN79kSrwxMRD1BBsYhUKzNnzsTpdP6zMe5dwLdAY77+eiL5+RYHJyKWUXIjIlXazJkzCQ4OBv4GriUg4AirV8PcuVZHJiJWUXIjIlVafHw8ubm5/yQ4qeTn3w7A9OnwxRfWxiYi1lByIyJV1rHbNOTk5HDVVVcBT2OzLSMvzz09/OBBq6MUEW9TciMiVVJJ+0+tWLGCs88+G2NuxW5PY8cOuPVW7SAuUt0ouRGRKsnlchVJbABsNhsff/wxwcFHcLmGEhDg4rXX4KmnLAxURLxOU8FFxO889dRT3HrrrQQGTiUvL4EaNdz7T3XoYHVkIlJemgpeAi3iJ1L9jB49mujoaPLyZnP66Z9z5Ahcey1kZ1sdmYh4g3puRMQvpaWlERERwaFDNQgN3cahQ7UZPVpDVCJVjXpuRET+0bx5c+bOnQvs5ciRYdhshsWL4dVXrY5MRCqbkhsR8VsFw1NHj66iSZPnABgzBnbutDgwEalUSm5ExG/NmjWLiIgIQkND2bnzFpo23c2BAzBqFDid2jlcxF8puRERv2W325k3bx79+vUD8tiz53JCQvJZvRpiY//QzuEifirQ6gBERCpLwRo4DoeDs88+m19//Zr69e8nJ2cmQUHzufbaIIsjFJHKoJ4bEfFrBTuH//rrrwDs3x9LePi3HD0axPDhkJdncYAi4nFKbkTE7/1v53AAQ1DQrYSFGTZtgoQES0MTkUrgk8nN4MGDqVu3LkOHDi1yPD09nb59+9KuXTs6duzIa6+9ZlGEIlKVFN05HH7//XMuvfQtAJxOSE62MjoR8TSfTG7Gjx/P0qVLix0PDAzk0UcfJTU1lQ8//JC77rqLbC05KiKlOH7n8GHDhgGwYsUwrrjiMHl58K9/wd9/WxyoiHiMTyY3/fr1IzQ0tNjxxo0b06lTJwDOOOMM6tWrx19//eXl6ESkqihp5/Bly5bRokULXK48fv01mkaN4IcfYNo0i4MVEY8pd3Kzbt06Bg4cSJMmTbDZbKxYsaJYm4ULF9KyZUtq1KhBVFQU69ev90SsRXz55Zfk5+fTrFkzj99bRPzDiXYOX7FiBTabjR9//Izx478BYP58+OgjqyIVEU8qd3KTnZ1NZGQkCxYsKPH8smXLmDBhAjNmzGDLli307t2bmJgY0tLSCttERUURERFR7LVr164yxbB//35uuukmnnzyyfKGLyLVSFxcXJHEpkCnTp0YM2YMAC+//C9Gj84HYORIOHjQiwGKSKU4pY0zbTYby5cvZ9CgQYXHevToQZcuXVi0aFHhsbZt2zJo0CASyjEt4ZNPPmHBggW8/vrrRY7n5ORw2WWXccsttzB8+PATXp+Tk0NOTk7h+8zMTJo1a6aNM0UEcP8lqXXr1vz111/MnfsECxfeyq+/uutvnn/e6uhEpIDlG2fm5uaSnJxMdHR0kePR0dFs2LDhlO9vjGHkyJFcfPHFpSY2AAkJCYSFhRW+NHwlIseqX78+9913HwD33TeVxx47SEAAvPACLF9ucXAicko8mtzs27cPl8tFeHh4kePh4eHs3r27zPfp378/w4YNY+XKlTRt2pTNmzcD8Nlnn7Fs2TJWrFhBp06d6NSpE1u3bi3xHtOmTSMjI6PwlZ6eXvEPJiJ+6dZbbyUyMpIDBw7wzjvTmDLFffyOO+DQIWtjE5GKq5TtF2w2W5H3xphix0qzevXqEo9feOGF5Ofnl+keISEhhISElPmZIlL92O12HnvsMS666CKeeOIJPvvsPyxb1ont293r3zz8sNURikhFeLTnpkGDBtjt9mK9NHv27CnWm+MtiYmJtGvXjm7dulnyfBHxbX369OG6667DGMOUKXfw+OPuMsRHHoFvv7U4OBGpEI8mN8HBwURFRZGUlFTkeFJSEr169fLko8ps3LhxpKamFg5tiYgc7+GHH6ZmzZp8+umnZGa+wqBB4HLB2LFQ8SkXImKVcic3WVlZpKSkkJKSAsD27dtJSUkpnOo9ceJEFi9ezDPPPMP333/PXXfdRVpaWuG0S29Tz42IlCYuLo5nn32WGTNmADB58mTuvz+bmjVh/Xq4+mpVF4tUNeVObr788ks6d+5M586dAXcy07lzZxwOBwDXXnstjz76KE6nk06dOrFu3TpWrlxJixYtPBt5GannRkRKY7fbcTgc5OTk0KpVK3bu3Mkbb8zjggvcK/olJUVz4IDFQYpIuZzSOjdVSUXmyYtI9VCwTcPQoUN5/fXXqVWrFtnZuTRsuJO9exty222wcKHVUYpUT5avc+OLNCwlIiczc+ZMnE5n4aKh2dnZDBp0Ja++2hCA//4XvvzSyghFpDzUcyMi8o+QkBByc3MBaN++PVu3bmX4cBsvvghdu8IXX4DdbnGQItWMem5ERCooPj6e3NxcgoODAfjuu+9ISkpizhyoU8fdc6Pt7ESqBiU3IlLtFdTcOJ1OcnJyOP/88wH3hIRGjeD++93tpk+HPXssDFREysTvkxvV3IhIaY5NbAp2EH/xxRex2Wz88ssv3HHHHdx2G3Tu7N4xvGCLBhHxXaq5EZFqLS4uDrvdXpjYFBg6dChvvPEGnTt35quvvmLjRujZ072o39q10KePRQGLVDMV+f1WciMiUoINGzZwwQUXEBISQlpaGmeccQb/+Y+77qZ9e9iyBYKCrI5SxP+poFhExEN69uxJjx49yMnJYdGiRQA88AA0aADffQePP25xgCJyQn6f3KjmRkQqwmazMXHiRMD958iRI0eoXx8efNB93umEvXstDFBETkjDUiIiJ5CXl8fZZ59NWloaixcvZtSoUbhc7jVvUlLcG2smJlodpYh/07CUiIgHBQYGcueddwIwb948jDHY7fDII+7zTzzhHqISEd+i5EZEpBSjRo0iNDSU1NRUPvjgAwD69oXBg8HlgrvvtjY+ESnO75Mb1dyIyKkICwtj1KhRgLv3psBDD7lnS61eDe+/b1V0IlIS1dyIiJzE9u3bOeecc8jPz2fr1q1EREQAMHkyzJkDbdrAN99oarhIZVDNjYhIJWjZsiVXX301AI8UFNwA997rnhr+ww/u+hsR8Q1KbkREyqBgWvgLL7zAn3/+CUBYmHtKOEBsLBw4YFV0InIsJTciImXQs2dPzj//fHJzcwsX9QO45Rb3isV//fW/REdErKXkRkSkjAp6bxYuXMjff/8NQGAgFNQZL1gAP/1kVXQiUkDJjYhIGQ0ePJgWLVqwd+9eXnrppcLj0dFw5ZWQl+cuMhYRa/l9cqOp4CLiKYGBgdx6660AvPXWW0XOzZnj7sV5+2346CMrohORApoKLiJyEnFxcdjtdmbOnMmWLVvo0qULtWvXZv/+/cyePRuXy0VcXBzjx7s31OzYEb76Cux2qyMXqfo0FVxEpBLY7XYcDgfx8fFERkbSsGFDsrKyuO2223A4HNj/yWJiY6FuXfeaN888Y3HQItWYem5ERMogPj4eh8OB0+nkhx9+KKy5cTqdzJw5s7Dd/PkwYQKccQb8/DPojxuRU1OR328lNyIiZVSQ4AQGBpKXl8eZZ57J77//XqTN0aPQoQP8+CPMmAH33WdRsCJ+QslNKZTciIgnhISEkJubC4DNZmPfvn3Uq1evSJvly+HqqyE0FHbsgONOi0g5qOZGRKQSxcfHk5ubS3BwMADGGD4qYWrUoEEQGQmHDsExuzWIiJcouRERKYNja25ycnLo2bMnALNnzy7W1mZzFxeDuwbnr7+8GamIKLkRETmJYxObguJhh8MBQHJyMs4S9l34v/9zTwk/dAgefdSb0YqI3yc3WsRPRE6Vy+UqNiuqT58+hcNTe/bsKXZNQEDR3httqiniPSooFhGpoEsvvZSPPvqIxx9/nNtvv73Y+fx86NQJtm4FhwNmzfJ+jCJVnQqKRUS86LLLLgPggw8+KPF8QIA7qQF3783Bg14KTKSaU3IjIlJB0dHRAKxZs6Zwevjxrr4aIiIgI0O1NyLeouRGRKSCjt2K4YsvviixzbG9N48+qt4bEW9QciMiUkEBAQEnHZoCGDIE2rd3997Mn++t6ESqLyU3IiKnoGBoqrTk5vjem4wMLwQmUo0puREROQUFPTdffvkl+/fvP2G7oUOhXTv3sNRjj3kpOJFqSsmNiMgpaNKkCe3bt8cYw8cff3zCdgEBULBMzrx56r0RqUw+mdwMHjyYunXrMnTo0CLHDx06RLdu3ejUqRMdOnTgqaeesihCEZH/KcvQFMCwYdC2rbv35vHHvRCYSDXlk8nN+PHjWbp0abHjNWvWZO3ataSkpLBx40YSEhJK7QYWEfGGY5Ob0tZFtduL9t5kZnojOpHqxyeTm379+hEaGlrsuN1up2bNmgAcOXIEl8tV6h8kIiLeULAVQ1paGj/99FOpba+5Btq0cW/HcKLem7i4OOLj40s8Fx8fT1xcXJnaiFRX5U5u1q1bx8CBA2nSpAk2m40VK1YUa7Nw4UJatmxJjRo1iIqKYv369Z6IFYCDBw8SGRlJ06ZNmTJlCg0aNPDYvUVEKqJmzZr07t0bOPnQVFl6b+x2Ow6Ho1jyUrCBp91uL1MbkWrLlNPKlSvNjBkzzBtvvGEAs3z58iLnX3nlFRMUFGSeeuopk5qaau68805Tq1Yt89tvvxW26dKli2nfvn2x186dOwvbrFmzxgwZMuSEcezevdv06tXL7N69u0xxZ2RkGMBkZGSU7wOLiJTB7NmzDWAGDBhw0rZ5ecacd54xYMy995bcxul0GsA4nc4S35elTWxsbJH2x98/Nja2HJ9QxBoV+f0ud3JT5OISkpvu3bubMWPGFDnWpk0bM3Xq1HLd+2TJjTHGjBkzxrz66qslnjty5IjJyMgofKWnpyu5EZFKs2XLFgOYWrVqmZycnJO2f+01d3ITGGhMcnLJbQqSleDg4GKJTVnalJQQlXZcxBdZntzk5OQYu91u3nzzzSLtxo8fb/r06VOue5eU3Ozevbvww2VkZJh27dqZr7/+usTrY2NjDVDspeRGRCqDy+UyDRs2NIBZu3btSdvn5xszZIg7wenQwZgjR0puV5C0BAcHn/BepbUpSw+QiC+rSHLj0YLiffv24XK5CA8PL3I8PDyc3bt3l/k+/fv3Z9iwYaxcuZKmTZuyefNmAH7//Xf69OlDZGQkF154IbfffjsdO3Ys8R7Tpk0jIyOj8JWenl7xDyYichJl3YqhgM0GixZBw4awdSs4ncXbxMfHk5ubS3BwMLm5uSUWEJ+szcyZM3E6nTgcDkJCQnA4HDidTmb+U/ijwmTxS6eSTXFcz83OnTsNYDZs2FCk3X333WfOO++8U3lUhS1YsMC0bdvWtG7dWj03IlKplixZYgDTrVu3E7Y5vg7m9dfdvTcBAcbceuviwjqY8tTcjBkzxvz8889m5syZJ+yVOVHvjoauxNf51bCUp6mgWEQqW8Ff8Gw2m9m3b1+JbUpKGq6/3p3gQKqZOfOBEyYWs2bNKjxe0MZmsxUZeq9Ro4YBTKtWrczPP/9c5Jknqt3R0JX4MsuHpYKDg4mKiiIpKanI8aSkJHr16uXJR4mI+JwmTZoQERGBMYaPPvqoxDbHDhMVDAe1aPEw8AfQlsOHp+FyuXA6ndx777189NFHxMXFcfnll/PII49w55134nK5cLlc9O/fH2MMDRo04LTTTgPca4ABbNu2jaCgoMKp4RdeeCG5ubnUrVsXh8PBkCFDSEtLKxaThq7EL5Q3gzp06JDZsmVL4cyAefPmmS1bthRO9S6YCv7000+b1NRUM2HCBFOrVi2zY8eO8j7KIzQsJSLeNHHiRAOYUaNGldru+N6UG2982YAxNpsxn35qTH5+vhk1alSxSRFLly4tvMfevXtNenq6Mcbd/uDBgyY1NdUkJSWZ5557rnBihdPpNNOnTzcBAQHF7nfWWWeZkSNHmgMHDmjoSnySV4al1qxZU+IspBEjRhS2SUxMNC1atDDBwcGmS5cuZZo5UNk0LCUi3rBq1SoDmGbNmpn8/PxS2x6fTIwc6R6eOuccY+65x504BAQEmBtvvNEsWLDAfPnllyY3N7fMsRxf35ORkWHee+89M2XKFHPmmWcWDmmFhYWZuLi4IsnWmDFjisSvoSuxitdrbqoC9dyIiDdlZ2cXJgjff//9CduVVAdz4IAxZ55ZUH8z3wDmiSeeqLRYMzMzzfvvv28GDx5cGEN+fr4JCwsrTHomTJhg1q5da/Ly8sq07o6Ipym5KYV6bkTEWy655BIDmMcee6zE86X1gqxaVZDcGDN8+DOVHuvxsWRkZJghQ4aYoKCgIr3zDRs2NKNHjy48fvzQlVZDlspieUGxiIhAp06dAAoLdo9VUOB7bMHusQW9mzbFc+ut7rbr14/k0KHKjbWgeLkgljp16vD666+TkZHB9ddfT8eOHalbty579+5l8eLFHD16tHBNHafTSV5eHlC2/bBEvEXJjYiIhwUEnPiP1uOTiQLXXHMNd999Ny6Xizlz4KyzYMcOG5MnV26scXFxxWIBOO2003jppZf4+uuv+fPPPxkxYgQAo0ePJicnB6fTSWxsLGFhYYwdO5Y+ffowa9asIglOSYmciFdUYk+ST1DNjYh42+TJkw1gJk2aVKb2O3fuNM2bNzdnnnmmSU1NNcYY8/HH/xueWrWqMqM9uRMVD19wwQVFhq6aNGlizj//fNXliEdpWKoE48aNIzU1tXALBxERX5KRkUFMTAxpaWnUqlWLhg0bAtCvH9x+u7vNzTfDvn3WxXii3qY1a9YwfPhwOnXqRFhYGLt27eKLL74AKNwS4thrtF6OeIvfJzciIr4qJyeHwYMH88033xAeHs6qVato0KBB4fnZs6FtW/jjDxg1yt2PY4UTDV0FBQWxdOlStmzZwp9//snbb79dZL+/Y/e6Sk1NBVBdjnhFoNUBiIhUR/n5+YwYMYI1a9ZQu3Zt3n//fVq2bFmkTc2a8NJL0KMHvP02PPEEjBljUcAnERISQkpKCt988w333nsvQ4YM4Z133sHhcOByuVi0aBE2m40LLrgAh8MBuAupVZcjlaISh8l8gmpuRMTbylJz8+yzzxrABAYGmg8++KDU+82d6669Oe00Y/4pyfE5J1vFuFatWsUWfw0MDFRdjpxURWpu/L7nZty4cYwbN47MzEzCwsKsDkdEBICffvoJgFGjRnHZZZeV2nbCBFi9Gj74AK6/HjZuhJAQLwRZDieqyyl4n5ubS1RUFE8//TQrV64kPz+/cBp5mzZtCtvHxcVht9tL7MWJj4/H5XKpNkdOyu+TGxERX3TTTTdxwQUXcNZZZ520bUAALFkCHTvC11/D9Okwd26lh1gupSUcxyYqgwYNYvLkycyZM6fw2CeffMKwYcMA+Pvvv3nooYeKXXfs8JXIyaigWETEAm3atOHKK6+kffv2ZWrfuDE8/bT73+fNc/fiVEXx8fHMmTMHp9OJy+Vi9OjRLFy4sLDIeNeuXdSoUQOHw8Edd9xReI3qcqQ81HMjIlJFXHUV3HYbLFoEI0bAN9/APzPHq4SSkpSnnnqK5s2b43A4yM/PJzk5mSNHjgCwYMECEhMTMcYQGxurxEbKzO97bhITE2nXrh3dunWzOhQRkULJycksWbKETZs2leu6OXPc08N377Z2enhFlFaX43Q6Mcbw7bffsmrVKgYPHgyA+ecDPvHEEyxYsADQejlycn6f3GgRPxHxRW+88QY333wzL7/8crmuq1kTXn4ZgoPhnXfgv/+tpAArwYnWywF3ghMXF0dAQAD9+/enc+fOAIVr3+zevZusrKzCYyXV32i9HCng98mNiIi/iYyEBx90//vEifDP+nh+49jhq7y8PGJjYwE4ePAg4E6ErrnmGmJjYxkwYACHDh1SXY4UoZobEZEq6M47YdUqd2HxDTf45vTwiigpSSmYHu5wOKhVqxYzZ84kIyMDgPfee486deoAMGHCBCU2Aii5ERGpko6fHj5hgrvQuKo72Xo5LpcLgGXLlrF06VLuvPPOwrqc+fPns23bNu68807Wrl1LYGCg1supppTciIhUUY0bw3PPwYAB7tqb9u3/t9lmVVXW9XLCwsI4ePAgxhgCAwPJy8vDGMPbb7/N3r17iYmJKbLNQwGtl1M9+H3NjWZLiYg/u+KK/9Xf3Hln1V3/pryOTVKOHj1amKx0796dyZMnF87AcjgcXHTRRezcuVN1OdWI3yc3mi0lIv5u8mT3ujf5+XDNNfD991ZHVLlKSlIKkplNmzbx7bffFh677LLLWLduHU2bNsXhcHDbbbcpsakGNCwlImKB66+/ns6dO3Puueee8r1sNveO4b/+Cp9+CgMHuguM69f3QKA+qKx1OQCTJk3iww8/LKzLWbRoET/++CN33303GzduVF2On1JyIyLiIQWzekpy/I9lhw4d6NChg8eeHRICb74J3bu7k5whQ9xDVMHBHnuEzyhrXQ7Axo0bMcYQFBTE0aNHsdlsfPzxx3z88cc0atSI3bt3F7tOdTlVn98PS4mIeEvBdOUNGzYUOe6txeUaNoR334XQUFi7FsaOrVorGHvasUlKbm5u4SrIvXr1IjQ0lEsvvbSwLic+Pl7r5fgTU01kZGQYwGRkZFgdioj4MafTaQADmLvvvrvwvdPpLNLum2++Ma+++qrZsmWLx2NYudKYgABjwJi5cz1++yrhRP/dC47PmDHD7Ny5s8ixgteECROsCFlOoCK/30puREQ87MILLzSAsdvtJf7AGmPM1KlTK/WH9JFH3MmNzWbMO+9UyiN8WmxsbIn/3Y1xJzOxsbFFjhX8bwWYgIAAc80115jNmzd7IVI5mYr8fqvmRkTEw3r16sWnn36Ky+UiODjYkuGNO+90b8vw1FNw/fXw+ecQEeH1MCxTnrqcgnqogrqc/Px8Xn31VV599VXOOussevfuzXPPPYfNZivxOhUd+x7V3IiIeFhBzY3dbic3N/eEO1hXJpsNEhOhb1/IynIv9Ldnj9fD8Hkl1eUAREZGEhgYyI4dO3j++ecL97c6/jpt0umb/L7nJjExkcTExCJTA0VEKkt8fDyffvop4N7rKCwsrMSVcr0hKAhefx3OPx9++QUuvRSSkiA83Kth+KwTrZcD4HA4uPvuuzHG8N133xEfH09QUBAzZszgmmuu4Y033lDRsQ/z++Rm3LhxjBs3jszMTMLCwqwOR0T8WMGP5YUXXliY4Bz7Y3nse2+pXx/ee8/dg7N1K/TpAx99BE2bejUMn1SW9XIKhpyO36m8du3a1KxZk6ysLGrXru3t0OUk/D65ERHxloIfy6ysrMLkBkpeXM6bWreGdevgkkvgp5+gd293gtOqlSXh+Izy1OXMnDmT++67j9zcXACysrKYNGkSCQkJdOzYkZ49e3L//fcXu4/qcqyhmhsREQ+Ji4s7Yc/MzJkzLf2BO+ccWL/e/c8dO9w9OD/8YFk4VU58fDy5ubkE/7Mq4qBBgzjnnHPYv38/a9as4YEHHuCiiy4iJyenyDWqy7GGkhsREQsMHTqUZ599luuvv95rz2ze3N2D064d7NzpTnC++cZrj6+yjh2SysnJwel0smLFCv71r3/x4osv0r59ewDWrVvH7Nmzi12juhzv07CUiIgFoqKiiIqK8vpzGzd2r14cHQ1btrhrcVatcm/bIMWdrOjY6XTyzTff8NZbb/H2228TGxvL/fffT25uLr179+aWW26xMvxqSz03IiLVTIMG8PHH0LMnHDjgnkW1fr3VUfmm0oqOnU4nLpeLgIAABg8ezLPPPktwcDC5ubnY7XbWr19Py5YtmTBhArt27bLoE1RPSm5ERCzwww8/8N5775GammrJ808/3b2xZr9+cOgQ9O/vniYuRZWnjurYuhyXy0WzZs04cuQI8+fPp3nz5px//vns3Lmz2H3i4+NVcOxhSm5ERCzw3HPPMWDAABYvXmxZDLVru6eJx8TA33+7F/p77TXLwqnSSqrLSU9P56abbuKCCy7A5XKxceNGWrRowfjx4zH/7GiqouPK4ZPJzeDBg6lbty5Dhw4t8fzhw4dp0aIFkyZN8nJkIiL+5bTTYPlyuPpqyM2Fa66ByZMhL8/qyKqOE9XlOJ1Oli5dSnR0NB9++CEtWrTA5XKxZs0abDabio4rkU8mN+PHj2fp0qUnPH///ffTo0cPL0YkIuK/QkJg2TJ3UgMwZw5cdhn8+ae1cVUVJ6vLyc/P55JLLmH79u3cfPPNfPvtt4SEhOBwOJgwYQIZGRns0d4YHuWTyU2/fv0IDQ0t8dzPP//MDz/8wBVXXOHlqERE/FdgIDz0kHtYqnZt+OQTiIpyb7gppStrXY7NZuOZZ54pLDoODg7m4MGDzJ07l5YtWzJ16lT27dvnxcj9V7mTm3Xr1jFw4ECaNGmCzWZjxYoVxdosXLiQli1bUqNGDaKioljvwTL8ghUhRUTE84YOhc2boU0b91o4F10ECxfCPyUicoqOLTrOzc0lPz+fbt26cfjwYWbPns2ZZ57JzJkzOXDgQLHrVHRcduVObrKzs4mMjGTBggUlnl+2bBkTJkxgxowZbNmyhd69exMTE0NaWlphm6ioKCIiIoq9TjZV7q233qJ169a0bt26vGGLiEgZtWkDmza5E52jR2HcOBgxAg4ftjqyqq2kouOlS5cyYMAA3nnnHRo3bkxubi733XcfLVu2ZO7cuUWuU9FxOZhTAJjly5cXOda9e3czZsyYIsfatGljpk6dWq57r1mzxgwZMqTIsalTp5qmTZuaFi1amPr165s6deqYWbNmlXj9kSNHTEZGRuErPT3dACYjI6NccYiIlNeUKVMMYO6+++4Ttpk6daoBzF133eXFyMonP9+YOXOMsduNAWM6djTml1+sjqpqcjqdBjBOp/OEx/Pz8811111nAAMYh8Nxwuuqk4yMjHL/fnt0heLc3FySk5OZOnVqkePR0dFs2LDhlO+fkJBQOCS1ZMkSvv3228KddktqO2vWrFN+pohIZbjqqqto2rQpkZGRVodyQjYb3H23u/bm2mvdWzVERcFLL4HKHsunLDuQ22w2Xn75Zdq2bUtsbCwPPvggubm5OJ1O+vbty3//+19GjRpFUFCQFR+hajmVbIrjem527txpAPPZZ58VaXf//feb1q1bl/m+0dHRpkGDBua0004zZ555ptm0aVOxNs8++2ypfytSz42IWKUsPTdVze+/G9Ozp7sHx2Yz5oEH3D07UjmCg4MNYIKDg01+fr7p0aOHAUyrVq3M888/b/Ly8qwO0Wsq0nNTKbOlbDbb8QlUsWOlWb16NXv37uXw4cP8/vvvdOvWrVibkSNHMmfOnBPeIyQkhDp16vD8889z/vnnc8kll5T9A4iISBFnnumeQTVmjLu4ePp0uO46yM62OjL/c3zRsdPp5MYbbyQ8PJxt27YxfPhwIiMjWb58eeFigAXXqejYzaPJTYMGDbDb7ezevbvI8T179hAeHu7JR5XZuHHjSE1NZfPmzZY8X0SqL1PKFKNt27bxySef8Msvv3gxolMTHAyLFsETT0BQELz6KlxwAezYYXVk/qOkouO4uDgOHjzIr7/+yqWXXgrAd999x9VXX02PHj349NNPVXR8HI8mN8HBwURFRZF03AYlSUlJ9OrVy5OPEhGp0p588kn69evHwoULrQ6l3G691b3x5hlnwNdfQ9eu7l4dOTWlrXTscDiYN28eSUlJTJ8+HYCgoCA2b97M/PnztdLxccpdUJyVlVXkbxrbt28nJSWFevXq0bx5cyZOnMjw4cPp2rUrPXv25MknnyQtLY0xY8Z4NPCySkxMJDExEZfLZcnzRaT6Kc8wfFV14YXw5ZcweDAkJ7t3Fn/kEbj9dnchspRfWYqOwb1Kf40aNQp7al5//fXC65KSkjjrrLM499xzvR6/TylvYc+aNWsKp6kd+xoxYkRhm8TERNOiRQsTHBxsunTpYtauXVvex3hcRQqSREQq4p577jGAmThx4knb+PJU8LI4fNiYf/3LXWgMxtx8szFHjlgdVfVwbNGxMcZkZ2ebRo0aGbvdbv7zn/+YnTt3WhyhZ3iloLhv374YY4q9lixZUthm7Nix7Nixg5ycHJKTk+nTp8+pZ2EiIuJzTjsNli6FuXMhIACefda9qvHvv1sdmX87vug4Pj6egwcP0qVLF1wuF0888QTnnHMO06ZNK7LacXUpOvbJvaU8KTExkXbt2pU440pERE6dzQYTJ8KqVVC3LmzcCO3aQWIiqCLA80oqOnY4HDz99NO89957/Pvf/wbg77//5sEHH+Tss89m7ty5xMbGVpuiY48u4ueLxo0bx7hx48jMzCQsLMzqcERE/NZll7n3pbrxRneCc/vt8Pzz8OST0LGj1dH5hxMVHQOFi9o+/fTTtGjRgtjYWBo2bMjevXuZNGkSQLUpOvb7nhsREfGes8+Gzz5z99rUqeNOcrp0gXvu0d5UnlBa0bHT6SwsOi5IgPbu3UtgoLsf49jrqtISBBVhM8a/93o9drbUTz/9REZGBnXq1LE6LBHxY1OnTmX27NlMnDixcPPD461du5ZNmzbRtWtX+vXr5+UIvWPXLhg/Ht54w/3+rLPc6+RcfrmlYVUrISEhhbU5OTk5AOzYsYPzzjuPvn37Mnv2bDp16mRtkCdRMPJSnt9vv++50SJ+IuKLLrroIiZPnuy3iQ1Akybw+uvw9tvQrJl7sb+YGLj+ejhurVepBCUVHQNs2LCBvLw8PvjgA7p06cLw4cP57bffilxX1YuO/T65ERERaw0cCKmpcNdd7hlVr7wCbdvC4sXuCeTieScqOo6Pj+eGG25g/PjxgHsV7RdeeIHzzjuPe+65hxkzZvhF0bHfFxSLiPii9PR0du/eTXh4OM2bN7c6nEpXuzbMm+cuNr71VvjqK7jlFli2DJ56yj1kJZ5RlqLjRx55hHr16uFwOGjZsiXbt2/noYceAiA2NrbKFx37fc+NpoKLiC9KTEyke/fuzJ8/3+pQvCoqyl1k/PDDUKMGfPghdOgACxdCfr7V0fmHshYdF7zfvn17YdHxRRddVGRIqqqW5fp9cqOaGxER3xIYCJMmufeluvBCyMqCcePg4ovBzyfxeEVcXNwJe15mzpxZJHmZOXMmwcHB5OXlERwczKpVqwrPrV+/nj59+rBx48bKDtnj/D65ERER39S6NaxdC489BjVruv+9Y0f3HlVa/M87ji86fvjhhwvPDR8+nE8//ZTzzz+f6667jm3bthW5zpeLjpXciIiIZQIC4I47YOtWd8/N33+7Vzvu3Rt++MHq6PxbaUXHAMOGDStsu2zZMtq2bcuUKVOqRNGxCopFRMRyrVq562+eeso9ZPX559CpE8yaBXff7R7KEs8pS9Hxww8/TJ06dXA4HLRq1Ypt27YV9uwMGTLEp4uO/b7nRgXFIiJVg83mnkn17bfQvz/k5MDUqXDBBfD991ZH51/KW3S8bdu2wqJjgP/7v//zarzl5ffJjQqKRUSqlubN4f334ZlnICwMNm2Czp1h9mzIy7M6Ov9wKkXHL774IjfeeGPh+Q8++ICtW7dWdsjl4vfJjYiIL7r00kuJjY0lOjra6lB8ks0GN9/s7sWJiSnai5OaanV01cvxRce//vorAQHu9GH69OkMGTKETp06ceutt7J3794i11lVdKzkRkTEApdeeilxcXH079/f6lB8WtOm8N578OyzRXtxHnxQvTjecLKi4/z8fLKyssjPz+fFF1/k6NGjRa6zquhYJVoiIuLTbDYYORIuu8xdk7NyJUybBm++CUuWQLt2Vkfon8pSdPzggw9Sq1YtHA4H/fv3p0mTJiVe521KbkRELPDnn3+yb98+6tevT6NGjawOp0o480x4911YuhTuvBM2b3b34jid7hlWPjwzuUoqrei44Pyx7x0OR+Eu5FYmNlANhqU0W0pErFLa0vXz5s0jIiKiyKJpcnI2G4wYAd99B1deCbm57lqciy6CX3+1Ojr/UpGi44LaHKunift9cqPZUiLibTabzeoQ/N6ZZ8I777hnVIWGwmefQWQk/Pe/2mncCscXHRfU5FjF75MbERHxTwUzqrZuhb59ITsbbrsNrrgCdu60Orrq42RFx1ZQciMiIlVaixbw0UfuPalq1IBVq9w7jb/yitWR+b8TFR1bneAouRERkSovIAAmTICvvoKoKDhwAK6/Hq67Dvbvtzo6/1XWlY69TbOlRETEb7Rt696X6oEHID4eli1z7zb+zDPuxQDFs0pbpE+zpURERDwkKAhiY+GLL6BNG9i9212Hc/vtcPiw1dGJNyi5ERGxQJ8+fZg0aRJ9+/a1OhS/1bWre5hq/Hj3+8RE95DVV19ZG5dUPiU3IiIWuPLKK3n44YcZOHCg1aH4tdNOg/nzYfVqaNwYfvgBevRwb99gUTmIeIHfJzdaxE9ERKKj3VPGr77avSfVtGnQrx/s2GF1ZFIZ/D650SJ+IuKL/vrrL3bs2MFff/1ldSjVRv368Prr7k04a9eG9euhY0d4/nkt/Odv/D65ERHxRbNnz6Zly5Y88MADVodSrRRswvn119CrFxw6BDfd5J4yrjzTfyi5ERGRaqdVK/cU8fvug8BAePVVd7Hx119bHZl4gpIbERGplgIDYcYM2LDBnezs2AE9e2plY3+g5EZERKq1bt1g82bo3x/+/tu9svHkye7CY6malNyIiEi1V68evPceTJ3qfj9njntFY23dUDUpuREREQHsdkhIcNff1KwJH37o7tVRHU7Vo+RGRETkGMOGubduaNUKtm9XHU5V5JPJzeDBg6lbty5Dhw4tdi4wMJBOnTrRqVMnRo8ebUF0IiKnrmfPnowdO5ZevXpZHYqUoEOH4nU4U6aoDqeq8MldwcePH8+///1vnnvuuWLnTj/9dFJSUrwflIiIBw0aNIhBgwZZHYaUoqAO59573ds1PPwwbNzoXgSwVSuro5PS+GTPTb9+/QgNDbU6DBERqeaOrcOpVQvWrXOvavzf/2pVY19W7uRm3bp1DBw4kCZNmmCz2VixYkWxNgsXLqRly5bUqFGDqKgo1q9f74lYAcjMzCQqKooLL7yQtWvXeuy+IiLelJWVxd69e8nKyrI6FCmDYcPgm2+gTx/IzobbboPLL4f0dKsjk5KUO7nJzs4mMjKSBQsWlHh+2bJlTJgwgRkzZrBlyxZ69+5NTEwMaWlphW2ioqKIiIgo9tq1a9dJn79jxw6Sk5P573//y0033URmZmZ5P4KIiFeYUv5q73Q6OeOMM4iLi/NeQHJKWrWCNWvgkUegRg344AOIiIAlS9SL42vKXXMTExNDTEzMCc/PmzePUaNGFRb7Pvroo6xevZpFixaRkJAAQHJycgXDhSZNmgAQERFBu3bt+Omnn+jatWuxdjk5OeTk5BS+VxIkIt5is9msDkEqSUAATJjgXgNn5Ej3rKqbb4Y334Qnn4RGjayOUMDDNTe5ubkkJycTHR1d5Hh0dDQbNmw45fsfOHCgMGH5/fffSU1NpdUJqroSEhIICwsrfDVr1uyUny8iIgJw3nnuXcUTEiA4GN55B9q3h2XLrI5MwMPJzb59+3C5XISHhxc5Hh4ezu7du8t8n/79+zNs2DBWrlxJ06ZN2bx5MwDff/89Xbt2JTIykgEDBjB//nzq1atX4j2mTZtGRkZG4StdA6MiIuJBgYHuFY2Tk6FzZ/eu4tddB9deqx3GrVYpU8GP75I1xpSrm3b16tUlHu/Vqxdbt24t0z1CQkIICQkhMTGRxMREXC5XmZ8vIiJSVhER7ini99/v3mX81Vfhs89g6VK4+GKro6uePNpz06BBA+x2e7Femj179hTrzfGWcePGkZqaWtj7IyIi4mlBQRAXB59/DueeCzt3wiWXwKRJcEz5p3iJR5Ob4OBgoqKiSEpKKnI8KSlJq3CKiIjf69YNtmyBW291v587F7p3h+++szau6qbcyU1WVhYpKSmFqwRv376dlJSUwqneEydOZPHixTzzzDN8//333HXXXaSlpTFmzBiPBl5WiYmJtGvXjm7dulnyfBGRkkRFRTFixIgSZ3tK1VarFjzxBLz1FjRo4F4fJyoKHnsM8vOtjq6aMOW0Zs0aAxR7jRgxorBNYmKiadGihQkODjZdunQxa9euLe9jPC4jI8MAJiMjw+pQRMTPTZs2zQBmwoQJVociFvvjD2NiYoxxr4RjTP/+xuzaZXVUVUtFfr/LXVDct2/fUhemAhg7dixjx44td6IlIiLiTxo1cu9PtXChu/5m9Wr3ppyLF4O2Fqs8Prm3lCdpWEpEfFFubi7Z2dnk5uZaHYpUMpsNxo1zTxnv1An274fBg2HMGDh82Oro/JPfJzeaLSUivujee++ldu3azJgxw+pQxEvatXNPGZ882f3+iSega1f4+mtr4/JHfp/ciIiI+IrgYHjoIUhKgsaN4fvv3bOp5s/X/lSe5PfJjYalRETE11x6qbvHZuBAyM1171d15ZWwZ4/VkfkHv09uNCwlIiK+qGFD93TxxET3LuPvvw8dO7qLjuXU+H1yIyIi4qtsNhg7FjZvdm/j8OefcPnlcPfdWtn4VCi5ERERsVhEBGzaBLff7n4/bx6cfz78+KO1cVVVfp/cqOZGRESqgtNOg8cfh3feca9snJLiXtl4yRIVG5eX3yc3qrkREV/UoUMHhg4dSseOHa0ORXzMgAHuYuOLL4bsbLj5ZrjxRsjMtDqyqsPvkxsREV80fPhwXnvtNYYPH251KOKDmjSBDz6ABx4Aux1efhk6d3YPXcnJKbkRERHxQXY7TJsG69dDixawbRtccIF7nRxtwFk6JTciIiI+rGdPd/3NNddAXh7cc497RtXu3VZH5rv8PrlRQbGI+KLJkydjs9mYXLAWv0gpTj8dXnnFveHmaae5VziOjIRVq6yOzDf5fXKjgmIRsYrRFBfxIJsNRo1yb8DZsaN7NeOYGHdPztGjVkfnW/w+uRER8TabzWZ1COLH2rZ1b8A5bpz7/UMPQZ8+sGOHpWH5FCU3IiIiVUyNGrBgAbzxBoSFwRdfuGdTvfmm1ZH5BiU3IiIiVdTVV7uLjXv0gIMHYcgQ9yrHR45YHZm1lNyIiIhUYWed5Z4uPmWK+31ionvrhp9+sjQsS/l9cqPZUiIi4u+CgmD2bPfO4g0bulc47tIFnn/e6sis4ffJjWZLiYgvatOmDVdccQVt2rSxOhTxI5df7h6m6tvXvXXDTTe5t2/IzrY6Mu/y++RGRMQXjRo1ivfee49Ro0ZZHYr4mSZN4MMPYdYsCAhwb7zZvTukplodmfcouREREfEzdjs4HPDRR9CokTux6dYNli61OjLvUHIjIiLip/r2dQ9TXXopHD4MI0a4FwI8fNjqyCqXkhsREQtMnz6dWrVqMWPGDKtDET8XHu7epsHpdA9TPfOMe5jq+++tjqzyKLkREbHA0aNHOXz4MEe1br54gd0OM2e6a3HCw+G776BrV/+dTaXkRkREpJro1889THXJJe6hqZtugtGj/W+YSsmNiIhINdKoEaxe7Z5NZbPB00/736J/fp/caBE/ERGRogpmUxUMU23d6h6meu01qyPzDL9PbrSIn4iISMkuvhi2bIGLLoJDh+Caa+DOOyE31+rITo3fJzciIiJyYo0bu3twpk51v3/sMejTB9LSrI3rVCi5ERGxQKtWrbjoooto1aqV1aGIEBgICQnwzjtw+umwcSN07uyeQl4VKbkREbHAbbfdxieffMKYMWOsDkWk0IAB8NVXEBUFf/0FV1zhnkLuclkdWfkouREREZFCLVvCZ5/B2LFgDNx3H/TvD3v2WB1Z2Sm5ERERkSJCQiAxEV58EWrVcu9R1bkzfP651ZGVjZIbERELxMbG0rBhQ2bNmmV1KCIndMMNsHkztG0Lu3a5Z1UtWODu0fFlSm5ERCqJKeUXIDs7m3379pGdne3FiETKr21b2LTJPU386FG44w4YPhx8+avrk8nN4MGDqVu3LkOHDi12bvv27fTr14927drRoUMH/cEgIj7HZrNZHYKIR9WuDa+8Ao884l4A8MUX3asa//yz1ZGVzCeTm/Hjx7N06dISz40cORKn00lqaipr164lJCTEy9GJiIhUPzYbTJgAa9a4t3D49lv3qsZvvWV1ZMX5ZHLTr18/QkNDix3/7rvvCAoKonfv3gDUq1ePwMBAb4cnIiJSbfXu7Z4ufuGFkJkJgwbB9Om+NV283MnNunXrGDhwIE2aNMFms7FixYpibRYuXEjLli2pUaMGUVFRrF+/3hOx8vPPP1O7dm2uuuoqunTpwgMPPOCR+4qIiEjZNW4MH3/s7skB9wKAl18Oe/daGlahcic32dnZREZGsmDBghLPL1u2jAkTJjBjxgy2bNlC7969iYmJIe2YdZyjoqKIiIgo9tq1a1epzz569Cjr168nMTGRzz//nKSkJJKSksr7EUREROQUBQW5a3Befhlq1nRv4RAV5RvbNpR7TCcmJoaYmJgTnp83bx6jRo1i9OjRADz66KOsXr2aRYsWkZCQAEBycnKFgm3atCndunWjWbNmAFxxxRWkpKRw2WWXFWubk5NDTk5O4fvMzMwKPVNEpDI0bdqUrl270rRpU6tDETkl110HHTrA1VfDOeeAL3ylPVpzk5ubS3JyMtHR0UWOR0dHs2HDhlO+f7du3fjzzz85cOAA+fn5rFu3jrZt25bYNiEhgbCwsMJXQUIkIuILJkyYwObNmxk/frzVoYicsvbt3evhvPACBPhANa9HQ9i3bx8ul4vw8PAix8PDw9m9e3eZ79O/f3+GDRvGypUradq0KZs3bwYgMDCQBx54gD59+tCxY0fOPfdcBgwYUOI9pk2bRkZGRuErPT294h9MRERESlWnDtSta3UUbpUy1ej4NR6MMeVa92H16tUnPHeyYbECISEhmiYuIiJSDXm056ZBgwbY7fZivTR79uwp1pvjLYmJibRr145u3bpZ8nwRkZLcf//9nHXWWTz44INWhyLidzya3AQHBxMVFVVsBlNSUhK9evXy5KPKbNy4caSmphYObYmI+IIDBw7w22+/ceDAAatDEfE75R6WysrK4pdffil8v337dlJSUqhXrx7Nmzdn4sSJDB8+nK5du9KzZ0+efPJJ0tLSGDNmjEcDL6vExEQSExNx+dLqQiIiIlJpyp3cfPnll/Tr16/w/cSJEwEYMWIES5Ys4dprr2X//v04nU7++OMPIiIiWLlyJS1atPBc1OUwbtw4xo0bR2ZmJmFhYZbEICIiIt5T7uSmb9++pe50CzB27FjGjh1b4aBEREREKsoHZqNXLhUUi4iIVC9+n9yooFhERKR68fvkRkTEF4WHh9OuXTvLlskQ8WeVsoifL9FsKRHxRZMnT2by5MlWhyHil/y+50bDUiIiItWL3yc3IiIiUr0ouRERscDDDz9M+/btmTdvntWhiPidalNzk5eXB0BmZqbFEYmIv8vJySn854n+zPntt99ITU0lLS1Nfy6JlKLg/x8nW2PvWDZTntZV2O+//06zZs2sDkNEREQqID09naZNm5apbbVJbvLz89m1axehoaHYbDaP379bt26WFC1XxnM9cc+K3KM815S1bVnaldYmMzOTZs2akZ6eTp06dcoUmy/T9/TU76HvaeXT9/TU7+FP31NjDIcOHaJJkyYEBJStmsbvh6UKBAQElDnjqwi73W7JHyqV8VxP3LMi9yjPNWVtW5Z2ZWlTp04dv/jR0Pf01O+h72nl0/f01O/hb9/T8u4NqYJiDxk3bpzfPNcT96zIPcpzTVnblqWdVf/bWUHf01O/h76nlU/f01O/R3X/nlabYSmRiijYTT4jI8Mv/kYs/knfU6kKvPk9Vc+NSClCQkKIjY0lJCTE6lBETkjfU6kKvPk9Vc+NiIiI+BX13IiIiIhfUXIjIiIifkXJjYiIiPgVJTciIiLiV5TciHjQ4cOHadGiBZMmTbI6FJFiDh06RLdu3ejUqRMdOnTgqaeesjokkWLS09Pp27cv7dq1o2PHjrz22mvlvodmS4l40IwZM/j5559p3rw5c+bMsTockSJcLhc5OTnUrFmTw4cPExERwebNm6lfv77VoYkU+uOPP/jzzz/p1KkTe/bsoUuXLvz444/UqlWrzPdQz42Ih/z888/88MMPXHHFFVaHIlIiu91OzZo1AThy5Agul6tcOy2LeEPjxo3p1KkTAGeccQb16tXjr7/+Ktc9lNxItbBu3ToGDhxIkyZNsNlsrFixolibhQsX0rJlS2rUqEFUVBTr168v1zMmTZpEQkKChyKW6sgb39ODBw8SGRlJ06ZNmTJlCg0aNPBQ9FJdeON7WuDLL78kPz+fZs2ales6JTdSLWRnZxMZGcmCBQtKPL9s2TImTJjAjBkz2LJlC7179yYmJoa0tLTCNlFRUURERBR77dq1i7feeovWrVvTunVrb30k8UOV/T0FOP300/n666/Zvn07L730En/++adXPpv4D298TwH279/PTTfdxJNPPln+II1INQOY5cuXFznWvXt3M2bMmCLH2rRpY6ZOnVqme06dOtU0bdrUtGjRwtSvX9/UqVPHzJo1y1MhSzVUGd/T440ZM8a8+uqrFQ1RpNK+p0eOHDG9e/c2S5curVBc6rmRai83N5fk5GSio6OLHI+OjmbDhg1lukdCQgLp6ens2LGDOXPmcMstt+BwOCojXKmmPPE9/fPPP8nMzATcmxiuW7eO8847z+OxSvXlie+pMYaRI0dy8cUXM3z48ArFEVihq0T8yL59+3C5XISHhxc5Hh4ezu7duy2KSqQoT3xPf//9d0aNGoUxBmMMt99+Ox07dqyMcKWa8sT39LPPPmPZsmV07NixsJ7n+eefp0OHDmWOQ8mNyD9sNluR98aYYsfKYuTIkR6KSKS4U/meRkVFkZKSUglRiRR1Kt/TCy+8kPz8/FN6voalpNpr0KABdru92N8q9uzZU+xvHyJW0fdUqgJf+Z4quZFqLzg4mKioKJKSkoocT0pKolevXhZFJVKUvqdSFfjK91TDUlItZGVl8csvvxS+3759OykpKdSrV4/mzZszceJEhg8fTteuXenZsydPPvkkaWlpjBkzxsKopbrR91SqgirxPa3QHCuRKmbNmjUGKPYaMWJEYZvExETTokULExwcbLp06WLWrl1rXcBSLel7KlVBVfieam8pERER8SuquRERERG/ouRGRERE/IqSGxEREfErSm5ERETEryi5EREREb+i5EZERET8ipIbERER8StKbkRERMSvKLkRERERv6LkRkRERPyKkhsRERHxK0puRERExK8ouRERERG/8v9bd4+hniKr4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(times, -data, 'k')\n",
    "plt.loglog(times, data, '--k')\n",
    "plt.loglog(times, -data_emip, 'kx')\n",
    "plt.loglog(times, data_emip, 'kx')\n",
    "plt.loglog(times, -data_em, 'b')"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}