notebooks/TDEM_vmd_sounding_over_sphere.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TDEM Magnetic Dipole Sounding over a Sphere\n",
    "\n",
    "In this notebook, we will simulate an TDEM sounding over a conductive sphere. A cylindrical mesh and the [SimPEG](http://simpeg.xyz) Electromanetics module is used to perform the simulation. \n",
    "\n",
    "For more on SimPEG and SimPEG.EM see:\n",
    "\n",
    "- [(Cockett et al., 2015)](http://www.sciencedirect.com/science/article/pii/S009830041530056X): *SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications*\n",
    "- [(Heagy et al., 2016)](https://arxiv.org/abs/1610.00804): *A framework for simulation and inversion in electromagnetics*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Setup\n",
    "\n",
    "An inductive loop transmitter is centered over an electrically conductive sphere (radius 30m). The transmitter is 20m above the surface, and the center of the sphere is 50m below the surface. A coil receiver is offset by 8m horizontally from the transmitter. \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/simpeg/em_notebooks/master/images/cylsetup.png\" width=60%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting Started: Package Imports\n",
    "\n",
    "- [NumPy](http://www.numpy.org/), [SciPy](https://www.scipy.org/) and [Matplotlib](http://matplotlib.org/) are core Python packages that can be installed when you install python with a distribution such as [Anaconda](https://www.continuum.io/downloads). \n",
    "\n",
    "- [SimPEG](http://simpeg.xyz) is a Simulation and Inversion package for geophysics\n",
    "\n",
    "If you would like to set up and run SimPEG on your machine, see http://simpeg.xyz and http://tutorials.simpeg.xyz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic python packages\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.constants import mu_0\n",
    "\n",
    "# Modules of SimPEG we will use for forward modelling\n",
    "from SimPEG import Mesh, Utils, Maps\n",
    "from SimPEG.EM import TDEM\n",
    "from SimPEG import SolverLU as Solver\n",
    "\n",
    "# Set a nice colormap! \n",
    "plt.set_cmap(plt.get_cmap('viridis'))\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Parameters\n",
    "\n",
    "We define a \n",
    "- resistive halfspace and \n",
    "- conductive sphere \n",
    "   - radius of 30m \n",
    "   - center is 50m below the surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# electrical conductivities in S/m\n",
    "sig_halfspace = 1e-6\n",
    "sig_sphere = 1e-2\n",
    "sig_air = 1e-8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# depth to center, radius in m\n",
    "sphere_z = -50.\n",
    "sphere_radius = 30. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Survey Parameters\n",
    "\n",
    "- Transmitter and receiver 20m above the surface \n",
    "- Receiver offset from transmitter by 8m horizontally\n",
    "- We will sample 30 times from $10^{-7}$s to $10^{-4}$s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "boom_height = 20. \n",
    "rx_offset = 8. \n",
    "\n",
    "times = np.logspace(-7, -4, 30)\n",
    "\n",
    "# source and receiver location in 3D space\n",
    "src_loc = np.r_[0., 0., boom_height]\n",
    "rx_loc = np.atleast_2d(np.r_[rx_offset, 0., boom_height])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min diffusion distance: 3.61e+00 m\n",
      "max diffusion distance: 1.14e+04 m\n"
     ]
    }
   ],
   "source": [
    "# print the min and max diffusion distances to make sure mesh is \n",
    "# fine enough and extends far enough \n",
    "\n",
    "def diffusion_distance(sigma, time):\n",
    "    return 1.28*np.sqrt(time/(sigma * mu_0))\n",
    "\n",
    "print(\n",
    "    'min diffusion distance: {:.2e} m'.format(diffusion_distance(sig_sphere, times.min()))\n",
    ")\n",
    "print(\n",
    "    'max diffusion distance: {:.2e} m'.format(diffusion_distance(sig_halfspace, times.max()))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mesh\n",
    "\n",
    "Here, we define a cylindrically symmetric tensor mesh. \n",
    "\n",
    "The figure below shows a cell in a cartesian mesh (a) and a cylindrically symmetric mesh (b). Note that edges are rotational and faces are radial and vertical in the cylindrically symmetric mesh. \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/simpeg/em_notebooks/master/images/finiteVolume-02.png\" width=60%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh Parameters\n",
    "\n",
    "For the mesh, we will use a cylindrically symmetric tensor mesh. To construct a tensor mesh, all that is needed is a vector of cell widths in the x and z-directions. We will define a core mesh region of uniform cell widths and a padding region where the cell widths expand \"to infinity\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x-direction\n",
    "csx = 2  # core mesh cell width in the x-direction\n",
    "ncx = np.ceil(1.2*sphere_radius/csx)  # number of core x-cells (uniform mesh slightly beyond sphere radius)\n",
    "npadx = 30  # number of x padding cells\n",
    "\n",
    "# z-direction\n",
    "csz = 1  # core mesh cell width in the z-direction\n",
    "ncz = np.ceil(1.2*(boom_height - (sphere_z-sphere_radius))/csz) # number of core z-cells (uniform cells slightly below bottom of sphere)\n",
    "npadz = 32  # number of z padding cells\n",
    "\n",
    "# padding factor (expand cells to infinity)\n",
    "pf = 1.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cell spacings in the x and z directions\n",
    "hx = Utils.meshTensor([(csx, ncx), (csx, npadx, pf)])\n",
    "hz = Utils.meshTensor([(csz, npadz, -pf), (csz, ncz), (csz, npadz, pf)])\n",
    "\n",
    "# define a SimPEG mesh\n",
    "mesh = Mesh.CylMesh([hx, 1, hz], x0 = np.r_[0.,0., -hz.sum()/2.-boom_height])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the mesh\n",
    "\n",
    "Below, we plot the mesh. The cyl mesh is rotated around x=0. Ensure that each dimension extends beyond the maximum skin depth. \n",
    "\n",
    "Zoom in by changing the xlim and zlim. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum diffusion distance (in background) is: 1.14e+04 m. Does the mesh go sufficiently past that?\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEWCAYAAADcsGj7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAG4hJREFUeJzt3X20ZWV92PHvT8igEXFAEOhAZAgjDdouReWlsfYusDCQtGO6sGLTMCrtWJWupNpVMXRVo9JK0kRLNOi4oA7WCNSXxSRBAZG71AYQSBUhOM6VF5lAIHQYmAF8AX/9Yz+Xuzlz7vs5Zz937vez1ll3n2c/e+9n/2bf+5v97GfvHZmJJEk1ek7XDZAkaTomKUlStUxSkqRqmaQkSdUySUmSqmWSkiRVyyQlARHxmxFxzZDW/ZmI+PAilt8VEUcOsk2jEBH3RMTru26HljaTlJaNiHhtRPxlRDwaEdsj4v9ExGsAMvNzmXlKBW0cj4h/0y7LzH0z864hbOueiPhpRBzYU/6diMiIOGLQ25TmyySlZSEi9gP+HPhj4ABgFfB7wE+6bFcF7gbePPklIv4B8LzumiM9m0lKy8VLATLz85n5dGY+mZnXZOZtABHxloj41mTlcibxzojYGhE7I+JDEfHLEXFDRDwWEVdExIp+y7aWP6q3ERGxf0T8eUT8XUQ8UqYPK/POB/4x8PHSxffx3nVFxAsj4tKy/L0R8Z8j4jntdkTEfy/rvjsiTpslLp8Fzmp9Xw9c2tPmfco6fxQRD0bEJyPieWXegWUfdpSz029Otqd4RUTcVs5eL4+I587SHulZTFJaLn4APB0RmyLitIjYfw7LrAVeBZwA/CdgI/CbwOHAy2mdgczDc4D/CbwE+CXgSeDjAJl5HvBN4JzSxXdOn+X/GHghcCTwT2gSzFtb848HtgAHAr8PXBwRMUN7bgT2i4hfiYi9gDcB/6unzgU0Sf4VwFE0Z6H/pcx7D7ANOAg4GPhdoP2stX9JE8fVwD8E3jJDW6TdmKS0LGTmY8Braf6Afhr4u4jYHBEHz7DYBZn5WGbeAdwOXJOZd2Xmo8BXgFcuoB3/LzO/mJlPZOZO4HyaZDOrVhJ5X2buzMx7gD8EfqtV7d7M/HRmPg1sAg6lSR4zmTyb+qfA94G/aW0zgH8L/IfM3F7a/F+BM0uVn5VtvCQzf5aZ38xnPxD0wsy8PzO3A39Gk+ikOTNJadnIzDsz8y2ZeRjNmdDfAz42wyIPtqaf7PN93/m2ISJ+MSI+VbrqHgO+AawsCWg2BwIrgHtbZffSnNlM+tvJicx8okzO1s7PAv+K5izn0p55BwG/CNxauvR2AF8t5QB/AEwA10TEXRFxbs/yf9uafmIObZGexSSlZSkzvw98hiZZLdbjNH/IAYiIQ2ao+x7gaOD4zNwPeN3kYpNNm2HZh2nOXF7SKvslWmc+C5GZ99IMoDgd+FKfbT4JvCwzV5bPCzNz37Lszsx8T2YeCfwz4N0RcfJi2iO1maS0LETE34+I97QGKRxOc03pxgGs/rvAyyLiFWVgwAdmqPsCmj/6OyLiAOD9PfMfpLnetJvShXcFcH5EvCAiXgK8m92vIS3E2cBJmfl4zzZ/TtM9+tGIeDFARKyKiFPL9K9HxFGlW/Ax4OnykQbCJKXlYifNoIKbIuJxmuR0O82ZzaJk5g+ADwJfA7YC35qh+sdohng/XNrw1Z75/wM4o4zOu7DP8v+e5sztrrKdPwUuWdQOAJn5w8y8ZZrZ76Xp0ruxdFF+jeZsEGBN+b4LuAH4k8wcX2x7pEnhSw8lSbXyTEqSVC2TlCSpWiYpSVK1TFKSpGrt3XUDarNy5co86qjdHrm2LD3++OM8//nP77oZVTAWU4zFFGMx5dZbb304Mw+aveb8mKR6HHzwwdxyy3QjcZeX8fFxxsbGum5GFYzFFGMxxVhMiYh7Z681f3b3SZKqZZKSJFXLJCVJqpZJSpJULZOUJKlaJilJUrVMUpKkapmkJEnVMklJkqplkpIkVcskJUmqlklKklQtk5QkqVomKUlStUxSkqRqmaQkSdUySUmSqtVZkoqIwyPi+oi4MyLuiIjfLuUHRMS1EbG1/Ny/lEdEXBgRExFxW0Qc21rX+lJ/a0Ssb5W/KiK+V5a5MCJi9HsqSVqoLs+kngLek5m/ApwAvCsijgHOBa7LzDXAdeU7wGnAmvLZAFwETVID3g8cDxwHvH8ysZU6G1rLrR3BfkmSBqSzJJWZD2TmX5XpncCdwCpgHbCpVNsEvKFMrwMuzcaNwMqIOBQ4Fbg2M7dn5iPAtcDaMm+/zLwhMxO4tLUuSdISsHfXDQCIiCOAVwI3AQdn5gPQJLKIeHGptgq4r7XYtlI2U/m2PuX9tr+B5oyLgw46iPHx8UXtz55i165dxqIwFlOMxRRjMXydJ6mI2Bf4IvA7mfnYDJeN+s3IBZTvXpi5EdgIcPTRR+fY2NgsrV4exsfHMRYNYzHFWEwxFsPXaZKKiF+gSVCfy8wvleIHI+LQchZ1KPBQKd8GHN5a/DDg/lI+1lM+XsoP61N/Rvc89nOOOPcv5tT+41cfwE13b39mGtjt+3Rlc5k3nzoLqTvbcjt2PMlFW26Y13oW04battHWG4tRb38hhtXGhR4Xg1JT7LuORdtkXO75yK913JLB6nJ0XwAXA3dm5h+1Zm0GJkforQeubJWfVUb5nQA8WroFrwZOiYj9y4CJU4Cry7ydEXFC2dZZrXVJkpaALs+kfhX4LeB7EfGdUva7wEeAKyLibOBHwBvLvKuA04EJ4AngrQCZuT0iPgTcXOp9MDO3l+l3AJ8Bngd8pXwWpN//2i5/+4m86VM3PDMN7PZ9urK5zJtPnYXUnW25pitjfutZTBtq20ZbbyxGvf2FGFYbF3pcDEpNse86Fm2TcdnTdJakMvNb9L9uBHByn/oJvGuadV0CXNKn/Bbg5YtopiSpQz5xQpJULZOUJKlaJilJUrVMUpKkapmkJEnVMklJkqplkpIkVcskJUmqlklKklQtk5QkqVomKUlStUxSkqRqmaQkSdUySUmSqmWSkiRVyyQlSapWNO8S1KR9Dl2Th67/2JzqHr/6AG66e/sz08Bu36crm8u8+dRZSN3ZltuxYwcrV66c13oW04battHWG4tRb38hhtXGhR4Xg1JT7LuORdtkXO75yK91sv2IuDUzXz3o9XomJUmqVmevj19q+v2v7fK3n8ibPnXDM9PAbt+nK5vLvPnUWUjd2ZYbHx9nbGx+61lMG2rbRltvLEa9/YUYVhsXelwMSk2x7zoWbZNx2dN4JiVJqpZJSpJULZOUJKlaJilJUrVMUpKkapmkJEnVMklJkqplkpIkVcskJUmqlklKklQtk5QkqVomKUlStXxVRw9f1eGrOvrxVR1Tun49RU2x7zoWbb6qQ5KkEfNVHXPkqzoWt65B81Uds/NVHcPXdSzafFXHEETEJRHxUETc3io7ICKujYit5ef+pTwi4sKImIiI2yLi2NYy60v9rRGxvlX+qoj4XlnmwoiI0e6hJGkxuu7u+wywtqfsXOC6zFwDXFe+A5wGrCmfDcBF0CQ14P3A8cBxwPsnE1ups6G1XO+2JEkV6zRJZeY3gO09xeuATWV6E/CGVvml2bgRWBkRhwKnAtdm5vbMfAS4Flhb5u2XmTdkMzrk0ta6JElLQNdnUv0cnJkPAJSfLy7lq4D7WvW2lbKZyrf1KZckLRFLaeBEv+tJuYDy3VccsYGmW5AVhxzVd+M7duzYrWx8fJwdO558Zrqp9+zv05XNZd586iyk7mzL7dq1a97rWUwbattGW28sRr39hRhWGxd6XAxKTbHvOhZtNcVlkGpMUg9GxKGZ+UDpsnuolG8DDm/VOwy4v5SP9ZSPl/LD+tTfTWZuBDZCc59Uvzr97oUYGzuRi7bc8Mw0sNv36crmMm8+dRZSd7blmpFLY/Naz2LaUNs22npjMertL8Sw2rjQ42JQaop917Foqykug1Rjd99mYHKE3nrgylb5WWWU3wnAo6U78GrglIjYvwyYOAW4uszbGREnlFF9Z7XWJUlaAjo9k4qIz9OcBR0YEdtoRul9BLgiIs4GfgS8sVS/CjgdmACeAN4KkJnbI+JDwM2l3gczc3IwxjtoRhA+D/hK+UiSlggfi9TDxyL5WKR+fCzSlK4fBVRT7LuORZuPRZIkacRqHDhRJR+LtLh1DZqPRZqdj0Uavq5j0eZjkSRJGjGTlCSpWiYpSVK1TFKSpGqZpCRJ1TJJSZKq5c28PbyZ15t5+/Fm3ild38BaU+y7jkWbN/NKkjRi3sw7R97Mu7h1DZo3887Om3mHr+tYtHkzryRJI2aSkiRVyyQlSaqWo/t6OLrP0X39OLpvStcj2mqKfdexaHN0nyRJI+bovjlydN/i1jVoju6bnaP7hq/rWLQ5uk+SpBEzSUmSqmWSkiRVy9F9PRzd5+i+fhzdN6XrEW01xb7rWLQ5uk+SpBFzdN8cObpvcesaNEf3zc7RfcPXdSzaHN0nSdKIeU2qh9ekvCbVj9ekpnR9Haam2HcdizavSUmSNGJek5ojr0ktbl2D5jWp2XlNavi6jkXbnnpNyu6+Hnb32d3Xj919U7ru4qop9l3Hos3uPkmSRswkJUmqlklKklQtk5QkqVoOnOgxn4ETklQbB05IkjQiJilJUrX2+CQVEWsjYktETETEuV23R5I0d3t0koqIvYBPAKcBxwBvjohjum2VJGmu5pSkIuK6iDi9p2zjcJo0UMcBE5l5V2b+FLgMWNdxmyRJczTXZ/etBt4bEa/JzN8rZQMfxTEEq4D7Wt+3Acf3VoqIDcAGgBWHHDWalknSEIyPj3fdhIGaa5LaAZwMXBgRfwb86+E1aaCiT9luY+4zcyOwEZoh6MNulCQNy9jYWNdNGKi5XpOKzHwqM98JfBH4FvDi4TVrYLYBh7e+Hwbc31FbJEnzNNck9cnJicz8DPAW4JohtGfQbgbWRMTqiFgBnAls7rhNkqQ5mlN3X2Z+quf7rcDbhtKiAcrMpyLiHOBqYC/gksy8o+NmSZLmaI9/6WFmXgVc1XU7JEnzt8cnqWHypYeDbUNt22jzpYdTun7RX02x7zoWbZNx2dPs0TfzSpKWNpOUJKladvfNUb+uhcvffiJv+tQNz0wDu32frmwu8+ZTZyF1Z1tufHycsbH5rWcxbahtG229sRj19hdiWG1c6HExKDXFvutYtE3GZU/j+6R6zOd9Ul6TGmwbattGm9ekpnR9Haam2Hcdi7bJuPg+KUmSRsTuvjmyu29x6xo0u/tmZ3ff8HUdiza7+5YJu/vs7uvH7r4pXXdx1RT7rmPRZnefJEkjZnffHNndt7h1DZrdfbOzu2/4uo5F257a3eeZlCSpWl6T6uE1Ka9J9eM1qSldX4epKfZdx6LNa1KSJI2Y16TmyGtSi1vXoHlNanZekxq+rmPR5jUpSZJGzCQlSaqWAyd6OHDCgRP9OHBiSteDBWqKfdexaHPghCRJI+bAiTly4MTi1jVoDpyYnQMnhq/rWLQ5cEKSpBEzSUmSqmWSkiRVyyQlSaqWQ9B7OATdIej9OAR9StfDrmuKfdexaHMIuiRJI+YQ9DlyCPri1jVoDkGfnUPQh6/rWLQ5BF2SpBEzSUmSqmWSkiRVyyQlSaqWSUqSVC3vk+rhfVLeJ9WP90lN6freoJpi33Us2rxPSpKkEfM+qTnyPqnFrWvQvE9qdt4nNXxdx6LN+6QGKCLeGBF3RMTPI+LVPfPeFxETEbElIk5tla8tZRMRcW6rfHVE3BQRWyPi8ohYUcr3Kd8nyvwjRrV/kqTB6Kq773bgXwDfaBdGxDHAmcDLgLXAn0TEXhGxF/AJ4DTgGODNpS7ABcBHM3MN8Ahwdik/G3gkM48CPlrqSZKWkE6SVGbemZlb+sxaB1yWmT/JzLuBCeC48pnIzLsy86fAZcC6iAjgJOALZflNwBta69pUpr8AnFzqS5KWiNoGTqwC7mt931bKpit/EbAjM5/qKX/Wusr8R0t9SdISMbSBExHxNeCQPrPOy8wrp1usT1nSP5nmDPVnWtfuG43YAGwAWHHIUX0btmPHjt3KxsfH2bHjyWemm3rP/j5d2VzmzafOQurOttyuXbvmvZ7FtKG2bbT1xmLU21+IYbVxocfFoNQU+65j0VZTXAZpaEkqM1+/gMW2AYe3vh8G3F+m+5U/DKyMiL3L2VK7/uS6tkXE3sALge3TtHUjsBGa+6T61el3L8TY2IlctOWGZ6aB3b5PVzaXefOps5C6sy3XjFwam9d6FtOG2rbR1huLUW9/IYbVxoUeF4NSU+y7jkVbTXEZpNq6+zYDZ5aReauBNcC3gZuBNWUk3wqawRWbs7kT+XrgjLL8euDK1rrWl+kzgK+ndy5L0pLS1RD034iIbcCJwF9ExNUAmXkHcAXw18BXgXdl5tPlLOkc4GrgTuCKUhfgvcC7I2KC5prTxaX8YuBFpfzdwDPD1iVJS0MnN/Nm5peBL08z73zg/D7lVwFX9Sm/i2b0X2/5j4E3LrqxkqTO+Oy+Hj67z2f39eOz+6Z0/by6mmLfdSzafHafJEkj5rP75shn9y1uXYPms/tm57P7hq/rWLT57D5JkkbMJCVJqpZJSpJULZOUJKlaJilJUrVMUpKkapmkJEnVMklJkqplkpIkVcskJUmqlklKklQtk5QkqVq+qqOHr+rwVR39+KqOKV2/nqKm2HcdizZf1SFJ0oj5qo458lUdi1vXoPmqjtn5qo7h6zoWbb6qQ5KkETNJSZKqZZKSJFXLJCVJqpZJSpJULZOUJKlaJilJUrVMUpKkapmkJEnVMklJkqplkpIkVcskJUmqlklKklQtk5QkqVomKUlStUxSkqRq+fr4Hr4+3tfH9+Pr46d0/cr0mmLfdSzafH38AEXEH0TE9yPitoj4ckSsbM17X0RMRMSWiDi1Vb62lE1ExLmt8tURcVNEbI2IyyNiRSnfp3yfKPOPGOU+SpIGIDNH/gFOAfYu0xcAF5TpY4DvAvsAq4EfAnuVzw+BI4EVpc4xZZkrgDPL9CeBd5TpdwKfLNNnApfPpW0vfelLU43rr7++6yZUw1hMMRZTjMUU4JYcQr7o5EwqM6/JzKfK1xuBw8r0OuCyzPxJZt4NTADHlc9EZt6VmT8FLgPWRUQAJwFfKMtvAt7QWtemMv0F4ORSX5K0ROzddQOAtwGXl+lVNElr0rZSBnBfT/nxwIuAHa2E166/anKZzHwqIh4t9R/ubUBEbAA2ABx00EGMj48vbo/2ELt27TIWhbGYYiymGIvhG1qSioivAYf0mXVeZl5Z6pwHPAV8bnKxPvWT/tfOcob6M61r98LMjcBGgKOPPjrHxsb6VVt2xsfHMRYNYzHFWEwxFsM3tCSVma+faX5ErAd+HTi59GdCcyZ0eKvaYcD9Zbpf+cPAyojYu5xNtetPrmtbROwNvBDYvvA9kiSNWlej+9YC7wX+eWY+0Zq1GTizjMxbDawBvg3cDKwpI/lW0AyE2FyS2/XAGWX59cCVrXWtL9NnAF9vJUNJ0hLQ1TWpj9OM4Lu2jGW4MTP/XWbeERFXAH9N0w34rsx8GiAizgGuphnpd0lm3lHW9V7gsoj4MPB/gYtL+cXAZyNiguYM6szR7JokaVA6SVKZedQM884Hzu9TfhVwVZ/yu2hG//WW/xh44+JaKknqko9FkiRVyyQlSaqWSUqSVC2TlCSpWiYpSVK1TFKSpGqZpCRJ1TJJSZKqZZKSJFXLJCVJqpZJSpJULZOUJKlaJilJUrVMUpKkapmkJEnVMklJkqoVvlH92SJiJ7Cl63ZU4kDg4a4bUQljMcVYTDEWU47OzBcMeqVdvT6+Zlsy89VdN6IGEXGLsWgYiynGYoqxmBIRtwxjvXb3SZKqZZKSJFXLJLW7jV03oCLGYoqxmGIsphiLKUOJhQMnJEnV8kxKklQtk5QkqVomqSIi1kbEloiYiIhzu27PsETEPRHxvYj4zuSQ0Yg4ICKujYit5ef+pTwi4sISk9si4tjWetaX+lsjYn1X+zMfEXFJRDwUEbe3yga27xHxqhLbibJsjHYP526aWHwgIv6mHBvfiYjTW/PeV/ZrS0Sc2irv+3sTEasj4qYSo8sjYsXo9m5+IuLwiLg+Iu6MiDsi4rdL+bI7NmaIRXfHRmYu+w+wF/BD4EhgBfBd4Jiu2zWkfb0HOLCn7PeBc8v0ucAFZfp04CtAACcAN5XyA4C7ys/9y/T+Xe/bHPb9dcCxwO3D2Hfg28CJZZmvAKd1vc/zjMUHgP/Yp+4x5XdiH2B1+V3Za6bfG+AK4Mwy/UngHV3v8wyxOBQ4tky/APhB2edld2zMEIvOjg3PpBrHAROZeVdm/hS4DFjXcZtGaR2wqUxvAt7QKr80GzcCKyPiUOBU4NrM3J6ZjwDXAmtH3ej5ysxvANt7igey72Xefpl5Qza/fZe21lWdaWIxnXXAZZn5k8y8G5ig+Z3p+3tTzhJOAr5Qlm/HtTqZ+UBm/lWZ3gncCaxiGR4bM8RiOkM/NkxSjVXAfa3v25j5H2YpS+CaiLg1IjaUsoMz8wFoDlLgxaV8urjsSfEa1L6vKtO95UvNOaUL65LJ7i3mH4sXATsy86me8upFxBHAK4GbWObHRk8soKNjwyTV6Nc/vKeOzf/VzDwWOA14V0S8boa608VlOcRrvvu+J8TkIuCXgVcADwB/WMqXRSwiYl/gi8DvZOZjM1XtU7ZHxaNPLDo7NkxSjW3A4a3vhwH3d9SWocrM+8vPh4Av05yWP1i6JCg/HyrVp4vLnhSvQe37tjLdW75kZOaDmfl0Zv4c+DTNsQHzj8XDNF1ge/eUVysifoHmj/LnMvNLpXhZHhv9YtHlsWGSatwMrCmjTlYAZwKbO27TwEXE8yPiBZPTwCnA7TT7OjkSaT1wZZneDJxVRjOdADxauj2uBk6JiP3Laf8ppWwpGsi+l3k7I+KE0u9+VmtdS8LkH+TiN2iODWhicWZE7BMRq4E1NAMB+v7elOsu1wNnlOXbca1O+fe6GLgzM/+oNWvZHRvTxaLTY6Pr0SS1fGhG7PyAZkTKeV23Z0j7eCTNKJvvAndM7idNP/F1wNby84BSHsAnSky+B7y6ta630VwknQDe2vW+zXH/P0/TVfEzmv/pnT3IfQdeXX55fwh8nPJElxo/08Tis2Vfbyt/fA5t1T+v7NcWWiPTpvu9Kcfat0uM/jewT9f7PEMsXkvT5XQb8J3yOX05HhszxKKzY8PHIkmSqmV3nySpWiYpSVK1TFKSpGqZpCRJ1TJJSZKqZZKSJFXLJCVJqpZJSqpcRLymPNjzueWpIXdExMu7bpc0Ct7MKy0BEfFh4LnA84BtmfnfOm6SNBImKWkJKM8/uxn4MfCPMvPpjpskjYTdfdLScACwL83bUp/bcVukkfFMSloCImIzzdtNV9M83POcjpskjcTes1eR1KWIOAt4KjP/NCL2Av4yIk7KzK933TZp2DyTkiRVy2tSkqRqmaQkSdUySUmSqmWSkiRVyyQlSaqWSUqSVC2TlCSpWv8fIqa9N9hw8dIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd15007cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# X and Z limits we want to plot to. Try \n",
    "xlim = np.r_[0., 2.5e4]\n",
    "zlim = np.r_[-2.5e4, 2.5e4]\n",
    "\n",
    "fig, ax = plt.subplots(1,1)\n",
    "mesh.plotGrid(ax=ax)\n",
    "\n",
    "ax.set_title('Simulation Mesh')\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(zlim)\n",
    "\n",
    "print(\n",
    "    'The maximum diffusion distance (in background) is: {:.2e} m. '\n",
    "    'Does the mesh go sufficiently past that?'.format(\n",
    "        diffusion_distance(sig_halfspace, times.max())\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Put Model on Mesh\n",
    "\n",
    "Now that the model parameters and mesh are defined, we can define electrical conductivity on the mesh. \n",
    "\n",
    "The electrical conductivity is defined at cell centers when using the finite volume method. So here, we define a vector that contains an electrical conductivity value for every cell center. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a vector that has one entry for every cell center\n",
    "sigma = sig_air*np.ones(mesh.nC)  # start by defining the conductivity of the air everwhere\n",
    "sigma[mesh.gridCC[:,2] < 0.] = sig_halfspace  # assign halfspace cells below the earth\n",
    "\n",
    "sigma_background = sigma.copy()\n",
    "\n",
    "# indices of the sphere (where (x-x0)**2 + (z-z0)**2 <= R**2)\n",
    "sphere_ind = (mesh.gridCC[:,0]**2 + (mesh.gridCC[:,2] - sphere_z)**2) <= sphere_radius**2 \n",
    "sigma[sphere_ind] = sig_sphere  # assign the conductivity of the sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Conductivity Model')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHwFJREFUeJzt3XmYJFWd7vHva4Mg4MIqWyOgrQ6o03dsGXdlFX1UFC8OqAhuLSp3xm1UxAUXXHAQx3FtvSg6IoKItMIINA7XZQahUQZpAWkQpReBBlkEaXt57x8RpdnZlVVZUVkZkZnv53ni6YwTGRG/qIT85TknzgnZJiIioooH1B1AREQMriSRiIioLEkkIiIqSxKJiIjKkkQiIqKyJJGIiKgsSSRqJeloST+ZgePuJumPkmZ18d4vSHpvr2PoJUmXSHptl++1pEfNdEwRkCQSHUh6maTF5RfxSkn/IenpdcfViaSbJB0wtm77d7a3sr1usn1tH2P7Q+Vxni1p2TTiOKH8Ev/HtvI3l+UnVD12RBMlicRGJL0V+BTwEeDhwG7A54BD6oxrgPwaOKqt7JVlecRQSRKJDUh6KPBB4E22v2P7XttrbH/P9j+X79lM0qckrSiXT0narNz2bEnLJL1N0q1lLeZVLcffVtJCSXdLugx4ZMu23ctf65u0lG3QjCPpdZKukXSPpF9J+jtJX6dIdN8ra07vaD2WpMMlLW67zrdIWli+/qqkD0vaEvgPYOfyOH+UtLOk+yRt27LvEyXdJmnTDn/Gy4EtJO1dvn9v4EFleWsMr5O0VNId5d9k55ZtB0q6VtJdkj4DqG3fV5d/hz9IukDSIzp9phEzKUkk2j0F2Bw4Z4L3HA88GZgL/C2wD/Celu07Ag8FdgFeA3xW0tblts8C9wM7Aa8ul65IOgw4geJX/UOAFwK32z4S+B3wgrIJ66S2XRcCj5E0p6XsZcDprW+yfS/wXGBFeZytbK8ALgFe2vLWVwBn2F4zQbhfL+OEolbytbZr2Q/4aHncnYDfAmeU27YDzqb4m24H3AA8rWXfFwHvBg4Ftgd+DHxzglgiZkySSLTbFlhle+0E73k58EHbt9q+DfgAcGTL9jXl9jW2zwf+SPElPgt4CfC+soZzNXDaFGJ7LXCS7ctdWGr7t5PtZPs+4FzgCIAymTyWIrl04zSKxEF5DUdQJImJ/DtwRFlbObxcb/Vy4FTbP7e9GjgOeIqk3YHnAb+y/e0yUX0K+H3Lvq8HPmr7mvJz+ggwN7WRqEOSSLS7HdiutUlpHDtT/HIe89uy7C/HaEtC9wFbUfxq3gS4uW3fbs2m+FVexemUSYSiFvLdMrl041xgL0l7AgcCd9m+bKIdbP8OWErxBX+97Zvb3rLB39D2Hyn+9ruU225u2WY2/Js9AvhXSXdKuhO4g6K5a5curyeiZ5JEot1/UzQ3vWiC96yg+CIbs1tZNpnbgLUUyaB13zH3lv9u0VK2Y8vrm2npQ2kz2XTUF1Ikx7kUyeT0Du/b6Di27wfOpKg9HMnktZAxXwPeRltTVmmDv2HZH7MtsBxYScvfSJLY8G92M/B62w9rWR5k+7+6jCuiZ5JEYgO27wLeR9GP8SJJW0jaVNJzJY31NXwTeI+k7cv2+/excXPNeMdeB3wHOKE87l603MVUNo0tB14haZakV7Nh0vgy8PayY1uSHtXShHMLsOcE514LfBv4BLANcFGHt94CbFveYNDqa8DRFP0wk15r6VvAQRQJqN3pwKskzS1vSvgI8DPbNwHnAXtLOrSsEf4jGybTLwDHtXTcP7TsL4rouySR2IjtTwJvpejYvY3il++xwHfLt3wYWAxcBfwS+HlZ1o1jKZq2fg98FfhK2/bXAf9M0bSzN/CXX9e2zwJOpPgCvqeMZ5ty80cpEtudkt7e4dynAwcAZ3Xq87F9LUWSvLE81s5l+U+B9cDPyy/6Sdn+k+1Ftv80zraLgfdSdKCvpEiWh5fbVgGHAR8r/w5zgJ+27HsO8HHgDEl3A1dT3BAQ0XfKQ6kiuiPph8Dptr9cdywRTZEkEtEFSU+iaAKbbfueuuOJaIo0Z0VMQtJpwCLgzUkg0QSSPlEORr1K0jmSHlZbLKmJREQMFkkHAT+0vVbSxwFsv7OOWFITiYgYMLYvbLk55FJg17pimWhA2VB4oDbz5mxZdxgRMQDu4Q+rbG9fdf/n7Lulb79j0omjAbjiqtVLKMZkjVlge0GF076a4nbyWgx9EtmcLfl77V93GBExABb521OZQWEjt9+xjssu2G3yNwKzdrr+ftvzOm2XtIgNxweNOd72ueV7jqcYwPuNCuH2xNAnkYiIfjGwnvW9OZZ9wETbJR0FPB/Y3zV2bieJRET0iDFrJn8O2rRJOhh4J/CsKcwBNyOSRCIieqhXNZFJfAbYDLiomFqNS20f048Tt0sSiYjoEWPW9aFlyfajZvwkXUoSiYjoofWTTig9XJJEIiJ6xMC6JJGIiKgqNZGIiKjEwJoRm0oqSSQiokeM05zVFJJuonjw0Dpgre15krahGN6/O3AT8FLbf6grxoiIDRjWjVYOafwEjPvantsyNcC7gIttzwEuLtcjIhqhGLHe3TIsmp5E2h0CnFa+Pg14UY2xRES0Eeu6XIZFY5uzKJL6hZIMfLGc3fLhtlcC2F4paYfxdpQ0H5gPsDlb9CveiBhxRcf68CSIbjQ5iTzN9ooyUVwk6dpudywTzgKAh2ibEWuhjIi6FONEkkQawfaK8t9bJZ0D7APcImmnshayE3BrrUFGRLRZP2I1kUb2iUjaUtKDx14DBwFXAwuBo8q3HQWcW0+EEREbG6uJpE+kfg8Hzilnp9wEON32DyRdDpwp6TXA74DDaowxImIDRqxr5m/zGdPIJGL7RuBvxym/HchjCiOisUatOauRSSQiYhAZ8WfPqjuMvkoSiYjokWKwYZqzIiKiomHqNO9GkkhERI/YYp1TE4mIiIrWpyYSERFVFB3ro/W1OlpXGxExg9KxHhER07Iu40QiIqKKjFiPiIhpWZ+7syIioopiAsYkkYiIqMCINZn2JCIiqrDJYMOIiKhKfRlsKOlDwCHAeoqH8x099iC/fhutlBkRMYNMURPpZpmmT9h+gu25wPeB9007+IpSE4mI6KF+dKzbvrtldUuK/FWLJJGIiB4x6ttDqSSdCLwSuAvYty8nHUeasyIiesTAGm/S1QJsJ2lxyzK/9ViSFkm6epzlEADbx9ueDXwDOLbvF1tKTSQiomc0leeJrLI9r9NG2wd0eZzTgfOA93d74l5qZE1E0mxJ/ynpGklLJP1TWX6CpOWSriyX59Uda0TEGFOMWO9mmQ5Jc1pWXwhcO60DTkNTayJrgbfZ/rmkBwNXSLqo3HaK7X+pMbaIiI769GTDj0l6DMUtvr8FjunHScfTyCRieyWwsnx9j6RrgF3qjSoiYmK2+jJ3lu2XzPhJutTI5qxWknYH/hfws7LoWElXSTpV0tYd9pk/1lm1htV9ijQiRl3RsT6rq2VYNDqJSNoKOBt4c3lf9OeBRwJzKWoqJ4+3n+0FtufZnrcpm/Ut3ogYderXYMPGaGRzFoCkTSkSyDdsfwfA9i0t279EMVIzIqIRio71PJSqdpIE/F/gGtufbCnfqewvAXgxcHUd8UVEdJKp4JvhacCRwC8lXVmWvRs4QtJcioR/E/D6esKLiNhYP0esN0Ujk4jtn8C498md3+9YIiKmYn1qIhERUYUNa9YniURERAVFc1aSSEREVNSnEeuNkSQSEdEjucU3IiKmIc1ZERExDf14xnqTJIlERPRIcXfW8MyL1Y0kkYiIHslgw4iImJY0Z0VERCW5OysiIqYld2dFREQltlibJBIREVWlOSsiIipJn0hERExLkkhERFSScSIRETEtGScSERGV2LA2D6WKiIiqRq05a+BSpqSDJV0naamkd9UdT0TEmLE+kW6WYTFQSUTSLOCzwHOBvYAjJO1Vb1QREX9lq6tlWAxUEgH2AZbavtH2n4EzgENqjiki4i/Wo66WYTFoSWQX4OaW9WVl2QYkzZe0WNLiNazuW3ARMdps+tqcJentkixpu54csIJB61gf7y/vjQrsBcACgC12mO3bX/LUmY4rIobBF749zQOIdX26O0vSbOBA4Hd9OWEHg1YTWQbMblnfFVhRUywRERvpY5/IKcA7GOeHdD8NWk3kcmCOpD2A5cDhwMvqDSkiojDFubO2k7S4ZX1B2YoyKUkvBJbb/h+p3v6VgUoittdKOha4AJgFnGp7Sc1hRUQUXPSLdGmV7XmdNkpaBOw4zqbjgXcDB005vhkwUEkEwPb5wPl1xxERMZ5e3Xll+4DxyiU9HtgDGKuF7Ar8XNI+tn/fk5NPwcAlkYiIpnIfOtZt/xLYYWxd0k3APNurZvTEHSSJRET00BSas4ZCkkhERA/1ezS67d37esI2SSIRET1i9z+J1C1JJCKih4ZpcsVuJIlERPRQ+kQiIqISI9bnoVQREVHViFVEkkQiInomHesRETEtI1YVSRKJiOih1EQiIqISA+vXJ4lEREQVBlITiYiIqjJOJCIiqhuxJDJao2IiImZUd4/GrbPzXdI7JT24V8dLTSQiopeaXxP5CHCVpE2BlbYvn87BkkRiKC1+/+cr7TfvA2/ocSQxUgxu/t1ZAr4P3A1sIen3wBvKp8ZOWZJIDIWqSWOy4ySpxNQ1PokAHG77LEmbA4cBX5H0MtsXT/VAjUsikj4BvAD4M3AD8Crbd0raHbgGuK5866W2j6klyGiMXiWPbo6fhBJdaX5z1r3AtwFs3w98XdIK4APA4CcR4CLgONtrJX0cOA54Z7ntBttz6wstmmCmE8dk500yiQk1P4ncBDwZ+O+Wsv8EzqpysMYlEdsXtqxeCvzvumKJZqkrebRLMomOBmOw4aeBsyS9ETjP9jrgUOCOKgdrXBJp82rgWy3re0j6BUWH0Hts/3i8nSTNB+YDbLrV1jMeZMyspiSPdkkmMZ6mDza0/SVJDwC+DGwp6V5gW8rvzKmqJYlIWgTsOM6m422fW77neGAt8I1y20pgN9u3S3oi8F1Je9u+u/0gthcACwC22GF2wz/SiBgqzb87C9tflPRl4JnAzsDVtv+nyrFqSSK2D5hou6SjgOcD+9tFXre9Glhdvr5C0g3Ao4HFMxxuRETXNCA/W22vk3QJ8EhgW0mPouh3ntIVNK45S9LBFB3pz7J9X0v59sAd5YXvCcwBbqwpzOiDpjZjtUuzVvyFGYSOdQAk7QecAvwGuBPYGthd0lts/7Db4zQuiQCfATYDLpIEf72V95nAByWtBdYBx9iu1BEUETEzNAgd62M+BOzb+j0qaRvge8DgJhHbj+pQfjZwdp/DiRoMSg2kXcaUBDAwNRGKUZHtfcr3MMXRko1LIhERA2193QF07XPAzyT9CPgDsA3wdOBTUzlIkkg0wqDWPjpJP8mIGoxxIgDY/ndJ3wP+nuIW39uB99u+ayrHSRKJiOihftydJekE4HXAbWXRu6tMoFgmjAsnfeMEkkSiVsNWA2mXGskI6l+fyCm2/6XXB5V0ru1Dun1/V0lE0sXAya2ZTtIC25VGOEZERL0kvXu8YuBvpnKcbp9suAfwTknvbymbN5UTRUSMArm7BdhO0uKWZao/yo+VdJWkUyVVmd/pLcAyYHnLsoxyUHe3um3OuhPYH/h02RHziqmcJCJiJJipTHuyynbHH+MTTQ8FfJ5inIfLf0+mmGtwKq4EFtle0XbeZ07lIN0mEdleC7xR0tHATyhGN0ZERKse9YlMNj3UGElfonhS4VSPf2CH8tdM5TjdJpEvtJzgq5J+CbxpKieKaDXsHert0sE+Ovp0d9ZOtleWqy8Grp7Cvrt12LTa9i1TjaWrJGL7i23rVzD1qlNExPDrz91ZJ0maW57tJuD1U9j3JjpEKWk1cAbw5vFmSB9Ptx3rERHRDXe5TOcU9pG2H2/7CbZf2FIr6cZrgR8BB1LMhH4gcAnwRopazeOArm8dzjiRiIgeabnzqsneATzd9qpy/QZJVwM/sv1YSddTJJmuJIlERPRS8x9KtSPwp7ay+8pybN8o6aHdHixJJPpq1DrU22Wm3+E3ADWRHwNflfR24GZgN+BjZTmSHg/8vtuDpU8kIqKX+tAnMk2vpZix9zfAGuAGYPuyHIpR66/r9mCpiURE9MoA9ImUt/HuL2kXYBdgue3lLduvmsrxkkQiInqp4UmkxVqKp5+snc5B0pwVEdFDWt/dUlt80taSvg+sBC4DVkj6Xvlo3ClLEomIGC2nlP8+FtiUYtZeA5+scrDGJRFJJ0haLunKcnley7bjJC2VdJ2k59QZZ0TEuJrfsX4Q8HLbv7a9zvavgaOASt+pTe0T2ehhK5L2Ag4H9gZ2BhZJerTtdXUEGBGxkQHoWC+1R1m5ga1xNZEJHAKcYXu17d8AS4F9ao4pImJDza+JLAK+LmlPSQ+QtCfwVeCiKgdrahIZ72Eru1AMjBmzrCzbiKT5Yw96Wfune2c61oiIv2p+EnkzsBnFD/E15b8PAt5a5WC1NGdVfNjKeHMJjPtR2F4ALADYYofZg1G5jIiBJ+q986obtu8ADpa0EzAbuHmKEzhuoJYkUvFhK8soLnjMrsCKjXaKiKhLQ/tEOjxPvXU7ALY/MtVjN65jfYKHrSwETpf0SYqO9TkU9zhHRDRHA5MIxXTvkzEw+EmEDg9bsb1E0pnAryhGWL4pd2ZFROM0MInY3nemjt24JGL7yAm2nQic2MdwIiKmpInNWTOpcUkkImKgJYlEREQlbv7dWb2WJBIR0UupiUTMnLGn+Y3qEw7zNMPhlz6RiIioLkkkIiIqqX9Kk75LEomI6BGR5qyIiJiGJJGIiKguSSRi5o3aXVq5K2uEjFgSaerzRCIiBk85i283y3RJ+j/lo8KXSDpp+kesJjWRiIhe6kNNRNK+FE97fYLt1ZJ2mPmzji9JJCKih/o07ckbgI/ZXg1g+9a+nHUcSSJRq2HvG0lfyOiZQlPVdpIWt6wvKJ/K2o1HA8+QdCJwP/B225d3H2XvJIlERPTK1AYbrrI9r9PGSR4jvgmwNfBk4EnAmZL2tN33bv0kkYiIXurR1/hEjxGX9AbgO2XSuEzSemA74LbenL17SSLRCK3NPsPQtJVmrNHUxxHr3wX2Ay6R9GjggcCqvpy5TZJIREQPaX1fssipwKmSrgb+DBxVR1MWJIlEAw1qZ3tqH9GvCRht/xl4xcyfaXKNSyKSvgU8plx9GHCn7bmSdgeuAa4rt11q+5j+RxgR0VnmzqqZ7X8Yey3pZOCuls032J7b/6iiDoNSI0kNJDaQJNIMkgS8lKLzKCJiIKQm0hzPAG6xfX1L2R6SfgHcDbzH9o/H21HSfGA+wKZbbT3jgcbMamqNJDWQGFeSyMybaBCN7XPL10cA32zZthLYzfbtkp4IfFfS3rbvbj9IOepzAcAWO8wesY90eDUlmSR5REfu27QnjVFLEploEA2ApE2AQ4EntuyzGhibJ+YKSTdQDP1fPO5BYmjVNaYkySMmM4pPNmzqVPAHANfaXjZWIGl7SbPK13sCc4Aba4ovImJ8dnfLkGhqn8jhbNiUBfBM4IOS1gLrgGNs39H3yKJRZrqJK7WPmKpRq4k0MonYPnqcsrOBs/sfTQyC9i/7qkklSSOmpU+DDZukkUkkYrqSDKIu6ViPiIjKkkQiIqIaM1Sd5t1IEomI6KF0rEdERHVJIhERUcUoDjZMEomI6BW7Xw+laowkkYiIXhqtHJIkEhHRS2nOioiIagykOSsiIiobrRySJBIR0UtpzoqIiMpyd1ZERFSTWXwjIqKqYrDhaGWRJJGIiF7KLL4REVFVaiIREVHNCPaJPKCuE0s6TNISSeslzWvbdpykpZKuk/SclvKDy7Klkt7V/6gjIiZSzJ3VzTIdkr4l6cpyuUnSlT26gCmrsyZyNXAo8MXWQkl7AYcDewM7A4skPbrc/FngQGAZcLmkhbZ/1b+QIyIm0YfmLNv/MPZa0snAXTN+0g5qSyK2rwGQ1L7pEOAM26uB30haCuxTbltq+8ZyvzPK9yaJREQzuL+Px1XxBfpSYL/+nXVDtTVnTWAX4OaW9WVlWafyiIjmsLtbeuMZwC22r+/VAadqRmsikhYBO46z6Xjb53babZwyM37CG/eTkDQfmA+w6VZbdxFpRESPdJ8ftpO0uGV9ge0FYytdfn8eAXyzSpi9MqNJxPYBFXZbBsxuWd8VWFG+7lTeft4FwAKALXaYPWL3SkREnbS+6/asVbbnddo42fenpE0o+pWf2H10vdfE5qyFwOGSNpO0BzAHuAy4HJgjaQ9JD6TofF9YY5wRERsyxWDDbpbpOwC41vaynhytoto61iW9GPg3YHvgPElX2n6O7SWSzqToMF8LvMn2unKfY4ELgFnAqbaX1BR+RMRGhPs52PBwam7KgnrvzjoHOKfDthOBE8cpPx84f4ZDi4iork9JxPbRfTnRJDJiPSKilzLtSUREVDLWJzJCkkQiInpoCndnDYUkkYiInunpQMKBkCQSEdErJkkkIiKmYbRas5JEIiJ6KQ+lioiI6pJEIiKiEhvWjVZ7VpJIREQvpSYSERGVJYlEREQlBqb5/PRBkyQSEdEzBqdPJCIiqjDpWI+IiGlIn0hERFSWJBIREdVkAsaIiKjKQKaCj4iIykasJvKAOk4q6TBJSyStlzSvpfxASVdI+mX5734t2y6RdJ2kK8tlhzpij4jorJz2pJtlSNRVE7kaOBT4Ylv5KuAFtldIehxwAbBLy/aX217cpxgjIqbG4IwTmXm2rwGQ1F7+i5bVJcDmkjazvbqP4UVEVDdiI9Zrac7q0kuAX7QlkK+UTVnvVXsGiohoAru7ZUjMWE1E0iJgx3E2HW/73En23Rv4OHBQS/HLbS+X9GDgbOBI4Gsd9p8PzAfYdKutK0QfEVGBnbuzesX2AVX2k7QrcA7wSts3tBxvefnvPZJOB/ahQxKxvQBYALDFDrOHJ+VHRPMNUS2jG426xVfSw4DzgONs/7SlfBPgYbZXSdoUeD6wqKYwIyI6MF63ru4g+qquW3xfLGkZ8BTgPEkXlJuOBR4FvLftVt7NgAskXQVcCSwHvlRH7BERHY1NBd/NMiTqujvrHIomq/byDwMf7rDbE2c0qIiIXhixW3ybfHdWRMRAMeD17mqZDklzJV1attYslrRPb65g6pJEIiJ6xeVDqbpZpuck4AO25wLvK9dr0aiO9YiIQdenjnUDDylfPxRY0Y+Tjkce8tvRJN0G/HYGDr0dxTQtwyrXN7iG+dpgZq/vEba3r7qzpB9QxNeNzYH7W9YXlMMTujnP31BMCyWKFqWn2p6J77nJYxn2JDJTJC22PW/ydw6mXN/gGuZrg+G/vjETDdgG9gf+n+2zJb0UmF91bN50pTkrIqKBJkoKkr4G/FO5ehbw5b4ENY50rEdEDJ4VwLPK1/sB19cVSGoi1XXVdjnAcn2Da5ivDYb/+rrxOuBfy9k87qecK7AO6ROJiIjK0pwVERGVJYlERERlSSJd6PRM+HLbcZKWls9/f05L+cFl2VJJ7+p/1FMn6QRJy1smv3xey7Zxr3PQDOLnMhlJN0n65dgUGGXZNpIuknR9+e/APFhH0qmSbpV0dUvZuNejwqfLz/MqSX9XX+SjKUmkO2PPhP9Ra6GkvYDDgb2Bg4HPSZolaRbwWeC5wF7AEeV7B8EptueWy/nQ+TrrDLKKAf9cJrNv+ZmN/ch5F3Cx7TnAxeX6oPgqxX9nrTpdz3OBOeUyH/h8n2KMUpJIF2xfY/u6cTYdApxhe7Xt3wBLKR6WtQ+w1PaNtv8MnFG+d1B1us5BM2yfy0QOAU4rX58GvKjGWKbE9o+AO9qKO13PIcDXXLgUeJiknfoTaUCSyHTtAtzcsr6sLOtUPgiOLZsFTm1pAhnk62k1LNfRzsCFkq4oHw0N8HDbKwHKf3eoLbre6HQ9w/qZDoyMEylVfCa8xikz4yfnRtxLPclUCp8HPkQR64eAk4FX0/k6B82wXEe7p9leUT7A7SJJ19YdUB8N62c6MJJEShXnnVkGzG5Z35W/zqbZqbxW3V6npC8B3y9XJ7rOQTIs17EB2yvKf2+VdA5Fs90tknayvbJs3rm11iCnr9P1DOVnOkjSnDU9C4HDJW0maQ+Kzr3LgMuBOZL2kPRAik7phTXG2ZW2tuQXU9xQAJ2vc9AM5OcyEUlbSnrw2GvgIIrPbSFwVPm2o4BOtelB0el6FgKvLO/SejJw11izV/RHaiJdkPRi4N+A7SmeCX+l7efYXiLpTOBXwFrgTbbXlfscSzFV8yzgVNtLagp/Kk6SNJeiOeAm4PUAE13nILG9dkA/l4k8HDhHEhT/P59u+weSLgfOlPQa4HfAYTXGOCWSvgk8G9hO0jLg/cDHGP96zgeeR3Gzx33Aq/oe8IjLtCcREVFZmrMiIqKyJJGIiKgsSSQiIipLEomIiMqSRCIiorIkkYiIqCxJJCIiKksSiZEh6Unl5JKblyO9l0h6XN1xRQyyDDaMkSLpw8DmwIOAZbY/WnNIEQMtSSRGSjln1uXA/cBTB3H6logmSXNWjJptgK2AB1PUSCJiGlITiZEiaSHFEw33AHayfWzNIUUMtMziGyND0iuBtbZPL5+3/l+S9rP9w7pjixhUqYlERERl6ROJiIjKkkQiIqKyJJGIiKgsSSQiIipLEomIiMqSRCIiorIkkYiIqOz/A6NcgxQIZr2iAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd150220b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot a cross section of the conductivity model\n",
    "fig, ax = plt.subplots(1,1)\n",
    "cb = plt.colorbar(mesh.plotImage(np.log10(sigma), ax=ax, mirror=True)[0])\n",
    "\n",
    "# plot formatting and titles\n",
    "cb.set_label('$\\log_{10}\\sigma$', fontsize=13)\n",
    "ax.axis('equal')\n",
    "ax.set_xlim([-120., 120.])\n",
    "ax.set_ylim([-100., 30.])\n",
    "ax.set_title('Conductivity Model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the Survey\n",
    "\n",
    "Here, we define sources and receivers. For this example, the receivers sample $db/dt$, the time-derivative of the magnetic flux. The source is a vertical magnetic dipole with unit moment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the receivers, we will sample the real secondary magnetic flux density as well as the \n",
    "# imaginary magnetic flux density \n",
    "\n",
    "dbdt_z = TDEM.Rx.Point_dbdt(locs=rx_loc, times=times, orientation='z')  # vertical db_dt\n",
    "rxList = [dbdt_z]  # list of receivers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There is 1 source. Each source has 1 receivers sampling the resulting db/dt-field at 30 times\n"
     ]
    }
   ],
   "source": [
    "# Define the list of sources - one source for each frequency. The source is a point dipole oriented\n",
    "# in the z-direction\n",
    "\n",
    "srcList = [\n",
    "    TDEM.Src.MagDipole(\n",
    "        rxList, loc=src_loc, orientation='z', waveform=TDEM.Src.StepOffWaveform()\n",
    "    )\n",
    "]\n",
    "\n",
    "print(\n",
    "    'There is {nsrc} source. Each source has {nrx} receivers '\n",
    "    'sampling the resulting db/dt-field at {ntimes} times'.format(\n",
    "        nsrc = len(srcList), \n",
    "        nrx = len(rxList),\n",
    "        ntimes = len(times)\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up Forward Simulation \n",
    "\n",
    "A forward simulation consists of a paired SimPEG problem and Survey. For this example, we use the E-formulation of Maxwell's equations, solving the second-order system for the electric field, which is defined on the cell edges of the mesh. This is the `prob` variable below. The `survey` takes the source list which is used to construct the RHS for the problem. The source list also contains the receiver information, so the `survey` knows how to sample fields and fluxes that are produced by solving the `prob`.\n",
    "\n",
    "\n",
    "\n",
    "**An Aside: Mappings** \n",
    "\n",
    "The `sigmaMap` defines a [mapping](http://dev-docs.simpeg.xyz/content/api_core/api_Maps.html) which translates an input vector to a physical property - specifically, $\\sigma$. This comes in handy when you want to invert for $\\log(\\sigma)$, in that case, the model is defined in terms of log-conductivity, and you would need an `ExpMap` to translate $\\log(\\sigma) \\to \\sigma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum time is 1.3e-04. \n",
      " There are 80 timesteps, 6 of them distinct (which is the same as the number of matrices that need to be factored)\n"
     ]
    }
   ],
   "source": [
    "# solve the problem at these times\n",
    "timeSteps = [(1e-8, 20), (1e-7, 20), (2e-7, 10), (1e-6, 10), (2e-6, 10), (1e-5, 10)] \n",
    "\n",
    "# define a problem - the statement of which discrete pde system we want to solve\n",
    "prob = TDEM.Problem3D_e(mesh, timeSteps = timeSteps, sigmaMap=Maps.IdentityMap(mesh))  \n",
    "prob.solver = Solver\n",
    "\n",
    "survey = TDEM.Survey(srcList)\n",
    "\n",
    "# tell the problem and survey about each other - so the RHS can be constructed for \n",
    "# the problem and the resulting fields and fluxes can be sampled by the receiver. \n",
    "prob.pair(survey)  \n",
    "\n",
    "print(\n",
    "    'The maximum time is {:1.1e}. \\n There are {} timesteps, '\n",
    "    '{} of them distinct (which is the same as the number of matrices that need to be factored)'.format(\n",
    "        prob.times[-1], prob.nT, (len(timeSteps))\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve the forward simulation \n",
    "\n",
    "Here, we solve the problem for the fields everywhere on the mesh, for with and without the sphere. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solving with sphere ... \n",
      "... done \n",
      "CPU times: user 6.31 s, sys: 415 ms, total: 6.72 s\n",
      "Wall time: 4.93 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "print('solving with sphere ... ')\n",
    "fields = prob.fields(sigma)\n",
    "print('... done ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solving without sphere ... \n",
      "... done \n",
      "CPU times: user 6.51 s, sys: 487 ms, total: 7 s\n",
      "Wall time: 6.62 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "print('solving without sphere ... ')\n",
    "fields_background = prob.fields(sigma_background)\n",
    "print('... done ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the fields\n",
    "\n",
    "Lets look at the physics!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# log-scale the colorbar\n",
    "from matplotlib.colors import LogNorm \n",
    "import ipywidgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bb711db917ac4490b952a8724589df31",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='time_ind', max=79, min=1), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_dbdtSphere>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def plot_dbdtSphere(\n",
    "    time_ind=0 # which frequency would you like to look at?\n",
    "#     ax=ax\n",
    "):\n",
    "    fig, ax = plt.subplots(1,1, figsize=(6,5))\n",
    "    \n",
    "    # Plot the magnetic flux\n",
    "    field_to_plot = fields[srcList, 'dbdt', time_ind]\n",
    "    max_field = np.abs(field_to_plot).max() #use to set colorbar limits\n",
    "    cb_range = 5e2  # dynamic range of colorbar\n",
    "    \n",
    "    cb = plt.colorbar(mesh.plotImage(\n",
    "        field_to_plot, \n",
    "        vType='F', view='vec', \n",
    "        range_x=[-120., 120.], range_y=[-180., 80.],\n",
    "        pcolorOpts={\n",
    "                'norm': LogNorm(), \n",
    "                'cmap': plt.get_cmap('viridis')\n",
    "            },\n",
    "        streamOpts={'color': 'w'}, mirror=True, ax=ax, \n",
    "        clim=[max_field/cb_range, max_field]\n",
    "    )[0], ax=ax)\n",
    "\n",
    "    ax.set_xlim([-120., 120.])\n",
    "    ax.set_ylim([-180., 70.])\n",
    "    cb.set_label('|db/dt|')\n",
    "\n",
    "    # plot the outline of the sphere\n",
    "    x = np.linspace(-sphere_radius, sphere_radius, 100)\n",
    "    ax.plot(x, np.sqrt(sphere_radius**2 - x**2) + sphere_z, color='k')\n",
    "    ax.plot(x, -np.sqrt(sphere_radius**2 - x**2) + sphere_z, color='k')\n",
    "\n",
    "    # plot the source and receiver locations\n",
    "    ax.plot(src_loc[0],src_loc[2],'co', markersize=6)\n",
    "    # ax.plot(rx_loc[0,0],rx_loc[0,2],'co', markersize=6)\n",
    "\n",
    "    # plot the surface of the earth\n",
    "    ax.plot(np.r_[-200, 200], np.r_[0., 0.], 'w:')\n",
    "\n",
    "    # give it a title\n",
    "    ax.set_title(\n",
    "        'db/dt, {time:10.2e} s'.format( \n",
    "            time=prob.times[time_ind]\n",
    "        )\n",
    "    )\n",
    "    plt.show()\n",
    "    return ax\n",
    "\n",
    "ipywidgets.interact(\n",
    "    plot_dbdtSphere, \n",
    "    time_ind=ipywidgets.IntSlider(min=1, max=len(prob.timeSteps)-1, value=1) \n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot data with and without sphere\n",
    "\n",
    "In what follows, we plot db/dt data with and without the sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpred = survey.dpred(sigma, f=fields)\n",
    "dpred_background = survey.dpred(sigma_background, f=fields_background)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'time (ms)')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEaCAYAAAAyinE1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd8zWf7wPHPlUGIRBA7RsggQ4m9Q0ujtFpVpWiVPjx2q4sOPPoU1aE6VFXVVkqVoqW/1mztvWKv2DFiBknu3x8n8qSpkXGS70lyvV+v83LOd145t3Ouc3/v+3vfYoxBKaWUshcnqwNQSimVs2hiUUopZVeaWJRSStmVJhallFJ2pYlFKaWUXWliUUopZVeaWJRKQUQmich/E583FJG9VsfkSETEiIif1XEox6WJRTkMEWkgIn+JSIyIXBCRP0WkppUxGWNWGWMCrYwhM4hINxGJFJErInJGRBaJiIfVcamcwcXqAJQCEBFPYCHQE5gN5AEaAjetjCu7ExEXY0xcimWNgeFAhDFmi4gUBh53hNhUzqA1FuUoAgCMMTONMfHGmBvGmKXGmO0AIuIkIu+IyFEROSsiU0SkYOK6cBGJSn4wETkiIo8kPh8qIrMT97kiIrtEpEaybauJyObEdbMAt2Tr/nbsxOO+JiLbE2tWs0Qk+fZviMgpETkpIi/d77KRiJQSkQWJtbMDIvKvZMtvJH7hJ48xWkRcE193FZE9InJRRJaISLlk2xoR6S0i+4H9dzl1TWCNMWZL4nt+wRgz2RhzJXH/SSIyTkR+S3xPViQ/fqJHRGR/4vm/FBFJdv40xSYilRLPdUFE9opIu7u9Xyr70MSiHMU+IF5EJotICxEplGJ9l8RHE6ACUAD4Ig3HfwL4HvACFtzZV0TyAD8BU4HCwA/A0w84VjsgAvAFqiTGhYhEAAOARwA/oPEDjjMTiAJKAW2B4SLysDHmJLAmRRzPAXOMMbdF5EngLaANUBRYlXis5J4EagNBdznvOuBREfmPiNQXkbx32aYj8B7gDWwFpqdY3wpbgnoI2/vxaOJ7kKbYRMQd+A2YARQDOgBjRST4LjGp7MIYow99OMQDqAxMwvZlG4ctARRPXPc70CvZtoHAbWyXc8OBqBTHOgI8kvh8KPB/ydYFATcSnzcCTgKSbP1fwH8Tn//t2InH7ZTs9ShgXOLzicCIZOv8AAP43eVvLQPEAx7Jlo0AJiU+fwn4I/G5AMeBRomvfwG6JdvPCbgOlEt8bYCmD3ivWwA/A5eAq8AngHPiuknA98m2LZAYa5lkx2+QbP1sYGB6YgOeBValiO1rYIjV/x/1kf6H1liUwzDG7DHGdDHG+AAh2H7Jf5q4uhRwNNnmR7ElleKpPPzpZM+vA24i4pJ43BMm8Rst2bHTcqwCyWI8nmxd8ucplQIumMTLT8nOWzrx+RygroiUwpb8DLZf/wDlgDEicklELgEXsCWf0smOdb9zY4z5xRjzOLZaWmtsta6X7ra/MeZq4jlKJVt/r/cgrbGVA2rf2T5xn45AifvFrxybNt4rh2SMiRSRSUCPxEUnsX0J3VEWW63mDLYvvPx3VoiIM7bLMKlxCigtIpIsuZQFDqYj7FOAT7LXZe6z7UmgsIh4JEsuZYETAMaYSyKyFNtlpsrAzGTxHQfeN8akvDyVXKqGLTfGJAC/i8gf2JL5P2IXkQLYEtDJVBwyrbEdB1YYY5qlJl6VPWiNRTmExAbcV0XEJ/F1GWzX29cmbjITeEVEfBO/6IYDs4ytV9E+bDWQlomN2+8Ad2s3uJs12BJUPxFxEZE2QK10/hmzgRdFpLKI5AcG32tDY8xxbJfcRoiIm4hUAbrx97aMGcDz2NpaZiRbPg4YdKcdQkQKisgzqQ1SRFqLSHsRKSQ2tbC1B61NttljYuv+nQdbW8u6xJgfJK2xLQQCRKSziLgmPmqKSOXU/j3K8WhiUY7iCrYG3XUicg3bl9xO4NXE9ROxNbCvBA4DsUBfAGNMDNALmIDtF/81bO00D2SMuYWtobkLcBHbNf8f0/MHGGN+AT4DlgEHsCUtuHeX6Q5AeWw1gXnY2hV+S7Z+AeAPnDHGbEt2nnnAB8D3InIZ2/vUIg2hXgT+ha1X1mVgGvBhilrGDGAItktZ1bFdnnqgtMaWWFtrDrTH9j6cTtw/tT8MlAOSv19aVkrZS+Kv7p1AXpON7tdIvAQZZYx5x+pYVPakNRal7EhEnhKRPIndpT8Afs5OSUUpe9DEopR99QDOYWv8j8c2koBSuYpeClNKKWVXWmNRSillV5pYlFJK2VWuvEHS29vblC9fPs37JSQkAODkpPnYkWi5OB4tE8eU0XLZtGlTtDHmgTcf58rEUr58eTZu3Jjm/a5csd0g7eGh01Y4Ei0Xx6Nl4pgyWi4i8qDhjgC9FKaUUsrONLEopZSyK00sSiml7CpXtrEopbKH27dvExUVRWxsrNWh5Aipbbx3c3PDx8cHV1fXdJ1HE4tSymFFRUXh4eFB+fLlSTb7sUqn+Ph4AJydne+5jTGG8+fPExUVha+vb7rOo5fClFIOKzY2liJFimhSyUIiQpEiRTJUS9TEkgZHL9xg9cELnLtyr1HQlVL2pkkl62X0PdfEkga/7DpLz+93UvP9/6P28P+j66QNfLx0L7/uPM3xC9fRcdeUyn0ee+wxLl26xKVLlxg7dmzS8uXLl9OqVSu7niszjpkZckQbi4g0xDYRkQsQZIyplxnn6VizNNXLFuTQxdvsPnmZnSdjWL73LAmJ+aRgPleCS3kmPgoSXMqTCkUL4Oykv7iUyqkWL14MwJEjRxg7diy9evWyOKJ7i4uLw8Ul87/2LU8sIjIRaAWcNcaEJFseAYwBnIEJxpiR9zqGMWYVsEpEngQ2ZFasHm4u1CznRdOQ/921Gns7nsjTV9h5IoZdJy+z+2QMk9cc5VacrfeFm6sTlUr8PdkElvDAzfXejWdKKccwatQo3Nzc6NevH6+88grbtm3jjz/+4Pfff+e7775j2rRpSSN5DBw4kIMHD1K1alWaNWtGy5YtuXr1Km3btmXnzp1Ur16dadOm/eMy02effca4ceNwcXEhKCiI77//nqFDh3Lw4EFOnDjB8ePHeeONN/jXv/4FcM9jbtq0iQEDBnD16lW8vb2ZNGkSJUuWJDw8nHr16vHnn3/SqlUrOnfuTO/evTl27BgAn376KfXr17fr+2Z5YgEmAV8AU+4sEBFn4EugGbYpZjeIyAJsSWZEiv27GmPOJj5/DnjpQSdMSEhIGtogLe61T0UvZyp6FaZ1cGEAbscncPj8DSJPX2XPmatEnr7Kgq0nmL7OVpDOAr7e+alUvACVShSgcvECBBZ3p2C+9HXty+3SU5Yqc9mrTBISEpJ6Mr23aA97Tl22y3HvqFzSk3dbVr7n+vr16zN69Gh69+7Nhg0buHnzJrGxsaxcuZL69esnxRYfH8/777/Pzp072bRpE2C7bLVlyxa2b99OqVKlaNiwIStXrqRBgwZ/O8fIkSM5cOAAefPm5dKlS8THx5OQkMD27dv5888/uXbtGjVq1CAiIoL4+Pi7HrN27dr06dOHefPmUbRoUWbPns1bb73FhAkTMMZw8eJF/vjjD+Lj43n++efp168fDRo04NixYzz22GPs3LnzH397er8nwQESizFmpYiUT7G4FnDAGHMIQES+B1obY0Zgq938g4iUBWKMMXf9nyci3YHuAGXKlLFP8Pfg6uxEQDF3Aoq58wTFAVsXvqhLsUSeucqe01eJPH2NdUcusXDn2aT9SnjmJbCYOwHF3QksXoDAYu6ULZwPJ228VMoS1atXZ/PmzVy5coW8efMSFhbGxo0bWb16NZ9++ukD969ZsyY+Pj4AVK1alaNHj/4jsYSGhtK5c2dat25N69atk5Y//vjj5MuXj3z58hEeHs6GDRsoWLDgXY/p5eXFrl27iIiIAGyJrkSJEknHeuaZZ5Ke//HHH0RGRia9vnz5MleuXLHruG6WJ5Z7KA0cT/Y6Cqj9gH26Ad/da6UxZjwwHqBGjRomI29ievcN8vQkqGyxvy07d+Umu09dZk/iY/fJy6w+dJH4xIab/HmcCSzhQeWSnlQu6Yl/sQL4ertTzCOv9pZJQQc8dDwZLRMnJ6ekey6GPhHygK3tz9nZmfLlyzNlyhTq169PlSpVWLlyJYcOHSIkJCTpM+js7JwUZ/J/3dzckl67uLiQkJDwj3tIFi9ezMqVK1mwYAHvv/8+u3btwsnJCRFJ2vbO83sd08nJieDgYNasWfOPv0FE8PT0TNonISGBNWvWkC9fvvv+7U5OTukuP0dNLHf7xrxvlytjzJBMiiVTFfXIS2OPojQO+N9I1LG349l/5qot0SQ+ft52khmJl9LAlnDKF3HHt6g7vkXcKe/tjm/io1B+V006StlJo0aN+Oijj5g4cSKhoaEMGDCA6tWr/+Mz5uHhkeZLRwkJCRw/fpwmTZrQoEEDZsyYwdWrVwGYP38+gwYN4tq1ayxfvpyRI0eyb9++ux4nMDCQc+fOsWbNGurWrcvt27fZt28fwcHB/9i2WbNmfPHFF7z++usAbN26lapVq6Yp7gdx1MQSBSS/XuUDnLQolizn5upMqE9BQn0KJi0zxnAyJpZD565yOPpa0mPXiRh+3Xk6qYYDtt5pFYq6E1TSk6BSngSV9KRSCU/y5dEOA0qlVcOGDXn//fepW7cu7u7uuLm50bBhw39sV6RIEerXr09ISAgtWrSgZcuWDzx2fHw8nTp1IiYmBmMMr7zyCl5eXgDUqlWLli1bcuzYMd59911KlSp1z8SSJ08e5syZQ79+/YiJiSEuLo6XX375ronl008/pV+/flSpUoW4uDgaNWrEuHHj0viu3J9DzHmf2May8E6vMBFxAfYBDwMnsPX0es4Ys8se56tRo4bJSfOx3I5P4PiF6xw5f43D0dc5HH01qcZzOTYOACcBX293gkoV/FvCKeqR1+LoM85RyyU3s1eZ7Nmzh8qV7924nlMNHTqUAgUK8Nprr9n1uKkZ0uWOu733IrLJGFPjQftaXmMRkZlAOOAtIlHAEGPMtyLSB1iCrSfYRHsllZzI1dmJCkULUKFogb8tN8Zw4tINdp9MvKR28jJbjl3k523/q/wV9chLSClPqvh48VCZgoSW9soRyUYpZR3LE4sxpsM9li8GFmdxODmKiOBTKD8+hfLTPPh/PURibtxO6iiw6+Rldp6IYcW+/Uk3epYq6EaoT0FbsvHxIrR0QQrm167QSmWVoUOHWh1ChlieWFTWK5jPlToVilCnQpGkZdduxrHr5GW2R11ie1QM26MusWTXmaT15Yvk56EyXlQr40W1soUIKuWJq7OOCKSU+idNLAoA97wu1PItTC3fwknLYq7fZseJGLZFXWJ71CXWHjrP/K22y2h5XZyo4lOQamULEVbWi7CyhSjm6WZV+EopB6KJRd1TwfyuNPD3poG/d9Kyk5dusPnYRbYcu8TmYxeZ9OcRxq+0DV9T2isf1RKTTM3yhQkq5anjpCmVC2liUWlSyisfpbzy0apKKQBuxsWz6+RlNh9NTDZHL7Jw+ykAPPK6UKN8IWpXKEIt38KEli6ol8+UygU0sagMyeviTFjZQoSVLZS07FTMDdYfvsC6wxdYd+g8y/aeA2w3dVYvV4javoWpXaEIVXwKktdF761RjuvIkSO0atXqrmNppdby5cv56KOPWLhwoR0js48uXbrQqlUr2rZta9fjamJRdleyYD5aVy1N66qlAduwNesPX2D94fOsO3yBj5babvLK6+JEWNlCtsttft6ElC6ol86USsYYgzHmgXPUO5rsFa3Klop65KVllZL8p3UIv77ciC3vNuPrztXpWLscl27c5sMle2n95Z9U/+9v9Jq+iRnrjnH8wnWrw1YKsM1h8sILL1ClShXatm3L9evXGTZsGDVr1iQkJITu3bsnTfJ34MABHnnkER566CHCwsI4ePDg3461YcMGqlWrxqFDhzh37hzNmjUjLCyMHj16UK5cOaKjozly5AiVK1emV69ehIWFcfz4cWbOnEloaCghISG8+eabSccrUOB/967NmTOHLl26ALaaSL9+/ahXrx4VKlRgzpw5gC1R9evXj6CgIFq2bMnZs2fJDFpjUVmukHseHg0uwaOJ99acu3KTvw5Gs2p/NKv3R7N4x2kAyhXJT30/bxr6eVOvorfeS6MIDw//x7J27drRq1cvrl+/zmOPPfaP9V26dKFLly5ER0f/45LP8uXLH3jOvXv38u2331K/fn26du3K2LFj6dOnD4MHDwagc+fOLFy4kMcff5yOHTsycOBAnnrqKWJjY5PGAgP466+/6Nu3L/Pnz6ds2bL06dOHpk2bMmjQIH799VfGjx//t3N+9913jB07lpMnT/Lmm2+yadMmChUqRPPmzfnpp5948skn7xv3qVOnWL16NZGRkTzxxBO0bduWn376ib1797Jjxw7OnDlDUFAQXbt2feB7kFaaWJTlinrkTbp0Zozh4LlrrN5/jtUHolmw1Tb4ppNAtbKFeLhyMR6uVJyA4gV0oE2VJcqUKZM0EVanTp347LPP8PX1ZdSoUVy/fp0LFy4QHBxMeHg4J06c4KmnngLAze1/3e/37NlD9+7dWbp0KaVK2Tq+rF69mnnz5gEQERFBoUL/a6csV64cderUAWy1nPDwcIoWtQ1U27FjR1auXPnAxPLkk0/i5OREUFAQZ87Y7klbtWoV7du3x9nZmVKlStG0aVN7vEX/oIlFORQRwa9YAfyKFaBLfV9uxyew7fglVu6PZlnkWUb9updRv+7Fp1A+HqlcnKaVihFcNA95XPSqbm5wvxpG/vz577ve29s7VTWUlFL+gBERevXqxcaNGylTpgxDhw4lNjaW+427WLJkSWJjY9myZUtSYrnf9u7u7knP77dd8thiY2P/ti5v3v8NzZT8GFnxg0w/jcqhuTo7UaN8YQY0C+Dnvg1YO+hhRrQJpVIJD77fcIznJ66n4eg1vDJnNz9sPE701ZtWh6xymGPHjiXNczJz5sykibq8vb25evVqUvuFp6cnPj4+/PTTTwDcvHmT69dtbYVeXl4sWrSIt956Kym5NWjQgNmzZwOwdOlSLl68eNfz165dmxUrVhAdHU18fDwzZ86kcePGABQvXpw9e/aQkJCQVPu5n4YNGzJr1izi4+M5deoUy5YtS+e7cn9aY1HZSomCbnSoVZYOtcoSezuevw5G88u2E6w4cJ7/2xuNCFQr48VjoSVpEVqS0l73n8xIqQepXLkykydPpkePHvj7+9OzZ08uXrxIaGgo5cuXp2bNmknbTp06lR49ejB48GBcXV354YcfktYVL16cn3/+mRYtWjBx4kSGDBlChw4dmDVrFo0bN6ZkyZJ4eHgkzcdyR8mSJRkxYgRNmjTBGMNjjz2WNNPkyJEjadWqFWXKlCEkJOQf+6b05JNPsmzZMkJDQwkICEhKUPbmEMPmZ7WcNmx+bnflyhWMMRy7Yvgj8ixLdp1m10nbDNVVy3jRMrQkLUJL4FMov8WR5h46bP6D3bx5E2dnZ1xcXFizZg09e/Zk69atmXrOXDNsvlL2ICKElPYkpHRB+j3sz5HoayzeeYrFO07x/uI9vL94Dw+V8aJlaAlahJSkTGFNMspax44do127diQkJJAnTx6++eYbq0OyG00sKkcq7+1Or3A/eoX7cfT8NRbtsCWZ4YsjGb44kod8CvJYaEkef6gUpfRymbKAv78/W7ZssTqMTKGJReV45Yr8L8kcO389KcmM+CWSkb9G0sDPm7bVfXg0uARurjrEjFIZpYlF5Spli+SnZ3hFeoZX5Oj5a8zdfIK5m6Lo//1WPNxcePyhUjxT3YeqZbz0PhkHYYzRsshiGW1718Sicq1yRdwZ0CyAlx/2Z+2h8/ywKYofN0cxY90x/IoVoG11H9pUK63zzFjIzc2N8+fPU6RIEU0uWcQYw/nz5/92g2daaa+wNLhy5Qpz584lNDSU6tWrZ0JkKj3s2VvvcuxtFm8/xQ+both09CJOAo0DivJszTI8Urk4Ljrsf6rYq0xu375NVFTUP27+U+mTkGCbO+lBg1q6ubnh4+ODq+vfh1FKba8wTSxpcPHiRWrWrMnhw4fp168fw4YN067HDiCzuoEfOneVOZui+HHzCU5fjqW0Vz5eqFeOZ2uWpWA+HbfsfrRrvmPKaLmkNrHoz680cHFxYdmyZfTo0YMxY8YQFBTE/PnzrQ5LZZIKRQvwRkQl/hzYlPGdq+NTKB/DF0dSd8TvDJ6/k0Pn7n8zmlK5ldZY0iB5tl+zZg09evRgx44d7Nixg5CQEHuHqVIpK38d7zwRw3d/HuHnbSe5FZ/Aw5WK0bWBL/UqahtAclpjcUxZVWPRxJIGKQvl9u3b/Pbbb0lDda9YsYIGDRqk6q5WZT9WfImdvRLL9LXHmL7uKNFXbxFY3IOuDcrTumpp7bKMJhZHpZfCsgFXV9ekpLJ3716aNm1K7dq12bx5s8WRqcxWzMONV5oFsPrNpnzYtgpOTsKbc3dQb+QffPb7fi7H3rY6RKUso4nFTgICApg5cyYnTpygZs2aDBgw4IEDwqnsz83VmWdqlGFxvwbM/FcdqpXx4pPf9tFg5B98+n/7iLmhCUblPnopLA1SU428dOkSgwYNYty4cfj7+7Nr165/dNlT9uVol112nohhzO/7+W33GTzcXHixvi/d6vvmqhkwHa1MlI22sWSirBjdeM2aNezZs4euXbtijOHs2bMUL148zedUD+aoX2K7Tsbw2e/7WbLrDB55XehSvzzdGvjilT+P1aFlOkctk9xO21iyubp16ybNJT1//nwqVKjAqFGjuH1bL43kFsGlCvJ15xr80r8hDQO8+fyPAzT4YBkfLonk4rVbVoenVKbRxJIFwsLCaN68OW+++SZhYWGsXr3a6pBUFqpc0pOxHauz5OVGNA4sytjlB2nwwR98uCSSazfjrA5PKbvTxJIFypYty7x585g/fz4xMTE0bNiQN9980+qwVBYLLOHBl8+FseTlRjStXJwvlx3k4Y9XsHD7yQwP+qeUI9HEkoWeeOIJdu/ezeuvv85DDz0E2GZ00y+V3CWguAefd6jG3J71KOyehz4zttDp23UcOKu9CFXOoIklixUoUIBRo0bx3HPPAfD5558THh6e6VOSKsdTvVwhfu7bgGGtg9kRFUOLMSsZ+YteHlPZnyYWixUpUoTdu3cTFhbGSy+9xOnTp60OSWUhZyfh+brl+eO1cJ6sWppxKw7yyCcrWLT9lNZkVbalicVinTt3Zv/+/QwYMIApU6bg7+/PlClTrA5LZTHvAnn58JmHmNuzHoXy56H3jM10/na9Xh5T2ZImFgfg5eXFRx99xO7du3nkkUcoV64cALGxsfqrNZdJfnlsW9SlpMtj12/p5TGVfWT7xCIiQSIyW0S+EpG2VseTEX5+fsybN4/GjRsD8Oabb1K/fn3WrVtncWQqK925PLbstXBaJ14eaz56Jav3R1sdmlKpYmliEZGJInJWRHamWB4hIntF5ICIDHzAYVoAnxtjegLPZ1qwFggLC+Pw4cPUqVOHTp06cfz4catDUlnIu0BePnrmIWb3qIursxOdvl3HG3O26fhjyuFZOqSLiDQCrgJTjDEhicucgX1AMyAK2AB0AJyBESkO0TXx3yHAdaCeMab+g84bFhZmVqxYkeZ4rRim4sqVK4wePZrPP/8cJycnJkyYwOOPP55l588OcsPwIbG34xm3+hiT1hynsHse3o7w4+FAb6vDuqfcUCbZUUbLxdPT0/GHdDHGrAQupFhcCzhgjDlkjLkFfA+0NsbsMMa0SvE4m/joDQwE7nmtQES6i8hGEdkYHZ19Lil4eHgwePBgNm3axNNPP42vry9gm/ulb9++/P777zpMTC7g5urMy018mf5iNQrnd+XlObt59cfdRF/VoWGU43GxOoC7KA0kv+YTBdS+18YiUh54C3AHPrzXdsaY8cB4sA1CmZFfUlb8CgsODmbq1KlJr0+ePMncuXOZPHkyhQsX5qmnnqJdu3Y88sgjODll+6azdMkNv47reHiwsH8Jxq88xJj/28/6o5sY8ngQT1Yt7ZAzWOaGMsmOMrtcHPEb6G6fjnterzPGHDHGdDfGdDTG5JpBuLp37865c+f46aefiIiIYNasWbz44otJ6w8ePKg9ynIoV2cnejfxY3H/BlTwdueVWdt4cdIGTly6YXVoSgGOmViigDLJXvsAJy2KxaG5ubnRunVrpk+fzrlz51iyZAlOTk4kJCTQsGFD6tSpw6JFizTB5FB+xTz44d/1GPp4EOsPX6D5JyuYtvaolreynCMmlg2Av4j4ikgeoD2wwOKYHJ6bmxshISEAGGMYNmwYZ8+epVWrVtSqVYuff/5Zv3ByIGcnoUt9X5a83IiwcoV456eddJ20gXNXblodmsrFrO5uPBNYAwSKSJSIdDPGxAF9gCXAHmC2MWaXlXFmN87Ozrz00kvs27ePCRMmcP78eZ544glWrlxpdWgqk5QpnJ8pXWvxnyeC+fPgeVqMWcmyyLNWh6VyKZ1BMg2yaxfK27dvs3DhQp588klEhLFjx1KyZElat26dIxr6s2u5ZJZ9Z67Qb+YWIk9f4YW65Rj0WGXcXJ2zNAYtE8ekM0gqu3F1deWpp55CREhISGDChAm0adOGatWqMWfOHBISEqwOUdlRQHEPfupdn24NfJm85iiPf76a3ScvWx2WykU0seQyTk5OrF+/nqlTpxIbG8szzzxD9erVWb9+vdWhKTtyc3Xm3VZBTOlai0s3bvPkl38yYdUhEhJy3xUKlfU0seRCLi4udOrUid27dzN16lSuXLmCm5sbAHFxOthhTtIooGjSlMj/XbSHF75bz5nLsVaHpXI4TSy5mLOzM506dWLv3r1UqVIFgBdeeIG2bduyd+9ei6NT9lLYPQ/jO1dn+FOhbDhygYhPV7Jkl877ozKPJhaFs7OtYdcYQ6VKlViyZAnBwcF0796dEydOWBydsgcR4bnaZVnYtyGlC+Wjx9RNjPo1Urugq0yhiUUlERHeffddDh48SK9evZg0aRKFFVW2AAAgAElEQVR+fn4sWKC3EeUUfsUK8GPP+nSoVYaxyw8yeP4ubXdRdqeJRf1DsWLF+Oyzz4iMjKR9+/bUrm0bqu3YsWPExur1+ewuj4sTw58KpUfjCkxde5TXfthGXLz2DFT2o4lF3VOFChX47rvvKF68OMYY2rdvT6VKlZgxY4Z2Uc7mRISBEZV4/dFAftxygl7TN3MzLt7qsFQOoYlFpYqI8N///pfChQvTsWNHatWqxfLly60OS2WAiNC7iR9DHw9i6e4zvDR5o06BrOxCE4tKtaZNm7Jx40amTJnC2bNnadKkCT/88IPVYakM6lLfl4+eeYg/D0TT+dv1OkOlyjBNLCpNnJyc6Ny5M3v37uXTTz9Nms1y3bp1nDlzxuLoVHq1re7Dl8+FsT3qEh3GryX6qg5iqdJPE4tKl3z58tG/f3/c3NxISEigU6dO+Pn58d5773Ht2jWrw1Pp0CK0JBNeqMmh6Ku0+3oNJ3V+F5VOmlhUhjk5ObFo0SKaN2/O4MGDCQgIYMqUKdrAnw01DijK1G61OXf5Js+MW8ORaP2RoNJOE4uyi4CAAObOncuqVavw8fHhhRde4Ndff7U6LJUONcsXZmb3Oty4Hc8zX68h8rQOYKnSRhOLsqsGDRqwZs0a5s+fT4sWLQD47bffOH/+vMWRqbQIKV2Q2T3q4CTQYfxaDmvNRaWBJhZld05OTjzxxBOICNeuXePZZ5/F39+fsWPHEh+v90pkF37FPJjVvS4iwovfrefCtVtWh6SyCU0sKlO5u7uzatUqqlatSu/evalevTqrVq2yOiyVSuW93fnm+RqcionlpckbiL2tPwzUg2liUZkuODiY33//ndmzZ3PhwgUaN27M/v37rQ5LpVL1coX49NmqbDl+iVdnb9OxxdQDaWJRWUJEeOaZZ9izZw+zZ8/G398fsLW/3Lqll1gcXYvQkrzVojKLdpzig18jrQ5HOThNLCpLubu707ZtWwAOHjxIREQEVatW1eFhsoGXGvryfN1yfL3yEFPXHrU6HOXANLEoy1SsWJH58+dz48YNmjRpQseOHTl16pTVYal7EBEGtwri4UrFGDJ/J39E6kgL6u40sShLtWrVit27d/Puu+8yZ84cqlWrxo0bese3o3JxduKzDtUIKuVJnxlb2HkixuqQlAPSxKIsly9fPoYNG8aOHTsYPXo0+fLlwxjD7t27rQ5N3YV7XhcmvlCTQvnz0HXSBk7o0C8qBU0symEEBATQoUMHABYtWkRwcDAvvfQS0dHRFkemUirm6cZ3L9bkxq14un63gcuxOiKy+h9NLMohhYeH8/rrrzN58mQCAwP55ptvdOwxBxNQ3INxnatz8NxVek7bxK04LR9lo4lFOaQCBQowatQotm7dSkhICN27d6dNmzZWh6VSqO/nzcinq/DngfO8NW8Hxug9LgpcrA5AqfsJDg5m+fLlTJs2jfz58wMQFxdHbGwsBQoUsDg6Bba5XI5fuM6Y3/cTWNyDfzWqYHVIymJaY1EOT0To3LkzTz/9NABffPEFlStXZt68efoL2UG8/Ig/jwYXZ9SSSO0pplKfWESkXGYGolRq1a5dm8KFC9OmTRsef/xxjhw5YnVIuZ6IMLJNFQq75+HlWVu5oWOK5WppqbHMS7lAROrYMRalUqVu3bps2rSJjz/+mOXLl1O7dm2mTZtmdVi5XiH3PHz8TFUOnL3Kx78fsjocZaEHJhYRaSciIwEPEaksIs7JVo/PvNCUujcXFxcGDBjAnj17aNasGb6+vgDac8xiDfy9+VdDX2ZtOsWK/ToHT26VmhrLn8BuoBDwCbBfRDaLyEJA74xSlipTpgzTpk2jfv36AAwYMIAuXbpw9uxZiyPLvV57NJDAYu4MXriPs1dirQ5HWSA1iaUK8H/AE8aYFsaYCsAjwBCgaWYGp1RaGGPw8PBg+vTpBAYG8tVXX+nEYhbI6+LMyCcrce1WPK//sF07WORCqUksTwMLgFki8puIfAg8BtwE9OeIchgiwnvvvcf27dupVq0avXr1ok6dOuzcudPq0HIdv6LuvPpwBVbsO8fkv45YHY7KYg9MLMaYl4wxNYCPgX3AYaAJsB7QsbOVw6lcuTK///47M2bM4Pz587i5uVkdUq7UvnpJmgQWZfgvkew7c8XqcFQWSkuvsBeNMb2NMWONMd2AhsDqTIrrnkSkgoh8KyJz7rdM5W4iQocOHdi/fz9+fn4AdOvWje+++04b+LOIiDCq7UN4urnQb+YWndY4F0lLYrksItXvvDDGbAIC0nIyEZkoImdFZGeK5REisldEDojIwPsdwxhzKDGx3XeZUgDOzrZOjFevXmXv3r107dqVhg0bsm3bNosjyx2KeuTlw7YPEXn6Ch8u2Wt1OCqLpGVIl27AVBHZDWwCQoG0Dmk6CfgCmHJnQWL35S+BZkAUsEFEFgDOwIgU+3c1xmS4u09CQgJXrqS9ap6efVTmS225LFq0iJkzZ/LOO+8QFhZGjx49ePvtt/H09MzkCHOf5GVSo3Q+2lcvxberD1OrjDv1KhS2MLLcLau+w1JzH0tdERFjzD6gHrAYKA7swdaIn2rGmJXAhRSLawEHEmsdt4DvgdbGmB3GmFYpHulOKiLSXUQ2ishGHYY9d3JycqJjx45s2rSJF198kblz52qvsSzy6sO+VPDOz9sL9nHxug6xn9PJg7oCisg4bF/++4BfgV+NMafTfUKR8sBCY0xI4uu2QIQx5qXE152B2saYPvfYvwjwPrYazgRjzIi7LbtfDDVq1DAbN25Mc+x3sr2Hh0ea91WZJ73lcvnyZTw9PYmLi6NXr1707duX0NDQzAgx17lbmew6GcNTX/5F48CijO9cHRGxKrxcK6PfYSKyKbEz132lplfYv40xYcBQbDdJThKRNSIyXEQapbgTP12x3u2094nnfGJMFe8kkLstU+pB7lwC27dvHz/++CPVqlXjlVdeISZGB1HMDMGlCvJGRCC/7T7D9xuOWx2OykSpbrw3xkQaY0YbYyKw3Ri5GngGWJfBGKKAMsle+wAnM3hMpVItKCiIvXv38tJLLzFmzBgCAwOZNm2a3tiXCbrW96WBnzfvLdzNqRgduCOnSu+w+V8Dy4wxfbFdfsqIDYC/iPiKSB6gPbYbMpXKMkWKFGHcuHGsX7+esmXLMmLECOLi4qwOK8dxchKGPxVKXIJh5C+RVoejMkl6E4sT8JWI5ANeSe1OIjITWAMEikiUiHQzxsQBfYAl2DoEzDbG7EpnXEplSI0aNVi7di1Lly7F1dWVy5cv88477+jlMTsqWyQ/PRpVYP7Wk6w/nLIvj8oJ0ptYDmNrc/kKcE/tTsaYDsaYksYYV2OMjzHm28Tli40xAYltJO+nMyal7MLJyYnSpUsDsHTpUoYPH05gYCBTp07Vy2N20ivcj1IF3RiyYBfxCfqe5jTpTSzfGGOOYEsuEXaLRikH07ZtW9atW0e5cuV4/vnnadiwIVu3brU6rGwvXx5n3m4ZxJ5Tl5mx/pjV4Sg7S1diMcYcS/z3iDEm2L4hKeVYatasyZo1a/j222/Zt28fb731ltUh5QiPhZagboUifLx0Lxev3bI6HGVHaZmauJGIrBSRXSIyQ0RqZWZgSjkSJycnunbtyr59+xg/3ja/3ZEjRxg/frzeZJlOIsKQJ4K4EhvHx7/pcC85SVpqLBOB94BwbEOyfCoi7TIjKKUclZeXFz4+PgBMnDiRHj16ULt2bdauXWtxZNlTpRKedK5TjhnrjrHrpHaQyCnSkliijTG/GWPOGWN+BZoDgzMpLqUc3n/+8x+mT5/OqVOnqFu3Ll27duXMmTNWh5XtvPJIAF758/CfBbu1c0QOkZqxwqaIyMvAahEZLCJ3Bq7Uib5UriYiPPfcc0RGRvLGG28wdepUhg8fbnVY2U7B/K68/mgg649cYME2vTc6J0hNjeVbbEOsFAaeBA6IyP8BkdgGpFQqV/Pw8OCDDz5gx44dDB5sq8SvW7eOH374QX+Bp1K7GmUILV2Q4Yv3cO2m3pia3aUmseQHfjDGdE0cM6witpsihwD5MjM4pbKTSpUqUaRIEQDGjh1Lu3btqF+/PmvWrLE4Msfn7CQMfSKYM5dv8uWyA1aHozIo1XPei8hxEfkNGAk8BGwF7jspl1K51cSJE5kwYQKHDx+mXr16PPvssxw6dMjqsBxa9XKFaBNWmgmrDnMk+prV4agM0DnvlcoEzs7OdOvWjf379zNkyBAWLlzIvHnzrA7L4Q2MqEQeFyfeW7jb6lBUBmS7Oe+Vyk4KFCjA0KFD2b9/P3362KYYmjNnDmPGjOHWLb0pMKVinm70e9iP3yPPsiwyw5PFKotk6Zz3SuVWpUqVIm/evIBtiuSXX36Z4OBgfvrpJ23gT6FLPV8qFHVn2MLd3IzTm0+zo7Qklq7AWBH5TkT6iMjXpH3Oe6VyvYkTJ7J48WLy5MnDU089RdOmTdm2bZvVYTmMPC5ODHk8mMPR15i4+ojV4ah0SMtEX/v535z3VUnHnPdKKdv9Ly1atGDbtm2MHTuWnTt3cuCA9oRKrnFAUZoFFefzP/ZzOkZvl8tu0jQIpTEm3hjzAxBmjPnUGHM+k+JSKsdzcXGhZ8+eHDx4kDZt2gDwySef8N5773H9+nWLo7Peuy2DiEswfLRUxxHLbtI7bP7d5qlXSqWDp6cnIraP1Pbt2xk8eHDS9MgJCQkWR2edskXy82K98szdHMXuk5etDkelQbrnY7FrFEopACZNmsTKlSspXrw4nTt3pk6dOmzevNnqsCzTq4kfBfO5MuKXPVaHotIgvfOxjLV3IEopm4YNG7J+/XomT57MqVOnkmotuXF4/oL5XOnb1J9V+6NZue+c1eGoVEpvjUUplYmcnJx4/vnnOXToEDVq1ACge/futGvXjt27c9fNg53rlKNs4fwMX7xHpzHOJjSxKOXAXF1dATDGULZsWX755RdCQkLo3LlzrulJlsfFiTciAok8fYUfN0dZHY5KBU0sSmUDIsKQIUM4fPgwr732GnPnzqVSpUpMmzbN6tCyRMvQklQt48XHS/dx41buuySY3WhiUSob8fb2ZtSoURw6dIi+ffsSHh4OwN69ezlx4oS1wWUiEeHtlpU5fTmWb1frYJ6OThOLUtlQiRIlGD16dNI0yf3796dixYq8/PLLnDyZMyfLqlm+MM2DijNuxSGir960Ohx1H5pYlMoBxo0bR6dOnfjiiy/w9fWlV69eHD2a8wYff7NFJW7cjmfM/+23OhR1H5pYlMoBypcvz4QJE9i/fz9dunRhwoQJzJgxw+qw7K5i0QI8V6ssM9Yf4+C5q1aHo+5BE4tSOYivry9ff/01Bw8eTBqmf9asWXTu3JnIyEiLo7OP/o/4k8/VmQ9+yRl/T06kiUWpHKhMmTJ4eHgAcPr0aX788UeCgoJ49tln2bFjh8XRZYx3gbz8u3EFlu4+w/rDF6wOR92FJhalcrj+/ftz5MgRBg4cyC+//EKVKlV49dVXrQ4rQ7o1qEAJTzeGL96j89k4IE0sSuUCRYsWZfjw4Rw5coQhQ4ZQp04dAG7cuEFUVPa76TBfHmcGNA9g6/FLLNpxyupwVAqaWJTKRQoXLszQoUN55plnABgzZgz+/v4MGjSImJgYi6NLm6fDfKhUwoNRv+7VmSYdjCYWpXKx9u3b07ZtW0aOHEnFihUZM2YMN29mj3tEnJ2EQY9V5tiF60xbe8zqcFQymliUysXKly/P1KlT2bx5M9WqVePll1/mxRdftDqsVGscUJSG/t58/sd+Ym7oTOmOQhOLUopq1arx22+/sWTJEl577TUATp06xbJlyyyO7MEGtahMzI3bjF2WOwblzA40sSilkjRv3pywsDAARo8eTdOmTWnZsqVDD9UfVMqTNtV8+O6vI5y5HGt1OApNLEqpexg2bBgffPABf/75J1WqVKFnz56cOXPG6rDuqv/D/sQnGMatOGh1KIpsmFhEpIKIfCsic5Itqywi40Rkjoj0tDI+pXIKNzc33njjDQ4cOECvXr2YMGECb731ltVh3VXZIvlpU600M9Yd46zWWiyXpYlFRCaKyFkR2ZlieYSI7BWRAyIy8H7HMMYcMsZ0S7FsjzHm30A7oIb9I1cq9/L29uazzz5j586dDBs2DIAdO3Ywffr0pGmTHUGfpn7EJRjGrdBh9a3mksXnmwR8AUy5s0BEnIEvgWZAFLBBRBYAzsCIFPt3NcacvduBReQJYGDi8e8rISGBK1eupDn49OyjMp+WS9YoVaoUYHu/v/jiC8aPH88nn3zC8OHDqVev3t+2taJMCueBViHFmL7uKJ1rFMe7QJ4sj8HRZVW5ZGmNxRizEkg5uE8t4EBiTeQW8D3Q2hizwxjTKsXjrkkl8dgLjDH1gI53Wy8i3UVko4hsjI6OttefpFSuNGrUKL7++mvOnDlDREQEHTt2dIipkrvXL0tcfAIT1xy3OpRcLatrLHdTGkj+vyAKqH2vjUWkCPA+UE1EBhljRohIONAGyAssvtt+xpjxwHiAGjVqmDsD9KVHRvZVmUfLJWt1796dTp06MXr0aEaMGEGVKlX473//+7dtsrpMgj08eLKaDz9sOUnfZpUo5uGWpefPLjK7XBwhschdlt1zVDljzHng3ymWLQeW2zUqpdQD5c+fn7fffptu3brh7u4OwNKlS9m2bRvdunV7wN6Zo09TP+ZtiWL8ikO80yrIkhhyO0foFRYFlEn22gfImXOrKpVDlShRIulX8Ny5c3njjTeoW7cuixYtyvLRh3293XmyammmrTvKuSvZY3ianMYREssGwF9EfEUkD9AeWGBxTEqpdBo3bhyzZs0iISGBVq1aERERwa5du7I0hj5N/bgVl8A3q7SHmBWyurvxTGANECgiUSLSzRgTB/QBlgB7gNnGmKz9X6iUshsRoUWLFqxdu5bRo0ezfv161q5dm6UxVChagNZVSzNlzRGir2qtJatJbpwkp0aNGmbjxo1p3u9OVz1tJHYsWi6OJ3mZnD9/Hi8vL5ydnfn222+5ePEiffv2JW/evJkaw8FzV2n2yQr+1bACgx6rnKnnyi4y+lkRkU3GmAfeK+gIl8KUUjlYkSJFcHZ2BmDFihW8/vrrBAcH8+OPP2Zq+0vFogV44qFSTFlzlPNaa8lSmliUUllmypQpLFmyBDc3N55++mnCw8PZvn17pp2vT1N/YuPiGa9tLVlKE4tSKks1b96crVu38tVXX7Fnzx4y84Zlv2K2WsvUNUe5cO1Wpp1H/Z0mFqVUlnNxceHf//43R44coWnTpgC8++67DBs2jOvXr9v1XH2b+nHjdrz2EMtCmliUUpbJnz8/AMYYDh06xJAhQwgMDGTatGl2G+DSr5gHraqUYvJfR7TWkkU0sSilLCciTJ8+nRUrVlC8eHE6d+5M3bp17db+0i+x1jJBay1ZQhOLUsphNGrUiPXr1zNp0iTOnDmDm5ttrK+M9h7zL+5By9CSTP7rCBe11pLpNLEopRyKk5MTL7zwAgcPHiQgIACA9u3b079//ww19Pd72J/rt+OZsFprLZlNE4tSyiHdufclLi4OLy8vvvjiCypWrMgHH3zAjRs30ny8gOIePBZaksl/HdVaSybTxKKUcmguLi58/fXX7Nixg0aNGjFw4EACAwNZv359mo/Vr6k/127FMXzxnkyIVN2hiUUplS0EBQXx888/s2zZMgIDA/Hz8wMgJiYm1ccILOFB3yZ+/LApirmbojIr1FxPE4tSKlsJDw/nt99+o3DhwiQkJNCkSRMiIiLYsWNHqvbv/0gAdSoU5p2fdnLgrE5rnRk0sSilsq34+Hg6duzIunXrqFq1Kj169ODs2XvOYA6As5PwWftquOd1ptf0zdy4FZ9F0eYemliUUtmWq6srr776KgcPHqRv375MnDgRf39/Nm3adN/9inm6MfrZquw/e5XB83dmUbS5hyYWpVS2V7hwYT799FN27NhBhw4dCA0NBeDEiRP3vAemoX/RpPaWOdreYleaWJRSOUalSpUYN24cefLk4erVq9SsWZOHH374nnfw32lvefennew/o+0t9qKJRSmVI7m5ufH222+zbds2qlWrdtf2l5TtLddvxVkUbc6iiUUplSO5uLjQu3dvDhw4QL9+/ZLaXw4d+vud93faWw6cu8qQ+Toruj1oYlFK5WiFChVi9OjR7Ny5k/79++Pr6wvAgQMHktpftL3FvjSxKKVyhcDAQIYNG4aIEBUVRZUqVYiIiGD37t1A8vtbdrBP21syRBOLUirXKV68OB988AHr16+nSpUq9OvXj5hLF/msfTUK5HWht7a3ZIgmFqVUruPq6krfvn3Zv38/PXr04Msvv6RSpUrkSYhNam8ZrO0t6aaJRSmVa3l7e/Pll1+ydetW3n77bby8vGjoX5R2FYU5m6KYvfG41SFmSy5WB6CUUlYLDQ1Nuqly48aNjPpXS0pXbcwrlzrz14EwXns0EJ9C+S2OMvvQGotSSiUTEhLC8OHDubR/I6cn/JvvRrxBg0FTGflLJJdjb1sdXragiUUppZJxc3Nj0KBB7Nu3j359+3Br/5+cnvoaY/9vF+EfLmfyX0e4HZ9gdZgOTROLUkrdRalSpRg9ejTHjh5lwby5LB7wCAHFCtCnT29q9R7DrztP3XMcstxOE4tSSt1H0aJFad68OSGlCzLy0ZK4Ht/A1q8H0PrRJjTs/RFbjl6wOkSHo4lFKaVSqXz58pw4dpTPv/iCguYaf371BrVrhPHCpws4dv661eE5DE0sSimVBvny5aNP796cOHqIbyZOwrecD2vPQpOPl/OvzxZw+Oxlq0O0nCYWpZRKB1dXV1568QX2blzNqoHN6VC9JJPe6UZAYCAte/+Ho+dyb4LRxKKUUhlUzNON99pUZcL4ryjqXYTFY4fi5x9Aq15DOJYLE4wmFqWUsgMR4YX2bTmxbweTZ/1IsWLFWPTVMOr2G8N/ft7F2SuxVoeYZTSxKKWUHYkIz7d7iqi92/hh4VKea/M4U9YcJaRNX1r9+22iomOsDjHTaWJRSqlMICK0bdmMD5+pyv+90gj36D0s+no4vhUq0n7A+8Rcy7k1mGyXWESkgoh8KyJzki0LF5FVIjJORMItDE8ppf7Bt2gBDm1exdS5CylcrCSzRr9DCd9KDB7/Y468iz9LE4uITBSRsyKyM8XyCBHZKyIHRGTg/Y5hjDlkjOmWcjFwFXADdPo3pZTDERE6tWnJ6f3b+XDcZFydYcLaEzw6eiULtxzPUXfxS1b+MSLSCFsCmGKMCUlc5gzsA5phSwobgA6AMzAixSG6GmPOJu43xxjTNvG5kzEmQUSKA58YYzreL46wsDCzYsWKNMd/5YptVjkPD48076syj5aL49EyebD4+HhWHrzEp8sOs37Se3g4x/HBf//DU+E1Mu2cGS0XT0/PTcaYBwaYpTUWY8xKIOX4B7WAA4k1kVvA90BrY8wOY0yrFI+z9zjunbrkRSDv3bYRke4islFENkZHR9vpL1JKqfRxdnamSUAR5rwURkS9apzdu4kXWj9CjZadWLPzoNXhZYgjzMdSGkg+m04UUPteG4tIEeB9oJqIDDLGjBCRNsCjgBfwxd32M8aMB8YD1KhRw2Tkl5T+CnNMWi6OR8skdeaM/5ijQ16jS79BrPhpOhGNlvDUq6P46u3uFPNws/v5MrtcHCGxyF2W3fP6nDHmPPDvFMt+BH60c1xKKZVlypUuybK5k9i8cyDd+g9i843ChH+4nGeD3OnfqjpeBbLPRGOO0CssCiiT7LUPcNKiWJRSylJhIZXY8vs8lr37FI0DvBnxWg9KlA/gtVHfEJdNepA5QmLZAPiLiK+I5AHaAwssjkkppSxV3tudsR2r88F/h5DH1ZWP3+xOUf+qjPvhV6tDe6Cs7m48E1gDBIpIlIh0M8bEAX2AJcAeYLYxZldWxqWUUo5IROjf5VnOH9tHn8Efcu1cFD3btSC8zygiTzvuGGRZ2t3YUdSoUcNs3LgxzftpF0rHpOXieLRMMsf5SzF0Hzicfd4NuBbvRMNCVxnUriGVypdK1f4ZLRcRcbzuxkoppdKviFdB5o77gFVvPUqXOmWZOfIVgisH8nj3N7l45ZrV4SXRxKKUUtmMV/48DG4dysL5P+ETUIWF34yiZDk/Xv1gHPEO0MCviUUppbKp5g1rcXTbn3w++QfyuOXnk4E9qdP7Y9YcPG9pXJpYlFIqm+vzfFvOH43k9VHjcC0XRodv1vLogNGs3hppSTyaWJRSKgdwdXVl1Os9+OO1cF5uUp7fxw+jUc2HaPxsD46fydoajCYWpZTKQdxcnXn50WDWrVtL5bqPsHL2eCpU9OeFN0YQe+t2lsSgiUUppXKg6sGB7Fq5iJkLf8erhA9TPnyLiP/MYt2Ri5l+bk0sSimVg7Vv2ZQz+7by2YyFeJYJJD4h8+9ddIRBKJVSSmUiJycn+nZoSaeYy7g43W3cXzufL9PPoJRSyiFkRVIBTSxKKaXsTBOLUkopu9LEopRSyq40sSillLIrTSxKKaXsShOLUkopu9LEopRSyq5y5QySIhID7E+xuCAQc5fNUy73BqIzKbQHuVeMmX2c1G7/oO3utz617/+9lllVLlaVSVr2SW+5ZHS5flbSv52jflbKGWOKPnArY0yuewDjU7PsbsuBjY4Ud1YcJ7XbP2i7+61P7ft/n2WWlItVZZIV5ZLR5fpZsX+ZpLVcrPqs5NZLYT+nctn9llvBXrGk9Tip3f5B291vfVrefy2TtO2T3nKx13Ir6GcldefJFLnyUlhGiMhGY0wNq+NQf6fl4ni0TBxTVpRLbq2xZMR4qwNQd6Xl4ni0TBxTppeL1liUUkrZldZYlFJK2ZUmFqWUUnaliUUppZRdaWKxMxFxF5FNItLK6lgUiEhlERknInNEpKfV8SgbEXlSRBLBpAIAAASrSURBVL4Rkfki0tzqeBSISAUR+VZE5mT0WJpYEonIRBE5KyI7UyyPEJG9InJARAam4lBvArMzJ8rcxR5lYozZY4z5N9AO0K6vdmCncvnJGPMvoAvwbCaGmyvYqUwOGWO62SUe7RVmIyKNgKvAFGNMSOIyZ2Af0AyIAjYAHQBnYESKQ3QFqmAbLsENiDbGLMya6HMme5SJMeasiDwBDAS+MMbMyKr4cyp7lUvifh8D040xm7Mo/BzJzmUyxxjTNiPxuGRk55zEGPP/7d1NiFVlHMfx78/epBe0dFFO0ZSmEZKOi4RaaBG00wqKIAjJFpK2qE1J2K4SLCqDCt3MriiEmEl7IXAyegGjlHGYRaJEg9AbvcxChJh/i/PMcLjN3Obc+1zvofv7rObcee5znjt/Dj+ee+7c/xFJ/Q0P3wacjIhTAJLeATZHxIvAv97qknQncBlwC3BW0qGImOrowv/HctQkzTMEDEk6CDhY2pTpWhGwG/jQodK+XNdKLg6W5vqAH0vHE8D6uQZHxLMAkrZQ7FgcKvlVqomkjcD9wCXAoY6urLdVqgvwBHA3sEjSioh4q5OL61FVr5UlwPPAgKSdKYBa4mBpTrM89p/vHUbEYP6lWFKpJhExAox0ajE2o2pd9gJ7O7cco3pNfgO25Tixb943NwFcVzq+FjjTpbVYwTWpJ9elfrpWEwdLc0eBmyTdIOli4CFgqMtr6nWuST25LvXTtZo4WBJJbwNfAaskTUjaGhF/AzuAj4Fx4N2IGOvmOnuJa1JPrkv91K0m/rixmZll5R2LmZll5WAxM7OsHCxmZpaVg8XMzLJysJiZWVYOFjMzy8rBYmZmWTlYzCqQtFjS46XjZTkaI81xrnslPZdhnpck3ZVjTWbz4X+QNKsgfTX5B9M9Lzp8ri+BTRHxa5vzXA/sjwh3arTzwjsWs2p2A8slHZO0R1L/dNc+SVskvS9pWNJpSTskPSXpO0lfS7oqjVsu6aPUwvpzSTc3nkTSSuDcdKhIGpT0pqTDkk5J2pC6Bo5LGkxjLkjjTkgalfQkQET8ACyRdPX5+RNZr/PX5ptV8wywOiLWwswOpmw1MEDRRfQk8HREDEh6BXgEeBXYB2yLiO8lrQfeABrfqroDaGyAdWUatwkYTmMeA45KWkvRGbCv1EFwcem536bxB1p72Wbz52Axy+twREwCk5L+pAgAgFHgVkmXA7cD7xVNFIGiCVmja4BfGh4bjoiQNAr8FBGjAJLGgH7gM+BGSa8DB4FPSs/9GVjW7oszmw8Hi1le50o/T5WOpyiutwXAH9M7nibOAovmmLs878zcEfG7pDXAPcB24EHg0TRmYZrTrON8j8WsmkngilafHBF/AaclPQBF7/cUBo3GgRVV5pa0FFgQEQeAXcC60q9XAidaW7VZNQ4WswpS+9Yv0g3yPS1O8zCwVdJxYAzYPMuYIxS9x2drLzuXPmBE0jFgENgJIOkiipD6psX1mlXijxub1ZSk1yjuq3za5jz3AesiYleelZk15x2LWX29AFyaYZ4LgZczzGM2L96xmJlZVt6xmJlZVg4WMzPLysFiZmZZOVjMzCwrB4uZmWX1D38NV91h+pOAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd1a098fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "fig, ax = plt.subplots(1,1)\n",
    "\n",
    "ax.loglog(1e3*times, -dpred)\n",
    "ax.loglog(1e3*times, -dpred_background, '--k')\n",
    "ax.grid(True, color='k',linestyle=\"-\", linewidth=0.1)\n",
    "ax.legend(['with sphere', 'background'])\n",
    "\n",
    "ax.set_title('Sounding over Sphere')\n",
    "ax.set_ylabel('-$db_z/dt$')\n",
    "ax.set_xlabel('time (ms)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}