{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D TDEM inversion using the cylindrical mesh\n",
    "\n",
    "In this example, we perform a 1D time domain electromagnetic inverion using a cylindrical mesh. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SimPEG import (Mesh, Maps, Utils, DataMisfit, Regularization,\n",
    "                    Optimization, Inversion, InvProblem, Directives)\n",
    "import numpy as np\n",
    "from SimPEG.EM import FDEM, TDEM, mu_0\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "try:\n",
    "    from pymatsolver import Pardiso as Solver\n",
    "except ImportError:\n",
    "    from SimPEG import SolverLU as Solver\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up cylindrically symmeric mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cs, ncx, ncz, npad = 10., 15, 25, 13  # padded cyl mesh\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')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x104a73940>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAG+hJREFUeJzt3X+0ZWV52PHvUyhoIOkwIlcE0rl0jV2LuFaMXBkmabIu/uCXLtFUK6w04o90aAPpj6ysOpSuamJZxdQ2ia0dmUQqtupANdQRqVQId9muEmQmVQR1ZIRRrhCRwDEZpSZcn/5x9oHjmXPvOWffd9997uX7Weuse/a733efZ7937nlm7/fde0dmIklSSX+t7QAkSRuPyUWSVJzJRZJUnMlFklScyUWSVJzJRZJUnMlFklScyUWSVFyrySUirouIRyPi3r6yd0XEtyLiC9Xrwr51V0bEwYg4EBHn9ZWfX5UdjIida70fkqQfFW1eoR8RvwAcBj6cmS+uyt4FHM7M9w7UPQP4GHAW8ELgNuBF1eqvAa8CFoG7gUsy88srffaJJ56YW7ZsqRX39773PY477rhabdeasTbDWJuxXmJdL3FC+Vj379//WGY+f1S9o4t9Yg2Z+bmI2DJm9YuAPZn5A+DBiDhIN9EAHMzMBwAiYk9Vd8XksmXLFvbt21cr7oWFBebn52u1XWvG2gxjbcZ6iXW9xAnlY42Ib4xTb1rHXK6IiHuq02YnVGWnAA/11VmsypYrlyS1pNXTYgDVkcvNfafFZoDHgATeDZycmW+LiPcDd2bmf63qfRC4hW6CPC8zf6Uq/2XgrMz8tSGftQPYATAzM3Pmnj17asV8+PBhjj/++Fpt15qxNsNYm7FeYl0vcUL5WM8555z9mTk3ql6rp8WGycxv995HxO8DN1eLi8BpfVVPBR6u3i9XPrjt3cBugLm5uax7qPhsPiRukrE2w1jLWy9xQnuxTt1psYg4uW/x9UBvJtle4OKIODYiZoGtwOfpDuBvjYjZiDgGuLiqK0lqSatHLhHxMWAeODEiFoF3AvMR8RK6p8UOAZcBZOZ9EXEj3YH6p4DLM3Op2s4VwK3AUcB1mXnfGu+KJKlP27PFLhlS/MEV6l8NXD2k/Ba64y+SpCkwdafFJEnrX+uzxdoyNzeXda5z2bLz0yuu3za7mbsefLxWTJO07a+7bXYzwMjlnuXKx11ft+6oNp1Oh02bNq16O6s1znYnjbXU59ax2libimuYJvq1p+R+NBlnaf2x9vrg0DWvrr29iBhrtphHLpKk4qZuKvK02za7ecX/tdxw2XbedO2dtbY9Sdv+ujdcth1g6HKn03l6uWew3qBR6+vWHdWmO2Vy9dtZrXG2O2mspT63jtXG2lRcwzTRrz0l96PJOEvrj7Xud1MdHrlIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoqLzGw7hlbMzc3lvn37Jm63ZeenV1y/bXYzdz34eK2YJmnbX3fb7GaAkcs9y5WPu75u3VFtOp0OmzZtWvV2Vmuc7U4aa6nPrWO1sTYV1zBN9GtPyf1oMs7S+mPt9cGha15de3sRsT8z50bV88hFklTc0W0HsN5sm9284v9abrhsO2+69s5a256kbX/dGy7bDjB0udPpPL3cM1hv0Kj1deuOarOwsMD8/Oq3s1rjbHfSWEt9bh2rjbWpuIZpol97Su5Hk3GW1h9r3e+mOjxykSQV12pyiYjrIuLRiLi3r2xzRHw2Iu6vfp5QlUdEvC8iDkbEPRHx0r42l1b174+IS9vYF0nSM9o+cvkQcP5A2U7g9szcCtxeLQNcAGytXjuAXdBNRsA7gW3AWcA7ewlJktSOVpNLZn4OGJwedRFwffX+euB1feUfzq4/BjZFxMnAecBnM/PxzHwC+CxHJixJ0hqaxgH9mcx8BCAzH4mIk6ryU4CH+uotVmXLlR8hInbQPephZmaGhYWFiYPrdJ5kaWmJTqczdP3CwgKdzpMTb3fStv11e/sxbHlpaemI/RysN2jU+rp1R7U5fPhwke2s1jjbnTTWUp9bx2pjbSquYZro156S+9FknKX1x7qWv8tpTC7LiSFluUL5kYWZu4Hd0L3OZX5+fuIgdh24c8XZYvPz29l1oN6MjEna9tftzQQZttzpdBjcz8F6g0atr1t3VJvurJb5IS3KfXap7U4aa6nPrWO1sTYV1zBN9GtPyf1oMs7S+mNdy99l22Muw3y7Ot1F9fPRqnwROK2v3qnAwyuUS5JaMo3JZS/Qm/F1KfDJvvI3V7PGzga+W50+uxU4NyJOqAbyz63KJEktafW0WER8DJgHToyIRbqzvq4BboyItwPfBN5YVb8FuBA4CHwfeCtAZj4eEe8G7q7q/VZm1rv/iiSpiFaTS2ZessyqVwypm8Dly2znOuC6gqFJklZhGk+LSZLWOZOLJKk4k4skqTif5zIhn+dSv+6oNj7Pxee5+DyX8nyeiyRpw1hPV+hPBZ/nUr/uqDY+z8Xnufg8l/J8noskacMwuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkorzxpUT8saV9euOauONK71xpTeuLM8bV0qSNgxvXDkhb1xZv+6oNt640htXeuPK8rxxpSRpwzC5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSipvaK/Qj4hDwF8AS8FRmzkXEZuAGYAtwCPh7mflERATwe8CFwPeBt2Tmn6y0fa/Q9wr91WzXK/Sb4RX65XmF/nDnZOZL+nZkJ3B7Zm4Fbq+WAS4AtlavHcCuNY9UkvS09XaF/kXAfPX+emABeEdV/uHsHob9cURsioiTM/OR0gF4hX79uqPaeIW+V+h7hX55XqF/pAT+Z0Tsj4gdVdlML2FUP0+qyk8BHupru1iVSZJaMM1HLj+XmQ9HxEnAZyPiqyvUjSFlRwwmVUlqB8DMzAwLCwsTB9XpPMnS0hKdTmfo+oWFBTqdJyfe7qRt++v29mPY8tLS0hH7OVhv0Kj1deuOanP48OEi21mtcbY7aaylPreO1cbaVFzDNNGvPSX3o8k4S+uPdS1/l1ObXDLz4ernoxFxE3AW8O3e6a6IOBl4tKq+CJzW1/xU4OEh29wN7IbugP78/PzEce06cOeKp8Xm57ez60C9Q89J2vbX7R3yDlvudDoM7udgvUGj1tetO6pN9/B9fkiLcp9daruTxlrqc+tYbaxNxTVME/3aU3I/moyztP5Y1/J3OZWnxSLiuIj48d574FzgXmAvcGlV7VLgk9X7vcCbo+ts4LtNjLdIksYzlVORI+J04KZq8Wjgo5l5dUQ8D7gR+Engm8AbM/PxairyfwTOpzsV+a2ZueI8Y6ciOxV5Ndt1KnIznIpcXltTkafytFhmPgD89JDyPwNeMaQ8gcvXIDRJ0himMrlMM6ci1687qo1TkZ2K7FTk8pyKLEnaMEwukqTiTC6SpOKmcrbYWnC2mLPFVrNdZ4s1w9li5XnjSknShuFssQk5W6x+3VFtnC3mbDFni5XnbDFJ0obhmMuEHHOpX3dUG8dcHHNxzKU8x1wkSRuGYy4Tcsylft1RbRxzcczFMZfyHHORJG0YjrlMyDGX+nVHtXHMxTEXx1zKc8xFkrRhOOYyIcdc6tcd1cYxF8dcHHMpr60xF0+LTcjTYvXrjmrjaTFPi3larDxPi0mSNgyTiySpOE+LTWjUaTFJmnaeFpMkrUsmF0lScSYXSVJxJhdJUnEmF0lScSYXSVJxJhdJUnEbJrlExPkRcSAiDkbEzrbjkaRnsw2RXCLiKOD9wAXAGcAlEXFGu1FJ0rPXhkguwFnAwcx8IDP/EtgDXNRyTJL0rDVWcomI2yPiwoGy3c2EVMspwEN9y4tVmSSpBeM+z2UWeEdEvCwzf7MqG3lvmTUUQ8qOuGlaROwAdgDMzMywsLDQcFiSNH3W4rtv3OTSAV4BvC8iPgX8/eZCqmUROK1v+VTg4cFKmbkb2A3dG1fOz89P/kmf8caVkta3Wt99Exp3zCUy86nM/FXgE8D/Bk5qLqyJ3Q1sjYjZiDgGuBjY23JMkvSsNe6Rywd6bzLzQxHxJeDyZkKaXGY+FRFXALcCRwHXZeZ9LYclSc9aYyWXzLx2YHk/8LZGIqopM28Bbmk7DknSxpmKLEmaIiYXSVJxJhdJUnEmF0lScSYXSVJxJhdJUnHjXueiMW2b3cxdDz7eeNv+uttmNwOMXO5Zrnzc9XXrjmrT6XTYtGnTqrezWuNsd9JYS31uHauNtam4hmmiX3tK7keTcZbWH2vd76Y6PHKRJBVncpEkFWdykSQVZ3KRJBVncpEkFWdykSQV51TkCW2b3bziNMQbLtvOm669s9a2J2nbX/eGy7YDDF3udDpPL/cM1hs0an3duqPaLCwsMD+/+u2s1jjbnTTWUp9bx2pjbSquYZro156S+9FknKX1x1r3u6kOj1wkScVF5hGPmn9WmJuby3379k3cbsvOlR9z7EWUo3kRZf3PrcOLKLu8iPKZPjh0zatrby8i9mfm3Kh6HrlIkopzzGVCjrnUrzuqjWMujrk45lKeYy6SpA3DMZcJOeZSv+6oNo65OObimEt5jrlIkjYMx1wm5JhL/bqj2jjm4piLYy7lOeYiSdowTC6SpOJMLpKk4pwtNiFni9WvO6qNs8WcLeZssfKcLSZJ2jCmbrZYRLwL+AfAd6qif5GZt1TrrgTeDiwB/zgzb63Kzwd+DzgK+IPMvKap+JwtVr/uqDbOFnO2mLPFymtrttjUJZfK72Tme/sLIuIM4GLgp4AXArdFxIuq1e8HXgUsAndHxN7M/PJaBixJesa0JpdhLgL2ZOYPgAcj4iBwVrXuYGY+ABARe6q6JhdJasm0jrlcERH3RMR1EXFCVXYK8FBfncWqbLlySVJLWjlyiYjbgBcMWXUVsAt4N5DVz38HvA2IIfWT4Qly6BS4iNgB7ACYmZlhYWFh0tDpdJ5kaWmJTqczdP3CwgKdzpMTb3fStv11e/sxbHlpaemI/RysN2jU+rp1R7U5fPhwke2s1jjbnTTWUp9bx2pjbSquYZro156S+9FknKX1x7qWv8tWkktmvnKcehHx+8DN1eIicFrf6lOBh6v3y5UPfu5uYDd0pyLPz8+PH3TlLZ/5NN0898Oh63cdOJYDT3xv4u1O2ra/7q4DxwIssxxPL/cM1hs0an3duqPadDpPsmnT6rezWuNsd9JYS31uHauNtam4hmmiX3tK7keTcZbWH2uvD+p8901q6k6LRcTJfYuvB+6t3u8FLo6IYyNiFtgKfB64G9gaEbMRcQzdQf+9axmzJOlHTeOA/m9HxEvonto6BFwGkJn3RcSNdAfqnwIuz8wlgIi4AriV7lTk6zLzvqaCcypy/bqj2jgV2anITkUuz6nIlcz85RXWXQ1cPaT8FuCWJuOSJI1v6k6LSZLWP5OLJKk4k4skqTiTiySpOJOLJKk4k4skqTgfFjYhHxZWv+6oNj4szIeF+bCw8nxYmCRpw5i6iyinnVfo1687qo1X6HuFvlfol9fWFfoeuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSijO5SJKKM7lIkoozuUiSivPGlRPyxpX1645q440rvXGlN64szxtXSpI2DG9cOSFvXFm/7qg23rjSG1d648ryvHGlJGnDMLlIkoozuUiSijO5SJKKayW5RMQbI+K+iPhhRMwNrLsyIg5GxIGIOK+v/Pyq7GBE7Owrn42IuyLi/oi4ISKOWct9kSQdqa0jl3uBXwQ+118YEWcAFwM/BZwP/KeIOCoijgLeD1wAnAFcUtUFeA/wO5m5FXgCePva7IIkaTmtJJfM/EpmHhiy6iJgT2b+IDMfBA4CZ1Wvg5n5QGb+JbAHuCgiAng58PGq/fXA65rfA0nSSqZtzOUU4KG+5cWqbLny5wGdzHxqoFyS1KLGLqKMiNuAFwxZdVVmfnK5ZkPKkuFJMFeov1xMO4AdADMzMywsLCxXdVmdzpMsLS3R6XSGrl9YWKDTeXLi7U7atr9ubz+GLS8tLR2xn4P1Bo1aX7fuqDaHDx8usp3VGme7k8Za6nPrWG2sTcU1TBP92lNyP5qMs7T+WNfyd9lYcsnMV9Zotgic1rd8KvBw9X5Y+WPApog4ujp66a8/LKbdwG7o3ltsfn5+4gB3HbhzxSv05+e3s+tAvatgJ2nbX7d39e2w5U6nw+B+DtYbNGp93bqj2nSvJJ4f0qLcZ5fa7qSxlvrcOlYba1NxDdNEv/aU3I8m4yytP9a1/F1O22mxvcDFEXFsRMwCW4HPA3cDW6uZYcfQHfTfm927bt4BvKFqfymw3FGRJGmNtDUV+fURsQhsBz4dEbcCZOZ9wI3Al4HPAJdn5lJ1VHIFcCvwFeDGqi7AO4Bfj4iDdMdgPri2eyNJGtTKjSsz8ybgpmXWXQ1cPaT8FuCWIeUP0J1NJkmaEj7PZUI+z6V+3VFtfJ6Lz3PxeS7l+TwXSdKG4fNcJuTzXOrXHdXG57n4PBef51Kez3ORJG0YJhdJUnEmF0lScSYXSVJxJhdJUnEmF0lScSYXSVJxJhdJUnEmF0lScSYXSVJxJhdJUnEmF0lScSYXSVJxJhdJUnEmF0lScSYXSVJxJhdJUnEmF0lScSYXSVJxkZltx9CKubm53Ldv38Tttuz89Irrt81u5q4HH68V0yRt++tum90MMHK5Z7nycdfXrTuqTafTYdOmTavezmqNs91JYy31uXWsNtam4hqmiX7tKbkfTcZZWn+svT44dM2ra28vIvZn5tyoeh65SJKKO7rtANabQ9e8moWFBebn59sOZSzdWLe3HcZYjLUZxlreeokT2ovVIxdJUnEmF0lSca0kl4h4Y0TcFxE/jIi5vvItEfFkRHyhen2gb92ZEfGliDgYEe+LiKjKN0fEZyPi/urnCW3skyTpGW0dudwL/CLwuSHrvp6ZL6le/7CvfBewA9havc6vyncCt2fmVuD2almS1KJWkktmfiUzD4xbPyJOBn4iM+/M7tzpDwOvq1ZfBFxfvb++r1yS1JJWr3OJiAXgNzJzX7W8BbgP+Brw58C/zMz/VZ06uyYzX1nV+3ngHZn5mojoZOamvm0+kZlDT41FxA66Rz/MzMycuWfPnlpxHz58mOOPP75W27VmrM0w1masl1jXS5xQPtZzzjlnrOtcyMxGXsBtdE9/Db4u6quzAMz1LR8LPK96fybwEPATwMuA2/rq/Tzwqep9Z+BznxgnvjPPPDPruuOOO2q3XWvG2gxjbcZ6iXW9xJlZPlZgX47xHdvYdS5ZHWVM2OYHwA+q9/sj4uvAi4BF4NS+qqcCD1fvvx0RJ2fmI9Xps0dXF7kkabWm6iLKiHg+8HhmLkXE6XQH7h/IzMcj4i8i4mzgLuDNwH+omu0FLgWuqX5+cpzP2r9//2MR8Y2aoZ4IPFaz7Voz1mYYazPWS6zrJU4oH+vfHKdSK2MuEfF6usnh+UAH+EJmnhcRfxf4LeApYAl4Z2Z+qmozB3wIeC7wP4Bfy8yMiOcBNwI/CXwTeGNm1ru51/jx78txzjlOAWNthrE2Y73Eul7ihPZibeXIJTNvAm4aUv4J4BPLtNkHvHhI+Z8BrygdoySpPq/QlyQVZ3KpZ3fbAUzAWJthrM1YL7GulzihpViftc9zkSQ1xyMXSVJxJpcJRcT5EXGguoFmq/cxi4jTIuKOiPhKdSPQf1KVD72ZZ3S9r4r9noh4aQsxHxUR/zcibq6WZyPirirWGyLimKr82Gr5YLV+yxrHuSkiPh4RX636d/u09mtE/LPq939vRHwsIp4zLf0aEddFxKMRcW9f2cT9GBGXVvXvj4hL1zDWf1v9G7gnIm6KiP67gVxZxXogIs7rK2/8O2JYrH3rfiMiMiJOrJbb6ddxrrT09fTV/0cBXwdOB44Bvgic0WI8JwMvrd7/ON3b5pwB/DawsyrfCbynen8h3WncAZwN3NVCzL8OfBS4uVq+Ebi4ev8B4B9V738V+ED1/mLghjWO83rgV6r3xwCbprFfgVOAB4Hn9vXnW6alX4FfAF4K3NtXNlE/ApuBB6qfJ1TvT1ijWM8Fjq7ev6cv1jOqv/9jgdnqe+GotfqOGBZrVX4acCvwDeDENvt1Tf4ANsoL2A7c2rd8JXBl23H1xfNJ4FXAAeDkquxk4ED1/lrgkr76T9dbo/hOpXvn6pcDN1f/2B/r++N9un+rP5Dt1fujq3qxRnH+RPWFHQPlU9evdJPLQ9UXxNFVv543Tf0KbBn4wp6oH4FLgGv7yn+kXpOxDqx7PfCR6v2P/O33+nUtvyOGxQp8HPhp4BDPJJdW+tXTYpPp/SH3LFZlratOb/wM3TsYzGTmIwDVz5Oqam3H/7vAPwd+WC0/j+694Z4aEs/TsVbrv1vVXwunA98B/nN1Cu8PIuI4prBfM/NbwHvpXkD8CN1+2s909mvPpP3Y9r/bnrfRPQKAKYw1Il4LfCszvziwqpVYTS6TiSFlrU+3i4jj6V58+k8z889XqjqkbE3ij4jXAI9m5v4x42mzr4+me8phV2b+DPA9Vn5OUJv9egLdx07MAi8EjgMuWCGeqfw3XFkuttZjjoir6N455CO9oiHVWos1In4MuAr4V8NWDylrPFaTy2QW6Z7T7Om/gWYrIuKv000sH8nMP6yKvx3dm3j2noXTu5lnm/H/HPDaiDgE7KF7aux3gU0R0btTRH88T8darf8bQKO39emzCCxm5l3V8sfpJptp7NdXAg9m5ncy86+APwR+luns155J+7HVv7tqoPs1wC9ldf5ohZjaivVv0f0Pxherv7FTgT+JiBe0FavJZTJ3A1urmTjH0B0Q3dtWMBERwAeBr2Tmv+9b1buZJ/zozTz3Am+uZo+cDXy3d3qiaZl5ZWaemplb6PbbH2XmLwF3AG9YJtbePryhqr8m/1vNzD8FHoqIv10VvQL4MlPYr3RPh50dET9W/XvoxTp1/dpn0n68FTg3Ik6ojtTOrcoaFxHnA+8AXpuZ3x/Yh4ur2XezdG+y+3la+o7IzC9l5kmZuaX6G1ukO9nnT2mrX5sYaNrIL7ozL75Gd0bIVS3H8nfoHsbeA3yhel1I9xz67cD91c/NVf0A3l/F/iX6nqWzxnHP88xssdPp/lEeBP4bcGxV/pxq+WC1/vQ1jvElwL6qb/873dk0U9mvwG8CX6X7vKT/QncG01T0K/AxumNBf0X3C+/tdfqR7njHwer11jWM9SDdcYne39cH+upfVcV6ALigr7zx74hhsQ6sP8QzA/qt9KtX6EuSivO0mCSpOJOLJKk4k4skqTiTiySpOJOLJKk4k4skqTiTiySpOJOLNCUi4mXV8zaeExHHRfcZLS9uOy6pDi+ilKZIRPxrulfRP5fu/c3+TcshSbWYXKQpUt2P6m7g/wE/m5lLLYck1eJpMWm6bAaOp/tk0ee0HItUm0cu0hSJiL10H0kwS/dpjVe0HJJUy9Gjq0haCxHxZuCpzPxoRBwF/J+IeHlm/lHbsUmT8shFklScYy6SpOJMLpKk4kwukqTiTC6SpOJMLpKk4kwukqTiTC6SpOJMLpKk4v4/qcLcIu3NEL4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104b362b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh.plotGrid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conductivity model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "layerz = np.r_[-200., -100.]\n",
    "layer = (mesh.vectorCCz >= layerz[0]) & (mesh.vectorCCz <= layerz[1])\n",
    "active = mesh.vectorCCz < 0.\n",
    "sig_half = 1e-2  # Half-space conductivity\n",
    "sig_air = 1e-8  # Air conductivity\n",
    "sig_layer = 5e-2  # Layer conductivity\n",
    "sigma = np.ones(mesh.nCz)*sig_air\n",
    "sigma[active] = sig_half\n",
    "sigma[layer] = sig_layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'conductivity')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOkAAAEMCAYAAAAhyZFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAE1BJREFUeJzt3X+Q3HV9x/Hn6y73Kz+ESCCE5Gj4cbbGMipzDViqIFAIzNgIFQu0GpFpagfaTtuxhWKL1uJYbeuUwapRU4MWMGWkpJoRAy2COkgyihRQ5ABDjsSEAFGS7CV3t+/+sd9Lvlx2L3u53dvP3r0eMzu3+/l+v5/v6zb32u/u9za3igjMLF0tjQ5gZmNzSc0S55KaJc4lNUucS2qWOJfULHEu6RQn6RxJ/XWae7ekk6tY728kfaEeGaaDGY0OYM1B0v3AVyLiQNkiYnY120bEx3LzLAaeBdoiYqi2KacmH0nNEueSNpikbklfk/SCpBcl3SKpRdKHJG2WtEPSrZKOytZfLCkkrZD0nKSdkm7Izdcl6UuSXpb0BPAbo/YXkk7N3f6SpH/I3V4u6RFJv5T0tKRlkm4C3grckj3FvSU/l6QzJf1cUmtunkskPZpd/7Ckr2SLHsi+7srmOlvSS5JOy217nKSCpGNrdDc3NZe0gbIf6q8Dm4HFwELgDuB92eXtwMnAbOCWUZv/FvCrwHnA30l6fTZ+I3BKdrkQWDGOPEuBW4EPAkcDbwN+FhE3AA8C10bE7Ii4Nr9dRDwE7AHOzQ1fCdxWZjdvy74enc317ex7/oPcOlcA90bEC9Vmn8pc0sZaCpwAfDAi9kTEQER8B/h94F8i4pmI2A1cD1wuKX8O4SMRUYiIHwE/At6Yjb8buCkiXoqILcDN48hzNbA6IjZERDEino+In1S57e2UyoWkOcDF2Vg11gBXShr5eXwP8OVx5J7SXNLG6gY2lzmBcgKlo+uIzZRO8s3Pjf08d30vpaPtyLZbRm07njxPj2P9vNuASyV1AJcCP4iIqvYdEd+ndCQ+W9KvAacC644wx5TjkjbWFuDEUUdIgK3Ar+RunwgMAdurmHMbpbLlt83bC8zM3T5+VJ5TKsw75n+XiognKD0gXETlp7pjzbOG0lPe9wB3RsTAWPubTlzSxnqYUqk+LmmWpE5JZ1F6mvjnkk6SNBv4GPDVKn9lsRa4XtJcSYuAPxm1/BFKTy1bJS0Dzs4t+yJwlaTzspNXC7MjG5QeIA73O9HbgD+l9LrzPyus8wJQLDPXl4FLKBX11sPsZ1pxSRsoIoaBd1B6evcc0A/8HrCa0g/tA5R+pzjAoWWr5COUjmjPAt/i0Nd2f5btcxel177/lcvzMHAV8CngF8C3OXhE/1fgXdlZ40qvc28HzgH+JyJ2Vvie9wI3Ad+VtEvSmdl4P/ADSkfaB6v8XqcF+T99WyokrQa2RsSHGp0lJX7HkSUheyfSpcCbG5skPX66aw0n6aPAY8AnI+LZRudJjZ/umiXOR1KzxLmkZolr+hNH8+bNi8WLF0/qPovFIgAtLc3zGNeMmaE5c49k/uEPf7gzIib8nwSavqSLFy9m06ZNk7rPQqEAQFdX16TudyKaMTM0Z+6RzDNnzhzPWzIrap6HJ7NpyiU1S5xLapY4l9QscS6pWeJcUrPETdmSbn5xDxue2M7gcLHRUcwmZMqWdMMT2/nDWzdRGBxudBSzCZmyJTWbKlxSs8S5pGaJc0nNEueSmiXOJTVLnEtqljiX1CxxLqlZ4lxSs8S5pGaJc0nNEueSmiXOJTVLnEtqlrialFTSakk7JD2WG3utpA2Snsq+zs3GJelmSX2SHpV0em6bFdn6T0laUYtsZs2uVkfSLwHLRo1dB9wXET3AfdltKH1ce092WQl8BkqlBm4EzgCWAjeOFNtsOqtJSSPiAeClUcPLgTXZ9TXAO3Pjt0bJQ8DRkhYAFwIbIuKliHgZ2MChxTebdur5MRPzI2IbQERsk3RcNr4Q2JJbrz8bqzQ+pmKxeODP+ucNDg4CMFAYoC2GjiR/RQMDAzWdbzI0Y2Zozty1ztyIE0cqMxZjjB86gbRS0iZJm3bu3FnTcGapqeeRdLukBdlRdAGwIxvvB7pz6y0Ctmbj54wav7/cxBGxClgF0NvbG+U+zKetrQ2Azq5OujrbJvJ9VNRMHyI0ohkzQ/PmroV6HknXASNnaFcAd+fG35ud5T0T+EX2tPge4AJJc7MTRhdkY2bTWk2OpJJup3QUnCepn9JZ2o8DayVdDTwHXJatvh64GOgD9gJXAUTES5I+CmzM1vv7iBh9Msps2qlJSSPiigqLziuzbgDXVJhnNbC6FpnMpgq/48gscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S1zdSyrpZ5L+T9IjkjZlY6+VtEHSU9nXudm4JN0sqU/So5JOr3c+s9RN1pH07RHxpojozW5fB9wXET3AfdltgIuAnuyyEvjMJOUzS1ajnu4uB9Zk19cA78yN3xolDwFHS1rQiIBmqZgxCfsI4FuSAvhcRKwC5kfENoCI2CbpuGzdhcCW3Lb92di2SpMXi0UKhcIh44ODgwAMFAZoi6FafB8HDAwM1HS+ydCMmaE5c9c682SU9KyI2JoVcYOkn4yxrsqMxSErSSspPR2mu7u7NinNElX3kkbE1uzrDkl3AUuB7ZIWZEfRBcCObPV+IN+6RcDWMnOuAlYB9Pb2RldX1yH7bWtrA6Czq5OuzrbafUM55fabumbMDM2buxbq+ppU0ixJc0auAxcAjwHrgBXZaiuAu7Pr64D3Zmd5zwR+MfK02Gy6qveRdD5wl6SRfd0WEd+UtBFYK+lq4Dngsmz99cDFQB+wF7iqzvnMklfXkkbEM8Aby4y/CJxXZjyAa+qZyazZ+B1HZolzSc0S55KaJc4lNUucS2qWOJfULHEuqVniXFKzxLmkZolzSc0S55KaJc4lNUucS2qWOJfULHEuqVniXFKzxLmkZolzSc0S55KaJc4lNUucS2qWOJfULHEuqVniXFKzxLmkZolzSc0S55KaJc4lNUucS2qWOJfULHEuqVnikiuppGWSnpTUJ+m6Rucxa7SkSiqpFfg0cBGwBLhC0pLGpjJrrLp+0vcRWAr0ZZ8QjqQ7gOXAE5U2KBaLFAqFQ8YHBwcBuOXeJ3lN5wwkIUACSbSI7LaQoKXC8haVrowsbxEMDg7RImhvb6dFjFo+sj6IbJ5szpHlGhk7sDwbe9Xyg/s7suwH99MisW/fPloEnYXBV2U7fHbV5l/2CA0MDDR0/0ei1plTK+lCYEvudj9wxuiVJK0EVgJ0d3eXnWjR3C4AVn1nc60zTjsjDwAHHwxGHhwOfTAot3x08cs/MBx8UCL3QASBEK2tLa9+0Drcg1rF5Qfzds5o5S/PP4UFR3U24F6tXmolLfewHYcMRKwCVgH09vZGV1fXIRu9480ncsFpCykWIQiKARFBADFqrBil2xEQAcVsvWKxtOtixKvGC4UBAmhv7yjNU6ywfYyMHZxvZJzgwH4PZBu5XaTs9vm5I5+p7PjBuYlg3/5BIoLWGW2Hbj96X2W2f3XW8t/Tq8Zz88DorPn7a2Tu/P1ycP3h4WGKAS0tLbms5e6vg9mGX5Upymbqf3kvrwwMcd6S4zn5+LlH9tM6SVIraT+QPzQuArYe6WQdM1onHKicQqE0b7kHh1SNvCRopsxQv9ybX9zD2Z+8v6Zz1ktSJ46AjUCPpJMktQOXA+sanMmsoZI6kkbEkKRrgXuAVmB1RDze4FhmDZVUSQEiYj2wvtE5zFKR2tNdMxsluSOp2WT67Lef5peFQWZ1zGB2xwxmdsxgdkcrszpmMKu9NDarYwbtMxp3PHNJbVqaO6udebPb+en23Xz4vyu+V+aA9tYWZna05orbyrFzOvjE776Ro2a21TWrS2rT0ms629h4w/nsGyqyZ98Qe/cPs3vfEHv2DbF71O3S2DB79w8dGPtu34v84LldvP+sX3LGycfUNatLatOWJDrbWulsa2W8Nfte306u/ML365JrNJ84MkucS2qWOJfULHEuqdk4RQS7CoOTtj+fODKroFgMnt9V4Kkdr/DU9t08taN06dv+Cnv2DwMwq6P+FXJJbdobLgZbXtqblfAV+rbv5qc7XuHpHXsoDA4fWO+4OR30zJ/NZb3d9MyfzRtOOIpfX3hU3fO5pDYtDQ0X+du7H+eRLbt4+oXd7B8qHlh2wlGdnDp/DmeccQw9x82mZ/5sTj12Tt3ftFCJS2rT0vO7Ctz+8HOctvAo3vebizn1uNm8bv4cTjl2FnM6G1PGSlxSm9auOmsxl56+qNExxuSzu2aJc0nNEueSmiXOJTVLnEtqljiX1CxxLqlZ4lxSs8S5pGaJc0nNEueSmiXOJTVLnEtqljiX1CxxLqlZ4lxSs8S5pGaJq1tJJX1Y0vOSHskuF+eWXS+pT9KTki7MjS/LxvokXVevbGbNpN5/PuVTEfFP+QFJS4DLgTcAJwD3SnpdtvjTwG8D/cBGSesi4vAfeWU2hTXibxwtB+6IiH3As5L6gKXZsr6IeAZA0h3ZumOWtFgsUigU6pn3EAMDA5O6v1poxsxQv9wj8+7fv7/mPz+1zlzv16TXSnpU0mpJc7OxhcCW3Dr92Vil8UNIWilpk6RNO3furEdus2RM6Egq6V7g+DKLbgA+A3wUiOzrPwPvB1Rm/aD8A0aU229ErAJWAfT29kZXV9e4s9dCo/Y7Ec2YGWqfu7Oz9Hd229vbk79PJlTSiDi/mvUkfR74enazH+jOLV4EbM2uVxo3m7bqeXZ3Qe7mJcBj2fV1wOWSOiSdBPQADwMbgR5JJ0lqp3RyaV298pk1i3qeOPqEpDdResr6M+CPACLicUlrKZ0QGgKuiYhhAEnXAvcArcDqiHi8jvnMmkLdShoR7xlj2U3ATWXG1wPr65XJrBn5HUdmiXNJzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiXNJzRLnkpolziU1S5xLapY4l9QscRMqqaTLJD0uqSipd9Sy6yX1SXpS0oW58WXZWJ+k63LjJ0n6vqSnJH01+yBhs2lvokfSx4BLgQfyg5KWUPqk7jcAy4B/k9QqqRX4NHARsAS4IlsX4B+BT0VED/AycPUEs5lNCRMqaUT8OCKeLLNoOXBHROyLiGeBPmBpdumLiGciYj9wB7BckoBzgTuz7dcA75xINrOpol6f9L0QeCh3uz8bA9gyavwM4BhgV0QMlVl/TMVikUKhMLG04zQwMDCp+6uFZswM9cs9Mu/+/ftr/vNT68yHLamke4Hjyyy6ISLurrRZmbGg/JE7xli/UqaVwEqA7u7uSquZTQmHLWlEnH8E8/YD+fYsArZm18uN7wSOljQjO5rm1y+XaRWwCqC3tze6urqOIOLENWq/E9GMmaH2uTs7iwC0t7cnf5/U61cw64DLJXVIOgnoAR4GNgI92Zncdkonl9ZFRAD/C7wr234FUOkobTatTPRXMJdI6gfeAnxD0j0AEfE4sBZ4AvgmcE1EDGdHyWuBe4AfA2uzdQH+GvgLSX2UXqN+cSLZzKaKCZ04ioi7gLsqLLsJuKnM+HpgfZnxZyid/TWzHL/jyCxxLqlZ4lxSs8S5pGaJc0nNEueSmiXOJTVLnEtqljiX1CxxLqlZ4lxSs8S5pGaJc0nNEueSmiXOJbVpqbVFdL+2i5nt9fozX7WTfkKzOlg0dyYP/tW5jY5RFR9JzRLnkpolziU1S5xLapY4l9QscS6pWeJcUrPEuaRmiVPpEx6al6QXgM0N2PU8Sp9h00yaMTM0Z+55wKyIOHaiEzV9SRtF0qaI6D38muloxszQnLlrmdlPd80S55KaJc4lPXKrGh3gCDRjZmjO3DXL7NekZonzkdQscS6pWeJcUrPEuaR1IOmtkj4r6QuSvtfoPNWQdI6kB7Pc5zQ6T7UkvT7LfKekP250nmpIOlnSFyXdWc36LukoklZL2iHpsVHjyyQ9KalP0nVjzRERD0bEB4CvA2vqmTfLNuHMQAC7gU6gv15Z82p0X/84u6/fDdT9DQ81yvxMRFxd9U4jwpfcBXgbcDrwWG6sFXgaOBloB34ELAFOo1TE/OW43HZrgdc0Q2agJdtuPvAfzXRfA78DfA+4slkyZ9vdWdU+G12KFC/A4lH/CG8B7sndvh64/jBznAh8vpkyZ+u1V/vDk1LubN1vNFPmau9nP92tzkJgS+52fzY2lquBf69bosMbV2ZJl0r6HPBl4JY6ZxvLeHOfI+nmLPv6eoerYLyZj5H0WeDNkq4/3OT+k57VUZmxMd8FEhE31ilLtcaVOSK+BnytfnGqNt7c9wP31ytMlcab+UXgA9VO7iNpdfqB7tztRcDWBmWpVjNmhubMXdfMLml1NgI9kk6S1A5cDqxrcKbDacbM0Jy565t5sk4QNMsFuB3YBgxSeoS8Ohu/GPgppbN4NzQ6Z7NnbtbcjcjsN9ibJc5Pd80S55KaJc4lNUucS2qWOJfULHEuqVniXFKzxLmkZolzSc0S9//HJQXxs2pm0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd145adcf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize = (3, 4))\n",
    "ax.semilogx(sigma, mesh.vectorCCz)\n",
    "ax.grid(which = 'major', linestyle = '-', linewidth=0.2)\n",
    "ax.set_title('conductivity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)\n",
    "mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap\n",
    "mtrue = np.log(sigma[active])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'model')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOIAAAEICAYAAABCuCiHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFDVJREFUeJzt3X1wZXV9x/H35ybZZEGQXVYR1sXFcZkWxFFYgdLp+IAgorgWrQXrsEVHqoNa21oFmY6tsjPYOrWiyKC4ozitiM/rCMXFSm0rT4sKuFqGACILVF15ht0k995v/zi/m5wkN9lskpP87s3nNZPJub/zcL+5uZ9z7j25OV9FBGa2uGqLXYCZOYhmWXAQzTLgIJplwEE0y4CDaJYBB9GmJekLki6c4bK/lPSqqmvqRg6iWQYcRLMMOIhdIr0s/FtJt0t6StLnJR0k6RpJT0i6TtKKtOzrJW2X9Kik6yX9fmk7L5H047TOV4CBCffzOkk/Tev+SNKLFvhH7UoOYnd5I3AScDhwGnAN8CFgFcXv+r2SDge+DLwPeBZwNfAdScskLQO+BXwJWAl8NW0TAElHA5uBvwAOBC4DtkjqX5Cfros5iN3lUxHx64h4APgv4KaI+ElEDAHfBF4C/Cnw3YjYGhEjwMeB5cAJwPFAH/AvETESEV8Dbilt/x3AZRFxU0Q0IuKLwFBaz+agd7ELsHn169L0rja3nwEcAtzXGoyIpqT7gdVAA3ggxv8nwH2l6ecBGyW9pzS2LG3T5sBBXHoeBI5q3ZAkYA3wABDAakkqhfFQ4O40fT+wKSI2LWC9S4Jfmi49VwGvlXSipD7gbyheXv4IuAGoU7yX7JV0OnBsad3PAe+UdJwK+0p6raT9FvqH6DYO4hITEXcCbwU+BeykOKlzWkQMR8QwcDrw58AjFO8nv1FadxvF+8RPp/mDaVmbI/kfg80Wn4+IZhlwEM0y4CCaZcBBNMtAx/8dcdWqVbF27dpZr99sNgGo1Tpzn9TJ9Xd77bfeeuvOiHjWTLbX8UFcu3Yt27Ztm/X6u3btAmD58uXzVdKC6uT6u712SfdNOXOCztsVmXUhB9EsAw6iWQayC6KkUyTdKWlQ0nmLXY/ZQsgqiJJ6gEuA1wBHAGdKOmJxqzKrXlZBpPik/2BE3JM+gHwlsGGRazKrXG5/vlhN8T9vLTuA46Zbodlsjp5Kno3du3fPet0cdHL9rn1MbkdEtRmb9O8hks6RtE3Stp07d7bd0OBvn+LiH9zDb58cmu8azeZdbkfEHRT/Ld7yXIr/KB8nIj4LfBZg/fr10e6Pqg88/hiX/vCXnPbiNRw6gz8Yd+Iflcs6uX7Xnt8R8RZgnaTD0hXFzgC2LHJNZpXL6ogYEXVJ7wauBXqAzRGxfZHLMqtcVkEEiIirKa61abZk5PbS1GxJchDNMuAgmmXAQTTLgINolgEH0SwDDqJZBhxEsww4iGYZyO6TNfPtyaE6T+weoSbRU1PpOxQdycwWX9cGsa+nCNmbL7thymUk6JGK7zXRoyKotdpYaGtpXjFeLF9rBXp0WdKyYyEvB79WEz1pO0rrFdNj9ysV2xm73zY1pG2P3q9Eo1GnJjHQ3ze5hprG1zupzvHbHvv5J/ycE36mqX/Gdvfb7mfyTnCirg3iH75gFRedfhRPDTeICBrNoBFBsxk0AxrNoBnB8PAIjQhqPb00S8s0oliu2Yy0LDQjRtdrTTeaFNtPtyNtuxHBSKOZpkn3W16f0v0EzebYeqP1trYXk6c7XWsnOGnHMS6043ccozuv0R1Z2nmN7shaOwDa72za7BzH7UQmrZdqmLDjrQmajTrL+3p46wnPZ7+Bvjk/Hl0bxIG+Hs449tA9LtepF7ltBfupp3fRiKC/f6DtjmY0+E1KoW+/4yh2AKTxsWXGrd/a/uh0lO6LcTupyTubtP10X0PDIzSaQU9vb2nH07rf0s5ryvttv3NsBtQbzTY70wm1taYnbLvd41her+yQlc9gw4tXz/n32bVB7Ha1mqghlvUW59uW93fer7KTd4L/++DDnPrpG6k35ufVic+amu2lWk301ub3Pa6DaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJaByj4XJemfgNOAYeBu4OyIeDTNOx94O9AA3hsR16bxU4BPUlzl+/KIuKiq+szKms1gd73B7pEmu0ca6avJ7nqDodHvY/Pv/90T83r/VX5AcStwfrqM/seA84EPpsajZwBHAocA10k6PK1zCXASRTOaWyRtiYifV1ijZajZDIbqKRATwjE6PtJkqF4KzIR5u9O8oYnzxm2vWYSr3mBkFp8Z7ZE45ID5+ZxsZUGMiO+Vbt4IvClNbwCujIgh4F5JgxQNSiE1KQWQ1GpSOm0Q3R+x2vojUijqzdKRId0uT4802F1vpoAUT/ThtN5oaErbGK43eXq4ztBIk+FGjJs3m1C09NaKD8IP9Nbo76sx0NvDQF+N/t4aA309rNynb/y80nR/Wq61bLt5xVgPNEZY3ldjxX77zOn5N1r3nLcwM28DvpKmV1MEs2VHGoMZNimVdA5wDsCaNWvaLdKVIoLhxvgn++NP7WKo3qRZ21U88ScEYzQ8o9OtI8nYNobqremJ85oMN5qzrrcmRp/ErSfw2JO7xoH79LGst8a+/X2jT/7ie43+3vL0nue1tttbW5jTHvO9/5tTECVdBzynzawLIuLbaZkLgDrwr63V2iwftD9x1HbXOJP+iHtrb7cxLhQTXg6Vn9Cj7zda4Zj4Emr0SDI2r932WvNmqyZKT+aeSdP7LV9WelIXR5GBvp4iPOnJPzB6hBmbX57XP2FeX8/0oejUf4Mqm6/a5xTEiHjVdPMlbQReB5wYEa1QTdeMdI9NSveiNgZ/8ySP7hoZfZ/Q7v3GE08PMVxvUkdtwzFUCtC49x71JjHLV1Aad6SY8KTvrXHAPqVQ9I2FZfSoUgpRLeoM9Paw/77Lx7bTOoL0jQWrr0e+PEXGqjxregrwQeBlEfF0adYW4N8k/TPFyZp1wM0UR8p1kg4DHqA4ofOW2d7/Dff8jrd87qYZLTv2BB9/ROjv6+GZy/sY2K9/NCQTw9Ffesk1GpjSdsa9T0nzlvXU5i0U3XBUsWrfI34a6Ae2pifdjRHxzojYLukqipMwdeDciGgAzGeT0sd31QG48A0v5PCD9hv/8qoUkubIEJL8RLZFVeVZ0xdMM28TsKnN+Lw3KT360BUcccj+U87fVffLNVt8/mSNWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJaByoMo6f2SQtKqdFuSLpY0KOl2SUeXlt0o6a70tbHq2sxyUWk3KElrKPod/qo0/BqKy+yvo+j2dClwnKSVwIeB9RTNZ25N/REfqbJGsxxU3ZbtE8AHgG+XxjYAV6SmNDdKOkDSwcDLga0R8TCApK3AKcCXp7uDqfojDg8PAzA0tJtdu/qmXN/9ERePax9T2UtTSa8HHoiI2ybMWs3kPoirpxlvt+1zJG2TtG3nzp3zWLXZ4qisPyLwIeDkdqu1GYtpxicPzqA/4rJlywDo7x+YUYOZTm9C08n1u/aK+iNKOgo4DLgtdYJ6LvBjSccydX/EHRQvT8vj18+lPrNOUclL04i4IyKeHRFrI2ItRciOjoj/o+iPeFY6e3o88FhEPETRju1kSSskraA4ml5bRX1muan6ZE07VwOnAoPA08DZABHxsKSPArek5T7SOnFj1u0WJIjpqNiaDuDcKZbbDGxeiJrMcuJP1phlwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGag0iJLeI+lOSdsl/WNp/PzUH/FOSa8ujZ+SxgYlnVdlbWY5qewCw5JeQdGC7UURMSTp2Wn8COAM4EjgEOA6SYen1S6h6Ke4A7gl9Uf8eVU1muWiyit9vwu4KCKGACLiN2l8A3BlGr9X0iBwbJo3GBH3AEi6Mi07bRDdH7Fz63ftY6p8aXo48EeSbpL0n5JemsbdH9Fsgir7I/YCK4DjgZcCV0l6PlP3QWy3U3B/xBnq5Ppde0X9EQEkvQv4Rmo6c7OkJrCKqfsjMs24WVer8qXpt4BXAqSTMcuAnRT9Ec+Q1C/pMGAdcDNFO7Z1kg6TtIzihM6WCuszy0aVJ2s2A5sl/QwYBjamo+N2SVdRnISpA+dGRANA0rspmpP2AJsjYnuF9Zllo7IgRsQw8NYp5m0CNrUZv5qikanZkuJP1phlwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGagsiJJeLOlGST9NDWOOTeOSdHHqgXi7pKNL62yUdFf62lhVbWa5qfJK3/8I/ENEXCPp1HT75cBrKC6zvw44DrgUOE7SSuDDwHqK5jO3pv6Ij1RYo1kWqgxiAPun6Wcy1lBmA3BFuvz+jZIOkHQwRUi3RsTDAJK2AqcAX57uTtwfsXPrd+1jqgzi+4BrJX2c4iXwCWl8XvojAucArFmzpt0iZh2lyv6IJwJ/FRFfl/Rm4PPAq5i6P+JU45MH3R9xkk6u37VX2x/xCuAv082vApen6an6I+6geHlaHr9+LvWZdYoq/3zxIPCyNP1K4K40vQU4K509PR54LCIeomjHdrKkFZJWACenMbOuV+V7xHcAn5TUC+wmvaejaLt2KjAIPA2cDRARD0v6KEXDUoCPtE7cmHW7Kvsj/jdwTJvxAM6dYp3NFA1OzZYUf7LGLAMOolkGHESzDDiIZhlwEM0y4CCaZcBBNMuAg2iWAQfRLAMOolkGHESzDDiIZhlwEM0y4CCaZcBBNMuAg2iWAQfRLAMOolkGHESzDDiIZhlwEM0y4CCaZcBBNMvAnIIo6U8kbZfUlLR+wrzzUw/EOyW9ujR+ShoblHReafwwSTel3ohfkbRsLrWZdZK5HhF/BpwO/LA8KOkI4AzgSIrWap+R1COpB7iEokfiEcCZaVmAjwGfiIh1wCPA2+dYm1nHmGsTml8ASJMaOW0AroyIIeBeSYPAsWneYETck9a7Etgg6RcU/THekpb5IvD3FE1Mp+X+iJ1bv2sfU9V7xL3tgXgg8GhE1CeMtyXpnNQOfNvOnTvntXCzxbDHI+J0PRAj4ttTrdZmLGgf/L3qjQjuj9hOJ9fv2mcQxOl6IE5jqh6ITDG+EzhAUm86KpaXN+t6Vb003QKcIalf0mHAOuBmipZr69IZ0mUUJ3S2pA5RPwDelNbfCEx1tDXrOnP988UfS9oB/AHwXUnXAkTEduAq4OfAvwPnRkQjHe3eTdGA9BfAVWlZgA8Cf51O7BxI0erbbEmY61nTbwLfnGLeJmBTm/GrKZqVThy/h7Ezq2ZLij9ZY5YBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJplwEE0y4CDaJYBB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJploJL+iJJOknSrpDvS91eW5h2TxgclXazUSkrSSklbU3/ErZJWzKU2s05SSX9Eil4Wp0XEURSXz/9Sad6lwDkUl+FfR9E/EeA84PupP+L3022zJaGS/ogR8ZPSze3AgKR+YCWwf0TckNa7AngDcA1FT8WXp3W+CFxPcRn+abk/YufW79rHLMR7xDcCP0lNS1dTdIpqKfdBPCgiHgJI35891QbdH9G6TVX9EVvrHknRkvvk1lCbxabsgzgV90ecrJPrd+3V9UdE0nMpGtScFRF3p+EdFL0PW8p9EH8t6eCIeEjSwcBvZnO/Zp2okpemkg4AvgucHxH/0xpPLzmfkHR8Olt6FmN9ELdQnNgB90e0JaaS/ogUPRBfAPydpJ+mr9Z7vncBlwODwN0UJ2oALgJOknQXcFK6bbYkVNIfMSIuBC6cYp1twAvbjP8OOHEu9Zh1Kn+yxiwDXRvEg/bv59SjnsN+A3M66JstiK59lr7k0BV85s+OWewyzGaka4+IZp3EQTTLgINolgEH0SwDDqJZBhxEsww4iGYZcBDNMqCIvf53wKxI+i1w3xw3s4ri8h6dqpPr7+banxcRz5rJhjo+iPNB0raIWL/nJfPUyfW79oJfmpplwEE0y4CDWPjsYhcwR51cv2vH7xHNsuAjolkGHESzDCzpIEp6v6SQtGqK+Y3Sxa+2LHR9ezKD+jemXiJ3SdrYbpmFJumjkm5Pj+n3JB0yxXLZPfZ7UfteP+5L9j2ipDUUV5P7PeCYiJj0h1lJT0bEMxa8uBnYU/2SVgLbgPUUF3G+NS33yELXOqGu/SPi8TT9XuCIiHhnm+Wye+xnUvtsH/elfET8BPABZnGl8Uzsqf5XA1sj4uH0JNjKWMOfRdN6Iif70kGP/wxrn9Xj3rXXrJmOpNcDD0TEbRMb6EwwIGkbUAcuiohvLUiBezDD+lcD95dul/uMLCpJmyguLv0Y8IopFsv1sd9T7bN63Ls2iNP17AA+xFg/jukcGhEPSno+8B+S7ii1D6jUPNQ/L31GZmNP/VIi4gLgAknnU1yM+sNtll2Ux34eap/d4x4RS+oLOIqir8Yv01cd+BXwnD2s9wXgTZ1SP3AmcFnp9mXAmYtd/4Qanwf8bAbLZfHYz6T22T7ui/4DLfZXejKvajO+AuhP06uAuyjenC96zTOsfyVwb/o5VqTplRnUu640/R7ga53y2M+w9lk97ov+RFrsr/ITmeJM1+Vp+gTgDuC29P3ti13r3tSfbr+NosfIIHD2Yteaavo6Rafp24HvAKs75bGfSe2zfdyX7J8vzHKylP98YZYNB9EsAw6iWQYcRLMMOIhmGXAQzTLgIJpl4P8BBmLwOFIVxJkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd14aff5c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize = (3, 4))\n",
    "ax.plot(mtrue, mesh.vectorCCz[active])\n",
    "ax.grid(which = 'major', linestyle = '-', linewidth=0.2)\n",
    "ax.set_title('model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min diffusion distance  114.18394343377736 max diffusion distance  510.64611891383385\n"
     ]
    }
   ],
   "source": [
    "# TDEM survey\n",
    "rxlocs = Utils.ndgrid([np.r_[50.], np.r_[0], np.r_[0.]])\n",
    "srcLoc = np.r_[0., 0., 0.]\n",
    "times = np.logspace(-4, np.log10(2e-3), 10)\n",
    "\n",
    "print('min diffusion distance ', 1.28*np.sqrt(times.min()/(sig_half*mu_0)),\n",
    "      'max diffusion distance ', 1.28*np.sqrt(times.max()/(sig_half*mu_0)))\n",
    "rx = TDEM.Rx.Point_b(rxlocs, times, 'z')\n",
    "src = TDEM.Src.MagDipole(\n",
    "    [rx],\n",
    "    waveform=TDEM.Src.StepOffWaveform(),\n",
    "    loc=srcLoc  # same src location as FDEM problem\n",
    ")\n",
    "\n",
    "surveyTD = TDEM.Survey([src])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "prbTD = TDEM.Problem3D_b(mesh, sigmaMap=mapping, Solver=Solver)\n",
    "prbTD.timeSteps = [(5e-5, 10), (1e-4, 10), (5e-4, 10)]\n",
    "prbTD.pair(surveyTD)\n",
    "\n",
    "std = 0.03\n",
    "surveyTD.makeSyntheticData(mtrue, std)\n",
    "surveyTD.std = std\n",
    "surveyTD.eps = np.linalg.norm(surveyTD.dtrue)*1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEDCAYAAADQunSaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFWhJREFUeJzt3X+Q3Hd93/Hnu7Jwrhz0guWeJWEiCEQZAxPsUw1NSEdKppVxUuw6pjWTMSGF0UDxNE2LCgozkGGmY4gmv/jpcYJrnFAfpVZVjVGiZowcm0kdkCxj2biyFUgbnTQY7Mrm4EIk5d0/9nv26nJ7u3ffvf3ufe75mPmOdr/f737vpa9Wr/3e5/u970VmIkkq099rOoAkaflY8pJUMEtekgpmyUtSwSx5SSqYJS9JBWu05CPitoh4MiIe6dP2/jgiTkfE3R2WfzwipvvxtSRpJWj6SP524Ko+bm83cON8CyJiCzDWx68lSUOv0ZLPzPuAp9vnRcSPVkfkhyPi/oj48UVs7x7gu3PnR8QaWh8A/7FuZklaSS5oOsA8bgXelZlPRMTrgU8BP1NzmzcB+zLzVETUDihJK8VQlXxEjAI/CXyhrYwvrJZdB3x4npdNZeb2Bba5AXgLsLWvYSVpBRiqkqc1fHQ6M183d0Fm7gH2LGGblwOvBI5XHxx/PyKOZ+YrayWVpBWg6ROv58nMZ4FvRsRbAKLlJ2pu84uZeUlmbsrMTcD3LXhJq0XTl1DeCfwvYHNEnIiIdwC/CLwjIr4GPApcs4jt3Q98AfjZansdh3EkaTUIbzUsSeUaquEaSVJ/dT3xGhGXAncAlwB/C9yamb87Z52twP8AvlnN2pOZ810J85x169blpk2blhB5+X3ve9/jhS98YdMxOhr2fDD8Gc1Xj/nqqZPv8OHD38nMi3t+QWYuOAHrgSuqxy8CHgcum7POVuDubttqnyYmJnJYHTx4sOkICxr2fJnDn9F89Zivnjr5gEO5iK7tOlyTmacy88Hq8XeBx4CNi/rokSQ1YlEnXiNiE3Af8JpsXe44O38rcBdwAjgJvDczH53n9TuAHQDj4+MTk5OTNaIvn+npaUZHR5uO0dGw54Phz2i+esxXT51827ZtO5yZW3p+Qa+H/MAocBi4bp5lLwZGq8dXA090257DNUs37Pkyhz+j+eoxXz1DNVwDEBFraR2pfy5bP3k694Pi2cycrh7vB9ZGxLqeP2kkScuia8lH614AnwEey8zf6rDOJdV6RMSV1Xaf6mdQSdLi9XLvmp+idY/2oxHxUDXv14CXAWTmLcD1wLsj4iwwA9xQfVvRV3uPTLH7wDFOnp5hw9gIO7dv5trLPQcsSZ10LfnM/DKw4P15M/MTwCf6FWo+e49MsWvPUWbOnANg6vQMu/YcBbDoJamDFfMTr7sPHHuu4GfNnDnH7gPHGkokScNvxZT8ydMzi5ovSVpBJb9hbGRR8yVJK6jkd27fzMjaNefNG1m7hp3bNzeUSJKG37D9ZqiOZk+uenWNJPVuxZQ8tIreUpek3q2Y4RpJ0uJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SCWfKSVDBLXpIKZslLUsEseUkqmCUvSQWz5CWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgrWteQj4tKIOBgRj0XEoxHxK/OsExHxsYg4HhEPR8QVyxNXkrQYF/SwzlngP2TmgxHxIuBwRPxJZn69bZ03Aa+qptcDn67+lCQ1qOuRfGaeyswHq8ffBR4DNs5Z7Rrgjmx5ABiLiPV9TytJWpTIzN5XjtgE3Ae8JjOfbZt/N/CRzPxy9fwe4H2ZeWjO63cAOwDGx8cnJicn6+ZfFtPT04yOjjYdo6NhzwfDn9F89Zivnjr5tm3bdjgzt/T8gszsaQJGgcPAdfMs+yLwxrbn9wATC21vYmIih9XBgwebjrCgYc+XOfwZzVeP+eqpkw84lD32dmb2dnVNRKwF7gI+l5l75lnlBHBp2/OXAid7/qSRJC2LXq6uCeAzwGOZ+VsdVtsHvK26yuYNwDOZeaqPOSVJS9DL1TU/BdwIHI2Ih6p5vwa8DCAzbwH2A1cDx4HvA7/c/6iSpMXqWvLZOpkaXdZJ4D39CiVJ6g9/4lWSCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SCWfKSVDBLXpIKZslLUsEseUkqmCUvSQWz5CWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgrWteQj4raIeDIiHumwfGtEPBMRD1XTB/sfU5K0FBf0sM7twCeAOxZY5/7M/Pm+JJIk9U3XI/nMvA94egBZJEl9FpnZfaWITcDdmfmaeZZtBe4CTgAngfdm5qMdtrMD2AEwPj4+MTk5udTcy2p6eprR0dGmY3Q07Plg+DOarx7z1VMn37Zt2w5n5paeX5CZXSdgE/BIh2UvBkarx1cDT/SyzYmJiRxWBw8ebDrCgoY9X+bwZzRfPearp04+4FD20LGzUy9j8t0+JJ5te7w/Ij4VEesy8zt1t70Ye49MsfvAMU6enmHD2Ag7t2/m2ss3DjKCJA2d2iUfEZcA38rMjIgraY3zP1U72SLsPTLFrj1HmTlzDoCp0zPs2nMUwKKXtKp1LfmIuBPYCqyLiBPAh4C1AJl5C3A98O6IOAvMADdU31IMzO4Dx54r+FkzZ86x+8AxS17Sqta15DPzrV2Wf4LWJZaNOXl6ZlHzJWm1KOInXjeMjSxqviStFkWU/M7tmxlZu+a8eSNr17Bz++aGEknScKh94nUYzI67e3WNJJ2viJKHVtFb6pJ0viKGayRJ87PkJalglrwkFcySl6SCWfKSVDBLXpIKZslLUsEseUkqmCUvSQWz5CWpYJa8JBXMkpekglnyklQwS16SCmbJS1LBLHlJKpglL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SCdS35iLgtIp6MiEc6LI+I+FhEHI+IhyPiiv7HlCQtRS9H8rcDVy2w/E3Aq6ppB/Dp+rEkSf3QteQz8z7g6QVWuQa4I1seAMYiYn2/AkqSli4ys/tKEZuAuzPzNfMsuxv4SGZ+uXp+D/C+zDw0z7o7aB3tMz4+PjE5OVkr/HKZnp5mdHS06RgdDXs+GP6M5qvHfPXUybdt27bDmbml5xdkZtcJ2AQ80mHZF4E3tj2/B5jots2JiYkcVgcPHmw6woKGPV/m8Gc0Xz3mq6dOPuBQ9tDbs1M/rq45AVza9vylwMk+bFeSVFM/Sn4f8LbqKps3AM9k5qk+bFeSVNMF3VaIiDuBrcC6iDgBfAhYC5CZtwD7gauB48D3gV9errCSpMXpWvKZ+dYuyxN4T98SLaO9R6bYfeAYJ0/PsGFshJ3bN3Pt5RubjiVJy6ZryZdi75Epdu05ysyZcwBMnZ5h156jABa9pGKtmtsa7D5w7LmCnzVz5hy7DxxrKJEkLb9VU/InT88sar4klWDVlPyGsZFFzZekEqyakt+5fTMja9ecN29k7Rp2bt/cUCJJWn6r5sTr7MlVr66RtJqsmpKHVtFb6pJWk1UzXCNJq5ElL0kFs+QlqWCWvCQVzJKXpIJZ8pJUMEtekgpmyUtSwSx5SSqYJS9JBbPkJalglrwkFcySl6SCWfKSVDBLXpIKZslLUsEseUkqmCUvSQWz5CWpYJa8JBVsVf0i78Xae2SK3QeOcfL0DBvGRti5fbO/CFzSimLJd7D3yBS79hxl5sw5AKZOz7Brz1EAi17SiuFwTQe7Dxx7ruBnzZw5x+4DxxpKJEmLZ8l3cPL0zKLmS9IwsuQ72DA2sqj5kjSMLPkOdm7fzMjaNefNG1m7hp3bNzeUSJIWr6eSj4irIuJYRByPiPfPs/ztEfHtiHiomt7Z/6iDde3lG7n5uteycWyEADaOjXDzda/1pKukFaXr1TURsQb4JPBPgRPAVyNiX2Z+fc6qn8/Mm5YhY2OuvXyjpS5pRevlSP5K4HhmfiMz/waYBK5Z3liSpH6IzFx4hYjrgasy853V8xuB17cftUfE24GbgW8DjwO/mpl/Nc+2dgA7AMbHxycmJyf79Nfor+npaUZHR5uO0dGw54Phz2i+esxXT51827ZtO5yZW3p+QWYuOAFvAX6/7fmNwMfnrHMRcGH1+F3Al7ptd2JiIofVwYMHm46woGHPlzn8Gc1Xj/nqqZMPOJRd+rV96mW45gRwadvzlwIn53xQPJWZP6ie/h4w0fOnjCRp2fRS8l8FXhURL4+IFwA3APvaV4iI9W1P3ww81r+IkqSl6np1TWaejYibgAPAGuC2zHw0Ij5M69uGfcC/jYg3A2eBp4G3L2PmFcWbnElqUk83KMvM/cD+OfM+2PZ4F7Crv9FWPm9yJqlp/sTrMvImZ5KaZskvI29yJqlplvwy8iZnkppmyS8jb3ImqWn+ZqhlNHty1atrJDXFkl9m3uRMUpMcrpGkglnyklQwS16SCmbJS1LBPPG6wnlvHEkLseRXMO+NI6kbh2tWMO+NI6kbS34F8944krqx5Fcw740jqRtLfgXz3jiSuvHE6wrmvXEkdWPJr3DDcG+c+S7jHGs0kaRZDteoltnLOKdOz5A8fxnnn50803Q0SVjyqqnTZZx3PW7JS8PAklctnS7XfOqvc8BJJM3HklctnS7XvOiHYsBJJM3HklctnS7j/IUfW9tQIkntvLpGtXS6jHPsmScaTtYbb/Cm0lnyqm2+yzjvvXf4S94bvGk1cLhGq5Y3eNNq4JG8Vi1v8NbikFXZPJLXquUN3jr/MNveI1NNR1OfWPJatbzBm0NWq4HDNVq1vMGbQ1aDNjs0NnV6ho0PfGkg7zdLXqvaMNzgrUkbxkaYmqfQV9OQ1aA0dTWXwzXSKuaQ1eA0NTTmkby0ijlkNThNDY31VPIRcRXwu8Aa4Pcz8yNzll8I3AFMAE8B/yoz/7K/USUth9U+ZDUoTQ2NdR2uiYg1wCeBNwGXAW+NiMvmrPYO4P9l5iuB3wY+2u+gkrSSNTU01suY/JXA8cz8Rmb+DTAJXDNnnWuAz1aP/xvwsxHhbQglqXLt5Ru5+brXsrE6ct84NsLN17122b+LisyF7/sdEdcDV2XmO6vnNwKvz8yb2tZ5pFrnRPX8L6p1vjNnWzuAHQDj4+MTk5OT/fy79M309DSjo6NNx+ho2PPB8Gc0Xz3mq6dOvm3bth3OzC29rt/LmPx8R+RzPxl6WYfMvBW4FWDLli25devWHr784N17770MazYY/nww/BnNV4/56hlkvl6Ga04Al7Y9fylwstM6EXEB8A+Ap/sRUJK0dL2U/FeBV0XEyyPiBcANwL456+wDfql6fD3wpew2DiRJWnZdh2sy82xE3AQcoHUJ5W2Z+WhEfBg4lJn7gM8AfxARx2kdwd+wnKElSb3p6Tr5zNwP7J8z74Ntj/8aeEt/o0mS6up6dc2yfeGIbwP/p5Ev3t064Dtd12rOsOeD4c9ovnrMV0+dfD+SmRf3unJjJT/MIuLQYi5RGrRhzwfDn9F89ZivnkHm8wZlklQwS16SCmbJz+/WpgN0Mez5YPgzmq8e89UzsHyOyUtSwTySl6SCWfKSVLLMLGYCrgKOAceB98+z/ELg89XyPwc2tS3bVc0/Bmzvtk3gc9X8R4DbgLXV/K3AM8BD1fTBhvLdDnyzLcfrqvkBfKxa/2Hgioby3d+W7SSwt6H9dxvwJPDInG29BPgT4Inqzx9uaP91yrcb+N9Vhv8OjFXzNwEzbfvvloby/Tow1Zbj6qW+V5Yp3+fbsv0l8NCg9x+t+30dBB4DHgV+pc77r2MvLrVQh22idcuFvwBeAbwA+Bpw2Zx1/s3sPxqtWy98vnp8WbX+hcDLq+2sWWibwNXVDg/gTuDd1fytwN1DkO924Pp5clwN/FGV+w3AnzeRb8527wLeNuj9Vy37J8AV/N0S+A2q/8jA+4GPDnr/dcn3z4ALqscfbcu3ae66DeX7deC98+RY8nuln/nmbPc3qQ4mBrn/gPVUJQ28CHic5///Lur9t9BU0nBNnV9ucg0wmZk/yMxv0vqUvHKhbWbm/qwAX6F1d86hybeAa4A7qugPAGMRsb6pfBHxIuBngL1dci9HPjLzPua/Y2r7tj4LXNs2f1D7r2O+zPyfmXm2evoAzbz/Ftp/ndR5L/c9X/X6f0nrQG0hfc+Xmacy88Eq53dpHdFvnGdbvbz/Oiqp5DcCf9X2/ATP77C/s071H+QZ4KIFXtt1mxGxFrgR+OO22f84Ir4WEX8UEa9uMN9/ioiHI+K3q9/Du1CORvYf8C+AezLz2bZ5g9p/CxnPzFPVtk4B/7BLjkHna/evaR3dzXp5RByJiD+NiJ/ukns5891Uvf9ui4gf7pKjqf3308C3MvOJtnkD338RsQm4nNYwDyz+/ddRSSVf55ebLHZ+u08B92Xm/dXzB2ndW+IngI/z/BHqoPPtAn4c+Ee0xvfet8SvsVz5Zr2V84+iBrn/lmKQ+697mIgPAGdpnSMCOAW8LDMvB/498F8i4sUN5Ps08KPA66pMv9mHr7Ec/75z338D338RMUpryPLfzTnYmc+i/64llXydX27S6bULbjMiPgRcTOvNAEBmPpuZ09Xj/cDaiFg36HzVt4KZmT8A/jPVt69L/bv2O1+1jYuqXF+cnTfg/beQb81+G1z9+eTcr7GYv+sy5CMifgn4eeAXq2FDqiGBp6rHh2mN//7YoPNl5rcy81xm/i3wezTz/ltQtY3raJ0snc090P1XjQTcBXwuM/e0rbPY919n3QbtV8pE67bJ36B1YmP2xMir56zzHs4/MfJfq8ev5vwTI9+gdWKk4zaBdwJ/BozM+RqX8PwPmV0J/F9an76Dzre++jOA3wE+Uj3/Oc4/cfOVJvZf9bp3AZ9tav+1vW4T81+90n7i6zcGvf+65LsK+Dpw8Zz5F/P8ScdX0LrC5SUN5Fvf9vhXaY1JL/m90u98bfvwT5vaf9V76A7gd+bJtqj334LduNzlO8iJ1pnnx2l9+n6gmvdh4M3V4x8CvkDrxMdXgFe0vfYD1euOAW9aaJvV/LPVvPMu9QNuonU51NdonRD7yYbyfQk4SusSzz8ERqv5AXyyWv8osKWJfNWye2n9Avj2eYPef3fS+hb9DK2jpHdU8y8C7qF1Cds9wEsa2n+d8h2nNTZ73qV+wC+07b8HgX/eUL4/qPbPw7R+c9z6HrY1sHzVstuBd815/w1s/wFvpDXU8jBzLjVlCe+/TpO3NZCkgpU0Ji9JmsOSl6SCWfKSVDBLXpIKZslLUsEseUkqmCUvSQX7/7EKigcQSafcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd14f804a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(times, surveyTD.dobs, 'o')\n",
    "plt.grid(which = 'both')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# set up the inversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    " # Inversion Directives\n",
    "beta = Directives.BetaSchedule(coolingFactor=4, coolingRate=3)\n",
    "betaest = Directives.BetaEstimate_ByEig(beta0_ratio=2.)\n",
    "target = Directives.TargetMisfit()\n",
    "directiveList = [beta, betaest, target]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmisfit = DataMisfit.l2_DataMisfit(surveyTD)\n",
    "regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])\n",
    "reg = Regularization.Simple(regMesh)\n",
    "opt = Optimization.InexactGaussNewton(maxIterCG=10)\n",
    "invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the inversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SimPEG.InvProblem will set Regularization.mref to m0.\n",
      "\n",
      "    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
      "    ***Done using same Solver and solverOpts as the problem***\n",
      "model has any nan: 0\n",
      "============================ Inexact Gauss Newton ============================\n",
      "  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   \n",
      "-----------------------------------------------------------------------------\n",
      "x0 has any nan: 0\n",
      "   0  1.81e+02  5.83e+02  0.00e+00  5.83e+02    2.54e+02      0              \n",
      "   1  1.81e+02  2.97e+02  4.71e-01  3.83e+02    4.76e+01      0              \n",
      "   2  1.81e+02  2.51e+02  6.78e-01  3.74e+02    1.51e+01      0   Skip BFGS  \n",
      "   3  4.53e+01  2.38e+02  7.46e-01  2.71e+02    1.24e+02      0   Skip BFGS  \n",
      "   4  4.53e+01  4.15e+01  2.63e+00  1.61e+02    2.40e+01      0              \n",
      "   5  4.53e+01  3.08e+01  2.75e+00  1.56e+02    8.05e+00      0              \n",
      "   6  1.13e+01  2.78e+01  2.80e+00  5.96e+01    6.12e+01      0   Skip BFGS  \n",
      "------------------------- STOP! -------------------------\n",
      "1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.8378e+01\n",
      "1 : |xc-x_last| = 5.6540e-01 <= tolX*(1+|x0|) = 2.4026e+00\n",
      "0 : |proj(x-g)-x|    = 6.1194e+01 <= tolG          = 1.0000e-01\n",
      "0 : |proj(x-g)-x|    = 6.1194e+01 <= 1e3*eps       = 1.0000e-02\n",
      "0 : maxIter   =      20    <= iter          =      7\n",
      "------------------------- DONE! -------------------------\n"
     ]
    }
   ],
   "source": [
    "inv = Inversion.BaseInversion(invProb, directiveList=directiveList)\n",
    "m0 = np.log(np.ones(mtrue.size)*sig_half)\n",
    "reg.alpha_s = 5e-1\n",
    "reg.alpha_x = 1.\n",
    "prbTD.counter = opt.counter = Utils.Counter()\n",
    "opt.remember('xc')\n",
    "moptTD = inv.run(m0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAInCAYAAACC+3E0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4FOX2wPHvSSOE0CQiCkgHAaVfAUUgAj+kXcR7UVDQKCIiiihKEVAuihQR5VpBQUBQRL2CIFgCKkrRIAoWqhQFRWpiaElI3t8fM4HNskk2yWZnNzmf55lnd2femTm7Sd49mXnnjBhjUEoppZRSqigIcToApZRSSimlfEWTW6WUUkopVWRocquUUkoppYoMTW6VUkoppVSRocmtUkoppZQqMjS5VUoppZRSRYYmt8WQiHQWka/y0P4dERlQmDEVFSLSRkQcq68nIvEiMj4P7Y2ItCnEkJQKGO59n4h8ISJjc2jfRUTW+Cc6Z4nIXhHp53QcBZHX/i8YicguEYmzn18nIol+3n8/Ednrz33mhya3xYyICPAc8EQeVnsCeFpESuaw3fEiclZETohIsojstudJQWMuSuzPxIjIy27zI0XkmL2sujPRKVV05afvM8asBMJF5F/ZbPM2u8/LnIyInHZ5/aqIVLfnn7T7xuMistHuC8q6bMu1nes217m0MSKSJiKXucUx0l42N6+fiwpexpivjDHlvGkrInEisquwYwoUmtwWP/8HRACfe7uCMWYbsAvom0vTL4wx0UAZ4A5ghP0YdEQkvBA3vwPoIyJRLvP+DRwsxH0qVdzlue+zzQEe9LTAGLPQGBOdOQHpQBeXefe6NK9njCkNVASGAh2AjSJSwW2z9Vy3aYy5xm35TuDOzBd20n43sDWP78sxhdy/BjwRCRURzb8KkX64xc+NQLxxuTWdiPQRkc0i8reI/CkiM0WklNt6n9nr5spYvgJ+Blq4LhORgSLyk4gkicj3IvJ/bstvso9qJInIQRGZ6LLsX3acSfZjL3t+mB13T7dtzROROd7s2z6KslpEponIX8CH9vzLReQ9e/t/isgsESntsl4d+9Rmsohsdn+/2fgd2ADc7DJvIPCae0MRGSwi2+2YN4jIdS7LRERGi8h++6jvc4C4rX+liHwiIkdE5DcRmZTdF4t95OgTEUm0jy59JyL1vHg/SgWDC/o+W4yILLePkv4sIl3cln8GtPGQhOaLMSbNGLMO6AmUBR7O4yZeBwa4nBVrD6QC67JdAxCRCiIy3+7HDtr940VuzWqKyNf2Z7FRRP7hsn5Hu9/82+5P4l2WRdl95x67L/pYRGq7LP9CRJ4XkSUi8jcw0gd9dq79n9u237PbuM67U0R+tbeV7/5P7KOiYh1B/1NEDonIs5l9rZw/Kj9ARH4BTgEV7e+ux0Rkh73ftSLS3GW74SIy3d7eQREZ6bbf9iJy1u0zuUdEfrR/Tr+LyBARaQ28ivXzzTwj0N5eJ8fvCBG52v5dOCEiXwM1vflMHGeM0akYTcA3wFC3eV2Ahlj/7NQGfgEmubX5F7A/h+2Ox/riwN5OLHAaeMClzT1YR4Ab2226AieA2i5xJAPdgTCsI8Bt7GWtgTN2mzCgm/26pb18KrDEZV/R9rav83Lf44GzwHCsoztRQKS9zgSgJFAeWAHMsdcJA7YBL9nL69ivTW6fE9ALWGvPqwv8hdVpGKC6Pb8vcARoae9rAHASqGYv7w8cAprbMY8B0oDx9vKKwFFgkL28MrAReNwlHuPyGb+FlWCXAEKBRsAlTv/O6qSTL6Zs+r4v7D6nk/03dpvdr1R3a5cMdPRiH2eB9m7zqtt/Z1U8tF8IbMitnUt7A7QBNgOd7HlvAcOwkt65Oaz7MbDM7sfKAx8BH7ks3wv84dKfjAIOA2Xs5X9gHTEWu4+IdVn3LWA5cIm97n/svjDc5XP+G7jeXj+KgvfZOfZ/Ht5/V/v9hLvM+xIY5/Ie8tX/AXH2vjO/C2phnaEb7fazXQVUsuMNBZ7G+r2sab8egNXnl7fXG2dvp7a93Vfs/cTZy9sDZ13iGGz/nNrYn1kMcLVLjLvc4s7xOwLrn6+j9u9CBPAPrDOMe53+e871Z+J0ADr5+Qdu/aHE5dLmfuBbt3mdgFM5rDPe/qNLBFLsP+SXgTCXNj8Bt7uttwwYaz9fATyTzfZnAQvd5r0NzLSf18c6elHRfn0XsCMP+x4P7HZb/m/gV7d5ze33Fwpcaz+Pclk+EO+S2zDgT6x/KqYBzwBVyJrcfgpMdFt/Pec7zM+AJ12WhWAdFR5vv34EWO22/r9cOziyJrdz7c+kvtO/pzrp5OvJU9+HlXS96Tbva+Axt3kHgJu92Edek9spwE63dkl2P5o5vezSPjO5vR94F6hgt6lADsktcJm9bh2XefXseZfar/e69ScC/Abc6rL8KaCS27Zj7O1c7jIvxH4fmX3LF9gHBVzaFLTPzrH/8/AZhNo/x17261r2z6uq/Trf/R9W4uj+XXB35vtx+dm2dft8k13n2fN/BPrZz3cCA1yWlbI/szj7dXuyJre/AENyiNE9uc3xOwLrn73fAHFZPpEgSG51WELxcxzriOg5ItJJRL4SkcP2KaMpwMVu65UBjuWy7S+NNbi9NPAY1h+e67jSGsBL9umXRLGu8ozF+m8RrA5gRzbbrgrsdpv3qz0fY8xWYBOQebXvncAbedg3WJ03butc7rbOKqxOqhJWMnrIGHPKZZ092cSfhTHmLFZnOgS4HeuLyV2O79ne/7mYjTEZwD63+K91i3+OHbsnj9rxL7NPrb0gItHevB+lgsAFfZ9tr4fXVdzmedP/5UcVrCNjrhoaY8q5TPd5WG8B1gGHR4AVxhj3bbjL7DNc+6df3ZZB1v7EYCU2mZ9FT6yzUz+KyC8iMsyeX8N+3OLSzxwDwrPbtr39gvbZufV/WRhj0oE3OT9eOQ5YZYz53X5d0P7P/btgLxf+Hu11eR6DdbR6mdt7rOmynvt7PIl1tDo71cn+O9ST3L4jqgD77N+FTF59xzlNk9vi53ugQeYLEYkAlgCLsP7zLgOM5MKxS1fa6+bKGJNqjJmEdQroPy6L9gF3uXXc0caYwfbyvVidpye/c74TzVTTnp/pDSDOHuvVCpifh30DZLhtfx/Wf97l3KZIY8wBrKMAFSXrhWHuMebkNazTQb8YY7Z7WJ7bez6A1ZkB5y4sqeYWf7xb7GWNdeHLBYwxh40xQ40xtbGOSrfHuihQqaIgS9/norqH1/szX4hINawjZj/4MhgRKY+VoK7O67rGmERgKVZfPcuLVTL7jOou82q6Lcuy3O5PLsf+LIwxm40xt2Cdyh4ETBKR6zmfUNZx62uijDFvu2zbvX+FgvXZufV/nrwB3CAil2IdVDiXTPug/3P/LqiOy++RzfUzOII1zKyj23ssZYyZbLdxf4+lsD7/7Owl++9QT59/bt8RB4Bq9mebKS/fcY7R5Lb4WYJ1lW6mCKyxpceNMadFpAHWKS93nex182IscJ/95QBWGZ7xItLEHvheUqy6sFfYy18C7hWrtmSYiJQRkWvtZXOBf4lVpzJUrIs+biLrf/qLsMYm/Rf4zE5AM+W2b0+WY5UBekxEStvrVRb7Qjasi8L2AZPt7dUCHvL2wzHG7AbaYp2O82QuMMge0B8mVm3DJljDMcA6CnGPiDSzLwAYRdajsvOBFiJyl1ilxkJEpKaI3OBpZyJyi4jUsDuyJKzTX2c9tVUqCLn3fZluFJEOdr/SF2tc4SKX5Z2wxscf8UUQ9t9yK+ADrNPS0/O5qVF2bF/m1tAY8wfWMKdnRaScnVg/C6w0xvzp0vQul/7kUawzbx+JSISI3CEiMfZRvONYydJZY8whrPGqL4tIZfs9lhORXl4c+SxIn51b/+fpc9iONaZ0NtYZxg8yl/mg/wvh/HdBTayj6vNyiMUAM4BpIlLHjiHa/o7LLPX2JvCoiNQSqxTnVHK4aA7rO/QxEWlt9/cxcv6iwINYCbjr2YvcviOWYx1dflSsi9uakf33VWBxelyETv6dsP4wtuIyLgxr4P5+rMH6nwOP4zKmBmts1l9AyRy2Ox77gjK3+fG4jAPDKg32PdY4scPAJ8BVLstvtpf/jTUm9Um3ZT/ay34E/u1hfwuxhg14WpbtvnOIvyrWKcAD9n63Af9xWX4F1pdLMtZFHsPwYsxtNsuyjLm1592PdZopCfjW7ecmWP9AHMA6Dfgc1rCJ8S5tGmBVfjhob2MzcJ/Lctcxt5OxTkOetNu/hssYMp10CuYpm77vC+B5rIurTtjLu7mtt85Tf5LNPnIac3vS7icSsU7HPwmUy6bdCZdpv0ubc3+vHvad2wVlF9t92UGs/vxNIMZl+V6sGsBf2/v9jvMX7EZgXRNxxF62G3jEZd0orPG4O+33+DtWwlvK5XMem01c+e2zc+3/stnfQHt/L7rNz7b/wxp7eiKHbcZhXfw20l73sP17lXlBXebPtorbemFY1TJ+4fx33geZ7ezPfYa9vYP29neR/ZhbwRrq9ov9c/gNu7+39/U+1jCYRKCdPT+374jW9u/CCft3I0t+EKiT2MGrYsT+r+wxY0xbL9u/jTU2ydO4UKWUCgr56Ps6YyVl1+XaWBVb9lm1scYa0qACgCa3SimllFL5pMlt4AnqMbcicoNYBe53icgop+NRSimllFLOCtojtyISijUWsRPWeNEEoK8x5hdHA1NKKaWUUo4J5iO3V2MVGt5tjEnFuuqyZy7rKKWUUkqpIiyYk9vKZK3Pt5+sBfmVUkoppVQxE+Z0AAXgqdZbljEWInIPVpkroqKimterVy/fO0tPTwcgNDS0UNfLS3tv2uY37qIg0N67v+Lx9X58sb1g/vvZsmXLEWOM+x37VAHExMSY6tWrOx2G32UOA8xaE985/orH1/vxxfbyu428rpeX9t60DbTfIX8yxrBp0yav+uNgTm73k/XWflWAP1wbGGNmYd+9pVmzZmbTpk353llycjIApUuXLtT18tLem7b5jbsoCLT37q94fL0fX2wvmP9+ypQpk+0tPVX+VKtWjY0bNzodht+lpaUBEB4e7nAkFn/F4+v9+GJ7+d1GXtfLS3tv2gba75A/paWlERER4VV/HMzDEhKAOvYdRSKAPliFiINPUhL06mU9KqWUUkqpfAva5NYYcxbr7k2fYN1VZrEx5mdno8qnDz+EJUtg2TKnI1FKKaWUCmpBm9wCGGNWGGPqGmNqGWMmOh1Pvs2Zk/VRKaWKMPn1Vz1TpZQqNEGd3Aatjh1B5Py0bp01f+1aSpcpQ+kyZaz5HTsWeFfx8fGsXLmS1NTUAm9LKaV8IilJz1QppQpNMF9QFrzGjIH16+HUKet1ZuLpmoBGRcHYsQXe1d13382xY8c4fPgwMTExBd6eUkr5xJw50K+f01EoVegyMjL47bffOHnyZI7ttFoClCpViipVqhASUrBjr5rcOiE2FpYvh+7dzye4LkzJkshHH0H79v6PTSml/GHtWusMVaYOHSA+3rl4lCokx44dQ0SoV69ejklbRkYGQIHbBKuMjAwOHDjAkSNHqFixYoG2VfQ+nWARGwvvvAORkVlmm8hITs+dq4mtUqpoK4QzVUoFor///ptLLrmkSCakvhQSEsIll1xCkg/G4+sn7aTERAgLg5AQKFnSegwNRTJ/sFoiTClV1EVFgZ6pUkVYenp6saxLmx/h4eGcPXu2wNvR5NZJs2dbwxIaN4alS63HU6cIX7DAWq4lwpRSRVhKSAjJr7+uia0q8orqGFlf89XnpMmtk8qWhWeegY0boVMnSEgg5amnMJl3TNISYUqpokiEs0BqRgbjH3qIrVu3Oh2RUqoI0QvKnLRkSdbXnTsTuWqV9VwEIiKs53rhhVKqCDG1apFRogSlfv6Zbn/9RcuWLXn77bfp1q2b06EppXIRFxdHWFgYr7/+utOhZEuP3AaSMWMwJUuef12IJcKUUsoxZcoQsXkzaU8/TfRll5GcnEyPHj2YMmXKuVJHSilnrF+/nhtuuIGyZcsSHR1N8+bNmTdvntNh5Ykmt4EkNpbT776bNcF1pRdeKKWKitBQSowezT/27+fJJ5/EGMOoUaPo168fp0+fdjo6pYqlTz/9lNjYWFq3bs3u3bs5dOgQI0eOZNiwYTzxxBNOh+c1TW4DTHrbtlYpMLcSYURGWqXDNLFVShUhIsLYsWP54IMPiI6O5q233uK6665j//79ToemVLEzZMgQ+vbtyxNPPEGFChWIiori5ptv5rnnnmPixIns3bsXgNOnT9O/f3/KlClDrVq1mDt37rlt7N27l86dO1OuXDnKly9P8+bN2b59u1/fhya3AUiSki4sERYWZpUOy4mWDlNKBakbb7yR9evXU6NGDb777jv+8Y9/sGHDBqfDUsrnRMTjFBoaSmhoaLbLvW3jPnlrx44d7Nq1i34e7hx46623Yozhs88+A2Dx4sV07tyZY8eO8eqrrzJ48GDWrVsHwGOPPcbll1/OX3/9xZEjR3jjjTcoV66cbz48L+kFZQEo/M03z5cImzIFRo6EzZtzv12la+kwva2lUipAGWNIS0u7YH69evVYt24dffv25YsvvqBdu3a8/PLL3H777Q5E6Xu+qN/pS/6Kx9f78cX28ruNvK6X2T7zzmJO8Hbff/31FwCXXnrpBeuEhYURExPDX3/9hTGGVq1aceuttwLQoUMHbrrpJt544w1atWpFeHg4f/75J7t27aJ+/fpceeWVeYoju/4hL5+9HrkNQKZ06QtKhDF1KpQpk/OKWjpMKRXARKSHiMzK6Q5EFSpU4KOPPmLw4MGkpqZy99138+ijjwZcYqhUfqWnp/t18tbFF18MwIEDBy5YlpqaypEjR4iJiQGgWrVqWZZXr1793FCiqVOnUqNGDXr27EnlypUZOnQoJ06cyO/HlS965DYAnVm0iPDMWrcAoaEwfLg1uSj5z38S9sUX52d4KB12FPiscMNVSvmBiFwJNAMuAo4B3xtjfnQ2qrwxxiwDljVv3nxgTndsCg8P5+WXX6ZJkyYMGTKEGTNmsHXrVhYtWkT58uX9F3AhCbS7VfkrHl/vxxfby+828rqeN7fezTyymVNbb9rkV7169ahZsyaLFi2iU6dOWZYtXrwYEaFz585s2LCBffv2ZYlh3759VK1a9dwtdF944QUAdu/eTc+ePZk2bRoTJkzwKg4RKfDPVo/cBrHURx7JtXTYKWCif8NSSvmIiISJyAMisgv4DhgF9LEfN4rILnt5kTxQcc8997B69WouvvhiPv30U1q2bKk3fFCqkIgIL774IgsWLOCpp57i2LFjnD59mvfee49hw4YxcuRIatSoAcCGDRt4++23SU9PZ/Xq1bz//vvnhg+988477NmzB2MMZcuWJSIigrAw/3ZRmtwGsfS2bTn97rtWiTBPoqLoU7o0X/o3LKWU72wG2gD3A2WNMQ2MMa2MMQ2Asvb8Nna7Ium6664jISGBJk2asHPnTlq1asWKFSucDkupIqlLly6sWrWKNWvWUL16dWJiYpg4cSLTpk1j4sTzh8puvvlmVqxYQfny5RkwYAAvvfQSbdq0AeD777+nXbt2REdH07BhQ5o1a8Yjjzzi1/dRJP/bL07S27a1SoT17g1nzpxfYJcOW3vHHTlvICkJ4uJg7lzrdsBKqUDSJ7uhB8aYM8DHwMf2kIUiq1q1anz99dfExcXx3nvv0b17dyZNmsSIESN8di96pZSlTZs2fPrpp9kudy375cnkyZOZPHmyj6PKGz1yWxQkJuavdBhkrbCglAooxpgfRSTXQx7GmJ/8EY+TSpUqxeLFi5kwYYLe8EEplSNNbouC2bPPlw5butR6PHXKu6oJWmFBqUD3uNMBBAoRYdy4cfzvf/+jVKlSvPXWW7Rt29bj1d1KqeJLk9uioGxZ70uHdexoVVLInOyiy+cqLGROHTv69z0opbKj593d9OrV69wNHzZu3EiLFi30hg9KqXM0uS0KliyBhx+2hiPA+dJhS5Zc2HbMmKwXoHmosEBUFIwdW3jxKqXyREQuFZHLspucjs8JV111Fd9++y2xsbEcPHiQdu3aMW/evOxX0Ds4KlVsaHJb3MTGwvLlOVZY4KOPoH17v4allMpWKWA/8LuHKXN+sRQTE8Mnn3zCkCFDSE1NJS4ujuHDh3u+4YNeX6BUsaHJbXEUG2tVWIiMzDrfrrCgia1SAeUUUDObqYb9WGyFh4fz4osvMnPmTMLCwpg+fTrdunXj+PHjWRvq9QVKFRua3BZXBamwoKf3lPKnDGPMvpwmpwMMBPfccw+rVq0iJiaGTz/9lF8qV9brC5QqpjS5La4KUmFBT+8p5U96QZmX2rZty8aNG2ncuDFjT5/mlOtCvb5AqWJDk9viKi8VFtzp6T2l/GmK0wEEk2rVqrF27Vpi/v1vugEns2uo1xeoYFKEzph27NiR8ePHF+o+NLktrvJSYUHLhynlCBGpaYx5ypt2/ognWGTe8OH6CRO4BbjgNg96fYEKNn46Y9q+fXtKlChBdHQ0ZcuWpWnTprz//vuFus/CoMmtyp2WD1PKKZ+IyOsicrWnhSJytYi8hnUbXuUi84YP4x98kLPAWeCMCCYv1xcoFSjsM6XyxhuFvqtx48Zx4sQJjh49St++fbnlllvYsWNHljbGGM9VSQKEJrcqd1o+TCmnNAL2AB+ISKKIbBCRT+3H48AH9vLGjkYZwFps3ky0CNvCw+lhDFsAc/KkDqtSgS27M6br1hESGuqXM6ZhYWHcd999pKen8+OPPyIizJgxgxYtWhAVFcXGjRsBeO2117jyyivPHen99NNPz23DGMOkSZOoUqUKF110EQ899BDGmEKLOZMmt8o7Wj5MKb8zxpw2xkwELge6AguAr+zH7sDlxpinjTEXnHlXtrJlkWnTqPznn5S56SaaZWTwiDF8/+uvnDyZ7YhcpZyVzRlT8eMZ09TUVF566SXCw8Np3Nj6/3n27Nm88847nDhxgqZNmzJr1iymTJnCwoULOX78OBMnTuSmm25i165dACxYsIDnnnuOpUuXcvDgQWJiYlizZk2hxZxJk1vlvfyWDytCA+GVcoIxJt0Ys84Y86Ix5kn7ca0xJt3p2AKefX1B+QoVeO+993jplVd4OTKSZr/9RosWLdi8ebPTESp1IQfPmE6cOJFy5cpRpUoVli5dyvvvv0/t2rUBeOSRR6hVqxahoaGUKFGC//73vzz++OM0btyYkJAQunbtSmxsLIsWLQJg/vz5DBo0iObNmxMREcHo0aOpVKmSz2N2p8mt8l5+y4dp6TClVAAQEe69914SEhJo0KAB27Zto2XLlrz44ot+OVWqVJ5kc8bUFPIZ0zFjxpCYmMihQ4dYt24dPXr0OLesevXqWdru2bOHIUOGUK5cuXPT559/zoEDBwDYv39/lnVCQkKoVq1aocTtSpNb5b38lg/T0mFKqQBy5ZVXkpCQwD333ENKSgoPPPAAvXr14ujRo06HplRWbmdMnb4gMiQka9pYrVo15syZQ2Ji4rnpxIkTvPLKKwBUrlyZvXv3nmtvjGHfvsK/74wmt8p73pYP09JhSqkAFxUVxcyZM1m8eDFly5Zl6dKlNGnSxC/jAZXymvsZ00aNvL/hkh889NBDjB8/nh9++AFjDKdPn+brr79m27ZtAPTv359Zs2axadMm0tLSmDx5MgcPHiz0uDS5Vb6npcOUUkGid+/e/PDDD7Ru3Zr9+/cTGxvL+PHjA7rMkSpG3M6Ymm+/xUyZ4t0Nl/xg4MCBjBgxgjvvvJPy5ctz+eWX8+STT5KWlgbA7bffzgMPPECPHj245JJLOHToEG3bti30uKS4jDNq1qyZ2bRpU77XT05OBqB06dKFul5e2nvTtkKFChw7dozDhw8TExPjVQw+8fnn0L279R+mOz+VDsvvz6yw+CseX+/HF9sL1r+f5ORkypQp850xpoVXAfiYiFzmTTtjzB+FHYsvNW/e3Hz33XdOh5FFWloa48ePZ9KkSRhjuO6661i4cCFVq1b16T4AwsPDfbbNgvBXPL7ejy+2l99t5HW9tLQ0du7cSYMGDXJtm5GRAVx46j+vbYLd1q1bqV+//gXz09LSiIiI8Ko/LrqfjnKWlg5Tyhf2A797MakCCg8PZ+LEicTHx3PppZfy1Vdf0aRJE5YuXep0aEqpPNLkVhWevJYO05JhSrmrAdS0p0HAN0A3oAFWndv19nzlI9dffz2bN2+ma9euHDt2jBtvvJH777+fM2fOOB2aUspLmtyqwpPX0mFaMkypLIwx+zInYDhwkzHmY2PMdmPMSuDfwCPORln0XHzxxSxfvpzp06cTHh7OSy+9RMuWLdm6davToSmlvKDJrSo8eS0dpiXDlMrJpUCy27wT9nzlYyLCQw89xPr166lduzZbtmyhRYsWzJ49W2viKhXgNLlVhSe30mFaMkypvPgaeENEqgCISFXgdXu+KiTNmzdn06ZN9O/fn1OnTnH33XfTt29fknT4lFIBS5Nb5RwtGaZUXtwNVAR+E5E0YC9QyZ6vClHp0qWZP38+8+fPJzo6mnfeeYemTZvyzTffOB2aChJ6tN87vvqcNLlVznHw3tlKBRtjzEFjTHusi8zaATWNMe2NMX/6Mw4RqSkis0XkPbf5pUTkOxHp7s94/Kl///5s2rSJZs2asWfPHtq0acOUKVPOlWdSypMSJUpw9OhRTXBzYYzh6NGjRLpXWcqHMB/Eo1T+ZZYM690bXK9G1pJhSnlkX1yWr/tXisgcrCoLh4wxV7rMvwGYAYQCrxtjJuew/93AAPfkFhgJLM5PXMGkTp06rFu3jscee4zp06czatQo4uPjefPNN6lUqZLT4akAdNlll3Hw4EEOHz6cY7vM5FdECtQmmEVGRlKlSpUCb0eTW+U815JhJUpASornkmFJSRAXB3PnWherKVXEichKINfDPcaYrl5uci7wIjDfZR+hwEtAJ6y6ugki8iFWojvJbf27jDGHPMTZEfgF8OqQizHmXDH8YBQSEsLkyZNp164dAwYMID4+nkaNGjFnzhw6d+6c7XqBdtczf8Xj6/34Ynv53UZe18ts703Cltlbn+w8AAAgAElEQVQ2LCz71MybNkWBp/4hL5+9DktQzvO2ZJiWClPFz0bgOy8mrxhj1gDH3GZfDewyxuw2xqQCi4CexpgfjTHd3aYLEltbLNAKuBUYKCIXfLeIyD0islFENh45csTbkANaly5d2LhxI7GxsRw+fJgePXowatQoUl2vG0hKIvTf/9b63Ur5UdFO/VVwyCwZNmyYdfT2+uvh+efhq6+ytnMtFdavn//jVMrPjDHj/LCbymS9y9l+oGV2jUWkAjARaCoio40xk4wxY+xlccARY8wFg1CNMbOAWWDdfjdQbkFbUNWqVeOzzz5j6tSpjBs3junTp7NmzRrefvttateuDStXwocfEvLxx3DLLUDg3H43k7/i8fV+fLG9/G4jr+vlpb03bQPtdyjQ6JFb5bzsSoadOKGlwpRyISKlReRmEXnYfl1RRAo60NPT4L1sh0IYY44aY+41xtQyxkxyWzbXGLO8gPEEndDQUEaPHs1XX31FtWrV2LhxI82aNWPhwoVav1spB2hyqwKXlgpT6hwRaQrsAsYDE+zZTbHG0BbEfqCqy+sqwB8F3Gax1Lp1a3ZWr44B/k5O5rZ+/Uhbs8ZauHYt4RERhEdE6D/lShUyHZagAldmqbDu3a0xuO60VJgqXp4HxhpjXhOR4/a8tcDsAm43AagjIjWAA0AfrLGzhSbYLyjLiYwZg0lIQOw+KzyzTJjLP+UmKor0UaMwDn8GekFZ4a+Xl/betA20ixL9SS8oU0VHZqkw97p3WipMFT9XYd2RDOxhA8aYE0C0txsQkbeB9UA9EdkvIgOMMWeB+4FPgK3AYmPMzz6N/Pz+e4jIrKJ8dy/Tvj3pS5ZgsqnfbUqWJH3pUky7dn6OTKniQ4/cqsCnpcKUAjiMNXzgt8wZIlKTPAwhMMb0zWb+CmBFQQP0Yv/LgGXNmzcfWKQviOnUyWP97tPA49WqcV+tWtQIoPevF5QV/np6QZl/6ZFbFfi0VJhSAG8Cb4tIK0BEpDHwGgUflqAKg+s/5SVLkhESQroIf27bRuPGjXnzzTf1jlVKFRJNblXgyywVtnGjdUQkIQGmToUyZbK206uSVdE2GWuM7SqgrP38W+A5J4NS2XD7pzykcWNKASNiYkhOTub222+nT58+HDvmXnZYKVVQmtyqwJfPUmGly5ShdJkyelWyKhKMMWeNMSOwxtheBpQxxoz2VFNWBQAP/5RnTJ5Mw9ateeONN4iOjmbx4sU0atSIVatWOR2tUkWKjrlVwWvMGFi//nwlBQ+lwkzJkoiWClNFiLHOZR90Oo6CKMrVEs55913rMT3dmoCzDzwADzzAbWFhtGrVijvvvJMNGzbQsWNHhg0bxoQJE4h0v3i2EGm1hMJfT6sl+I5WS1DFQ2apsByuSj793ntaUUEFLRFJcHm+VUR+8TQ5GWNeFIdqCd6qVasWq1ev5oknniA0NJTnn3+ea665hh9//NHp0JQKenrkVgW3zFJhblclExnJ6blzSb/uOudiU6rgXnB5Ps2xKHyk2FRLyEXmew8PD2f8+PF07dqV2267jZ9++olrrrmGyZMnM3ToUEJC/HP8SaslFP56Wi3Bv/TIrQp+blclExICYWGI+9GhpCTo1ct6VCo4XOXyfKcxZranybHolE9cffXVfP/99wwcOJCUlBQeeughOnfuzIEDB5wOTamg5HhyKyK9ReRnEckQkRZuy0aLyC4R2S4inV3m32DP2yUio/wftQoo2ZQKC1+wIGs7LRWmgs89Ls+XOxaFKnTR0dHMmjWLJUuWEBMTQ3x8PFdddRXvvfee06EpFXQcT26Bn4CbgDWuM0WkAdZtIBsCNwAvi0ioiIQCLwFdgAZAX7utKq6yKRVmSpfO2k5Lhangs09EnhSRm4Fw+2DAze6T00Eq3+nZsyc//vgjXbp04fjx4/Tu3Zs777yTv//+2+nQlAoajo+5NcZsBRAR90U9gUXGmBRgj4jsAq62l+0yxuy211tktw2aiyqc8OeffxIeHk6pUqUIC3P8x+5bS5ZkfW2XCpPly61SYJkiIqxHu1TYOR06QHx84cepVN4NwKpv2w8IB6Z7aGOAxf4MqqCKRbUED7y92rtChQosWbKEmTNnMmLECObOncuXX37JG2+8wTXXXOP3eAJtP1otoXjKy3sP5CynMrDB5fV+ex7A727zW+a2sfT0dJKTk/MdTH7Xzet6eWmfl7aNGjU69zwyMpJSpUoRHR1NdHQ0pUqVOvfa9Xnma0+P7m0DMWFOue8+LvrmG0JOn7ZmZFMq7PTDD5NegN8NbxXk98/J/fhie8H69+Ovn1l2jDEJQAcAEUk2xlR1NKACEpEeQI+aNWs6HUrAExHuvfde2rVrR1xcHN9//z3XX389I0eOZOzYsXpBkVI58EtGIiLxQCUPi8YYY5Zmt5qHeQbPQyk83sNQRO7BHrNWuXJlT02KvLvuuov333+fM2fOcOLECU6cOMGZM2c4c+YMR48e9dl+SpQoQalSpShdurTXyXJu8wvaeadeey3H5s2jwh13IJkJrovMUmFaUUEFA2NM6dxbBTatlmDJy3tv1KgRGzZs4IknnmDKlClMmjSJ+Ph4FixYQN26df0eTyDtR6slqOz4Jbk1xuTnFlH7AdejFFWAP+zn2c133+8sYBZAs2bNTGn3MZj5kN9t5HW9vLTPqe3jjz/O448/fq6NMYbTp0+fS3Qzp+Tk5Avm5TTffVlKSgopKSk+vZVkiRIlziW/rlPp0qU9zndfJiJEV69OxksvcfHgwUhKyvmNR0YiixcT1bWrz+L1li9+D53YT3H8+1EqEERERDBp0iS6dOlC//79SUhIoGnTpjz33HMMHDjQ07A+pYq1wDuXfN6HwFsiMh3rVpN1sO6jLkAdEakBHMC66OxWx6IMMiJCVFQUUVFRVKxY0SfbNMZkOTKcl6Q4u/nJycnnEuaCHmG+DXgFKAmkAhHAmZQUnhw0iC+qVLkgWY4JD+e2zz4jvl8/Ii6+ONckOiJzLK9SShWitm3bsmXLFu6//34WLFjAoEGDWL58Oa+//vqF/XlSEsTFwdy51kW3ShUjjie3ItILq1D5xcBHIvKDMaazMeZnEVmMdaHYWWCIMSbdXud+4BMgFJhjjPnZofAVVsJcsmRJSpYsycUXX+yTbRpjSElJKVCynJiYyMmTJ7lv715KpaayGRgBTAEaG0PnP/5g6h8XHvTvB9QHJo4bx0IvYg0PD8/1qHLp0qXPtatQoUKuR6IjIiL0aIxS6gJly5blzTffpFu3bgwePJhly5Zx1VVXMWfOHLp163a+oWvpw379nAtYKQc4ntwaYz4APshm2URgoof5K4AVhRyacpCIEBkZSWRkJDExMfnaRubFQKX794e2bWny4IMsP3uWE0lJJM2YQYv161n/9NMXJMVdn3kGdu/miapVKdWlS46JdHJyMmlpaRw/fpzjx4/77P2HhYV5PfzCfX5ISAilSpXikksuybKsRIkSmjArVUT06dOHa6+9ljvuuIPPP/+c7t27M3jwYKZNm0ZUVFTW0oea3KpixvHkVqlCZ5cKE6BEaCglKlaEidb/TK0AOnaEVavOt7eHGdT56y9mzpp1fn6HDvDttxdsPjU11asjy0ePHuXkyZOkpqbmetQ5NTW1UBLm/IxhDg0NpWLFinTo0EGTYweJyCxggTFmjcu8dkBfY8y9zkWmnFK1alXi4+OZPn06TUeMoMMrr8Arr1gLtfShKsY0uVVqzBhYv966yxl4LBlGVBSMHetx9YiICCpUqECFChVy3M25I8leXMDkmgDn9SK/pKQkTpw4ccGFgykpKSQmJpKYmJj7Z+LB66+/zoABA/K1rvKJnsBQt3nfYtW4DarkVuvc+taDDz7I7uhoTt93HyWNXTzIU+nDqCjSR43C2J+91rkt/PW0zq3vFJU6t0r5R2wsLF8O3bufT3BdRUXBRx9B+/Z+CykiIoKLLrqIiy66KM/rZpdEp6amcvLkyTxf5Pftt9+yY8cOvc+988KBdLd5Z4ESDsSSL1rntvDUHDCA1CpVSO3ViwgPSYCJiiJ96VJMu3YORKeUf2lyqxRYCe4770Dv3nDmzPn5kZHWfD8mtoUlIiKCiIgIypcvn6f1Ro4cydSpUwspKpUHvwC9gbdc5v0L2OpMOHmndW4thfXew7t3hw8+IP1f/yLU5YhtakgIh559liodPVfl1Dq3hb+e1rn1L083RFCqeEpMhLAwCAmBkiWtx7Awa743kpKgVy/rUSnfGwfMFpGFIvKEiCwA5gCex8uo4ikxkdCICExICKlhYZwFUjIyeOy++7j77rv5/fffc92EUsFOk1ulMs2ebQ1LaNwYli61Hk+dOn/VcW5cS+8o5WPGmM+Ba4AT9uNJ4FpjzKocV1TFi92PSePGRKxYQUbDhpQC7jSG2bNnU6dOHYYPH86RI0ecjlSpQqPJrVKZypaFZ56BjRuhUydISICpU6FMGe/Wdy29o1QhMMZ8b4wZZNcCH2SM+d7pmFSAcevHIjZvJmTaNFp06MAtt9xCSkoK06dPp2bNmjz55JPnxugrVZRocqtUpiVL4OGHreEIAKGhMHz4uVJiF+jY0SqxkzmtW2fNzyy9kzllM85NqdyISHOX51dnNzkZowow2fRjpePjWbRoEZs2baJLly4kJyfz5JNPUq9ePZ5//nnOuF5roFSQ0+RWqfwaM8aqpJApjyXElPLCFy7PN2Qzrfd/WCpYNW3alBUrVvDll1/SunVrjhw5wkMPPUTdunWZM2dOsS41pYoOrZagVH4FYAkxVeSUdXleZC6P1jq3zmvdujXx8fGsXLmS//znP/z4448MGDCAqVOn8p///IdevXr57KYtWufWd20D6XfI3/Ly3vXIrVIFkVlCLDIy6/wiVEJMOccYk+Hy8mZjTLr7hFUeLCiISA8RmZWkFUUCgojQtWtXEhISmDdvHjVr1mT79u306dOHa665hvj4eEzmTSGUCiJ65FapgnItIVaiBKSk5K2EmFLemQm87WH+y8AiP8eSL1rn1hJo7z08PJzbb7+dPn36MHv2bCZMmMB3331H165diY2NZdKkSbRs2dIn+/ElrXOrsqNHbpUqqIKWEMuJ1s5V511wjlhEqnDhXcuUypeIiAgGDx7Mr7/+yuTJkylXrhyff/45rVq1olevXvz8889Oh6iUVzS5VaqgClpCLCdaO7fYE5HTInIKiBKRU64TsA941+EQVRETFRXFyJEj2b17N6NHjyYqKoolS5bQqFEj4uLi2Lt3r9MhKpUjTW6VKqi8lhDLC62dq6A78E8gBejhMnUDrjTG3OdgbKoIK1++PE8//TS//vorQ4YMITQ0lHnz5lG3bl2GDh3KX3/95XSISnmkya1SgURr5yo3xphVxph4oJ79PHP63Biz1en4VNFXqVIlXnzxRbZt20a/fv04e/YsL7zwArVq1WLcuHHoBYIq0Ghyq1Qg0dq5KhvGmN9FpLWIvCwiHwCISFMRaeN0bKp4qFmzJm+++SabN2/mn//8JydPnuSpp56iZs2aPPPMM5w+fdrpEJUCNLlVKrBk1s51TXBdae3cYktE+gAr7Jcd7McQYIIzEani6qqrrmLp0qWsW7eOdu3acezYMUaMGEGdOnV4Y8YMMnr21ItglaM0uVUq0GjtXOXZWKCzPcY2s0LCj8CVzoWkirPWrVvz+eef8/HHH9OsWTMOHDjA6mHDCPnwQ7ZMnOh0eKoY0+RWqUDkWju3ZEnrUWvnFneVjTHf2s8zK+ufBUIdikcpRITOnTuTkJDA4sWLub9kSQCOPPMM/fr149ChQw5HqIojvYmDUoHItXbulCkwciRs3mxVTejXz+nolDN2i0grY8wGl3mtgJ1OBZRfevvdwOCLeEJvuIGQ1asJwbpVnomIAOBa4PqFC2HhQgDSY2NJ++STAu8P9Pa7xZXeflepYFeYtXNVsHoaWCoio4BwEXkQ645lTzkblvf09rtFT8aoURiXawTEvvi1hEubk8DDR4+yY8cO/wanii09cqtUIHKvkZtZO3f4cN9sPykJ4uJg7lwrkVYBzxjzvn3jhqHAAaArcK8xZqWzkXlPb79rCbT3XqB4OnWyLoLt3t062+TmbEQEt5UsydItW5jZvDljx45lxIgRRNhHeAtCb7+rsqNHbpUqjvTOZ0HJGLPSGNPFGHOFMaZzMCW2qgjL4SLYsPff59Vt24iLiyMlJYVx48bRtGlT1q5d60ysqljQ5Fap4kjvfBaURCRCROqLyNWuk9NxKZXTRbAVKlRg1qxZrF69mjp16vDLL7/Qpk0b7r33XhL1IllVCDS5Vao40DufBT0R6Q78CfwMbHCZ1jsZl1JA1otgly61Hk+dyvIPdGxsLFu2bGHcuHGEh4czc+ZM6tevz7vvvosxJoeNK5U3mtwqVRzonc+KgunAJKAcEO4yFXzwolIF5eVFsJGRkUyYMIEffviBa6+9loMHD3LzzTfzz3/+k99++82h4FVRo8mtUsWB3vmsKKhkjJlmjPnbGJPuOjkdmFIsWQIPP2wNR4DzF8G6Xxxra9CgAWvWrOHVV1+lbNmyLF++nAYNGvD888+Tnq6/0qpgNLlVqrjQO58Fu3gR+YfTQSjlKyEhIQwaNIitW7fSu3dvTp48yUMPPUTLli35/vvvnQ5PBTFNbpUqTvTOZ8FsJ7BcRJ4TkRGuk9OBKVUQl156KYsXL2bZsmVUrVqV7777jn/84x888sgjnDx50unwVBDS5Fap4sSLiz5UwLoG2AG0AHq4TN2dDEopX+nevTu//PILw4YNwxjDs88+S8OGDVm5UiveqbzR5Fap4kTvfBa0jDHXZTO1dTo2pXwlOjqa5557jm+++YYmTZqwb98+unbtSt++ffnrr7+cDk8FCU1ulSpO8njRh1JKOaFFixYkJCQwbdo0oqKiWLRoEVdccQWvv/46GRkZToenApzeflcppYKAiKQBnoqBpgD7gLeAZ4wxaX4NLB+MMaSlBXyYPnf27FmnQ8jCX/EUZD9Dhw6lR48eDB06lE8++YSBAwcyd+5cXnzxRRo2bOj3mPK6Xl7ae9M20H6H/Ckv712P3CqlVHB4FOuisvuxxtrejzUG9z/AK8BA+3nAEpEeIjIrKSnJ6VBUEKlRowYffvghb775JhUrVmTt2rW0bNmSCRMmkJKS4nR4KgDpkVul1HlJSRAXB3PnWuNzVSDpD/QwxuzOnCEiq4DFxpjmIrIW+AB4zKkAc2OMWQYsa968+cDw8HCnw3FMoL13f8VT0P3069ePrl278uijjzJnzhyeeuop3nvvPWbOnEnbxo3z1XflN6a8rpeX9t60DbTfoUCjR26VUud9+KE1/nbZMqcjUReqA/zuNu93ez7GmB+Ai/0dlFL+dNFFF/Hqq6+yatUq6tWrx7Zt22jXrh2ze/bUvkudo8mtUuq8zJJgWhosEH0PPC0iEQAiEg48BWy2X1cHjjkVnFL+dN1117F582aeeOIJwsPDqfnllwAcnDQJYzwNTVfFiQ5LUKo469gRVq06/zoiwnpcuxZEAJgCdATW+z045eYeYDkwWET+Ai4B/sQafwtQBRjtUGxK+VXoDTcQvno144HxQKoIGEP5X35BQlyO23XoAPHxzgSpHKNHbpUqzsaMgaio869TU7M+AqlhYUz0c1jqQsaY7UB9rJs2jAO6AfWNMdvs5V8bYxY6GKJSfpMxalSWvivCPlpbwqVNang4aSNH+jkyFQg0uVWqOIuNheXLsya4rqKieOPf/+ZL/0alsmGMOWuM+cIY85Yx5ktjTPGtC6SKNdO+fY5910ng/9LSaPzgg6xZs8avsSnnaXKrVHEXGwvvvAORkVnnR0bCO++w+/LLnYlLXUBE7hSRBSLyiYh8mjk5HZdSjsih79r11FP8WbcuW7dupV27dtx1110cPXrUmTiV32lyq5SCxEQIC7PuXFaypPUYFmbNVwFBRJ7EGgJ9GLgO+BloAvziZFxKOSqbvqtxtWps3ryZ8ePHExERwRtvvMEVV1zBvHnz9IKzYkCTW6UUzJ4Np05B48awdKn1eOqUVk0ILP2A/zPGPASk2I83Yl1IplTxlEPfFRkZyRNPPMGWLVuIjY3lyJEjxMXFcf3117Nt2zanI1eFSJNbpZRV9PyZZ2DjRujUCRISYOpUKFPG6cjUeRXsWrYAGSISYoxZB3RwMiilHOVF31WvXj1WrVrF/PnziYmJ4YsvvqBFixaMHz+eM2fOOBi8Kiya3CqlrOLnDz9sndIDCA2F4cOt+SpQ/CEimQOgdwOdReQfgF5UpoovL/suEaF///5s376du+++m9TUVJ5++mmuuuoq4rVUWJGjya1SSgWHV4EW9vPnsWrebgBeciwipYLMRRddxGuvvcbq1aupX78+u3btolOnTvTr149Dhw45HZ7yEU1ulVIqCBhjnjfG/M9+vhCoBTQ2xox3NDClglCbNm1ISEjg6aefJjIykoULF1KvXj1mzZpFRkaG0+GpAtLkVimlgpAxZq8x5ien41AqWEVERDB69Gh++uknOnfuTGJiIoMGDeK6667jp5/0TyuYaXKrlFJBQEQqisiLIvKNiPziOjkdm1LBrFatWqxcuZJFixZRqVIl1q1bR9OmTRk1ahSnTp1yOjyVD5rcKqVUcFgIXAksAJ51m5RSBSAi3HLLLWzdupXBgweTnp7OlClTaNiwIStXrnQ6PJVHmtwqpVRwuBroYox5wRgz23VyOjCliopy5crx8ssvs379eho3bszevXvp2rUrN998M3/88YfT4SkvaXKrlFLBYQdQ1ukglCoOWrZsycaNG5k2bRpRUVG8++671K9fn1deeYX09HSnw1O5CHM6AKWUjyQlQVwcvPCCVdhcFTX3AK+KyBzgoOsCY8y3zoSUP8YY0tLSnA7D786eDaySxP6Kx9f78cX2vN3G0KFD6dmzJ8OGDeOjjz7iwQcfZP78+bz88ss0bdoUkpIIHTCA9NmzPfa7eYnVm7aB9jvkT3l573rkVqmi4sMPYckSwnR8WFFVG4gFlmDVt82c1jsZVF6ISA8RmZWUlOR0KEp5rVq1avzvf/9j8eLFXHbZZXz33Xe0bt2aRx99lNT33iPkww+R5cudDlO50CO3ShUVc+YAEL5gAWf79HE4GFUIpgOPAfOBoLyE2xizDFjWvHnzgeHh4U6H45hAe+/+isfX+/HF9vKyjd69e9OhQwfGjx/PSy+9xIwZM7i5RAmuAcLmz7fOnPlgP960DbTfoUCjR26VClYdO4LI+WndOgBCN2ygdJky5+d37JjvXRhjOH36tK8iVgVT2hjzkjEm2RiT7jo5HZhSxUW53r357wsvkJ6RgQGap6QAkPbll1n74wL0u6rg9MitUsFqzBhYvx4y6zCmpgIg9iMAUVEwdmyumzpz5gw7d+5k27ZtWabt27dz8uRJAEJDQ33+FlSefCAinYwxnzkdiFLFVcaoUciGDYjd75aw54e73tXMy35XFR5NbpUKVrGxsHw5dO9+PsF1FRUFH30E7dsD1lHYw4cPX5DAbtu2jb1792KM8bibmJgYGjZsSO/evQvxzSgvfSAinwF/us40xtznUDxKFSumfXvSlywh7MYbPfa7J4G3e/XijmuvRQcOOEeTW6WCWWwsvPMO9O4NZ86cm50REUHCww/z5bffsm3+/HNJ7PHjxz1uJjQ0lJo1a3LFFVdkmerVq0dERAQApUuX9stbUtkKAz6wn+sPQymHmPbtPfa7aWFh9Dl7luULFzJrxw4WLlxI9erVHYuzOHM8uRWRZ4AeQCrwK3CnMSbRXjYaGACkA0ONMZ/Y828AZgChwOvGmMlOxK6UkxITE9m+fTupy5ZxdXo6YUCqCOHGcDo1lReeeoqFbuuUKVPmguT1iiuuoFatWpQoUcLTbkhOTi7096JyZ4zp73QMSilbYiKEhUFICJQoASkphEdG8sywYWyZP5+EhASaNGnCc889R1wOF5qpwuF4cgt8Bow2xpwVkSnAaGCkiDQA+gANgcuAeBGpa6/zEtAJ2A8kiMiHxhi9v7oqcjIyMvjtt9/Yvn37BUMJDh60Sp2uxvpD3gyMNIYpQBPgkfLlqdC/f5ZktlKlSoiIc29IKaWKgtmzrWEJjRvDlCkwciRs3swVa9eyefNmhgwZwltvvcWgQYNYsWIFs2fPpkKFCk5HXWw4ntwaYz51ebkB+Lf9vCewyBiTAuwRkV1Yt58E2GWM2Q0gIovstprcqqB16tQpduzYcUECu2PHjmyrFURGRlKvXj2iDx/msyuvJPGOO3imQQMqXXwxqfPm0eTbb5kxY4af34lSShUDZcvCM8/AsGHW0dvrr4fnn4evvqJcuXIsXLiQrl27ct9997F06VISEhKYN28eHbWKgl84nty6uQt4x35eGSvZzbTfngfwu9v8lrltOD09vUCnV/O7bl7Xy0t7b9oW51PKgfbe//77bw4dOsSff/7Jjh07zk07d+7kt99+y3a9ihUrUrduXerWrUudOnXOPa9atSohIRdW80tOTuZYXBylH3gAfPAZ+OJz1L8fpVSRsmRJ1tehoTB8uDXZbrvtNq6++mri4uJYt24dnTp14uGHH+bpp5/OdhiY8g2/JLciEg9U8rBojDFmqd1mDHAWzg0T9HTu1OC5Nq/Hy7xF5B6sW1ZSuXJlT02U8rnU1FT27NmTJYHNnLJLlsLCwqhZs+YFSWydOnUoV66cn9+BUkopX6hevTrx8fE8++yzjB8/nunTpxMfH89bb71Fw4YNnQ6vyPJLcmuMyfE4vIjcAXQHOpjz9Yj2A1VdmlUB/rCfZzfffb+zgFkAzZo1M7642ju/28jrenlp703b4nyle2G992PHjmWpB5v5/NdffyU93XNd/XLlylG/fv0sF3NdccUV1KxZ0+d3nPH1+y7Ofz9KKZVfYWFhjB07lubZ6QwAACAASURBVE6dOnHbbbexZcsWWrRowdSpU7n//vv1OohC4PiwBLvywUignTHGtWjch8BbIjId64KyOsC3WEd064hIDeAA1kVnt/o3alVcpKens2/fPo+1YQ8fPuxxHRGhRo0aF5TVqlKlCjExMZQpU8bP70IFKxHZSjZnplwZYxr4IRylVAG0bNmSH374gWHDhjF79myGDh3KypUrmTNnDpUqeTq5rfLL8eQWeBHrJh+f2f+9bDDG3GuM+VlEFmNdKHYWGJJ5m0kRuR/4BKsU2BxjzM/OhK6KihMnTpw7+up6FHbHjh2k2LdXdBcVFeWxLmydOnUoWbLkBe11/KbKh2lOB6CU8p3o6Ghef/11unTpwsCBA1m5ciWNGjVizpw5dO7c2enwigzHk1tjTO0clk0EJnqYvwJYUZhxqaLHGMMff/zh8Sjs/v37s13vsssuuyCJveKKK6hcubLHC7qU8hVjzGynY1BK+d6//vUvWrVqxR133MGqVavo0aMHgwYNYsqUKZQtW9bp8IKe48mtUr6WkpLCzp07+f7779mxYwd79+49l8SeOHHC4zoRERHUqVPnggS2bt26OoxABQwRCQFqAhfjctGtMWadY0EppfKlcuXKfPrppzz33HOMHj2amTNn8sUXX/D222/TtGlTp8MLaprcqqB15MgRj0dh9+zZQ0ZGhsd1YmJislzIlTlVr16dsDD9c1CBS0SaAO8DNbDG4Qrnx+OGOhWXUir/QkJCGD58OB06dKBv375s27aNli1bMnHiRIYPH65nB/NJv81VQDt79ix79uzJkrxmjok9evSox3VCQkKoXbs2tWvXpm7dujRq1OjceNiYmBg/vwOlfGYG8BHwOLAHqA5MBT53MCallA80adKEDRs2MHr0aF555RVGjBjBxx9/zLx586hSpYrT4QUdTW5VQPj777893mJ2586dpKWleVyndOnSWS7kynxeu3ZtSpQoce4CLi31pIqIRsD/GWNSRESMMUki8giwCVjkcGxKqQKKiopixowZdOvWjbvuuovVq1fTqFEjZs6cSe/evZ0OL6hocqv8JiMjg/3791+QwG7fvp0//vBYqhiAqlWreryg69JLL9X6gKo4cf0vL0lELgaSgEsdikcpVQi6devGli1buOuuu1ixYgU333wzcXFx/Pe//yUyMtLp8IKCJrfK506fPs3OnTs9JrGnTp3yuE6JEiWyHH3NfF63bl2io6P9/A6UCkibgA5YlWLWAPOAU8BPTgallPK9Sy65hOXLl/Pyyy/zyCOPMHfuXNasWcO8efNo2bKl1SgpCeLiYO5c0AoLWWhyq/LFGMOhQ4c8XtC1b98+zt9oLquKFSt6PAp7+eWXExpaRK6J0Q5HFY6BnL/9+MPAFKAMcKdjESmlCo2IMGTIEGJjY7n11lvZvHkz7du3Z8yYMYwbN46wDz+EJUtg2TLo18/pcAOKJrcqR2lpaezevdtjEpuYmOhxndDQUGrXru3xBgfly5f38ztwgHY4qhAYY353eX4YuMvBcJRSftKgQQO++eYbxo4dy7Rp05gwYQLx8fHEp6dTEmDOHP2ucaPJrQIgMTHRYwL766+/cvbsWY/rlC1b1uNR2Jo1axIREeHndxBA5sz5f/buPc6qut7/+OszMzAwcpHAS3EbPaHmDRVFj7/SGQW10tDOkUSUzAtJksek8jKQJqImVtZJPY6BUHEEyuJuhCVqggV5OqkZhRdg1PRoMTKOwDB8fn+sPbBnmMteM3uvvfbe7+fjsR979nd9L5817D185jvf9V17n/UDR9LEzG4BfuXua5PKTgVGuftt2YtMRDKttLSUmf/zP8xsKlizhqZ7Z/ozzzS//uTMM+HxxyOOMF6U3BaQ3bt3s2nTplZ3JXjrrbdabWNmlJeXt5rEHnjggbqgC2DUKPj1r/e+bkrsn3kGkr4/PSsq+GDJkoiDkzwyEfhOi7LngZ8CSm5F8l1VFb52LZa4dqU0UWw7d+6tU1YGU6dGH1vMKLnNc4899hgLFixg48aN/PWvf2X79u2t1uvZs+c+W2odccQRDBs2jLKysoijzjFVVbB2LTRdLNf0g6bFD5ydX/ta9LFJPtkPeL9F2fuArrgUKQSVlTQuWkTx+efvSXCT7SwpYeeCBfSqqIg+tphRcpvnvvKVrzTbZuvDH/5wq7OwgwYN0p1QOquyEpYtg3PP3ZvgJisrg+XLaRwxIvrYJJ9sBEYDK5PKzgReyU44IhI1r6igcd48SsaNg6TJqg+AC3ftYv2VV/Kd73yHcePGFfRfVpXc5rmmmdpf/epXjBw5kr66ej8zKithwQK48MJmP3Do0SMor6iAxE0lRDrpLmChmf0A+CswDLgG+GJWoxKRaG3dCiUlUFQEpaWwYwfde/TgpIEDWf63vzF+/Hh++MMfcv/993PEEUdkO9qs0FRdgTj++OOV2GZa8g+cnj2D55KSoFyki9z9Z8AlwInAzcBJwOfdfWFWAxORSBU9/HDwV8Lhw2HxYhg+nOLt2/nGwIHMmjWL/v3788QTT3DsscdSVVXV5v7y+UzJrUi6zJq1zw8c6uv37p4g0kXuvtTdz3b3wxPPkV+haGaHmtksM/tZUlmFmT1tZv9lZhVRxyRSUPr2hZkzYf16GD0a1q2Du+/G+vbl8ssvZ8OGDVx11VU0NDRwxx13cOSRR7J06dJsRx0pJbci6dLGDxz69Ml2ZJKjzGxA0tcHtvUI0d9sM3vbzF5oUX6OmW0ws41mdmN7fbj7K+5+RctioA7oAdSkGo+IhNf46KNw/fXBXwcBiothypRgf3Wgf//+VFdXs2bNGoYPH86mTZv4zGc+w5gxY3jttdeyF3iEtOZWJF0SP1j2aPqBM2VKduKRfPAKwV3IAP5OkEQms0RZqrf3mwP8APjRng7MioH7CC5WqwHWmdmSRJ93tmh/ubu/3Uq/T7v7k2Z2EMF2ZePbC8LdaWhoSDHk/NHWnuHZElU86R4nHf11to+w7cLUT6VumP5OPPFE1q5dywMPPMCtt97KkiVLWLVqFVVVVVx33XU5tx99mHPXzK2ISHwNT/p6GHBYi0dTWUrc/SngHy2KRwIbEzOyO4H5wBh3f97dz23xaC2xxd13J778J3u332zGzCaa2XozW//OO++kGrKIdEFJSQlf/vKXef7557nwwgv54IMPmDp1KiNGjGD16tXZDi9jNHMrIhJT7v5q0suD3f2ZlnUSdyl7uQvDDAS2JL2uAU5uq7KZ9QdmAMeb2U3ufqeZfRY4G9ifYGZ4H+5eDVQDjBgxwrt169aFkHNb3M49qnjSPU46+utsH2HbhamfSt2w4w8dOpSFCxeyatUqrrnmGjZs2MBZZ53F+PHjueeeezj44IND9Rd3mrkVEckNj7VRvqyL/ba2GWbL5Q97D7i/6+5Xu/u/uPudibKfu/sX3f1z7r66i/GISIaMHj2a559/nunTp9OjRw/mzZvHEUccwX333UdjY2O2w0sbJbciIrlhnyTUzPYDdrdSN4waYHDS60HAG23UFZEcV1paytSpU3nxxRf51Kc+RW1tLZMnT2bkyJGsW7cu2+GlhZYliIjEmJm9RDCT2tPM/tzi8MHAE10cYh0wzMwOAV4HLgIu7mKf7dIFZfGgC8oy3y6bF5R1ZPDgwfziF79g8eLFXH/99Tz33HOcfPLJTJw4kdtuu41+/fqlbax00AVlIiL54x6CHQgagG8nPWYS7EpwUaodmdkjwFrgcDOrMbMr3H0XMJngtr4vAQvd/cX0nsKe8c8zs+ra2tpMdC8iIZkZ559/Pn/605+YMmUKxcXFPPjggxx99NH8+Mc/xr3NFUqxpplbEZEYc/dZZlYCHAT8xN13dKGvcW2UrwBWdLbfEOMvBZaOGDHiqrhdVBWluJ27LijLfLtsX1DWkX79+nHPPffwhS98gS996Us89dRTXHHFFcydO5f777+fo446Kq3jZZpmbkVEYi4xu3pjVxJbEZGOHHXUUaxevZq5c+dywAEH8NRTT3Hcccfx9a9/nbq6umyHlzIltyIiueEPZnZ0toMQkfxmZkyYMIENGzYwadIkGhsbmTlzJkceeSS/+MUvcmKpgpJbEZHcsApYamY3mNlFZja26ZHtwEQk//Tr14/777+f3/3ud4wYMYItW7bw2c9+lrFnn0392WdDjNfOa82tiEhuuDrxPLlFuQMLI46lS7RbQjxot4TMt4vzbgmpOu644/jtb39LdXU13/jGN+ixahVlwEt3381Hb701sji0W4KISJ5x98FtPIZkO7ZUabcEkdxUXFzMpEmTeP7555mS2CLs73feydy5c7McWes0cysikkPM7EBgiLuvz3YsYWm3hEDczl27JWS+Xdx3S2jXqFHw618Dwd1eBnXvDsCp7lRedRVcdVVQ78wz4fHHo4urHZq5FRHJAWY2wMxWAn8HnkyUjTWz/8xuZCKS16qqoKxsz0vbuROA0qQq24uKqLvuuogDa5uSWxGR3PB9gjuIfRjYmShbDZyTrYBEpABUVsKyZc0S3GT1wDm7d3PS177Gxo0bo42tDUpuRURywxnANe7+FsFFZLj728ABWY1KRPJfZSUsWAA9ejQv79GD9x56iHeOOoq//OUvjBw5kt/85jfZiTGJ1tyKiOSGBlpMSJhZP+Cf2Qmn87RbQjxot4TMt8uH3RKa2DvvUFxSAkVFUFoKO3ZASQkDSkp48sknmTBhAitWrOCss87iu9/9LldffXXHnYag3RJERPLPKmBm4la8Tb4BPJaleELTbgkiuavo4Yfh/ffxY46h8dFH8WOOgfffp2juXPr06cOjjz7KV7/6VRobG7n22mv58pe/nLVfYjVzKyKSG74GLAX+AZSZ2bvAS8B5WY0qBO2WEIjbuWu3hMy3y+ndEpr06wf33INddx1FRUVw1llw773Y009T1K0b3bp1Y+bMmRx77LFceeWVPPjgg/z1r3/lpz/9Kf379480VM3cSv6rrYULLoj13VREOuLu77r7qcBo4BLgXOA0d8+5ZQkikoMWLYLrrw+WJQAUF8OUKUF5kksvvZQnn3ySgw46iCeeeIKTTz6Zl156KdJQldxK/luyJPjwLV2a7UhEQjOz75vZUU2v3f137j7f3de6++5sxiYi0ppTTjmFdevWcfzxx/Pyyy9zyimn8Nhj0a2gUnIr+W/27ObPIrnlGOB/zey3ZnapmfXosIWISJYNHjyYp59+mn//93/nvffe49xzz+U73/kO7p7xsZXcSv4ZNQrM6N2nD7379IE1a4LyZ54Bs72PUaOyG6dICty9EvgYsBa4B3jdzO41s49lNzIRkfbtt99+LFiwgFtvvZXdu3czZcoULr/8cnbs2JHRcZXcSv5pcTcVEndT2fMMwfGpU6ONS6ST3P1v7v41YBDwJeBo4AUze9rMLsludCIibSsqKuKWW25h4cKF9OzZkzlz5nDGGWfw9ttvZ2xM7ZYg+SdxNxX/9KexDz7Y93hZGSxfDhUVkYcm0hXu3gAsABaY2XHAo8Bc4CdZDSwk7XMbD9rnNvPt8mmf2646//zzWb16NZ/97GdZs2YNJ510Eo8++ijDhw9Pqb32uRWprOSDOXPwVu6mwoIFSmwlZ5nZSWb2EPAUwc/waVkOKWXa51aksB1//PGsWbOGkSNHsnnzZioqKljUYreFdNDMreQtq60NtippcTcVtm7NdmgioZhZH+BS4CrgSGA5MBZY6VFcnZEm2uc2ELdz1z63mW+XF/vcpsmQIUN48sknueqqq/jJT37C2LFjuf3227n55psxs7SMoZlbyVvdfvxjqK+H4cNh8eLgub5euyZITjGzucAbBDdx+Bkw1N0vcPdf5lJiKyLSpEePHvzoRz/irrvuwsyYOnUqF198MR+0tpSwE5TcSt7y3r3ZMWMGrF8Po0fDunVw993Qp0+2QxMJoy9wIXCIu9/u7m9mOyARka4yM2644QYWL15Mr169mD9/Pqeddhqvv/56l/tWcit5a/v8+TRMntzh3VRE4szdz3f3xzRLKyL56LzzzmPt2rWUl5ezfv16TjrpJNatW9elPpXcioiIiEjWHH300fz+97/ntNNO48033+S0007jkUceaV6ptpZh8C+p9KfkVkRERESy6oADDmDVqlVceeWVbN++nYsvvphp06axe3dwl3FbupQ+sH8qfWm3BBERiZT2uY0H7XOb+Xba5zYcM+O+++7jyCOP5Ktf/Sq33347zz//PA8//DB9Hn445X40cysiIpHQPrci0hEzY/Lkybx17LE4sGjxYvp96EPYs8+m3IdmbkVEJBLa5zYQt3PXPreZb6d9bsMb8J3vsPtTn6Jo+3YAikL8tUcztyIiIiISL5WVFK1YgffsGbqpklsRERERiZ/KSmzhQrxHj1DNlNyKiIiISDxt3YqVlOBFRTjsTqWJklsRERERiadZs6C+Hj/mGF6Bl1NpouRWREREROKpb1+YOZPG3/2OrfBeKk20W4KIiIiIxNOiRcGzdksQERERkUKk5FZERERE8kbWlyWY2XRgDMEVcG8Dl7n7G2ZmwPeATwH1ifLnEm0+D0xNdHG7u8+NPnIREekM3X43HnT73cy30+130yfMucdh5namux/r7scBy4BvJMo/CQxLPCYCDwCY2YeAW4CTgZHALWbWL/KoRUQkFN1+V0SikPWZW3dPvvJtP8ATX48BfuTuDjxrZvub2YeBCmCVu/8DwMxWAecAj0QXtYiIhNV0+10zG9e9e/cN2Y4nS/oCccruo4on3eOko7/O9hG2XZj6qdQdALwTYvx8MiyVSllPbgHMbAYwgeAftDJRPBDYklStJlHWVnm7Ghsb2bZtW6dj7GzbsO3C1E+lbvC7AdTV1VFaWhoqllzXlX/vTIgqnnSPk47+cvXzE7f3UB7Z4O4nZjuIbDCzanefmO04mkQVT7rHSUd/ne0jbLsw9VOpa2brC/nzk0q9SJYlmNnjZvZCK48xAO5e5e6DgXnA5KZmrXTl7ZS3Nu5EM1tvZuvffffddJyKiIhIVyzNdgAtRBVPusdJR3+d7SNsuzD14/b+iJuUvj+RzNy6+6gUq/43sJxgTW0NMDjp2CDgjUR5RYvy1W2MWw1UA5xwwgneu3fvMGG3qrN9hG0Xpn57dYPr8qBXr16djj3Xxe28o4on3eMU4udHJN0SSzNiI6p40j1OOvrrbB9h24WpH7f3R9yk+v3J+gVlZpa8fuIzwF8SXy8BJljgFKDW3d8EVgJnmVm/xIVkZyXKREQkN6T0p0URaZU+Px3ocObWzCqA84ETgA8B/wD+B1jk7k+kIYa7zOxwgq3ANgFXJ8pXEGwDtpFgK7AvALj7PxLbh61L1Lut6eIyERGJv8Rf1USkE/T56Vibya2ZVQL3Av2AXwOLCO7p2wc4GphjZluB67qS5Lr7v7VR7sA1bRybDczu7JgiIiIikp/am7mdAXyNYNutti7YOguYDnw8A7GJiIiIiITSZnLr7qd21NjdfwX8Kq0RiYiIiIh0UtYvKBMREUlmZvuZ2R/M7NxsxyKSS8zsY2b2X2b2MzOblO14siWl5Dbxg+ZGM3vUzH6V/Mh0gCIikhvMbLaZvW1mL7QoP8fMNpjZRjO7MYWubgAWZiZKkXhKx+fH3V9y96uBsUBB3ugBUt/n9kfAEcAygp0LREREWpoD/IDg/wwAzKwYuA8YTbBP+TozWwIUA3e2aH85cCzwZ6BHBPGKxMkcuvj5cfe3zewzwI2JvgpSqsntmUC5u2/NZDBS4Gpr4bLLYM4c6Ns329GISEju/pSZlbcoHglsdPdXAMxsPjDG3e8E9ll2kNipZz/gSOADM1vh7rszGrhIDKTj85PoZwmwxMyWE9wcq+CkmtxuoY1b3IqkzZIlsGgRLF0Kl1yS7WhEJD0GEvwf0qQGOLmtyu5eBWBmlwHvKLGVAhfq85O4N8FngVKC+wUUpFST2+uAB83sbuDvyQfc/Y20RyWFafbsvc9KbkXyhbVS1uFkibvPSX8oIjkn1OfH3VcDqzMVTK5IdbcEBz5BcFewLYlHDc1/mxAJZ9QoMNv7WLMmKH/mmeblo0ZlN04R6YoaYHDS60GAJkVEUqPPTyekmtw+SLDQ+Wjg0MTjkMSzSOdUVUFZ2d7XO3c2f4bg+NSp0cYlIum0DhhmZoeYWXfgImBJlmMSyRX6/HRCqsntQcDUxBYTm5IfmQxO8lxlJSxb1jzBTVZWBsuXQ0VFpGGJSOeY2SPAWuBwM6sxsyvcfRcwGVgJvAQsdPcXsxmnSBzp85M+qa65fZxgv7R1GYxFClFlJSxYABdeCNu37y3v0SMoV2IrkjPcfVwb5Sso4ItbRFKhz0/6pJrcvgosM7OFwJvJB9z9jrRHJYVl61YoKYGiIigthR07gtdbtfOciIiIhJPqsoQTCDbVPppgI+Gmh670ka6bNQvq62H4cFi8OHiur9+7e4KIiIhIilKauXX3ykwHIgWsb1+YOROuuy6YvT3jDLj3Xnj66WxHJiIiIjkm1WUJIpmzaFHz18XFMGVK8BAREREJoc1lCWb2pJmd3l5jMzvdzFanPSoRERERkU5ob+b2DuB+M+tGsFvCn4H3gD4E9/w+E9gFXJ/pIEVEREREUtHmzK27r3T3o4AvJ4ouAaoSzwZc6+5HufvKzIcpIiIicWVmj5nZ17Mw7l1mNj2FegeY2SYzGxBFXJJdHa65TSSvSmBFREQKkJnVJb0sTTzvaCpw917u/sloowIzGwJcSQp3S3X3/zOz/wZuYe+kneSpVLcCExERkQKUSF57uXsvYC4wr0VZtkwCFrv7eynWnw18wcz6ZDAmiQEltyIiItIlZrbazKYmvi43Mzezz5vZn83sfTNbYWb9EssI3jazv5vZNS36+ISZ/dbM/mFmL5vZFDOzdoY9H1iV1N7MbIaZvWFm28zsNTPbM0vr7n8D3kF79Oc9JbciIiKSCf8GfBwYApQDvwNeBj4CfAG4N7G0ADM7iuAWszOBA4BPA5OBS1vr2Mx6AkcQXOzeZDTweeBkd+8NnAw806Lp8wQ3ppI8puRWREREMmG6u//D3d8FlgEN7v6Qu+9y98eAfwLHJ+pOAn7q7ovdvdHd/wL8AJjQRt/9Es/JSxJ2Aj2Ao8ysh7u/5e7PtWj3HvChNJybxJiSWxEREcmEN5O+rm/xuqmsd+LrQ4BxZra16UFw8deH2+j7n4nnPetn3X01cDMwFXjbzFaa2Ykt2vUB/hH2RCS3pHSHMjPbD7gWOJG9b0QA3P2sDMQlIiIihWMTMNvdr+mwJuDuH5jZBoJ99/+UVF4NVJtZGXAr8HOCZRFNjgbmpClmialUb787m+BPB4uA9zMXjoiIiBSg+4EnzeyXwC8BBw4DDnD3J9tos4jg4rD5AGZ2EsFWZesItirbRnCzKRLHP0qwnvfxDJ2DxESqye1ZwGHu/n+ZDEZEREQKj7u/YGbnArcDDxMsm9wI3N1OsweAP5jZ9YntwHoD9wDDgEaCi8cuSqp/OTDH3WszcAoSI6kmt+8CdR3WEhERkbzl7le2UV6R9PVrBHcyTT5+ayttylu8XgucGSKWTWb2Q+BrwDR3/w1t7ISQuDPZeILllZLnUk1ubwa+b2Y3uLsWYouIiEjWufuNKdZ7Bxia4XAkJtpMbs2sgWDNS3Ldy82sMbmeu3fPUGwiIiIiIqG0N3OrO3iIiIiISE5pM7lNvjrRzE5z96da1jGzT2QqMBERERGRsFK9icOyNsoXpysQEREREZGuSjW5tX0KzHoDu9MbjoiIiIhI57W7W4KZ/Y3gorKeZvbXFocPBFZlKjARERERkbA62grsdoJZ2weAGUnlu4G/A7/JUFwiIiIiIqG1m9y6+1wAM/uLuz8bTUgiIiIiIp2T0k0c3P1ZMxsCXAwMAmqA+Ym7kIiIiIiIxEJKF5SZ2TnABuDTQN/E80uJchERERGRWEj19rszgSvc/b+bCsxsHPBt4JeZCExEREREJKxUtwIrB+a3KFsADElrNCIiIiIiXZBqcrsaqGhRdjrw5D41RURERESyJNVlCRuBX5jZIuA1gpnc84FZZnZzUyV3vyPdAYqIiIiIpCrV5PY44DmCZQhNSxGeA45PquOAklsRERERyZpUtwKrzHQgIiIiIiJdleqaW8ys2MxONbPPJV6XmVnPzIUmIiIiIhJOqvvc/gvwArACmJUoPgt4KENxiYiIiIiElurM7X8SbAX2IaAhUbYa+EQGYhIRERER6ZRULygbCXzG3XebmQO4+1Yz2z9zoYmIiIiIhJPqzO17QLNE1sw+AryV9ohERESkXWZ2tpk9nWLdSWb24xTqvWZmFe0c/3jTBFeK4642s6mp1k83MxtkZm5m5dkc18zKw3zfpOtSTW5/Dsw2s0EAZtYfuJd971omIiIiGWRmBnwXuCXFJg8Bp5vZiZmLSiQ+Uk1upwF1wGaCGdy3gR1oX1sREZGonQV0B55IpbK77wJ+DFybyaDyjZl1y3YM0jkpJbfu/oG7XwwcQLD+9mB3v9Tdt2c0OhEREWnpfOBxd9/zp24z62ZmN5vZBjPbZmYvm9m/JbVZBZxnZmG2AB2WWFqwzcz+FzixxfEyM/uemW0xs3fMbJGZDWnRzQAzW2ZmdWb2opl9Mqn98Wb2WzOrNbN/mNkaM+sX6jvRPJ6DzWxJor+/Aue0OD7HzOaZ2Y/N7L3E9+iypOOXmdlGM/uamdUAf0yU9zezWYnz/D8zW2hmB6U6rkQvzJv8MOAw4D13/7/MhSQiIiLtOAH4c4uy24FLgAuBPsDpwN+Sjj9P8JfXQ1MZwMxKgKXAi8CBwL8DV7eo9l3glMRjKPAOsNTMipPqXAF8LzH2HcAvktbA3gf8imAnpoOA64GdqcTXhnlAI8GdVE8DLmulzlhgZWLMq4EHzOzUpOPlwEeAYcBJiSUgiwjuwnp04jy36Fi8ggAAIABJREFUAf8dclyJUIfJrZmNMbPNwEvAM8BLZrbZzM7PeHSSW2pr4YILgmcREcmUfgQXegN71uBeA3zN3f/kgRp3/1NSm6b6H0pxjJOBQxJ9fuDufwO+nTRmETABmOrur7v7+8B1wMcI/sLbZJG7r3L3Xe4+D1gPXJw4tpMgIRzs7g3u/myin9DMbCBwBvBVd691978D32yl6rPu/pNEPKuAR2mejDYANybOuR4YkXhck+i3Hvg6cEbiwrFUx5UItZvcmtkJwE8JfrM6Gzgy8fwrYKGZjch4hJI7liyBRYtg6dJsRyIiks/+STA72+QAYD/gr+20aar/jxTHGAS8nUjmmrzaYswewCtNBe5eR3BNzuCkeq+16Pe1RN8AXyDIQ35rZq+a2fTEjHEzZjYksayh6dHaHvtNfW5qI95U4gF40913JL0+BCgF3jKzrWa2FXgZ2E6QmKc6rkSoo31urwO+5e7Tksr+AjxuZq8njl+aqeAkx8yevff5kkuyG4uISP76H4LJpib/B7xP8Kf0v7XaIviTei2pJ16vAweaWVlSgntIizF3JMpeBjCzXgRLGLYk1Stv0W85wd1OcfdXgcsTbY8hmDh7FZid3MDdNwO9UogXgmUDL7cSb3vx1CS93t3i+CaC7+2H3L3lMZp2kUphXIlQR8sSTgWq2zj2UOK4FKpRo8Bs72PNmqD8mWeal48ald04RUTyyyLgzKYXiQvLHgDuNrOjLTAwkTA2GQ0sdffGFMd4liCxu8vMeprZvwBfSRpzN/AjYLqZfcTMygiWLfwF+H1SP+eb2ZlmVmxm44CTSGwjamafT+yZD7AV2JV4hObuNQR3Tr3bzPokLvia1krVU8xsXCKeM4B/S5xHW9YTXFj2vcQ2qJjZAWZ2UchxJUIdJbf93X1LawcS/6D90x+S5IyqKigr2/t6587mzxAcn5q1PbxFRPLRSmCXNb/hQhWwkCDx3QY8STCT23Rx2KXA91MdILF92GeA4QRLDX7OvpNdXyFI/tYRbBX6YYK7mSYn0LMILhSrBb4BfNbdm5YynAH8wczqgLUEF2nNSzXGVlxMsIRgC/A0rSetC4FPESztmEWwlva3bXWYSOLPJ8iX/mBm24DfARUhx5UIdbQsoaPk19IViJl9FZgJHODu7yQWyH+P4E1YD1zm7s8l6n4eaMqYbnf3uemKQ0KorIRly+Dcc6G+ft/jZWWwfDlUVEQemohIvnJ3N7OvALcRXJ2Pu+9MvL6tlSZXAE+7+7qQ4/yFYNeFZPcmHX8f+HLi0Vr7ig76/3yYeDri7m8C57Yo/mGL1x+4+5VttJ8DzGml/B8EF+xd04VxJUIdJbelZnZzO8e7pyMIMxtM8CeTzUnFnyT4rXMYwVWbDwAnm9mHCO7KciLB1hx/MLMl7v7PdMQiIVVWwoIFcOGFsD1p2+MePYJyJbYiImnn7r8Efpli3QeBBzMbkUh8dJTcPkuQdLZ3PB2+S7C1xuKksjHAjxJriZ41s/3N7MMEfwpYlfhNCjNbRbBh8iNpikXC2roVSkqgqAhKS2HHjuD11q3ZjkxERFJ3L/vuJiBdtxVtDxapdpPbjv6kkA5m9hngdXf/32Alwh4DaX7FZU2irK3ydjU2NrJt27ZOx9nZtmHbhamfSt2mG9jU1dVRWloaKpZU9ayupvj999l9zDHsuO02Sr/xDYqef57Ghx7igzFjMjJmKrry750JUcWT7nHS0V+ufn7i9h4SySR3v7fjWrnL3S/L0rhbgVuzMXah6mjmNi3M7HHg4FYOVQE3E9wne59mrZR5O+WtjTsRmAgwcGCH+a90kvfuzY4ZM2j40pegqIj600+n2/33U9y0e4KIiIhIRCJJbt291b2gEtuUHAI0zdoOAp4zs5EEM7LJG0EPAt5IlFe0KF/dxrjVJK7uPOGEE7x3795dOQ0AOttH2HZh6rdXt2k2vFevXp2OvUPLl9ONYDfvPW4Olmp3y8yIoWTsvDspqnjSPU4hfn5ERCT3RJLctsXdnyfY8BkAM3sNODGxW8ISYLKZzSe4oKzW3d80s5XAHWbWL9HsLOCmiEMXEZFOGjBggJeXl2c7jMg1LRNrsQQva6KKJ93jpKO/zvYRtl2Y+qnUjdt7KEruznPPPfeOux/QUd2sJrcdWEGwDdhGgq3AvgDBlhxmNp1gXz2A25ouLhMRkfgbOnQo69evz3YYkWtoaACgW7c4/E0runjSPU46+utsH2HbhamfSt24vYei1NDQQPfu3Td1XDNmya27lyd97bS9p9xsWtyeT0RERESko5s0iIiIiIjkDCW3IiIiIpI3lNyKiIiISN5QcisiIikzs0PNbJaZ/SyprMLMnjaz/zKziiyG16p586C8PLiJYnl58FpE8peSWxGRAmFms83sbTN7oUX5OWa2wcw2mtmN7fXh7q+4+xUti4E6gu2ua9IbddfMmwcTJ8KmTeAePE+cqARXJJ8puRURKRxzgHOSC8ysGLgP+CRwJDDOzI40s2PMbFmLx4H7dgnA0+7+SeAG4JsZjD+0qiqor29eVl8flItIforVVmAiIpI57v6UmZW3KB4JbHT3VwASN84Z4+53Auem2O/uxJf/BEpTqL9nv85M27y5hNbu2r55s9PQsCuSGJrs2hXteB2JKp50j5OO/jrbR9h2YeqnUjdu76EohTl3zdyKiBS2gcCWpNc1ibJWmVl/M/sv4HgzuylR9lkzexD4MfCDNtpNNLP1Zrb+nXfeSV/0HRg8OFy5iOQ+zdyKiBS21u7j6W1Vdvd3gatblP0c+Hl7g7h7NVANMGLECI/qDkt33BGssU1emlBWBnfcYVm7y1Pc7i4VVTzpHicd/XW2j7DtwtRPpW7c3kNxo5lbEZHCVgMkz2MOAt7IUixpN348VFfD0KFgFjxXVwflIpKfNHMrIlLY1gHDzOwQ4HXgIuDiTA4Y5ZpbgLFjg0eyCIffI27rJbXmNvPttOY2fbTmVkRE9mFmjwBrgcPNrMbMrnD3XcBkYCXwErDQ3V/M0PjnmVl1bW1tJroXEQE0cysiUjDcfVwb5SuAFRGMvxRYOmLEiKsKec1g3M5da24z305rbqOlmVsRERERyRuauRURkUhFveY2LuK2XlJrbjPfTmtu00drbkVEJHa05lZEoqCZWxERiYTW3Abidu5ac5v5dlpzGy3N3IqIiIhI3tDMrYiIREprbuNBa24z305rbtNHa25FRCR2tOZWRKKgmVsREYmE1twG4nbuWnOb+XZacxstzdyKiIiISN5QcisiIiIieUPJrYiIiIjkDa25FRGRSGm3hHjQbgmZb6fdEtJHuyWIiEjsaLcEEYmCZm5FRCQS2i0hELdz124JmW+n3RKipZlbEREREckbSm5FREREJG8ouRURERGRvKE1tyIiEintlhAP2i0h8+20W0L6aLcEERGJHe2WICJR0MytiIhEQrslBOJ27totIfPttFtCtDRzKyIiIiJ5Q8mtiIiIiOQNJbciIiIikjeU3IqIiIhI3lByKyIiIiJ5Q7sliIhIpLTPbTxon9vMt9M+t+mjfW5FRCR2tM+tiERBM7ciIhIJ7XMbiNu5a5/bzLfTPrfR0sytiIiIiOQNJbciIiIikjeU3IqIiIhI3lByKyIiIiJ5Q8mtiIiIiOQNJbciIiIikjeU3IqIiIhI3tA+tyIiEindoSwedIeyzLfTHcrSR3coExGR2NEdykQkCpq5FRGRSOgOZYG4nbvuUJb5drpDWbQ0cyvtq62FCy4InkVERERiTsmttG/JEli0CJYuzXYkIiIiIh1Scivtmz27+bOIiIhIjCm5leZGjQKzvY81a4LyZ55pXj5qVHbjFBEREWmFkltprqoKysr2vt65s/kzBMenTo02LhEREZEUKLmV5iorYdmy5glusrIyWL4cKioiDUtEREQkFUpuZV+VlbBgAfTo0by8R4+gXImtiIiIxJSSW2nd1q1QUgJFRdCzZ/BcUhKUi4h0wXPPGeXlMG9etiMRkXyk5FZaN2sW1NfD8OGweHHwXF+vXRNEJC02bYKJE5Xg5pN586C8PJgL0S8vkk1ZT27N7FYze93M/ph4fCrp2E1mttHMNpjZ2Unl5yTKNprZjdmJPM/17QszZ8L69TB6NKxbB3ffDX36ZDsyEckT9fXBNayS++bNC35Z2bQJ3PXLi2RXXG6/+113vye5wMyOBC4CjgI+AjxuZoclDt8HjAZqgHVmtsTd/xxlwHlv0aLmr4uLYcqU4CEikiabN2c7AkmHqqrgl5VkTb+8jB2bnZikcMUluW3NGGC+u+8AXjWzjcDIxLGN7v4KgJnNT9RVcisikmMGD3YaGnZlO4xI7NoVr/NMZzybN5cA1kq5p/2809FfZ/sI2y5M/VTqxu09FKUw5x6X5HaymU0A1gNT3P2fwEDg2aQ6NYkygC0tyk/uaIDGxka2bdvW6QA72zZsuzD1U6nr7gDU1dVRWloaKpZc15V/70yIKp50j5OO/nL18xO391CuM7PzgPNgBABlZc706Y3ZDUrSYvDg1mfhBw+OPhaRSJJbM3scOLiVQ1XAA8B0wBPP3wYup7VfAYM6ra0T9jbGnQhMBBg4cGBrVUREJCLuvhRYanbiVUOHwowZxvjxcZljiU63bt2yHUIz6YjnjjuCNbbJSxPKyuCOO4ySkpK0jZMsHf11to+w7cLUT6Vu3N5DcRPJTxV3T+lerWb2ELAs8bIGSP6dbxDwRuLrtspbjlsNVAOccMIJ3rt37xBRt66zfYRtF6Z+e3XNgt8RevXq1enYc13czjuqeNI9TiF+fiQzTjjB+cMfsh2FpNP48cFzVVUwgztkCMyYEZQ3NGQ3Nik8Wf+V2cw+7O5vJl5eALyQ+HoJ8N9m9h2CC8qGAb8nmNEdZmaHAK8TXHR2cbRRi4iISLLx4/cmuSLZlPXkFrjbzI4jWFrwGvBFAHd/0cwWElwotgu4xt0bAcxsMrASKAZmu/uL2QhcREREROIl68mtu1/azrEZwIxWylcAKzIZl4iIiIjknqzfxEFEREREJF2U3IqIiIhI3lByKyIiIiJ5Q8mtiIiIiOQNJbciIiIikjeU3IqIiIhI3lByKyIiIiJ5Q8mtiIiIiOQNJbciIiIikjeU3IqIiIhI3lByKyIiIiJ5Q8mtiIiIiOQNJbciIiJ5ZN48KC+HoqLged68bEckEq2SbAcgIiIi6TFvHkycCPX1wetNm4LXAOPHZy8ukShp5lZERCRPVFXtTWyb1NcH5SKFQsmtiIikzMwONbNZZvazpLIiM5thZv9pZp/PZnyFbvPmcOX5QMswpCUltyIiBcLMZpvZ22b2Qovyc8xsg5ltNLMb2+vD3V9x9ytaFI8BBgINQE16o5YwhgwJV57rmpZhbNoE7nuXYSjBLWxacysiUjjmAD8AftRUYGbFwH3AaILEdJ2ZLQGKgTtbtL/c3d9upd/DgbXu/mBiRvfX7QXh7jQ0NHT6JHLVrl27Mj7GbbcZkyYVU19ve8rKypzbbmukocEjjycT4yT3d/PNJc3OFYJlGDff7Iwd2/a4nY0pbLsw9VOpG9W/WRyFOXcltyIiBcLdnzKz8hbFI4GN7v4KgJnNB8a4+53AuSl2XQPsTHzd2FoFM5sITAQYPHhwuMAlZePGOdDItGnFbNkCgwfD9OmNifL8s2VLuHIpDEpuRUQK20AgORWoAU5uq7KZ9QdmAMeb2U2JJPjnwH+a2SeAp1pr5+7VQDXAiBEjvFu3bmkKP/dk+twnTAgee7X/X31U/xbpHqdbt24MGRIsRWhpyBBLabzOxhS2XZj6mYy7UCi5FREpbNZKWZvTfO7+LnB1i7J6oOU6XJGMmzGj+dZnAGVlQbkULiW3IiKFrQZIXicwCHgjkwNqzW085MOa27FjobHR9lmGMXas095bTGtuc4/W3IqISKrWAcPM7BDgdeAi4OJMDGRm5wHnHXrooZnoXgrUuHHOuHGFm/TJvpTciogUCDN7BKgABphZDXCLu88ys8nASoIdEma7+4uZGN/dlwJLR4wYcVUhrxmM27nn8prbbPWhNbfxpuRWRKRAuPu4NspXACsiDkdEJCOU3IqISKS05jYe8mHNbdR9aM1t9oQ5d92hTEREImFm55lZdW1tbbZDEZE8pplbERGJhNbcBuJ27lpzm/l2WnMbLc3cioiIiEjeUHIrIiIiInlDyxJERCRSuqAsHnRBWebb6YKy9NEFZSIiEjuFdkHZI48YH/1oCaWlJXz0oyXMn1+c7ZBEuqTle/qRR1q7e3f2aeZWREQiUUgXlM2bB5MmQX198HrzZpg8uRvFxUVMmBCv/3p1QVnm2+XDBWWtvacnTSqhuBjGj488nHZp5lZERCTNqqr2JgFN6uuNadM0eyu5qfX3dFAeN0puRURE0mzz5tbLt2yJNg6RdGnrPd1WeTbF628jIiKS9wrhgrLBg0vYvHnf9YiDBjkNDfG4KEgXlGW+XT5dUNbWe3rw4Gje07qgTEREYqeQLiibPr2RsjJvVtazp/PNb8YjsZWuy5WLq9Kltfd0WZkzfXpjliJqm2ZuRUQkEoV0QdmECVBcHKxH3LwZhgyB225rZNw4i93dpXRBWXgLF3YLdXFVPlxQ1tp7esYMY/z4+KWS8YtIREQkD4wf3zzRaWjwtitLTmnv4qq47RyQTi3f03GlZQkiIiIiIeTSxVWFSMmtiIiISAhDhoQrl2gpuRUREREJYcYMKCtrXlZWFpRL9mnNrYiIRKoQtgJrTTa3cWqNtgLrfLuxY6GxMbgpx5YtMHhwsJvA2LFO8ls7n7YCy7Yw567kVkREImFm5wHnHXroodkORaTLxo1zxo0r3GQzzpTciohIJAppK7D2xO3ctRVY5tvlw1ZguURrbkVEREQkbyi5FREREZG8oeRWRERERPKGklsRERERyRu6oExERCKlrcDiQVuBZb6dtgJLnzDnrplbERGJhJmdZ2bVtbW12Q5FRPKYZm5FRCQS2gosELdz11ZgmW+nrcCipZlbEREREckbSm5FREREJG8ouRURERGRvKHkVkRERETyhpJbEREREckbSm5FREREJG8ouRURERGRvKHkVkRERETyRixu4mBmXwYmA7uA5e7+9UT5TcAVQCNwrbuvTJSfA3wPKAZ+6O53ZSVwEREJTbffjQfdfjfz7XT73fQJc+5ZT27NrBIYAxzr7jvM7MBE+ZHARcBRwEeAx83ssESz+4DRQA2wzsyWuPufo49eRERSZWbnAecdeuih2Q5FRPJY1pNbYBJwl7vvAHD3txPlY4D5ifJXzWwjMDJxbKO7vwJgZvMTdZXciojEmG6/G4jbuev2u5lvp9vvRisOye1hwCfMbAawHfiqu68DBgLPJtWrSZQBbGlRfnJHgzQ2NrJt27ZOB9nZtmHbhamfSl13B6Curo7S0tJQseS6rvx7Z0JU8aR7nHT0l6ufn7i9h0REpGORJLdm9jhwcCuHqhIx9ANOAU4CFprZoYC1Ut9p/SI4b2PcicBEgIEDB7ZWRURERETySCTJrbuPauuYmU0Cfu7BFOPvzWw3MIBgRnZwUtVBwBuJr9sqbzluNVANcMIJJ3jv3r07fQ5NOttH2HZh6rdX1yz4HaFXr16djj3Xxe28o4on3eMU4udHRERyTxy2AlsEnAGQuGCsO/AOsAS4yMxKzewQYBjwe2AdMMzMDjGz7gQXnS3JSuT5oLYWLrggeBYRERHJcXFIbmcDh5rZC8B84PMeeBFYSHCh2C+Ba9y90d13EWwbthJ4CViYqCudsWQJLFoES5dmOxIRERGRLsv6BWXuvhO4pI1jM4AZrZSvAFZkOLTCMHv23udLWv1nEBEREckZcZi5lSiNGgVmex9r1gTlzzzTvHxUm8ukRURERGJLyW2hqaqCsrK9r3fubP4MwfGpU6ONS0RERCQNlNwWmspKWLaseYKbrKwMli+HiopIwxIRERFJByW3haiyEhYsgB49mpf36BGUK7EVERGRHKXktlBt3QolJVBUBD17Bs8lJUG5iIiISI7K+m4JkiWzZkF9PQwfDt/6FtxwA/zv/2rXBBHJOHenoaEh22FEbteuXdkOoZmo4kn3OOnor7N9hG0Xpn4qdeP2HopSmHPXzG2h6tsXZs6E9eth9GhYtw7uvhv69Ml2ZCKSp8zsPDOrrtVNY0QkgzRzW6gWLWr+urgYpkwJHiIiGeDuS4GlI0aMuKpbt27ZDidr4nbuUcWT7nHS0V9n+wjbLkz9VOrG7T0UN5q5FREREZG8oeRWRERERPKGklsRERERyRtKbkVEREQkbyi5FREREZG8oeRWRERERPKGklsRERERyRva51YkRzU0NFBTU8P27duble/evRuAoqL0/O6ajv4620fYdmHqN9Xt1q0b+++/PwMGDEjb90xERLJHya1IjqqpqaF3796Ul5djZnvKGxsbASguLk7LOOnor7N9hG0Xpn5jYyPuzu7du3nrrbeoqalhyJAhoeITEZH40TSFSI7avn07/fv3b5bYSjhmRvfu3Rk4cCDvv/9+tsMREZE0UHIrksOU2KaHliOIiOQP/UQXERERkbyh5FZERERE8oYuKBORjOnVq9eer3fs2AFAaWnpnrK6urrIYxIRkfym5FZEMqYpeW1sbGTixIk0NjYyZ86cdts0NDTQrVu3CKITEZF8pGUJIpJVgwYN4vbbb+f0009nv/32Y/HixUydOpVzzjmnWb2Pf/zj3HXXXXte/+lPf2L06NEMGDCAIUOGUFVVRUNDQ9Thi4hIzCi5FckTZoaZUVJSQklJyZ7XXX201V86PfTQQ3z/+9+nrq6OT3/60x3W//vf/87pp5/O5z73Od58803WrFnDY489xj333JPWuEREJPcouRWRrPviF7/I8OHDMTN69uzZYf05c+Zw0kknceWVV9KtWzcGDRrEDTfcwE9+8pMIohURkTjTmluRPOHuQDzvUNaR8vLyUPVfffVVnnzySfbff/89Zbt379Z+tREws0OBKqCvu/97ouwTwHiC/1OOdPdTsxiiiBQ4/U8gIlnXMint1avXPncMe/PNN/d8PXToUM455xy2bt265/Hee+/x7rvvRhJvrjKz2Wb2tpm90KL8HDPbYGYbzezG9vpw91fc/YoWZU+7+9XAMmBu+iMXEUmdklsRiZ0TTzyRdevW8cc//pGGhgbuvfdeNm/evOf4ZZddxtq1a5k7dy7bt29n9+7dvPzyy6xcuTKLUeeEOUCzK/XMrBi4D/gkcCQwzsyONLNjzGxZi8eBHfR/MfBIJgIXEUmVliWISOyMGjWKa6+9lrPPPpuioiImTZrEySefvOf4Rz7yEX7zm99w00038fWvf50dO3ZQXl7OF7/4xSxGHX/u/pSZlbcoHglsdPdXAMxsPjDG3e8Ezk21bzMbAtS6+3spxFGQO1vs2rUr2yE0E1U86R4nHf11to+w7cLUT6Vu3N5DUQpz7kpuRSQS1dXVra7brampabX+3XffzZ133gkE632/+c1vNjt+9NFHs3Tp0mZlTeuDJZSBwJak1zXAyW3Uxcz6AzOA483spkQSDHAF8HA77SYCEwEGDx7c1ZhFRNqk5FZEpLC1tq+bt1XZ3d8Frm6l/Jb2BnH3aqAaYMSIEV7IN+qI27lHFU+6x0lHf53tI2y7MPVTqRu391DcaM2tiEhhqwGSp1IHAW9kKRYRScG8eVBeDkVFwfO8edmOKF40cysiUtjWAcPM7BDgdeAiggvDMkZrbuNBa24z3y4Ta27nzy9m8mSnvj74o8umTTBxotPY2Mi4cW3+0SXnhfleauZWRKRAmNkjwFrgcDOrMbMr3H0XMBlYCbwELHT3FzM0/nlmVl1bW5uJ7kUKwi23lOxJbJvU1xvTpmVuL/Jco5lbEZEC4e7j2ihfAayIYPylwNIRI0ZcVchrBuN27lpzm/l26VxzW1PT+u3Pt2yx2L23skUztyIiIiI5oq3NRoYMiTaOONPMrYiIREprbuNBa24z3y4Ta25vvdWZPLlbs6UJZWXObbc10tCgNbegmVsREYmI1tyKdN1FFzXywAONDBnimDlDhjgPPJDfF5OFpZlbERGJhNbcBuJ27lpzm/l26d7ndsKEEiZMSC5ROpdMM7cihay2Fi64IHgWERHJA0r1RQrZkiWwaBEsXQqXXJL27nv16rXn6x07dgBQWlq6p6yuro6KigrWrl1L9+7dKSoqon///vy///f/uO666zjuuOP21G2q13JWY+3atRxzzDFcdtllzJ07lzvvvJMbb7xxz/E33niDIUOG0NjYiLv+bBcHWnMbD1pzm/l2mVhzW6i05lZEUjN7dvPnNKurq6Ouro7a2lomTJjA+PHj95TV1dXtqTdt2jS2bdtGbW0tTzzxBEOHDuWUU05h0aJFzfqbNm1as/Z1dXUcc8wxe45/7GMfY3aLc5k9ezaHHXZYRs5PwtGaWxGJgmZuRQrJqFHw61/vfd29e/D8zDNgSXsnnnkmPP54tLElDB06lNtvv50333yT//iP/2DMmDEpt/3Xf/1Xnn32WVavXk1FRQXuzqxZs7j22mu5/vrrMxi1pEJrbgNxO3etuc18u3SvuY3beyhuNHMrUkiqqqCsbO/rnTubP0NwfOrUaONqxUUXXcTrr7/Ohg0bQrW7/PLLeeihhwBYtWoVffv25aSTTspEiCIiEkNKbkUKSWUlLFvWPMFNVlYGy5dDRUWkYbVm0KBBALz77rt7ymbMmMH+++/f7NHSpZdeyvLly3n33Xeprq7mqquuiixmERHJPiW3IoWmshIWLIAePZqX9+gRlMcgsQWoqakBoH///nvKqqqq2Lp1a7NHS/379+eTn/wk99xzD48//jjjx4+PLGYREck+rbkVKURbt0JJCRQVQWkp7NgRvG4lWcyWBQsWMHDgQA4//PDQbSdOnMiZZ57JhAkTWp1Us4TNAAAPSUlEQVTdlezSbgnxoN0SMt9OuyWkj3ZLEJH2zZoF9fUwfDgsXhw819dnbNeEMLZs2cItt9zCnDlz+O53v4slX+iWooqKClatWsWdd96ZgQils7RbgohEQTO3IoWob1+YOROuuy6YvT3jDLj3Xnj66ayEM336dL71rW9hZvTv359TTz2VNWvWMGLEiH3q3XXXXc3K5s+fz7nnntuszMw488wzMx63hKPdEgJxO3ftlpD5dtotIVpKbkUKUYv9YykuhilTgkeGVFdXU1xcvE/56tWr22zT2NiYUj2AOXPmNKuf7OMf/7hu4CAiUiC0LEFERERE8oaSWxERERGJrXnz4KMfLQFarFVrg5YliIhIpLRbQjxot4TMt9NuCV33yCPGpEnF1NenfnGxZm5FRCQS2i1BRMKaNi1cYguauRXJae7eqa2ypDldbBYN7ZYQiNu5a7eEzLfTbgmdt2VL+DaauRXJUcXFxQX5p91M+OCDDwrqPwsRkVwxZEj4NkpuRXLU/vvvz1tvvcXu3buzHUrOcnfq6+t5/fXXOfDAA7MdjoiItDBjBpSVhWujZQkiOWrAgAHU1NSwYcOGZuVNyW5RUXp+d01Hf53tI2y7MPWb6paWlnLQQQfRp0+fULGJiEjmjR8fPN98s7N5c2ptlNyK5KiioiKGtPL3mm3btgHQu3fvtIyTjv4620fYdmHqp/v7JCIimTF+PIwdu4vu3f/wh1TqZ31ZgpktMLM/Jh6vmdkfk47dZGYbzWyDmZ2dVH5Oomyjmd2YnchFREREJG6yPnPr7p9r+trMvg3UJr4+ErgIOAr4CPC4mR2WqHofMBqoAdaZ2RJ3/3OkgYuISKdon9t40D63mW+nfW7TJ8y5Zz25bWLBfkZjgTMSRWOA+e6+A3jVzDYCIxPHNrr7K4l28xN1ldyKiMSYmZ0HnHfooYdmOxQRyWOxSW6BTwBvufvfEq8HAs8mHa9JlAFsaVF+ckedNzY27llj1xmdbRu2XZj6qdRt2r+zrq6O0tLSULHkuq78e2dCVPGke5x09Jern5+4vYdynfa5DcTt3LXPbebbaZ/baEWS3JrZ48DBrRyqcvfFia/HAY8kN2ulvtP6OuFWd2A3s4nAxMTL7X369HkxtYjb1JfEsokMtwtTP5W6Aw455JB3QoyfTzr7b5YpUcWT7nHS0V+ufn4ODzG2pOC5556rM7MNHdfMS/qZFJ/+cvVn0gCgUP9PH5ZSLXfP+oMgyX4LGJRUdhNwU9LrlcC/Jh4r26rXzhjVaYizU32EbRemfip1gfXZ/jfO1iMd/+65GE+6x9HnJ/vvnXx6FPL3VD+T4tOffibl3iPV72XWd0tIGAX8xd1rksqWABeZWamZHUKQrf8eWAcMM7NDzKw7wUVnS1IYY2ka4uxsH2HbhamfjvPKZ3H7/kQVT7rH0edHJD3i9p7Tz6TMt9PPpPRJ6ftjiUw4q8xsDvCsu/9Xi/Iq4HJgF3Cduz+WKP8UcC9QDMx29xnRRpxbzGy9u5+Y7ThEcpE+P+mn76lI5+nz07FYJLeSWWY20d2rsx2HSC7S5yf99D0V6Tx9fjqm5FZERERE8kZc1tyKiIiIiHSZklsRERERyRtKbkVEREQkbyi5LXBmdr6ZPWRmi83srGzHI/+/vTuPsauswzj+fQABtQURShQRWlJEFo2SIhoB0Vjiwu4CRJayuKCAEhNBhYhpECIiiwRlK4UKFETDJoigVtwFVGQtFKTsDEuFFlAoffzjfS8ebmemnemdzszt80lu7txz3vOe3znJee9v3vO+58ZoImkjSedIunS4Y+kmkl4v6WZJOw53LBGjiaRNJf1I0qWSDh7ueIZLkttRTNI0ST2Sbmtb/hFJsyXNkXRkf3XYvsz2Z4EpwB5DGG7EiNKh6+c+2wcObaSjRyfOaXUEcMnQRBkxMnWoTbrT9heATwMr7OPC8rSEUUzSdsAC4HzbW9RlKwN3A5OBhyg/erEX5ZnAx7VVcYDtnrrdicAFtv+2nMKPGFYdvn4utf3J5RX7SNWJcwq8k/LzoqsDT9q+avlEHzG8OtUmSdoZOBI4zfaFyyv+kWSV4Q4gBs/2DZLGty1+DzDH9n0AkmYCu9g+DljsFp8kAccD1ySxjRVJJ66feLUOtUkfBF4PbAa8IOlq24uGNPCIEaBTbZLtK4ArJP0cSHIbXeEtwIONzw8BW/dT/lDKzx+vKWli+6/ERaxgBnT9SFobOBZ4t6Sv1y+ceLUBnVPb3wSQNIXSc5vENlZkA22Ttgd2B1YDrh7SyEawJLfdR70s63Psie1TgVOHLpyIUWWg189TwBeGLpyuMKBz+koBe3rnQ4kYdQbaJs0CZg1VMKNFJpR1n4eAtzY+rw88MkyxRIw2uX46L+c0YvBy/QxCktvucyOwsaQJklYF9gSuGOaYIkaLXD+dl3MaMXi5fgYhye0oJuki4E/AJpIeknSg7YXAIcC1wJ3AJbZvH844I0aiXD+dl3MaMXi5fjonjwKLiIiIiK6RntuIiIiI6BpJbiMiIiKiayS5jYiIiIiukeQ2IiIiIrpGktuIiIiI6BpJbiMiIiKiayS5jRFF0vWSjhnC+r8h6coBlF8g6X1DEMdMSQd2ut4l7PNgSTOW5z4jYsUg6RpJXxuG/R4vaepSlBsnaa6kdZZHXDG8ktxGnyRNknSZpCckPSvpbkknS3rzcMe2NCTNknRUc5nt79jeaWnrsD3G9p9qfdtLWtiBuN4LvAeY3lj2IUk3SJon6RlJsyUd28u2v5K0xyB3fRbwAUmTBrl9RKyA6j/5rddL9fXKMgDbH7X93eUc1wbAQcAJSypr+wngQuBbQx1XDL8kt9ErSZOB3wOzgXfZXgP4APBUfY/B+zJwru2XASRNAK6iJJ/rAmsDuwN3NTeStBawNXD1YHZaf+lmBnDYoCOPiBVO/Sd/jO0xwHnABW3LhsvBwOW2n13K8tOA/SWtMYQxxQiQ5Db6cjpwoe0jbD8MYPtR21NtzwSQ9DpJp0h6UNKTtZd3g1YFtef0REk/lTRf0r2Sdmmsl6Sv158ZfFrSSYAa6xfrKZV0jKTrG5/HSTpH0gO1d/lmSZtIOg3YFji69i7Mbt9e0iGS/t5W/wRJL0saXz9b0jaS1gOuAVZu9FjsJ+liSae01XGApDmSRBtJqwAfB65rLN4SmG97hu2XbC+0fbvt9iEEOwI32J4vaUrdx+H1/M2X9D1Ja9fz/aykuyRt01bHdcBOknLtR0THNO+USRpf2879JN0h6TlJV0taqw4j6JH0mKQvtdWxraTf1++DeyV9tbd2tGFXGm1p/U45VtIjtU28X9KhrfW27wGeBD7c2aOPkSZfcLEYSW8DJlJu4fTnJOC99bUhpdG4UtLKjTL7Ad8H1gROA86T9Lq6bm/gcGAX4E11++0GEOdKwOXAG4Ct6vv+lETxEOB3wNTau7BJL1VcAGwq6V2NZVOAWbbvbxa0/QjwUeDlRo/FecAZwN6SVmsUPwg4273/tvXGwFjgjsaym4AxkmZI2lXSW/s45N2AyxqfN6zHvBGwDXAoJQE/AVgL+Blwblsdtza2iYgYSp+gtE0bAOOBvwD3AutR2uqTWx0ikjan3JU6ARhH6QQ4BNint4olvRZ4O69uSydTvnO2tj2WcqfrD22b3krpUIguluQ2ejOuvj/cV4GaWO4LHGX7YdvPAV8BNqWMJ2252PYfbC8CzqQkuRvXdfsCZ9i+2faLwHHAYwOIcxIlqT3A9uO2F9n+Z01El8j2PEpyvH89JlEaxmkDiOE3lKEau9U6Nq1xTe+j/Fr1fX4jjrmURvi/wPeAubXXdddWGUmrU3obrmjU9QLwbdsv2r4FuAW40faf65CHHwMTJa3Z2KZ1++6NAzjGiIjBmGr7adtPUYZevWT7rHp36hpgHvDuWvZg4Ce2L7f9su27KB0i+/ZRd6stbQ5JeBFYHdhc0ur1e+Fvbds9S9q/rpfkNnrzRH1/Sz9lxlEakftaC2wvAHqAZs/jo431z9U/x9b39YH7G+sXAXMHEOd4oMf2MwPYpt25wGckrQp8iNKr+bOl3bj2zp5F6a2lvl9lu68kfV59H9tcaPs22wfZnkjpxf4F8JPaiw6wA3BrW7099Zy1PE/jfNfP7ftqjTV7uv8ji4hYZu3t0aNt65/n/+3TBGAvSf9uvSiTv/qawNxqS18ZP2t7FvAN4CigR9K1WnwC7Rqk/et6SW5jMbbvBuYAe/VT7AlKT+OE1gJJYygToh5cyl09TElQW9uLcqu9ZQFljGvzlv96jb/vB9ZV35MDFvWxvOmXwH8o41mnADNtvzDA+qYD75e0CeUW2ln97O8eynFt1lcB2z3A0cAqwBZ1cfuQhMHaAngG+FcH6oqI6JS5wDTbb2i81rC9eW+Fazs9m7a21PaZtrehdBLcwuKdFVsAfye6WpLb6MsXKT2a36mTqZC0bp0AtkftMTwfmCppvTqO9kTKDP+/LuU+ZgCfk7SlpNcAR1IapJbZlETwIEkr1clRn2ysvwm4GTi7xraSpHfo/48qe4wydrhPjeM4jPKEgv6GJDxGSbYnNBfWR8xcDlxEGSpwbT/7Wwj8nMaEhjqJ4jBJG9RjGAscUeu6qY5h3pHOJLeTgStbT2qIiBghTgf2lLSTpNdIWkXSZpL6ezrPZby6Ld2qTgBejdL5Mh9Y2Fg/kXLX8fr2iqK7JLmNXtm+jjIRYDPgVknzKQPz1wV+W4sdTkkwbwQeoNw+2nkAidP5wA+AK4HHa903NGKYTxkP+1VKb+OXKY+haa1fBOxMSQL/AfybMsygdZvrJGBSvcV1ez9xnEt5vNm/bPeZmNce7dOBv9Y6mxMdzqCMHZvWNlSgN6cAUxoT7+YB2wN/pIwHu48ySe9jth+gPPXh8TrTd9Dqkxr2AU5dlnoiIjrN9m2Uf+K/Qhm+0EO5Kzaun81+COzauHs3ltK+PUmZC7EDsGej/AHA9GUcyhajgHqf0B0RA1F7c+8BJthe4rAMSTOB62yfsxRlTwYW2D5qSWWXUM/ngW1t770s9UREjBSSjqdMVDt6CeXWodzpm1TvtkUXS3IbsYxqj+hpwNq2PzUE9X8O+LXtOZ2uOyIiotskuY1YBnUm7m8pQwl2rI/1ioiIiGGS5DYiIiIiukYmlEVERERE10hyGxERERFdI8ltRERERHSNJLcRERER0TWS3EZERERE10hyGxERERFd43/dB9WfkOFr5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd151f3668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "ax0 = plt.subplot2grid((2, 2), (0, 0), rowspan=2)\n",
    "ax1 = plt.subplot2grid((2, 2), (0, 1))\n",
    "ax2 = plt.subplot2grid((2, 2), (1, 1))\n",
    "\n",
    "fs = 13  # fontsize\n",
    "matplotlib.rcParams['font.size'] = fs\n",
    "\n",
    "# Plot the model\n",
    "ax0.semilogx(sigma[active], mesh.vectorCCz[active], 'k-', lw=2)\n",
    "ax0.semilogx(np.exp(moptTD), mesh.vectorCCz[active], 'r*', ms=10)\n",
    "ax0.set_ylim(-700, 0)\n",
    "ax0.set_xlim(5e-3, 1e-1)\n",
    "\n",
    "ax0.set_xlabel('Conductivity (S/m)', fontsize=fs)\n",
    "ax0.set_ylabel('Depth (m)', fontsize=fs)\n",
    "ax0.grid(\n",
    "    which='both', color='k', alpha=0.5, linestyle='-', linewidth=0.2\n",
    ")\n",
    "ax0.legend(['True', 'TDEM'], fontsize=fs, loc=4)\n",
    "\n",
    "# plot the data misfits - negative b/c we choose positive to be in the\n",
    "# direction of primary\n",
    "dpred = surveyTD.dpred(moptTD)\n",
    "ax1.loglog(times, surveyTD.dobs, 'k-', lw=2)\n",
    "ax1.loglog(times, dpred, 'r*', ms=10)\n",
    "ax1.set_xlim(times.min(), times.max())\n",
    "\n",
    "# plot the difference\n",
    "ax2.loglog(times, np.abs(dpred-surveyTD.dobs), 'bo')\n",
    "ax2.set_xlim(times.min(), times.max())\n",
    "ax2.grid(which='both', alpha=0.5, linestyle='-', linewidth=0.2)\n",
    "ax2.set_xlabel('Time (s)', fontsize=fs)\n",
    "ax2.set_title('(c) |dobs - dpred|', fontsize=fs)\n",
    "\n",
    "# Labels, gridlines, etc\n",
    "ax1.grid(which='both', alpha=0.5, linestyle='-', linewidth=0.2)\n",
    "ax1.set_xlabel('Time (s)', fontsize=fs)\n",
    "ax1.set_ylabel('Vertical magnetic field (T)', fontsize=fs)\n",
    "ax1.legend((\"Obs\", \"Pred\"), fontsize=fs)\n",
    "\n",
    "ax0.set_title(\"(a) Recovered Models\", fontsize=fs)\n",
    "ax1.set_title(\"(b) TDEM observed vs. predicted\", fontsize=fs)\n",
    "\n",
    "plt.tight_layout(pad=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}