1708_Nonlinear_inversion/1_MT1D_NumericalSetup.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Getting Started**\n",
    "\n",
    "This tutorial relies on standard python packages as well as [SimPEG](http://simpeg.xyz). If you do not have SimPEG installed, you can uncomment the next line and install it from [pypi](https://pypi.python.org/pypi/SimPEG)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# !pip install SimPEG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "from SimPEG import Mesh, Utils, Solver  \n",
    "from scipy.constants import mu_0, epsilon_0\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical simulation of the 1D Magnetotelluric (MT) problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Purpose\n",
    "\n",
    "With [SimPEG's](http://simpeg.xyz) mesh class, we discretize Maxwell's equations for a 1D magnetotelluric problem. We then solve for both electric and magnetic fields, and evaluate data at a receing location. There are some milestones to be accomplished:\n",
    "\n",
    "- Introduce differential operators and the terminology used in the SimPEG mesh class\n",
    "\n",
    "- Set up boundary conditions\n",
    "\n",
    "- Set up an linear system $\\mathbf{A}\\mathbf{u} = \\mathbf{rhs}$, compute the fields, $\\mathbf{u}$\n",
    "\n",
    "- Evaluate the data at a receiver location: apparent resistivity and phase\n",
    "\n",
    "- Recognize extensibility of this example to higher dimensions: 2D and 3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Physics: Maxwell's equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The governing equations for electromagnetic problems are Maxwell's equations. Here, we show them in the frequency domain. For more background on Maxwell's equations, we recommend http://em.geosci.xyz and [Ward & Hohmann, 1988](http://library.seg.org/doi/abs/10.1190/1.9781560802631.ch4).\n",
    "\n",
    "$$\\nabla \\times \\mathbf{E} + \\imath\\omega \\mu \\mathbf{H} = 0 $$\n",
    "\n",
    "$$\\nabla \\times \\mathbf{H} - (\\sigma + \\imath \\omega \\epsilon) \\mathbf{E} = 0$$\n",
    "\n",
    "where\n",
    "\n",
    "- $\\mathbf{E}$ is the electric field (V/m)\n",
    "- $\\mathbf{H}$ is the magnetic field (A/m)\n",
    "- $\\omega = 2\\pi f$ is the angular frequency\n",
    "- $\\mu$ is the magnetic permeability, often taken to be that of free spase ($\\mu_0 = 4\\pi\\times 10^{-7}$ H/m)\n",
    "- $\\sigma$ is the electrical conductivity (S/m). \n",
    "- $\\epsilon$ is the dielectric permittivity, often taken to be that of free space ($\\epsilon = 8.85 \\times 10^{-12}$ F/m)\n",
    "\n",
    "For convienence, we will make the substitution: $\\hat{\\sigma} = \\sigma + \\imath \\omega \\epsilon$ and write Maxwell's equations as\n",
    "\n",
    "$$\\nabla \\times \\mathbf{E} + \\imath\\omega \\mu \\mathbf{H} = 0$$\n",
    "\n",
    "$$\\nabla \\times \\mathbf{H} - \\hat{\\sigma} \\mathbf{E} = 0$$\n",
    "\n",
    "The first equation is [Faraday's Law](http://em.geosci.xyz/content/maxwell1_fundamentals/formative_laws/faraday.html), and the second is [Ampere's Law](http://em.geosci.xyz/content/maxwell1_fundamentals/formative_laws/ampere_maxwell.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the Magnetotelluric problem, we are interested in examining Maxwell's equations for a plane wave source. We consider a vertically propagating plane wave. For a 1D earth model, the fields and fluxes are defined by horizontal, orthogonal electric and magnetic fields, so we take\n",
    "\n",
    "$$\\mathbf{E} = E_x\\mathbf{\\hat{x}}$$\n",
    "$$\\mathbf{H} = H_y\\mathbf{\\hat{y}}$$\n",
    "\n",
    "<img src=\"http://em.geosci.xyz/_images/planewavedown1.png\" alt=\"plane_wave\" href=\"http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/index.html\" width=300px>\n",
    "\n",
    "The coordinate system we use is right-handed, and $z$ is positive up. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, our governing equations simplify to scalar equations\n",
    "\n",
    "$$ \\frac{\\partial E_x}{\\partial z}  + \\imath \\omega \\mu H_y = 0$$\n",
    "\n",
    "$$-\\frac{\\partial H_y}{\\partial z} +  \\hat{\\sigma} E_x = 0$$\n",
    "\n",
    "with the boundary conditions:\n",
    "\n",
    "$$E_x (z=0) = 1$$\n",
    "\n",
    "$$E_x (z=-\\infty) = 0$$\n",
    "\n",
    "To solve the forward problem, the \n",
    "- **knowns** are: $\\omega$, $\\mu$, $\\hat{\\sigma}$, boundary conditions\n",
    "- **unknowns** are: $E_x$, $H_y$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discretiation, the Short version. \n",
    "\n",
    "**TL;DR.** Here is the answer. \n",
    "If you want to see the full derivation, checkout <a href=\"#Discretization,-the-Gory-Details.\">The Gory Details</a> \n",
    "\n",
    "\n",
    "We define physical properties at cell centers, and stagger the electric and magnetic fields\n",
    "\n",
    "- $\\sigma$, $\\mu$, $\\epsilon$ : cell centers\n",
    "- $E_x$: cell centers\n",
    "- $H_y$: faces\n",
    "\n",
    "<img src=\"https://github.com/simpeg/tle-magnetotelluric_inversion/raw/master/images/1DMT_discretize.png\" width=200px> \n",
    "\n",
    "and use a finite difference approach to define the operators, this gives us the discrete system of equations\n",
    "\n",
    "$$\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "        \\mathbf{Grad} & \\imath \\omega \\mathbf{M}^{f}_{\\mu} \\\\[0.3em]\n",
    "        \\mathbf{M}^{cc}_{\\hat{\\sigma}} & \\mathbf{Div}           \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{A}}\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "       \\mathbf{e_x} \\\\[0.3em]\n",
    "       \\mathbf{h_y} \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{u}}\n",
    "=\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "       - \\mathbf{B}\\mathbf{e_x}^{BC} \\\\[0.3em]\n",
    "       \\boldsymbol{0} \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{rhs}}\n",
    "$$\n",
    "\n",
    "with \n",
    "\n",
    "- $\\mathbf{e_x}$: Discrete $E_x$, on cell centers $[\\text{nC} \\times 1]$\n",
    "\n",
    "- $\\mathbf{h_y}$: Dicrete $H_x$, on cell faces $[(\\text{nC}+1) \\times 1]$\n",
    "\n",
    "- $ \\mathbf{Grad}$: Discrete gradient operator $[\\text{nC} \\times (\\text{nC}+1)]$\n",
    "\n",
    "- $ \\mathbf{Div}$: Discrete divergence operator $[(\\text{nC}+1) \\times \\text{nC}]$\n",
    "\n",
    "- $\\mathbf{M}^{f}_{\\boldsymbol{\\mu}} = \\mathbf{diag}(\\mathbf{Av^{cc2f}}  \\boldsymbol{\\mu})$ $[(\\text{nC}+1) \\times (\\text{nC}+1)]$\n",
    "\n",
    "- $\\mathbf{M}^{cc}_{\\boldsymbol{\\hat{\\sigma}}} = \\mathbf{diag}(\\boldsymbol{\\hat{\\sigma}})$ $[\\text{nC} \\times \\text{nC}]$\n",
    "\n",
    "- $\\mathbf{B} \\mathbf{e_x}^{BC}$ handles the boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Designing a mesh\n",
    "\n",
    "When designing a mesh, we need to ensure that we can capture the physics (what should the thickness of the finest cells be?) and make sure that the boundary is far enough away so that the fields have decayed (how far do we need to go to approximate $\\infty$??). To address these, we look at the [skin depth equation](http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/frequency/analytic_solution.html#attenuation-and-skin-depth), which tells us over what distance we expect electromagnetic fields to have decayed by a factor of $1/e$ in a conductive medium:\n",
    "\n",
    "$$\n",
    "\\delta = \\frac{500}{\\sqrt{\\sigma f}}\n",
    "$$\n",
    "\n",
    "- The finest cells capture the behaviour of the highest frequency near the surface\n",
    "- The mesh needs to extend far enough so that the fields at the lowest frequency have sufficiently decayed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets start by considering:\n",
    "- a half-space with conductivity $\\sigma = 10^{-2}$ S/m \n",
    "- a maximum frequency of 1000 Hz\n",
    "- a minimum frequency of 0.01 Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sigma_halfspace = 1e-2\n",
    "fmax, fmin = 1e3, 1e-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def skin_depth(sigma, f):\n",
    "    return 500./np.sqrt(sigma*f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The minimum skin depth is 158.1 m\n",
      "The maximum skin depth is 50000.0 m\n"
     ]
    }
   ],
   "source": [
    "skin_depth_min = skin_depth(sigma_halfspace, fmax)\n",
    "skin_depth_max = skin_depth(sigma_halfspace, fmin)\n",
    "\n",
    "print(\"The minimum skin depth is {:2.1f} m\".format(skin_depth_min))\n",
    "print(\"The maximum skin depth is {:2.1f} m\".format(skin_depth_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ensure that we are capturing the physics and have a sufficiently far boundary, we will choose\n",
    "- a minimum cell size of $\\delta_{\\text{min}}/4$\n",
    "- a padding distance that extends to $2 \\delta_{\\text{max}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The smallest cell size should be 39.5 m\n",
      "The mesh should extend 1.0e+05 m\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    \"The smallest cell size should be {:2.1f} m\".format(\n",
    "        skin_depth_min / 4.\n",
    "    )\n",
    ")\n",
    "print(\n",
    "    \"The mesh should extend {:1.1e} m\".format(\n",
    "        skin_depth_max * 2.\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up a mesh\n",
    "\n",
    "Here, we use the [SimPEG Mesh class](http://docs.simpeg.xyz) to set up the mesh, differential operators, and handy properties and methods that handle counting and plotting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mesh extends 1.2e+05m, is that far enough?\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAADTCAYAAACss0uSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHTZJREFUeJzt3XuUFeWZ7/HvA32DeOHSmmG6hcYTw9UGpDFwNIphBm9Z\naoacOXFhIroMXllHV0wCYsYsOUYjOSZhKYPMDDIMJIg4jnq8cBPCingJqCDYjSAQaXTJ7diRoNJ0\nP+ePqm52N3v33r27+lLdv89ae+2q933rrfd9uuind1Wxy9wdERER6fi6tfcAREREJDNK2iIiIjGh\npC0iIhITStoiIiIxoaQtIiISE0raIiIiMaGkLSIiEhNK2iIiIjGhpC0iIhITOe09gMYKCwu9pKQk\n0j6PHTtGXl5epH3GkeIQUBwCikNAcQgoDoH2isOmTZsOuvsZ6dp1uKRdUlLCxo0bI+2zoqKCwYMH\nR9pnHCkOAcUhoDgEFIeA4hBorziY2Z8zaafT4yIiIjGhpC0iIhITStoiIiIx0eGuaYuISDxVV1dT\nWVnJF1980d5DyVp1dTXl5eWt1n9BQQHFxcXk5uZmtb2StoiIRKKyspJTTz2VkpISzKy9h5OVzz//\nnB49erRK3+7OoUOHqKysZODAgVn1odPjIiISiS+++IK+ffvGNmG3NjOjb9++LToToaQtIiKRUcJu\nWkvjo6QtIiISE0raIiLSPpYsgZIS6NYteF+ypL1HBMCnn37K3Llz23sYSSlpi4hI21uyBKZOhT//\nGdyD96lTO0TiziZpuzu1tbWtNKITlLRFRKTtzZwJR482LDt6NChvoUWLFlFaWsqIESP4/ve/z4ED\nB5g0aRJjxoxhzJgxvPrqqwD8/Oc/58Ybb2T8+PGcffbZzJkzB4Dp06fzwQcfMHLkSH784x8DMHv2\nbMaMGUNpaSn33XcfAHv27GHIkCHcdtttnHfeeezdu5cpU6YwfPhwzj33XH7961+3eC6N6b98iYhI\n2/vww+aVZ2jbtm088MADvPrqqxQWFnL48GHuuOMO7rrrLi688EI+/PBDLr300vr/i11RUcHatWv5\n7LPPGDRoEFOmTOGhhx5i69atvPPOOwCsXLmSHTt28Oabb+LuXHXVVaxfv57+/fuzfft2nnjiCebO\nncumTZvYt28fW7duBYJP7FFT0hYRkbbXv39wSjxZeQu88sorfPe736WwsBCAPn36sHr1at577736\nNn/5y1/47LPPALjyyivJz88nPz+fM888k08++eSkLz5ZuXIlK1euZNSoUQAcOXKEHTt20L9/fwYM\nGMDYsWMBOPvss9m1axfTpk3jyiuvZOLEiS2aSzJK2iIi0vYeeCC4hp14irxnz6C8Bdz9pP9WVVtb\ny2uvvZb0S1Py8/Prl7t3705NTc1JSdvdmTFjBjfffHOD8j179vCVr3ylfr13795s3ryZFStW8Nhj\nj7Fs2TIWLFjQovk0pmvaIiLS9iZPhvnzYcAAMAve588PyltgwoQJLFu2jEOHDgFw+PBhJk6cyKOP\nPlrfpu60dyqnnnpq/SdxgEsvvZQFCxZw5MgRAPbt28f+/ftP2u7gwYPU1tYyadIkZs2axVtvvdWi\nuSSjT9oiItI+Jk9ucZJubNiwYcycOZOLL76Y7t27M2rUKObMmcPtt99OaWkpx48f56KLLmLevHkp\n++jbty8XXHABw4cP5/LLL2f27NmUl5czbtw4AE455RQWL15M9+7dG2y3b98+brjhhvq7yB988MFI\n5wZg7h55py1RVlbmGzdujLRPPdw9oDgEFIeA4hBQHAJRxKG8vJwhQ4ZENKL20ZrfPV4nWZzMbJO7\nl6XbNu3pcTNbYGb7zWxrinozszlmttPMtpjZeY3qTzOzfWb2aLLtRUREJDOZXNNeCFzWRP3lwDnh\nayrwz43qZwF/yGZwIiIickLapO3u64HDTTS5GljkgdeBXmbWD8DMRgNfBVZGMVgREZGuLIob0YqA\nvQnrlUCRmX0C/B/g+8CEpjows6kEn9IpKiqioqIigmGdcPDgwcj7jCPFIaA4BBSHgOIQiCIO1dXV\nfP755xGNqH0cP3681edQXV2ddayjSNrJnjPmwG3Ai+6+N92jyNx9PjAfghvRor4pRDeaBBSHgOIQ\nUBwCikMgqhvRWvsmrtbWFjei5ebmZh3rKJJ2JXBWwnox8BEwDvimmd0GnALkmdkRd58ewT5FRES6\nnCi+XOU54AfhXeRjgSp3/9jdJ7t7f3cvAe4muO6thC0iIq1qzpw5DBkyhMkR/x/wjiDtJ20z+z0w\nHig0s0rgPiAXwN3nAS8CVwA7gaPADa01WBER6Vxeew3WrYPx4yH87pIWmzt3Li+99BIDBw6MpsMO\nJG3Sdvdr09Q7cHuaNgsJ/uuYiIh0AXfeCWm+LZSqKtiyBWproVs3KC2F009P3X7kSPjNb5ru85Zb\nbmHXrl1cddVVXHfddTz77LP116mfeOIJBg0aRE1NDT/96U9ZsWIFZsYPf/hDpk2bxqZNm7jzzjs5\nevQohYWFLFy4kH79+jFnzhzmzZtHTk4OQ4cOZenSpc0PSET0NaYiItIuqqqChA3Be1VV00k7E/Pm\nzePll19m7dq15OXl8aMf/YicnBxWr17NPffcw9NPP838+fPZvXs3b7/9Njk5ORw+fJjq6mqmTZvG\n0qVL6d+/P08++SQzZ85kwYIFPPTQQ+zevZv8/PxWedxmcyhpi4hI5NJ9Iobg1PiECXDsGOTlwZIl\n0Z0iB6iqquL6669nx44dmBnV1dUArF69mltuuYWcnCAF9unTh61bt7J161a+/e1v061bN2pqaujX\nrx8ApaWlTJ48mWuuuYZrrrkmugFmQUlbRETaxbhxsGZN9Ne06/zsZz/jkksu4ZlnnmHPnj2MHz8e\nSP74Tndn2LBhvPLKKyf9l68XXniB9evX89xzzzFr1iy2bdtWn/Dbmh7NKSIi7WbcOJgxI/qEDcEn\n7aKiIgAWLlxYXz5x4kTmzZvH8ePHgeDxnYMGDeLAgQO88cYbQPAFKNu2baO2tpa9e/dyySWX8PDD\nD/Ppp5/WP6KzPShpi4hIp/STn/yEGTNmcMEFF1BTU1NfftNNN9G/f39KS0sZMWIEv/vd78jLy2P5\n8uXce++9jBgxgpEjR7JhwwZqamq47rrrOPfccxk1ahR33XUXvXr1arc56fS4iIh0Knv27AGgsLCQ\n999/v7581qxZAOTk5PDII4/wyCOPNNhu5MiRrFq16qTT43/84x9bd8DNoE/aIiIiMaGkLSIiEhNK\n2iIiEpng+7YklZbGR0lbREQiUVBQwKFDh5S4U3B3Dh06REFBQdZ96EY0ERGJRHFxMZWVlRw4cKC9\nh5K16upqcnNzW63/goICiouLs95eSVtERCKRm5sb+4d0dPTnq+v0uIiISEwoaYuIiMSEkraIiEhM\nKGmLiIjEhJK2iIhITChpi4iIxISStoiISEwoaYuIiMSEkraIiEhMpE3aZrbAzPab2dYU9WZmc8xs\np5ltMbPzwvKRZvaamW0Ly/9n1IMXERHpSjL5pL0QuKyJ+suBc8LXVOCfw/KjwA/cfVi4/W/MrFf2\nQxUREena0n73uLuvN7OSJppcDSzy4LEur5tZLzPr5+7vJ/TxkZntB84APm3hmEVERLqkKB4YUgTs\nTVivDMs+risws/OBPOCDZB2Y2VSCT+kUFRVRUVERwbBOOHjwYOR9xpHiEFAcAopDQHEIKA6Bjh6H\nKJK2JSmrf5iqmfUD/gO43t1rk3Xg7vOB+QBlZWUe9RNWOvpTW9qK4hBQHAKKQ0BxCCgOgY4ehyju\nHq8EzkpYLwY+AjCz04AXgHvd/fUI9iUiItJlRZG0nwN+EN5FPhaocvePzSwPeIbgevdTEexHRESk\nS0t7etzMfg+MBwrNrBK4D8gFcPd5wIvAFcBOgjvGbwg3/UfgIqCvmU0Jy6a4+zsRjl9ERKTLyOTu\n8WvT1Dtwe5LyxcDi7IcmIiIiifSNaCIiIjGhpC0iIhITStoiIiIxoaQtIiISE0raIiIiMaGkLSIi\nEhNK2iIiIjGhpC0iIhITStoiIiIxoaQtIiISE0raIiIiMaGkLSIiEhNK2iIiIjGhpC0iIhITStoi\nIiIxoaQtIiISE0raIiIiMaGkLSIiEhNK2iIiIjGhpC0iIhITaZO2mS0ws/1mtjVFvZnZHDPbaWZb\nzOy8hLrrzWxH+Lo+yoFnpHdvMGPQkCFg1vJXt27Be05O8F5YGLy6dWu4XFICt90WvNetL1kSjGnJ\nkuTlddLVN7ddglOff77Z27R0n+3abwotjkPU2nj+kYv7+LPVVecdFcUvO+7e5Au4CDgP2Jqi/grg\nJcCAscAbYXkfYFf43jtc7p1uf6NHj/ZI9OrlDn4+rzrUZPGq9j584hex1q/hab+Gp/0W5vrj3OS/\nYHr9+wbG+gbG1i87nLTu4Bvyx/svvvFfviF/fH2Zg3vPnr7hvpf9F79w33Dfy+49e6au3xDObfHi\npO188WLfsMEbtvUT29QUFCTdpk7KbdPsM5km+2pBv83uP8n+0sUhE1nvP8l4spl/FGMsLy9v8T5a\nc/ypRBb7UFZxaId5p9PSuERyPGSqA8avTiZxiPoYdHcHNnqa/OjuWNC2aWZWAvxfdx+epO5xYJ27\n/z5c3w6Mr3u5+83J2qVSVlbmGzduTDumDAbNN3iVNxnX8r6a2g0OOE43ulHLf2MnH/A1asP1UjYD\nsIURDcpO5y8AVHFaWNedbtQ0qDupvhuUlsLp5a/Dl1+cNJaq3EK21AyntpYTbU8PK19Pvg35BTB2\nLFVVsGULybfNYPsG40jXV6Jm9JtV/xHsL9L9t8J4sh3j0aNH6dmzZ9b7AFpt/KlEGvtQVnFo43mn\nE0VcIjkeMtXB4pcoXRzqYu0OBQWwZg2MiyDNmNkmdy9L1y6n5buiCNibsF4ZlqUqP4mZTQWmAhQV\nFVFRUdHiQQ0C3mJ03R5a3N8JHvYXvHt9/0YtxkHOoJZu9etV9AI4qexE0u4V1nFS3Un1tc6BA9Wc\nnuxgB6qqe1IbjquubW7ucQB6ptiGL7/g6NGjHDiQQ21tbtJtM9k+Ubq+EjWn32z6j2J/Ue6/NcaT\n7Rirq6tbtA9ovfGnEmXs62QTh7aedzpRxCWK4yFTHS1+idLFITHWx445Tz11gN69D7fdADP5OA6U\nkPr0+AvAhQnra4DRwI+BexPKfwb8KN2+Ijs9Xn9qvDbyl1HtUOvdqPZ8Pvc8PvfuHPMe/NUf5ybv\nwV/r1+tOnzcuqz9tzljvYUe9e3f3Hna0Qd1J9T3C0zEDBjQ8rVTX9m++4z16eMO2dVJs4wMGuHvQ\nNuW2GWyfKG1fWfabVf8R7C/S/bfCeLIdYySnQ1tp/KlEGvtQVnFo43mnE0Vc2vT0eAeLX6J0cWiN\nY9A989PjUSTtx4FrE9a3A/2Aa4HHU7VL9dI1bV3T1jVtXdNuiq5pJ6dr2tHo6Ne0o0jaV9LwRrQ3\nw/I+wG6Cm9B6h8t90u0rsqTtXp+4a5P9RZfNyyx47949eO/bN3iZNVweMMD91luD97r1uoNx8eLk\n5XXS1Te3XYLKhx9u9jYt3We79ptCi+MQtTaef53Ifkm30/ijknUcYj7vxto0abt32Pi1eRxCmSbt\ntDeimdnvCW4qKwQ+Ae4DcsNT6/PMzIBHgcuAo8AN7r4x3PZG4J6wqwfc/Yl0p+sjuxEtQUVFBYMH\nD460zzhSHAKKQ0BxCCgOAcUh0F5xiOxGNHe/Nk29A7enqFsALEi3DxEREUlP34gmIiISE0raIiIi\nMaGkLSIiEhNK2iIiIjGhpC0iIhITStoiIiIxoaQtIiISE0raIiIiMaGkLSIiEhNK2iIiIjGhpC0i\nIhITStoiIiIxoaQtIiISE0raIiIiMaGkLSIiEhNK2iIiIjGhpC0iIhITStoiIiIxoaQtIiISE0ra\nIiIiMaGkLSIiEhMZJW0zu8zMtpvZTjObnqR+gJmtMbMtZrbOzIoT6h42s21mVm5mc8zMopyAiIhI\nV5E2aZtZd+Ax4HJgKHCtmQ1t1OxXwCJ3LwXuBx4Mt/3vwAVAKTAcGANcHNnoRUREupBMPmmfD+x0\n913ufgxYClzdqM1QYE24vDah3oECIA/IB3KBT1o6aBERka4oJ4M2RcDehPVK4BuN2mwGJgG/Bb4D\nnGpmfd39NTNbC3wMGPCou5c33oGZTQWmAhQVFVFRUdHsiTTl4MGDkfcZR4pDQHEIKA4BxSGgOAQ6\nehwySdrJrkF7o/W7gUfNbAqwHtgHHDezrwFDgLpr3KvM7CJ3X9+gM/f5wHyAsrIyHzx4cOYzyEBF\nRQVR9xlHikNAcQgoDgHFIaA4BDp6HDJJ2pXAWQnrxcBHiQ3c/SPgHwDM7BRgkrtXhZ+gX3f3I2Hd\nS8BYgsQuIiIizZDJNe0/AeeY2UAzywO+BzyX2MDMCs2srq8ZwIJw+UPgYjPLMbNcgpvQTjo9LiIi\nIumlTdrufhy4A1hBkHCXufs2M7vfzK4Km40HtpvZ+8BXgQfC8uXAB8C7BNe9N7v789FOQUREpGvI\n5PQ47v4i8GKjsn9KWF5OkKAbb1cD3NzCMYqIiAj6RjQREZHYUNIWERGJCSVtERGRmFDSFhERiQkl\nbRERkZhQ0hYREYkJJW0REZGYUNIWERGJCSVtERGRmFDSFhERiQklbRERkZhQ0hYREYkJJW0REZGY\nUNIWERGJCSVtERGRmFDSFhERiQklbRERkZhQ0hYREYkJJW0REZGYUNIWERGJCSVtERGRmMgoaZvZ\nZWa23cx2mtn0JPUDzGyNmW0xs3VmVpxQ19/MVppZuZm9Z2Yl0Q1fRESk60ibtM2sO/AYcDkwFLjW\nzIY2avYrYJG7lwL3Aw8m1C0CZrv7EOB8YH8UAxcREelqMvmkfT6w0913ufsxYClwdaM2Q4E14fLa\nuvowuee4+yoAdz/i7kcjGbmIiEgXk5NBmyJgb8J6JfCNRm02A5OA3wLfAU41s77A14FPzew/gYHA\namC6u9ckbmxmU4GpAEVFRVRUVGQxldQOHjwYeZ9xpDgEFIeA4hBQHAKKQ6CjxyGTpG1JyrzR+t3A\no2Y2BVgP7AOOh/1/ExgFfAg8CUwB/q1BZ+7zgfkAZWVlPnjw4IwnkImKigqi7jOOFIeA4hBQHAKK\nQ0BxCHT0OGRyerwSOCthvRj4KLGBu3/k7v/g7qOAmWFZVbjt2+Gp9ePAfwHnRTJyERGRLiaTpP0n\n4BwzG2hmecD3gOcSG5hZoZnV9TUDWJCwbW8zOyNc/xbwXsuHLSIi0vWkTdrhJ+Q7gBVAObDM3beZ\n2f1mdlXYbDyw3czeB74KPBBuW0Nw6nyNmb1LcKr9XyKfhYiISBeQyTVt3P1F4MVGZf+UsLwcWJ5i\n21VAaQvGKCIiIugb0URERGJDSVtERCQmlLRFRERiQklbREQkJpS0RUREYkJJW0REJCaUtEVERGJC\nSVtERCQmlLRFRERiQklbREQkJpS0RUREYkJJW0REJCaUtEVERGJCSVtERCQmlLRFRERiQklbREQk\nJszd23sMDZjZAeDPEXdbCByMuM84UhwCikNAcQgoDgHFIdBecRjg7meka9ThknZrMLON7l7W3uNo\nb4pDQHEIKA4BxSGgOAQ6ehx0elxERCQmlLRFRERioqsk7fntPYAOQnEIKA4BxSGgOAQUh0CHjkOX\nuKYtIiLSGXSVT9oiIiKxp6QtIiISE506aZvZZWa23cx2mtn09h5P1MxstplVmNkWM3vGzHol1M0I\n573dzC5NKE8aEzMbaGZvmNkOM3vSzPLC8vxwfWdYX9KWc2wOM7vbzNzMCsN1M7M54di3mNl5CW2v\nD+e6w8yuTygfbWbvhtvMMTMLy/uY2aqw/Soz6932M0zPzKaFP99tZvZwQnmXOR7MbKSZvW5m75jZ\nRjM7Pyzv9MeDmf2P8Gdfa2Zljepa/RhItY+2lioOZvb3ZrYp/JluMrNvJdQ162edzfEUCXfvlC+g\nO/ABcDaQB2wGhrb3uCKe40QgJ1z+JfDLcHloON98YGAYh+5NxQRYBnwvXJ4H3Bou3wbMC5e/BzzZ\n3vNOEYuzgBUEX8xTGJZdAbwEGDAWeCMs7wPsCt97h8u9w7o3gXHhNi8Bl4flDwPTw+XpdbHuSC/g\nEmA1kB+un9kVjwdgZcLP7QpgXVc5HoAhwCBgHVCWUN7qx0CqfXSwOIwC/jZcHg7sS6hr1s86m+Mp\nildn/qR9PrDT3Xe5+zFgKXB1O48pUu6+0t2Ph6uvA8Xh8tXAUnf/0t13AzsJ4pE0JuFflN8Clofb\n/ztwTUJf/x4uLwcm1P0F2sH8GvgJkHhn5dXAIg+8DvQys37ApcAqdz/s7v8PWAVcFtad5u6vefCv\nbxHJ45AYn47kVuAhd/8SwN33h+Vd7Xhw4LRw+XTgo3C50x8P7l7u7tuTVLXFMZBqH20uVRzc/W13\nrzsetgEF4ZmDbH7WzTqeoppbZ07aRcDehPXKsKyzupHgrz5IPfdU5X2BTxP+AEiMVf02YX1V2L7D\nMLOrCP5i3tyoqrlxKAqXG5cDfNXdPwYI38+MbALR+TrwzfCU5R/MbExY3qWOB+BOYLaZ7QV+BcwI\ny7va8ZCoLY6BuP3OnQS8Hf6Rm83PurkxjUROVB11QMn++o/d/28zs9XA3ySpmunuz4ZtZgLHgSV1\nmyVp7yT/I82baN9UX22qqTgA9xBcKjhpsyRlTc23Q8y1KWnikENwOm4sMAZYZmZn0/WOhwnAXe7+\ntJn9I/BvwN/RSY6HTH4nJNssSVnUx0CbxivLONRtO4zgkmLd741sxt4ucejMSbuS4DpnnWJOnCaL\nDXf/u6bqw5scvg1MCE/rQNNzT1Z+kODUTk74l3Ni+7q+Ks0sh+B04+HsZ5SdVHEws3MJrp9tDs/S\nFgNvhTcfpYpDJTC+Ufm6sLw4SXuAT8ysn7t/HJ4C2087aOp4MLNbgf8Mj4M3zayW4OEHXeZ4ADCz\nRcD/ClefAv41XO4Ux0O63wkptMUx0Ka/c7OMA2ZWDDwD/MDdPwiLs/lZN/d4ikZUF8c72ovgD5Jd\nBL/Q626wGNbe44p4jpcB7wFnNCofRsMbQnYR3HCSMiYEv9wSbzq5LVy+nYY3nSxr73mnickeTtyI\ndiUNbxR5MyzvA+wm+FTaO1zuE9b9KWxbdzPKFWH5bBrejPJwe881ydxvAe4Pl79OcIrOutrxAJQD\n48PlCcCmrnY8cPINWK1+DKTaRweLQ69wjJOStG3Wzzqb4ymSObX3wdXKP7ArgPcJ7mKc2d7jaYX5\n7ST4xfxO+JqXUDcznPd2wrsgm4oJwd2jb4Z9PsWJO5ALwvWdYf3Z7T3vNDHZw4mkbcBj4VzfbfSP\n98ZwTjuBGxLKy4Ct4TaPcuJbA/sCa4Ad4Xtk/wgjnHsesDgc/1vAt7ri8QBcCGwKfzm/AYzuKscD\n8B2CT3pfAp8AK9ryGEi1j44SB+Be4K+c+J35Dif+l0WzftbZHE9RvPQ1piIiIjHRme8eFxER6VSU\ntEVERGJCSVtERCQmlLRFRERiQklbREQkJpS0Rbo4M1sefnNapu3PNbOFrTgkEUlBSVukCwu/zrG7\nu+/KdBt3fxcoNrP+rTcyEUlGSVukkzKzWyx4pvQ7ZrbbzNYmaTYZeDZhmyNm9svwWcOrzex8M1tn\nZrvCB7PUeZ7g27BEpA0paYt0Uu4+z91HEjw8pBJ4JEmzCwi+PazOVwiePz0a+Az438DfE3zD1P0J\n7TYC32yNcYtIap35gSEiEvgt8Iq7P5+krh9wIGH9GPByuPwu8KW7V5vZu0BJQrv9wN+2wlhFpAlK\n2iKdmJlNAQYAd6Ro8jnBd0nXqfYT321cS/Ddzbh7bfhEpzoF4bYi0oZ0elykkzKz0cDdwHXuXpui\nWTnwtSy6/zrBwxVEpA0paYt0XncQPCZwbXgz2r8mafMCDZ/9m6lLwm1FpA3pKV8iXZiZ9QDWAhe4\ne02G2+QDfwAudPfjrTk+EWlISVukizOzS4Fyd/8ww/bnAEXuvq5VByYiJ1HSFhERiQld0xYREYkJ\nJW0REZGYUNIWERGJCSVtERGRmFDSFhERiYn/D4+s94xwmLVKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1156e16d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cs = 39.  # core cell size\n",
    "npad = 25  # number of padding cells\n",
    "ncz = 100   # number of core cells\n",
    "\n",
    "# define a tensor mesh\n",
    "hz = [(cs, npad, -1.3), (cs, ncz)] \n",
    "mesh = Mesh.TensorMesh([hz], x0='N')  # put the origin at the surface \n",
    "\n",
    "# plot the mesh\n",
    "fig, ax = plt.subplots(1,1, figsize=(8, 3))\n",
    "mesh.plotGrid(centers=True, faces=True, ax=ax)\n",
    "ax.legend([\"centers\", \"faces\"])\n",
    "ax.invert_xaxis()  # so that the surface is on our left hand side\n",
    "ax.set_xlabel('z (m)')\n",
    "ax.grid(which=\"both\", linewidth=0.5)\n",
    "\n",
    "print(\n",
    "    \"The mesh extends {:1.1e}m, is that far enough?\".format(\n",
    "        mesh.hx.sum()\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  ---- 1-D TensorMesh ----  \n",
      "   x0: -122984.33\n",
      "  nCx: 125\n",
      "   hx: 27520.00,  21169.23,  16284.02,  12526.17,  9635.52,  7411.94,  5701.49,  4385.76,  3373.66,  2595.12,  1996.25,  1535.58,  1181.21,  908.63,  698.94,  537.65,  413.58,  318.13,  244.72,  188.25,  144.80,  111.39,  85.68,  65.91,  50.70,  100*39.00,\n"
     ]
    }
   ],
   "source": [
    "print(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assemble the discrete system of equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Model parameters\n",
    "\n",
    "We start with a half space that has physical properties\n",
    "- $\\sigma = 10^{-2}$ S/m\n",
    "- $\\mu = \\mu_0$\n",
    "- $\\epsilon = \\epsilon_0$\n",
    "\n",
    "and define these on the mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 125 cell centers. \n",
      " sigma   is 125 elements long, all cells have a value of 1.00e-02 S/m \n",
      " mu      is 125 elements long, all cells have a value of 1.26e-06 H/m \n",
      " epsilon is 125 elements long, all cells have a value of 8.85e-12 F/m \n",
      "\n"
     ]
    }
   ],
   "source": [
    "sigma = np.ones(mesh.nC)*sigma_halfspace # conductivity values for all cells\n",
    "mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n",
    "epsilon = np.ones(mesh.nC)*epsilon_0 # dielectric constant values for all cells\n",
    "\n",
    "print(\n",
    "    \"There are {:1.0f} cell centers. \\n\"\n",
    "    \" sigma   is {:1.0f} elements long, all cells have a value of {:1.2e} S/m \\n\"\n",
    "    \" mu      is {:1.0f} elements long, all cells have a value of {:1.2e} H/m \\n\"\n",
    "    \" epsilon is {:1.0f} elements long, all cells have a value of {:1.2e} F/m \\n\".format(\n",
    "        mesh.nC, \n",
    "        len(sigma), sigma_halfspace,\n",
    "        len(mu), mu_0,\n",
    "        len(epsilon), epsilon_0\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will pick a single frequency to work with for now\n",
    "- f = 1000 Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "frequency = 1e3 # Frequency (Hz)\n",
    "omega = 2*np.pi*frequency # Angular frequency (rad/s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we will adopt the quasistatic assumption and ignore displacement current $(i \\epsilon \\omega)$. To explore the impacts of this assumption, uncomment the next second line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sigmahat = sigma  # quasi-static assumption\n",
    "# sigmahat = sigma + 1j*epsilon*omega  # includes displacement current"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system we want to solve is\n",
    "$$\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "        \\mathbf{Grad} & \\imath \\omega \\mathbf{M}^{f}_{\\mu} \\\\[0.3em]\n",
    "        \\mathbf{M}^{cc}_{\\hat{\\sigma}} & \\mathbf{Div}           \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{A}}\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "       \\mathbf{e_x} \\\\[0.3em]\n",
    "       \\mathbf{h_y} \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{u}}\n",
    "=\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "       - \\mathbf{B}\\mathbf{e_x}^{BC} \\\\[0.3em]\n",
    "       \\boldsymbol{0} \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{rhs}}\n",
    "$$\n",
    "\n",
    "so we need to construct each of the operators. For details, see: <a href=\"#Discretization,-the-Gory-Details.\">The Gory Details</a> \n",
    "\n",
    "We start by laying our the piece and will then assemble the full matrix system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Grad \n",
    "mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n",
    "Grad = mesh.cellGrad # Gradient matrix\n",
    "\n",
    "# MfMu\n",
    "Mmu = Utils.sdiag(mesh.aveCC2F * mu) \n",
    "\n",
    "# Mccsigma\n",
    "Msighat = Utils.sdiag(sigmahat) \n",
    "\n",
    "# Div\n",
    "Div = mesh.faceDiv # Divergence matrix\n",
    "\n",
    "# Right Hand Side\n",
    "B = mesh.cellGradBC  # a matrix for boundary conditions\n",
    "Exbc = np.r_[0., 1.] # boundary values for Ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Assemble the matrix\n",
    "\n",
    "# A-matrix\n",
    "A = sp.vstack([\n",
    "    sp.hstack([Grad, 1j*omega*Mmu]), # Top row of A matrix\n",
    "    sp.hstack((Msighat, Div)) # Bottom row of A matrix\n",
    "])\n",
    "\n",
    "# Right-hand side\n",
    "rhs = np.r_[\n",
    "    -B*Exbc, \n",
    "    np.zeros(mesh.nC)\n",
    "] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have all of the pieces, we can go ahead and solve the MT system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.97 ms, sys: 3.19 ms, total: 6.16 ms\n",
      "Wall time: 20 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Ainv = Solver(A) # Factorize A matrix\n",
    "sol = Ainv*rhs   # Solve A^-1 rhs = sol\n",
    "Ex = sol[:mesh.nC] # Extract Ex from solution vector u\n",
    "Hy = sol[mesh.nC:mesh.nC+mesh.nN] # Extract Hy from solution vector u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Impedance, Apparent Resistivity, and Phase\n",
    "\n",
    "MT data are natural source data, meaning that the source is free! but we don't know its amplitude. To account for this, the data we examine are typically transfer functions that involve ratios of the electric and magnetic fields. For the 1D problem, the Impedance is simply given by\n",
    "\n",
    "$$\n",
    "Z_{xy} = - \\frac{E_x}{H_y}\n",
    "$$\n",
    "\n",
    "(The negative is because we have defined a coordinate system such that z is positive up) \n",
    "\n",
    "$Z_{xy}$ is a complex number, so we can look at real and imaginary components or amplitude and phase. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Impedance: 6.2e-01 + 6.4e-01i\n",
      "or in terms of Amplidude: 8.9e-01 and phase: 45.9 degrees\n"
     ]
    }
   ],
   "source": [
    "Zxy = - 1./Hy[-1] # Impedance at the surface\n",
    "\n",
    "print(\"Impedance: {:1.1e} + {:1.1e}i\".format(Zxy.real, Zxy.imag))\n",
    "print(\"or in terms of Amplidude: {:1.1e} and phase: {:1.1f} degrees\".format(\n",
    "    np.absolute(Zxy), np.rad2deg(np.arctan(Zxy.imag / Zxy.real)))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often is useful to translate the impedance to an apparent resistivity ($\\rho_a$) and phase. \n",
    "\n",
    "$$\n",
    "\\rho_a = \\frac{1}{\\mu_0\\omega} \\big|Z_{xy}\\big|^2\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\phi = \\tan^{-1}\\left(\\frac{\\text{Im}(Z_{xy})}{\\text{Re}(Z_{xy})}\\right)\n",
    "$$\n",
    "\n",
    "For a  half-space, we expect the apparent resistivity to equal the true resistivity, and the phase to be $45^\\circ$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Apparent Resistivity: 100.0, Phase: 45.9\n"
     ]
    }
   ],
   "source": [
    "app_res = abs(Zxy)**2 / (mu_0*omega)\n",
    "app_phase = np.rad2deg(np.arctan(Zxy.imag / Zxy.real))\n",
    "\n",
    "print(\n",
    "    \"Apparent Resistivity: {:1.1f}, Phase: {:1.1f}\".format(\n",
    "        app_res, app_phase\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note that the apparent resistivity, 100.0 is the same as the true half-space 100.0\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    \"Note that the apparent resistivity, {:1.1f} \"\n",
    "    \"is the same as the true half-space {:1.1f}\".format(\n",
    "        app_res,\n",
    "        1./sigma_halfspace  \n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Put it all together\n",
    "\n",
    "Here, we define a function that performs an MT simulation so that we can readily compute the Magnetotelluric response at multiple frequencies and for a variety of models. We write this function to the file MTsimulation.py so that we can import it and use it in later notebooks. Uncomment the first three lines to write out the file again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %%writefile MTforward.py\n",
    "# import numpy as np\n",
    "# import scipy.sparse as sp\n",
    "# from scipy.constants import mu_0\n",
    "# from SimPEG import Utils, Solver\n",
    "\n",
    "\n",
    "def simulateMT(mesh, sigma, frequency, rtype=\"app_res\"):\n",
    "    \"\"\"\n",
    "       Compute apparent resistivity and phase at each frequency. \n",
    "       Return apparent resistivity and phase for rtype=\"app_res\",\n",
    "       or impedance for rtype=\"impedance\" \n",
    "    \"\"\"\n",
    "    \n",
    "    # Angular frequency (rad/s)\n",
    "    def omega(freq):\n",
    "        return 2*np.pi*freq\n",
    "    \n",
    "    # make sure we are working with numpy arrays\n",
    "    if type(frequency) is float:\n",
    "        frequency = np.r_[frequency]  # make it a list to loop over later if it is just a scalar\n",
    "    elif type(frequency) is list: \n",
    "        frequency = np.array(frequency)\n",
    "    \n",
    "    # Frequency independent pieces of the A matrix\n",
    "    # Grad \n",
    "    mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n",
    "    Grad = mesh.cellGrad # Gradient matrix\n",
    "\n",
    "    # MfMu\n",
    "    mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n",
    "    Mmu = Utils.sdiag(mesh.aveCC2F * mu) \n",
    "\n",
    "    # Mccsigma\n",
    "    sigmahat = sigma  # quasi-static assumption\n",
    "    Msighat = Utils.sdiag(sigmahat) \n",
    "\n",
    "    # Div\n",
    "    Div = mesh.faceDiv # Divergence matrix\n",
    "\n",
    "    # Right Hand Side\n",
    "    B = mesh.cellGradBC  # a matrix for boundary conditions\n",
    "    Exbc = np.r_[0., 1.] # boundary values for Ex\n",
    "    \n",
    "    # Right-hand side\n",
    "    rhs = np.r_[\n",
    "        -B*Exbc, \n",
    "        np.zeros(mesh.nC)\n",
    "    ] \n",
    "    \n",
    "    # loop over frequencies \n",
    "    Zxy = []\n",
    "    for freq in frequency: \n",
    "\n",
    "        # A-matrix\n",
    "        A = sp.vstack([\n",
    "            sp.hstack([Grad, 1j*omega(freq)*Mmu]), # Top row of A matrix\n",
    "            sp.hstack((Msighat, Div)) # Bottom row of A matrix\n",
    "        ])\n",
    "    \n",
    "        Ainv = Solver(A) # Factorize A matrix\n",
    "        sol = Ainv*rhs   # Solve A^-1 rhs = sol\n",
    "        Ex = sol[:mesh.nC] # Extract Ex from solution vector u\n",
    "        Hy = sol[mesh.nC:mesh.nC+mesh.nN] # Extract Hy from solution vector u\n",
    "\n",
    "        Zxy.append(- 1./Hy[-1]) # Impedance at the surface\n",
    "    \n",
    "    # turn it into an array\n",
    "    Zxy = np.array(Zxy)\n",
    "\n",
    "    # return impedance or apparent resistivity and phase \n",
    "    if rtype.lower() == \"impedance\":\n",
    "        return Zxy\n",
    "\n",
    "    elif rtype.lower() == \"app_res\":\n",
    "        app_res = abs(Zxy)**2 / (mu_0*omega(frequency))\n",
    "        app_phase = np.rad2deg(np.arctan(Zxy.imag / Zxy.real))\n",
    "        return app_res, app_phase\n",
    "    \n",
    "    else:\n",
    "        raise Exception, \"rtype must be 'impedance' or 'app_res', not {}\".format(rtype.lower())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 86.1 ms, sys: 2.2 ms, total: 88.3 ms\n",
      "Wall time: 93.7 ms\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "\n",
    "# Run the simulation over 25 frequencies from 1e-3 Hz to 100 Hz\n",
    "frequencies = np.logspace(-2, 3, 25)\n",
    "\n",
    "# for freq in frequencies:\n",
    "app_res, phase = simulateMT(mesh, sigma, frequencies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a half-space, the apparent resistivity should equal the true resistivity and the phase should be $45^\\circ$. How did we do??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYZVV56P/vW1XdQDdg44xghERU0ORqbDUmMZKgCSYi\nXq9G0USN3HDlUWPUGHHIT+PzM1GT/LxOUVvB9sYpSryKqNcBJcSoKBiuoh2FoMgoQ+ymR7qG9/fH\n3qf79OlTVWfXPuOu7+ehnqq99rv2WrVqUf3W2lNkJpIkSU0yNeoOSJIk9ZsJjiRJahwTHEmS1Dgm\nOJIkqXFMcCRJUuOY4EiSpMYxwZEkSY1jgiNJkhpnYhKciDgxIt4dEedHxNmj7o8kSRpfI01wIuK8\niLglIq7sKD81In4QEVdHxDkAmbklM58P/D6wcRT9lSRJk2HUKzibgVPbCyJiGngn8ATgJOCMiDip\n3Pck4KvARcPtpiRJmiQzo2w8My+JiOM6ih8JXJ2Z1wBExEeB04HvZ+YFwAUR8Rngw92OGRFnAWcB\nrF+//uEnnHBC3/s9Pz8PwPT0dNftKmVLlS+3r0pMnfh+16/TVt2xrrOvSsxKYgdRv257zu3BcW47\nt5eLqRPf7/p12ur33L7iiituy8x7LNePkSY4izgGuK5t+3rgURFxMvAU4BDgs4tVzsxNwCaAjRs3\n5mWXXdb3Dm7duhWADRs2dN2uUrZU+XL7qsTUie93/Tpt1R3rOvuqxKwkdhD167bn3B4c57Zze7mY\nOvH9rl+nrX7P7Yi4tpd+jGOCE13KMjMvBi4eblckSdIkGscE53rgvm3bxwI3ruRA8/Pz+7K+ftq+\nffuS21XKlipfbl+VmDrx/a5fp626Y11nX5WYlcQOon7d9pzbg+Pcdm5XMUnjPey5vZhRX2TczbeA\nEyLi+IhYCzwDuKDKASLitIjYtG3btoF0UJIkjbeRruBExEeAk4G7R8T1wGsz89yIeCHweWAaOC8z\nv1fluJn5aeDTGzdu/ONBnm/sPHa3tnotW6p8uX1VYurE97t+nbbqjnWdfVViVhI7iPp123NuD45z\n27ldxSSN97DndqdR30V1xiLln2WJC4klSZKWMo7X4PSN1+D0J77f9eu05XUKg+Xcdm73apLGult7\nzu3B8RqcAfIaHEmSVrdGruB4DU5vPJfrdQq9tufcHhzntnO7ikka71Ffg9PIFRxJkrS6NXIFp8Vr\ncPoT3+/6ddryOoXBcm47t3s1SWPdrT3n9uB4Dc4AeQ2OJEmrWyNXcLwGpzeey/U6hV7bc24PjnPb\nuV3FJI231+BIkiT1mQmOJElqHBMcSZLUOI28BqfFu6j6E9/v+nXa8k6TwXJuO7d7NUlj3a095/bg\neBfVAHkXlSRJq1sjV3C8i6o3Xo3vnSa9tufcHhzntnO7ikkab++ikiRJ6jMTHEmS1DgmOJIkqXEa\neQ1Oi3dR9Se+3/XrtOWdJoPl3HZu92qSxrpbe87twfEuqgHyLipJkla3Rq7geBdVb7wa3ztNem3P\nuT04zm3ndhWTNN7eRSVJktRnJjiSJKlxTHAkSVLjmOBIkqTGaeRFxi3eJt6f+H7Xr9OWt9IOlnPb\nud2rSRrrbu05twfH28QHyNvEJUla3Rq5guNt4r3xdkNvpe21Pef24Di3ndtVTNJ4e5u4JElSn5ng\nSJKkxjHBkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS1DiNfA5Oi08y7k98v+vXacunvQ6Wc9u5\n3atJGutu7Tm3B8cnGQ+QTzKWJGl1a+QKjk8y7o1PxPRpr72259weHOe2c7uKSRpvn2QsSZLUZyY4\nkiSpcUxwJElS45jgSJKkxjHBkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS1DiVn2QcEYcA9wEO\nA27NzFv73itJkqQaelrBiYgjIuLsiLgE2AZcDVwJ3BwR10XEeyPiEYPsqCRJUq+WXcGJiJcArwGu\nAS4A3gDcCOwG7go8BHgM8MWI+Abwosy8amA9rsC3ifcnvt/167TlG5cHy7nt3O7VJI11t/ac24Mz\nLm8T7+UU1a8Cj83MKxfZ/03gvIh4PnAm8FhgpAlORJwGnHb88cePshuSJGlElk1wMvNpvRwoM+8E\n/r52j/rAt4n3xrfS+sblXttzbg+Oc9u5XcUkjbdvE5ckSeqzldxFdW+K01b3pCNBysyxWMGRJEmr\nW6UEJyL+AHgfEMDPgGzbnYzJKSpJkrS6VV3BeQPwZuD1mTk3gP5IkiTVVvUanCOBzSY3kiRpnFVN\ncD4E/N4gOiJJktQvVU9RvRT4ZEScAnwXmG3fmZmv71fHJEmSVqpqgvM/gFOB24D7c/BFxiY4kiRp\n5KomOH8BvCwz3zKIzkiSJPVD1WtwpineRyVJkjS2qiY47weeNYiOSJIk9UvVU1TrgP8eEb8DfIeD\nLzL+k351TJIkaaWqJjgnAv9Wfv2gjn2JJEnSGKiU4GTmbw6qI5IkSf2y4reJR8SvRcQh/eyMJElS\nP6w4wQE+BxzTr45IkiT1S50EJ/rWC0mSpD6qk+BIkiSNpToJzv8AftqvjiwnIp4cEe+NiE9FxG8P\nq11JkjR5VpzgZOaHM3MnQERMRcTPVT1GRJwXEbdExJUd5adGxA8i4uqIOKds75OZ+cfAc4Gnr7Tf\nkiSp+Xq+Tby8Y+oc4AzgfsAdwL8A/y9wE/Ajilc5VLEZeAfwv9ramQbeCTweuB74VkRckJnfL0Ne\nU+5f1vz8PFu3bq3YpaXdsv1ONv3LjwBYs2YtALOze/dtR3ll0uze4hmIa9au2Vd3drYoW7tmf9kB\n5WvXHtTe3r17F93XLaaXC6N6OeYg69dpq1vbvZbV3VclZiWxg6hft71u7bfmey/jvdxYR8SyYx2L\n1O8U0Z/xzjy4freHfGUu/+iv9pAD3kycyd7ZWcj9vyNafW/9XgG4s/X9rFlLlkeYnZ0ly3qdXdh/\njDX72mvFzM7OkuW+VuFSMa2yVrsdn4r+AzMzMwd9f52x2aW9ubm5oq2ZmYPbautT9+PmATEHHLuj\nfqt0bm4egOnp4p+pufk5Mov+t2Lny7Lpmel9B5ubny/j9v/zllnEAkxNH1heHGeeBGbKfZ0/p2R/\n/enp3v4ZXip+KmB6KpiKYHqK4nMEU62vp4KFuTmmpuCQtWuZDphqxZdx0xFElN9bJpmwkMlCFj+/\n1uf5sjw7Pi+0fZ6dnWUh4eWPvz8/f/d1bN++/aA+91rWy77F9DSyEXEo8BXggcAHgB8CdwVOA74J\nvLpyy0BmXhIRx3UUPxK4OjOvKdv+KHB6RGwB3gh8LjO/vURfzwLOAjj22GNX0q0lbd09x/nfub3v\nx5UkqUm27Zkbafu9ruCcA9wDeGBm3tpW/oaIeC7w7j726Rjgurbt64FHAS8CHgfcJSLun5ld28zM\nTcAmgI0bN+aGDRv62DW4/5p1vPyU4wE47LDDANi9ezcAh5bb7WWtmMw8qGyx2F73dcYceujiMVWP\nOcj6ddrq1navZXX3VYlZSewg6tdt76B5TLJ7956y7NC2evvL2v9a3b2nrN9lbu7q49xu/UW/p+zH\noW19q2LPvu/j4Paiy/JoZ1H3mOi6v/W9rVu3bt92tLUd0a0s2uoddmAfIti9a9cBx2y1FwS7dhf7\n1rf2dfQpCHa16q/vjImD6uzatYsoY1vfY9fvv0tdYF9bh69f33U/baXtZYv16YDvpe17atm5cycA\n6w9fT7RtH3744fuOs2PnDgCOKMva67WXAezcsWN//QP6F+wo9x1xxOFt5QfaubNV/wh6sWPH9q7x\nxapKsrBA+TmZX8j9X2exvWPHTuaz+H+pVbZQfm59vZD7V3JaK0JTUYx16+upiHL1p/g6WmWt/VPB\n7l07mYrgYT9/bzas278i2e3f417LetnXqdcE5wzgnI7kBoDM3BwRRwF/13OrS+t2liUz823A2/rU\nxorddf1anvWI+wD7B7p1Gqx94HstW6p8uX1VYurE97t+nbbqjnWdfVViVhI7iPp123NuD45zexzn\n9pqeyoryma7lxb7pRfdViTkwfqpS/MH1RzC31w3n1Ppiek1w7sf+d1AdJDPfArylLz0qVmzu27Z9\nLHDjSg40iGtw4OBzgXXPL9Y971j13ORKzmX2s36dtgZ9Lrff4z1JY92tPef24Di3ndtVTNJ4D3tu\nL6bXu6i2A0cvtjMiHhoR51VuvbtvASdExPERsRZ4BnBBlQNExGkRsWnbtm196pIkSZokva7gfAV4\nAfDVzh0RcW/go8AJwPOqNB4RHwFOBu4eEdcDr83McyPihcDnKe7KOi8zv1fluJn5aeDTGzdu/ONB\nLsd1Hrvu+cW65x2rfq91x2ZYS8vd2hr0udx+j/ckjXW39pzbg+Pcdm5XMUnjPey53anXBOf1wDci\n4oPAm4Cr2H8X1WuAH1MkOJVk5hmLlH8W+GzV40mSJEGPCU5mXhkRpwLnAVe07ZoD3gq8Hbi2/92r\nx2tw+hPf7/p12vI6hcFybju3ezVJY92tPef24IzLNTg9P+gvM78aEQ8CHgEcT3Fdztcz8z8jYj3w\nl5VbH5CIOA047fjjjx91VyRJ0ggsm+BExPGZ+SOAzFwALi0/9ilf2fCXUTyY4NjMvO7gIw2P1+D0\nxnO5XqfQa3vO7cFxbju3q5ik8R71NTi93EX19Yg4NyIevVhARBwVEWcD3wdO77l1SZKkAejlFNWD\nKF7F8JmImAcup3j31B7gKOAk4ESKVzb8aWZ+fkB9lSRJ6smyCU5mbgVeHhH/D/B7wK9TPPjvMOA2\nindTfT4zr1z8KKPhRcb9ie93/TpteSHmYDm3ndu9mqSx7taec3twJvEi493A+eXHWPMiY0mSVree\nE5xJ4kXGvfFiNS/E7LU95/bgOLed21VM0niP+iLjWglORKyheBLxHuD7mXl7neNJkiT1Q90VnE9Q\nXHD8FOBnEbEO+G5mnlq7Z5IkSStUN8H5ucw8LSIemZkPjYgXAPfoR8f6wYuM+xPf7/p12vJCzMFy\nbju3ezVJY92tPef24IzLRca9vk18MbvLz3sjYm1mvhN4bM1j1ubbxCVJWt3qruC8PSLuCvwT8O6I\n+DZwXO1e1eRFxr3xYjUvxOy1Pef24Di3ndtVTNJ4T9xFxhHxUIoLi+8ELs/M/wTeFBHPBR6MTzKW\nJEkjVinBiYjXULxU86fALHBMRPwH8NzM3Nz/7kmSJFVX9RqclwKvysz7ZOb9gKOBjwJfiohT+t47\nSZKkFah6iupQiuttAMjMW4HXRsQO4M3Aw/vYt9q8i6o/8f2uX6ct7zQZLOe2c7tXkzTW3dpzbg/O\npN5FdQXwmC7ln6R46eZY8C4qSZJWt6orOC8BLizfKv6hzJwvyx8D/Htfe1aDd1H1xqvxvdOk1/ac\n24Pj3HZuVzFJ4z1Rd1Fl5qURcQawCXhLRHyH4rTVA4FnVjmWJEnSoFR+0F9mfokioXkG8HVgK8Ut\n45+NiNsj4isR8db+dlOSJKl3K3rQX2bOAl8sPwCIiKOBhwIPKz8kSZJGou6TjPfJzJsoXrz5uX4d\nU5IkaSXqvotKkiRp7PRtBWcc+Ryc/sT3u36dtnxWyGA5t53bvZqkse7WnnN7cCb1OTgTwefgSJK0\nujVyBcfn4PTG5yn4rJBe23NuD45z27ldxSSN96ifg9PIFRxJkrS6meBIkqTGMcGRJEmNY4IjSZIa\nxwRHkiQ1jgmOJElqnEbeJt7ig/76E9/v+nXa8mFog+Xcdm73apLGult7zu3B8UF/A+SD/iRJWt0a\nuYLjg/564wOjfBhar+05twfHue3crmKSxtsH/UmSJPWZCY4kSWocExxJktQ4JjiSJKlxTHAkSVLj\nmOBIkqTGMcGRJEmNY4IjSZIaxwRHkiQ1jgmOJElqHBMcSZLUOI18F1WLbxPvT3y/69dpyzcuD5Zz\n27ndq0ka627tObcHx7eJD5BvE5ckaXVr5AqObxPvjW+l9Y3Lvbbn3B4c57Zzu4pJGm/fJi5JktRn\nJjiSJKlxTHAkSVLjmOBIkqTGMcGRJEmNY4IjSZIaxwRHkiQ1jgmOJElqHBMcSZLUOCY4kiSpcUxw\nJElS45jgSJKkxjHBkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS1DgmOJIkqXEmJsGJiJ+PiHMj\n4vxR90WSJI23kSY4EXFeRNwSEVd2lJ8aET+IiKsj4hyAzLwmM88cTU8lSdIkGfUKzmbg1PaCiJgG\n3gk8ATgJOCMiThp+1yRJ0qSaGWXjmXlJRBzXUfxI4OrMvAYgIj4KnA58v+rx5+fn2bp1a91uHmT7\n9u1LblcpW6p8uX1VYurE97t+nbbqjnWdfVViVhI7iPp123NuD45z27ldxSSN97Dn9mJGmuAs4hjg\nurbt64FHRcTdgDcAD4uIV2bmX3erHBFnAWeVmzuOOuqoHwyon3cHbmvbvguwrSOm17KlypfbVyWm\nTny/61fR77Gus69KzEpiB1G/Kue2c9u5vXzZUuXL7asSUye+3/WrGOTcvl9PPcjMkX4AxwFXtm0/\nDXhf2/YfAm8fdT+79Puyju1NXWJ6KluqfLl9VWLqxPe7/ijHetjjPUljPYjxdm4Pb6yHPd6TNNaD\nGG/n9vDGeiX9H/U1ON1cD9y3bftY4MYR9aWKT9coW6p8uX1VYurE97t+v9uuMtZ19lWJWUnsIOrX\n5dweHuf2cDm3h2fQc/sgUWZFI1Neg3NhZj6k3J4BfgicAtwAfAt4ZmZ+b1R97CYiLsvMjaPux2rg\nWA+X4z08jvVwOd7DMw5jPerbxD8CfB14YERcHxFnZuYc8ELg88AW4GPjltyUNo26A6uIYz1cjvfw\nONbD5XgPz8jHeuQrOJIkSf02jtfgSJIk1WKCI0mSGscER5IkNY4JTh9ExIkR8e6IOD8izh51f5ou\nIp4cEe+NiE9FxG+Puj9N5ktuhysi1kfEB8r5/axR92c1ca4Pz7B+h5vgLKLii0C3ZObzgd8HvAVx\nBSqO9ycz84+B5wJPH0F3J5ovuR2uKuMNPAU4v5zfTxp6ZydUxTHuyrnemz6N9VB+h5vgLG4zFV4E\nGhFPAr4KXDTcbjbGZqq/ePU15X5VsxlfcjtMm+l9vI9l/6tq5ofYx0m3mR7HOCJ+MSIu7Pi45/C7\nPLE207+xHujv8HF8F9VYyIovAs3MC4ALIuIzwIeH2dcmqDLeEbEFeCPwucz89lA72gBV5/Zwe9c8\nFcf7eook5wr8A7RnVcY4i/cYPnG4PWyOfox1RARD+B3u/0DVdHsR6DERcXJEvC0i3gN8djRda6Su\n4w28CHgc8NSIeP4oOtZAi83tu0XEuylfcjuarjXSYnP7E8B/i4h3MfrXGEy6xca4K+d6LZXGmiH9\nDncFp5roUpaZeTFw8XC7siosNt5vA9427M403GJjfTtgEtl/i433TuCPht2Zhuo6xosFO9drqTrW\nQ/kd7gpONZP6ItBJ5XgPj2M9XI734DnGwzOWY22CU823gBMi4viIWAs8A7hgxH1qMsd7eBzr4XK8\nB88xHp6xHGsTnEXEZL8IdOI43sPjWA+X4z14jvHwTNJY+7JNSZLUOK7gSJKkxjHBkSRJjWOCI0mS\nGscER5IkNY4JjiRJahwTHEmS1DgmONIqFhFTEfGeiLg9IjIiTh51n8ZZRKyJiB9GxG/0+bi/GBE3\nRMT6fh5XWs1McKTV7Xcp3n10GnA08LXRdmfsnQXckJmXtArKxPCpnYER8Y6IuLiXg2bmd4FvAC/t\nV0el1c4ER1rd7g/clJlfy8ybM3NvZ0D56HUVXgScO6Bjvx84OyJ8CbLUByY40ioVEZuBtwA/V65C\n/Lgsvzgi3hURfxsRtwL/WpbfJSI2RcQtEbE9Iv45IjZ2HPPZEXFtROyKiAsj4gURkW37XxcRV3bU\neW5E7OgoOy0iLo+IPRHxo4h4Q3uiFRE/jojXlKfX7igfGf/yjmMcWX4fN5XH2RIRT4+I9WWdp3bE\nPz4iZiPiXouM10bgAcCFPQ5xe93jyjHu/PhxW9gXgLsCJ1c9vqSDmeBIq9eLgddTvAn4aOARbfv+\nAAjgMcCzIyKAzwDHAE8EHgZcAnw5Io4GiIhHAZuBTcBDgU+Xx68kIn4H+BDwDuDBwPOApwJ/1RH6\nEuC7wC8DbwLeHBGPLo8RwOeAx1KcgjuJ4vTP3szcCXykPG675wEXZuZPF+naY4CrM3Nr1e8JuI5i\njFsfDwCuBS5uBZSrZ1eUfZZUk0uh0iqVmdsiYjswn5k3d+z+UWa+rLUREb9FkbTcIzN3l8V/ERGn\nAX8IvJkiYbooM99Q7v9hRDwCOLNi114N/E1mvr/c/o+IeAXwwYh4ee5/gd4XMvMd5ddvj4g/AU6h\neBHg44BHAw/OzC1lzDVtbbwX+EZEHJOZN0TEUcCTgact0a/7ATctsu8fyhWxdmspr2nKzHngZigu\n7AbeV24/v6POjcBxS/RBUo9cwZHUzeUd2w8H1gG3RsSO1gfwEOAXypgTKZKLdp3bvXg48OqOdj4M\nrAfu3Rb3nY56NwL3LL9+GMW1RVvoIjMvo1j9eU5Z9EzgZxSrPos5DNizyL6XUySA7R//uEjsm4Bf\nAp6cmZ3H2122I6kmV3AkdbOzY3sK+CnFaZpOd5Sfo4fjLnSJW9Olrb8EPt6l/q1tX8927Ev2/9HW\nS1/eB/wpxamv5wGby5WWxdxGkTh1c3NmXt1eEBHbgPt2lD2HYtXm17usmkFxDc6Pe+i7pGWY4Ejq\nxbeBewELmXnNIjHfB36lo6xz+1bgXhERbaeaHtqlrQd1Jgwr6O/REXHiYqs4wAeBv4mIF1Jcx/OM\nZY75b8ALI2IqMxeqdigifhV4F3BGZv7fRcIeAnyi6rElHcwER1IvvkRxN9WnIuLPgX+nOF10KvCl\nzPwX4G3A1yLilcD5FHcD/deO41xMsUrxqoj4aBnT+QyZ1wMXRsS1wMeAOYp/+B+ZmX/eY38vAi4F\n/ikiXgL8kOKW+PWZ+UnYdw3Sx4G/Ay7JzKuWOeZXgEMpTi9d0WM/AIiIewP/G/h74NJyG4rrn24t\nY46juIj7C1WOLak7r8GRtKxyteV3gS9TXKD7A4rk44EU176Qmd+guKD4bIrrY54CvK7jOFvK/WeV\nMY+n4+6ozPw88HvAbwLfLD/OAX5Sob8LwBMokrIPAluAt1Jc+Nvu3LJs2WfbZObtFKsrz+q1H20e\nRHF90MsoLlRufXyrLeYMigunr13B8SV1iP2rxJLUX+WzZj6emb1cEzN0EfF04D3AfTJzVw/xD6ZY\nybl/Zt6xXHyFfhwCXEVx+upf+3VcaTVzBUfSqhMR6yLiJOBVwHt7SW4AMvN7wJ8Bx/e5S/cD3mBy\nI/WPKziSBmZcV3Ai4nUUz9v5KnB6P1djJI0HExxJktQ4nqKSJEmNY4IjSZIaxwRHkiQ1jgmOJElq\nHBMcSZLUOCY4kiSpcUxwJElS45jgSJKkxjHBkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS1Dgz\no+7AIN397nfP4447ru/HnZ+fB2B6errrdpWypcqX21clpk58v+vXaavuWNfZVyVmJbGDqF+3Pef2\n4Di3ndvLxdSJ73f9Om31e25ffvnlt2XmPZbrR6MTnOOOO47LLrus78fdunUrABs2bOi6XaVsqfLl\n9lWJqRPf7/p12qo71nX2VYlZSewg6tdtz7k9OM5t5/ZyMXXi+12/Tlv9ntsRcW0v/fAUlSRJahwT\nHEmS1DgmOJIkqXFMcCRJUuOY4EiSpMYxwZEkSY1jgiNJkhrHBEeSJDWOCY4kSWocExxJktQ4JjiS\nJKlxTHAkSVLjjGWCExEviYjvRcSVEfGRiDg0Io6PiEsj4qqI+MeIWDvqfkqSpPE0dglORBwD/Amw\nMTMfAkwDzwDeBLwlM08AfgacObpeSpKkcTYz6g4sYgY4LCJmgXXATcBvAc8s938AeB3wrqUOMj8/\nv+8V6/20ffv2JberlC1Vvty+KjF14vtdv05bdce6zr4qMSuJHUT9uu05twfHue3crmKSxnvYc3sx\nY7eCk5k3AH8L/IQisdkGXA5szcy5Mux64Jhu9SPirIi4LCIuu+2224bRZUmSNGbGbgUnIo4CTgeO\nB7YCHwee0CU0u9XPzE3AJoCNGzfmhg0bBtRT6Dx2t7Z6LVuqfLl9VWLqxPe7fp226o51nX1VYlYS\nO4j6ddtzbg+Oc9u5XcUkjfew53ansVvBAR4H/Cgzb83MWeATwK8CGyKilZAdC9w4qg5KkqTxNo4J\nzk+AX4mIdRERwCnA94GvAE8tY54DfGpE/ZMkSWNu7BKczLwUOB/4NvBdij5uAl4BvDQirgbuBpw7\nsk5KkqSxNnbX4ABk5muB13YUXwM8cgTdkSRJE2bsVnAkSZLqMsGRJEmNY4IjSZIaxwRHkiQ1jgmO\nJElqHBMcSZLUOCY4kiSpcUxwJElS45jgSJKkxjHBkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS\n1DgmOJIkqXFMcCRJUuOY4EiSpMYxwZEkSY1jgiNJkhrHBEeSJDWOCY4kSWocExxJktQ4JjiSJKlx\nTHAkSVLjmOBIkqTGMcGRJEmNM5YJTkRsiIjzI+LfI2JLRDw6Iu4aEV+MiKvKz0eNup+SJGk8jWWC\nA7wV+D+Z+SDgvwBbgHOAizLzBOCicluSJOkgM6PuQKeIOBL4DeC5AJm5F9gbEacDJ5dhHwAuBl6x\n1LHm5+fZunVr3/u4ffv2JberlC1Vvty+KjF14vtdv05bdce6zr4qMSuJHUT9uu05twfHue3crmKS\nxnvYc3sx47iC8/PArcD7I+LfIuJ9EbEeuFdm3gRQfr5nt8oRcVZEXBYRl912223D67UkSRobY7eC\nQ9GnXwZelJmXRsRbqXA6KjM3AZsANm7cmBs2bBhML4HOY3drq9eypcqX21clpk58v+vXaavuWNfZ\nVyVmJbGDqF+3Pef24Di3ndtVTNJ4D3tudxrHFZzrgesz89Jy+3yKhOenEXE0QPn5lhH1T5Ikjbmx\nS3Ay82aJFHxQAAAWVElEQVTguoh4YFl0CvB94ALgOWXZc4BPjaB7kiRpAozjKSqAFwEfioi1wDXA\nH1EkYx+LiDOBnwBPG2H/JEnSGBvLBCczrwA2dtl1yrD7IkmSJs/YnaKSJEmqywRHkiQ1jgmOJElq\nHBMcSZLUOCY4kiSpcUxwJElS4/TtNvGIOAS4D3AYcGtm3tqvY0uSJFVRawUnIo6IiLMj4hJgG3A1\ncCVwc0RcFxHvjYhH9KOjkiRJvVpxghMRLwF+DDwP+CJwOvBQ4AHAo4HXUqwQfTEi/k9EnFC7t5Ik\nST2oc4rqV4HHZuaVi+z/JnBeRDwfOBN4LHBVjfYkSZJ6suIEJzN7ehdUZt4J/P1K25EkSarKu6gk\nSVLjmOBIkqTGMcGRJEmNY4IjSZIaZ8UXGUfEGuDZ5eY/ZObe/nRJkiSpnjorOH8L/JTiAX9/05/u\nSJIk1VfnOThTFK9lmAKm+9Od8XfTtt287lM/YO30FEesP5S101MwP8va6eAuR9zGITNTrJ2ZYn7v\nnaydCY46cve+stk9u1g7M8Xddk9xyMwUh6yZZu30FLv3zDEzFRw6O8+a6Smmp2LU36YkSROtToLz\nMuCZQAB/1p/ujL//3LmXz2+5baBtRMCaqSLRmZ6CNVPBzMx08Xl6ipnpYGYqmJmaYs10QC4UCdIh\na5iemmI6YHoqiAimI5ieCqamgumAqalgfnaWqangsEMP2b8/iraKuNhXfypgqvwcEUTAnXv2MBXB\n+nX/CezfPzUVBJT19tel/Bzl91YURfl1W3mrbF9csGvnTiLg8PV3EgE7W9uH76WVBu7auRMCjjh8\nthg/gp07dwBwxOFztAJbx9+xYwcBHPGzhX3jXeyH7a192xbYX7o/BmDnjvLY29sjOCgOirYAjtgR\nB+3r9nM/YJtg+/adRf1dU4vGLVfeOlYv9rW3Z7rrdlG2C4Aj79xfdscdB5a1WrujSyxt+wI4cm/7\nvgP7uaOsf8TszEHfQXR8wwHcsX03AWybX7Pve+42Lvt+5uUXrZDtd9wJwJ7Yc2Bcl+51Hr+1q/OY\n0aUOAXfsmSOAqT3FvN1x5xwAM+XnAHbtnQdg7d65ffV3zxZlh+ydP+h7u3OumLd7ypj2/bPzxb69\ncwsH7euMaX3u+JYPMLeQxee22Cq61c8l4nORnblIrc74zrHp/AwcMLYHlOX+n0FL62ez5s795a0m\nd7b9LLNLx5Pi55+Z5Jq9JBRfl/1OkvK/fdvb7riTTNiRu/Z9b63Prd+rU22/e+nYbs23mTvLeVf+\nPm//3Po93xR1HvS3F9jcv65MhmM2HMZfP+kB3Dm3wMzaQ7lzboFt23eyd36BmFnL3rkF7pybZ/vO\nPeydXyBjhjvn5tk7v8DOPXvZO5fMZxRlcwvsmVtg79w8cwvJwgLMLiyQCXvnF2C+veW5xbokSVJt\n01PBhsPWcJd1a9hw2BqOWre2/HotR61bw4Z1a7jLurX79hXbazjikJmxTIz69jbx1WLDurU84aR7\nFF9v2ADA1q1bD9iuUtatfH4hmZ1fYH4huf1nW5mbX2D9EUcyO7/A3Hwyt7DA3EIyN1/Ebb1jO3Pz\nyaHr1pWJUjK/kCwkLGTr6+Lz/EKyY+cu5jM55NDD2mLL/dmqX9Rt/WWxkMXxMmH3nj1kJmvWHlLE\nZCumo055jIVs/+uk/a+Ucrv8CyWTA78G9u7dSybMrFlDZrJ3dg5IpmfWQFl/dnaOJFkzs2Zf3bm5\nOTJhema6jNv/19Xc7BwJzMxM7/9LqNw3O1f+5TU9va+s/Q+wBObnyvptMZ1xRWwyP1dkqVPTB65g\nHBzbsV0GLMx3r99psb9uW/3o1Xz51/T09FTX7UyYXyj6ND21v0/tZe2tLdX/+da+qfLYB3d833Gn\npjrHLztDD+hvlKd5O//S7XaM9p/zQhZ/YExNxUFz48Bj5IHHX+KYi7dXfBVE2/8buW+bfWWtmAMb\nO2jOtDcYcUBA6/+Lfbt6nBJLhrV9Pyv6x629Q22WOtKiq5eL1Wov7mhvsXbioHqtagfWKH5W3VcS\nWz+r/at83Vor6k9NTbWtXO9f0aZjRTuziJ+emtq/Shj7v7XWXGn/Xd363Z0U/2+0/+5txZGtr4t/\nA27fuZfbd1a7Z6gzMTp8TXDkYTO8+PEncsK9jqh0rH4ywRlDxamp4hf67KHFj2jDkYcuGr91azHL\nOxOnxeO7J1q9qlu/Tlt1k8k6+6rErCR2EPXrttfv5L3XfVVi6sT3u36dtpzbg+XcXr7+7PwCW3fN\nsm33XrbumuVnu2bZumsv23bP8rNdRdnW3UXZ1l2z5cdedu6d75oYPe83RvuO7Tq3iR+fmT/qMTaA\nYzPzupW2J0mSBmfN9BT3OOIQ7nHEIZXq7Z1bYFsr8dk9yw23/oytu+c4/m7rB9TT3tRZwfl6RHwG\neF9mfr1bQEQcBTwD+BPgncA7arQnSZLGzNqZAxOjEzYUp503rF87ym7VSnAeBLwa+ExEzAOXAzcB\ne4CjgJOAE4FvAn+amZ+v2VdJkqSerPhBf5m5NTNfDhwDnA38O7ABOJ7ilp8PAA/LzF9bSXITEdMR\n8W8RcWG5fXxEXBoRV0XEP0bEaFNDSZI0tmpfZJyZu4Hzy49+ejGwBTiy3H4T8JbM/GhEvBs4E3hX\nn9uUJEkN0LeXbUbEb0XEayPilRHxmJrHOhb4PeB95XYAv8X+JOoDwJPrtCFJkpqrL7eJR8RfAH8J\nXEdxeur1EXED8OzMvGQFh/yfwJ8DrRvo7wZszczW0+6upzg1tqT5+fl9t8b10/bt25fcrlK2VPly\n+6rE1Invd/06bdUd6zr7qsSsJHYQ9eu259weHOe2c7uKSRrvYc/txax4BSciXhIRJ5d3Sr0YeEFm\n3i8zfwG4F/Bu4MKIeFzF4z4RuCUzL28v7hLa9RlUEXFWRFwWEZfddttgX6kgSZLGU50VnGcBfwW0\nLvZ9YkTcDfg2cEVmvjEibgXeCGyscNxfA54UEb8LHEpxDc7/BDZExEy5inMscGO3ypm5CdgEsHHj\nxhzkQ6Q6j92trV7Llipfbl+VmDrx/a5fp626Y11nX5WYlcQOon7d9pzbg+Pcdm5XMUnjPey53anO\nXVQbgcOB/wLsBf4TeBzwIeC6iLgFeB7wSxHxtIg4MSKWbS8zX5mZx2bmcRTP0PlyZj4L+Arw1DLs\nOcCnVtp3SZLUbLUuMs7M+cy8EvgX4D8y8+TMPAp4AMWt45dTrBK9FfgesLNGc68AXhoRV1Nck3Nu\nnb5LkqTm6te7qP4M+OeI+AWKa2+uAL4M/AZwQ2beNyLuDjykykEz82Lg4vLra4BH9qm/kiSpwfqS\n4GTmdyLil4G/p0hIWitDc8AflTG3lfskSZIGqm9vEy9fvPmEcqXmV4BDgEsz8/p+tSFJktSLviU4\nLeVKzYX9Pq4kSVKv+vYkY0mSpHFhgiNJkhrHBEeSJDWOCY4kSWocExxJktQ4JjiSJKlxTHAkSVLj\nmOBIkqTGMcGRJEmNY4IjSZIaxwRHkiQ1jgmOJElqHBMcSZLUOCY4kiSpcUxwJElS45jgSJKkxjHB\nkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS1DgmOJIkqXFMcCRJUuOMXYITEfeNiK9ExJaI+F5E\nvLgsv2tEfDEirio/HzXqvkqSpPE0dgkOMAe8LDNPBH4FeEFEnAScA1yUmScAF5XbkiRJB5kZdQc6\nZeZNwE3l19sjYgtwDHA6cHIZ9gHgYuAVSx1rfn6erVu39r2P27dvX3K7StlS5cvtqxJTJ77f9eu0\nVXes6+yrErOS2EHUr9uec3twnNvO7SomabyHPbcXM44rOPtExHHAw4BLgXuVyU8rCbrnInXOiojL\nIuKy2267bVhdlSRJY2TsVnBaIuJw4J+AP83MOyKip3qZuQnYBLBx48bcsGHDwPrYeexubfVatlT5\ncvuqxNSJ73f9Om3VHes6+6rErCR2EPXrtufcHhzntnO7ikka72HP7U5juYITEWsokpsPZeYnyuKf\nRsTR5f6jgVtG1T9JkjTexi7BiWKp5lxgS2b+f227LgCeU379HOBTw+6bJEmaDON4iurXgD8EvhsR\nV5RlrwLeCHwsIs4EfgI8bUT9kyRJY27sEpzM/Cqw2AU3pwyzL5IkaTKN3SkqSZKkukxwJElS45jg\nSJKkxjHBkSRJjWOCI0mSGscER5IkNY4JjiRJahwTHEmS1DgmOJIkqXFMcCRJUuOY4EiSpMYxwZEk\nSY1jgiNJkhrHBEeSJDWOCY4kSWocExxJktQ4JjiSJKlxTHAkSVLjmOBIkqTGMcGRJEmNY4IjSZIa\nxwRHkiQ1jgmOJElqHBMcSZLUOCY4kiSpcUxwJElS40xUghMRp0bEDyLi6og4Z9T9kSRJ42liEpyI\nmAbeCTwBOAk4IyJOGm2vJEnSOJoZdQcqeCRwdWZeAxARHwVOB76/WIX5+Xm2bt3a945s3759ye0q\nZUuVL7evSkyd+H7Xr9NW3bGus69KzEpiB1G/bnvO7cFxbju3q5ik8R723F7MJCU4xwDXtW1fDzyq\nMygizgLOKjd3HHXUUT8YUH/uDtzWtn0XYFtHTK9lS5Uvt69KTJ34ftevot9jXWdflZiVxA6iflXO\nbee2c3v5sqXKl9tXJaZOfL/rVzHIuX2/nnqQmRPxATwNeF/b9h8Cbx9hfy7r2N7UJaansqXKl9tX\nJaZOfL/rj3Kshz3ekzTWgxhv5/bwxnrY4z1JYz2I8XZuD2+sV9L/ibkGh2LF5r5t28cCN46oL918\nukbZUuXL7asSUye+3/X73XaVsa6zr0rMSmIHUb8u5/bwOLeHy7k9PIOe2weJMisaexExA/wQOAW4\nAfgW8MzM/N6I+nNZZm4cRdurjWM9XI738DjWw+V4D884jPXEXIOTmXMR8ULg88A0cN6okpvSphG2\nvdo41sPleA+PYz1cjvfwjHysJ2YFR5IkqVeTdA2OJElST0xwJElS45jgSJKkxjHBkSRJjWOC0wcR\ncWJEvDsizo+Is0fdn6aLiCdHxHsj4lMR8duj7k+TRcTPR8S5EXH+qPuyGkTE+oj4QDm/nzXq/qwm\nzvXhGdbvcBOcRUTEeRFxS0Rc2VF+0BvNM3NLZj4f+H3AZyysQMXx/mRm/jHwXODpI+juRKs41tdk\n5pmj6WkzVBlv4CnA+eX8ftLQOzuhKo5xV8713vRprIfyO9wEZ3GbgVPbC5Z6o3lEPAn4KnDRcLvZ\nGJupMN6l15T7Vc1mqo+1Vm4zvY/3sex/5978EPs46TbT4xhHxC9GxIUdH/ccfpcn1mb6N9YD/R0+\nMQ/6G7bMvCQijusoXvSN5pl5AXBBRHwG+PAw+9oEVcY7IrYAbwQ+l5nfHmpHG6Dq3B5u75qn4nhf\nT5HkXIF/gPasyhhn5l8DTxxuD5ujH2MdEcEQfof7P1A13d5ofkxEnBwRb4uI9wCfHU3XGqnreAMv\nAh4HPDUinj+KjjXQYnP7bhHxbuBhEfHK0XStkRab258A/ltEvIvRv6dp0i02xl0512upNNYM6Xe4\nKzjVRJeyzMyLgYuH25VVYbHxfhvwtmF3puEWG+vbAZPI/ltsvHcCfzTszjRU1zFeLNi5XkvVsR7K\n73BXcKoZ9zeaN43jPTyO9XA53oPnGA/PWI61CU413wJOiIjjI2It8AzgghH3qckc7+FxrIfL8R48\nx3h4xnKsTXAWEREfAb4OPDAiro+IMzNzDmi90XwL8LERv9G8MRzv4XGsh8vxHjzHeHgmaax9m7gk\nSWocV3AkSVLjmOBIkqTGMcGRJEmNY4IjSZIaxwRHkiQ1jgmOJElqHBMcaRWLiKmIeE9E3B4RGREn\nj7pP4ywi1kTEDyPiN/p83F+MiBsiYn0/jyutZiY40ur2uxTvPjoNOBr42mi7M/bOAm7IzEtaBWVi\n+NTOwIh4R0Rc3MtBM/O7wDeAl/aro9JqZ4IjrW73B27KzK9l5s2ZubczoHz0ugovAs4d0LHfD5wd\nEb4EWeoDExxplYqIzcBbgJ8rVyF+XJZfHBHvioi/jYhbgX8ty+8SEZsi4paI2B4R/xwRGzuO+eyI\nuDYidkXEhRHxgojItv2vi4grO+o8NyJ2dJSdFhGXR8SeiPhRRLyhPdGKiB9HxGvK02t3lI+Mf3nH\nMY4sv4+byuNsiYinR8T6ss5TO+IfHxGzEXGvRcZrI/AA4MIeh7i97nHlGHd+/Lgt7AvAXYGTqx5f\n0sFMcKTV68XA6yneBHw08Ii2fX8ABPAY4NkREcBngGOAJwIPAy4BvhwRRwNExKOAzcAm4KHAp8vj\nVxIRvwN8CHgH8GDgecBTgb/qCH0J8F3gl4E3AW+OiEeXxwjgc8BjKU7BnURx+mdvZu4EPlIet93z\ngAsz86eLdO0xwNWZubXq9wRcRzHGrY8HANcCF7cCytWzK8o+S6rJpVBplcrMbRGxHZjPzJs7dv8o\nM1/W2oiI36JIWu6RmbvL4r+IiNOAPwTeTJEwXZSZbyj3/zAiHgGcWbFrrwb+JjPfX27/R0S8Avhg\nRLw8979A7wuZ+Y7y67dHxJ8Ap1C8CPBxwKOBB2fmljLmmrY23gt8IyKOycwbIuIo4MnA05bo1/2A\nmxbZ9w/lili7tZTXNGXmPHAzFBd2A+8rt5/fUedG4Lgl+iCpR67gSOrm8o7thwPrgFsjYkfrA3gI\n8AtlzIkUyUW7zu1ePBx4dUc7HwbWA/dui/tOR70bgXuWXz+M4tqiLXSRmZdRrP48pyx6JvAzilWf\nxRwG7Flk38spEsD2j39cJPZNwC8BT87MzuPtLtuRVJMrOJK62dmxPQX8lOI0Tac7ys/Rw3EXusSt\n6dLWXwIf71L/1ravZzv2Jfv/aOulL+8D/pTi1NfzgM3lSstibqNInLq5OTOvbi+IiG3AfTvKnkOx\navPrXVbNoLgG58c99F3SMkxwJPXi28C9gIXMvGaRmO8Dv9JR1rl9K3CviIi2U00P7dLWgzoThhX0\n9+iIOHGxVRzgg8DfRMQLKa7jecYyx/w34IURMZWZC1U7FBG/CrwLOCMz/+8iYQ8BPlH12JIOZoIj\nqRdforib6lMR8efAv1OcLjoV+FJm/gvwNuBrEfFK4HyKu4H+a8dxLqZYpXhVRHy0jOl8hszrgQsj\n4lrgY8AcxT/8j8zMP++xvxcBlwL/FBEvAX5IcUv8+sz8JOy7BunjwN8Bl2TmVcsc8yvAoRSnl67o\nsR8ARMS9gf8N/D1wabkNxfVPt5Yxx1FcxP2FKseW1J3X4EhaVrna8rvAlyku0P0BRfLxQIprX8jM\nb1BcUHw2xfUxTwFe13GcLeX+s8qYx9Nxd1Rmfh74PeA3gW+WH+cAP6nQ3wXgCRRJ2QeBLcBbKS78\nbXduWbbss20y83aK1ZVn9dqPNg+iuD7oZRQXKrc+vtUWcwbFhdPXruD4kjrE/lViSeqv8lkzH8/M\nXq6JGbqIeDrwHuA+mbmrh/gHU6zk3D8z71guvkI/DgGuojh99a/9Oq60mrmCI2nViYh1EXES8Crg\nvb0kNwCZ+T3gz4Dj+9yl+wFvMLmR+scVHEkDM64rOBHxOorn7XwVOL2fqzGSxoMJjiRJahxPUUmS\npMYxwZEkSY1jgiNJkhrHBEeSJDWOCY4kSWqc/x+WvcbSu4YIsQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11799b6d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2, 1, figsize=(8, 3*2))\n",
    "\n",
    "# plot apparent resistivity\n",
    "ax[0].loglog(frequencies, app_res, lw=2)\n",
    "ax[0].set_ylim(1./sigma_halfspace*np.r_[0.1, 10])\n",
    "ax[0].set_ylabel(\"$ \\\\rho_a (\\Omega$-m)\", fontsize=14)\n",
    "\n",
    "# plot phase\n",
    "ax[1].semilogx(frequencies, phase, lw=2)\n",
    "ax[1].set_ylim(np.r_[0., 90.])\n",
    "ax[1].grid(which=\"both\", linewidth=0.4)\n",
    "\n",
    "ax[1].set_xlabel(\"frequency (Hz)\", fontsize=14)\n",
    "ax[1].set_ylabel(\"$\\phi (^\\circ)$\", fontsize=14)\n",
    "\n",
    "for a in ax:\n",
    "    a.invert_xaxis()  # highest frequencies see the near surface, lower frequencies see deeper\n",
    "    a.set_xlabel(\"frequency (Hz)\", fontsize=14)\n",
    "    a.grid(which=\"both\", linewidth=0.4)\n",
    "\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discretization, the Gory Details. \n",
    "\n",
    "If you want to skip this section, we won't judge! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To numerically solve Maxwell's equations, we need to first discretize them so they are represented on a mesh. We will take a finite difference approach for this example.\n",
    "\n",
    "Since we are solving for a 1D model, we will use a 1D mesh and leverage the Mesh class in SimPEG to build the operators (see http://docs.simpeg.xyz for docs). \n",
    "\n",
    "We show a very small mesh in the derivation so that it is meaningful to print out the matrices. When we go to solve, we will use a larger mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This mesh has 4 cells and 5 faces. Each cell is 1.0m wide\n"
     ]
    }
   ],
   "source": [
    "cell_size = 1. # width of the cell in meters\n",
    "ncells = 4  # number of cells that make up our domain\n",
    "\n",
    "# define a Tensor Mesh\n",
    "dz = [(cell_size, ncells)]\n",
    "mesh = Mesh.TensorMesh([dz], x0='N')\n",
    "\n",
    "print(\n",
    "    \"This mesh has {nC} cells and {nF} faces. \"\n",
    "    \"Each cell is {h}m wide\".format(\n",
    "        nC=mesh.nC, nF=mesh.nF, h=mesh.hx.min()  # it is hx because SimPEG treats dimensions in the order (x, y, z), so if the mesh is 1D, we work with the first component\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two places where we can discretize variables on a 1D mesh for the electromagnetic problem: cell centers and cell faces. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x117dff4d0>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAADTCAYAAACss0uSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGm5JREFUeJzt3XuUlPWd5/H3h7saDJfWLEMrjWccDyCtaOPKkCiuGxTd\nqBmzu/Fg4mUTNF72jImJIpuY6HA0mk2yHDUMM4voSoyOWaNZLyCKY2JwtYmiIBoQLzSwchOEY4zQ\n/d0/6gGLppou6Oqu+lV/XufU6arn+dWvfl9+TX/6ufTzKCIwMzOzytej3AMwMzOz4ji0zczMEuHQ\nNjMzS4RD28zMLBEObTMzs0Q4tM3MzBLh0DYzM0uEQ9vMzCwRDm0zM7NE9Cr3AFqrqamJurq6kvXX\n3NwMQM+ePUvWZzlUSx3gWipVtdRSLXWAa6lUnVHL4sWLN0bEYe21q7jQrquro7GxsWT9bdmyBYAB\nAwaUrM9yqJY6wLVUqmqppVrqANdSqTqjFknvFtPOu8fNzMwS4dA2MzNLhEPbzMwsERV3TNvMzNK0\nY8cOmpqa+Pjjj/da19LSAsC6deu6elgl15Fa+vXrR21tLb179z6gz3Zom5lZSTQ1NdG/f3/q6uqQ\ntMe6nTt3AtCrV/qxc6C1RASbNm2iqamJ4cOHH9Bne/e4mZmVxMcff8zgwYP3CmzLkcTgwYML7oko\nlkPbzMxKxoG9bx3993Fom5mZJcKhbWZm5TF3LtTVQY8eua9z55Z7REDu4il33XVXuYdRkEPbzMy6\n3ty5MGUKvPsuROS+TplSEcF9IKEdEbvPKu9MDm0zM+t606bBRx/tueyjj3LLO+jee++lvr6e4447\njq997Wts2LCB888/n7FjxzJ27Fief/55AH74wx9y6aWXMmHCBI466ihmzJgBwPXXX89bb73F8ccf\nz3e/+10Abr/9dsaOHUt9fT0/+tGPAHjnnXcYMWIEV1xxBSeccAKrV6/m4osv5thjj2X06NH87Gc/\n63AtraV/7r2ZmaXnvff2b3mRli1bxvTp03n++eepqalh8+bNXHXVVVxzzTV8/vOf57333uOMM85g\n+fLlALzxxhssXLiQbdu2ccwxx/Ctb32LW2+9laVLl/LKK68AMH/+fFasWMGLL75IRPClL32J3/3u\ndwwfPpw333yTu+++m7vuuovFixezZs0ali5dCnx6jfJScmibmVnXO/LI3C7xQss74JlnnuErX/kK\nNTU1AAwaNIgFCxbw+uuv727z4Ycfsm3bNgDOPvts+vbtS9++fTn88MN5//339+pz/vz5zJ8/nzFj\nxgCwfft2VqxYwfDhwxk2bBgnn3wyAEcddRSrVq3i6quv5uyzz2bixIkdqqUQh7aZmXW96dNzx7Dz\nd5EffHBueQdExF5/VtXS0sKiRYs46KCD9mrft2/f3c979uy5+8IprfucOnUql112GfDpxVWampo4\n5JBDdrcbOHAgS5YsYd68edx55508+OCDzJ49u0P1tOZj2mZm1vUmT4ZZs2DYMJByX2fNyi3vgNNP\nP50HH3yQTZs2AbB582YmTpzIHXfcsbvNrt3ebenfv//uLXGAM844g9mzZ7N9+3YA1qxZw/r16/d6\n38aNG2lpaeH888/n5ptv5o9//GOHainEW9pmZlYekyd3OKRbGzVqFNOmTePUU0+lZ8+ejBkzhhkz\nZnDllVdSX1/Pzp07OeWUU5g5c2abfQwePJjx48dz7LHHMmnSJG6//XaWL1/OuHHjADjkkEO45557\n9thKh1yYX3LJJbvPIr/llltKWhuAIqLknXZEQ0NDNDY2lqy/arnxerXUAa6lUlVLLdVSB6RXy/Ll\nyxkxYkTBdb72+KcK/TtJWhwRDe29t93d45JmS1ovaWkb6yVphqSVkl6VdEKr9YdKWiPpjkLvNzMz\ns+IUc0x7DnDmPtZPAo7OHlOAX7RafzPwrwcyODMzM/tUu9v2EfGcpLp9NDkXuDdy+9lfkDRA0pCI\nWCfpROBzwJNAu5v9AM3NzSX927b8kwlSVi11gGupVNVSS7XUAenV0tLSUvDs613rgDbXp6SjtbS0\ntBxwzpXi7PGhwOq8103AUEk9gP8OfLe9DiRNkdQoqXHjxo0lGJKZmVn1KcUZAYXuMxbAFcDjEbG6\nvVuRRcQsYBbkTkTrjJMuUjmRoz3VUge4lkpVLbVUSx2QTi3r1q1r8+Qsn4j2qR49ehzwnJbiX68J\nOCLvdS2wFhgHfEHSFcBngD6StkfE9SX4TDMzs26nFLvHHwW+np1FfjKwNSLWRcTkiDgyIuqAa8kd\n93Zgm5lZp5oxYwYjRoxgcon/BrwStLulLel+YAJQI6kJuBHoDRARM4HHgbOAlcBHwCWdNVgzM6su\nixbBs8/ChAmQXbukw+666y6eeOIJhg8fXpoOK0gxZ49f0M76AK5sp80ccn86ZmZm3cDf/z3kXy00\noieQu2LpLlu3wquvQksL9OgB9fXw2c+23efxx8PPf77vz7388stZtWoV55xzDhdeeCGPPPIIf/7z\nnznooIO4++67OeaYY2hubua6665j3rx5SOKb3/wmV199NYsXL+bb3/4227dvp6amhjlz5jBkyBBm\nzJjBzJkz6dWrFyNHjuS+++7rwL9Mx6R/RoCZmSVp69ZcYEPu69at+w7tYsycOZMnn3yShQsX0qdP\nH77zne/Qq1cvFixYwA033MCvf/1rZs2axdtvv83LL79Mr1692Lx5Mzt27ODqq6/mkUce4bDDDuOB\nBx5g2rRpzJ49m1tvvZW3336bvn37dsrtNveHQ9vMzEqu9Rbxzp3NwJ5nXC9aBKefDp98An36wNy5\npdtFDrB161YuuugiVqxYgSR27NgBwIIFC7j88st3j2XQoEEsXbqUpUuX8sUvfhHIXTNkyJAhANTX\n1zN58mTOO+88zjvvvNIN8AA4tM3MrCzGjYOnny79Me1dvv/973Paaafx8MMP88477zBhwgSg8O07\nI4JRo0axaNGivfp57LHHeO6553j00Ue5+eabWbJkSdn+dM235jQzs7IZNw6mTi19YENuS3vo0KEA\nzJkzZ/fyiRMnMnPmzN1/b71582aOOeYYNmzYsDu0d+zYwbJly2hpaWH16tWcdtpp3HbbbWzZsmX3\nLTrLwaFtZmZV6Xvf+x5Tp05l/PjxNDc3717+jW98gyOPPJL6+nqOO+44fvnLX9KnTx8eeughrrvu\nOo477jiOP/54/vCHP9Dc3MyFF17I6NGjGTNmDNdcc01ZL3bjW3MmolrqANdSqaqllmqpA9Krxbfm\nLE6n3prTzMzMKoND28zMLBEObTMzK5lKO+RaaTr67+PQNjOzkujXrx+bNm1ycLchIti0aRP9+vU7\n4D7SPyPAzMwqQm1tLU1NTWzYsGGvdS3Zpc969Eh/W7EjtfTr14/a2toD/myHtpmZlUTv3r3bvElH\namfC70s5a0n/Vx4zM7NuwqFtZmaWCIe2mZlZIhzaZmZmiXBom5mZJcKhbWZmlgiHtpmZWSIc2mZm\nZolwaJuZmSWi3dCWNFvSeklL21gvSTMkrZT0qqQTsuXHS1okaVm2/D+XevBmZmbdSTFb2nOAM/ex\nfhJwdPaYAvwiW/4R8PWIGJW9/+eS0r9+nZmZWZm0e+3xiHhOUt0+mpwL3Bu527q8IGmApCER8ae8\nPtZKWg8cBmzZ1+c1Nzfvvq5rKWzbtq1kfZVTtdQBrqVSVUst1VIHuJZKVc5aSnFMeyiwOu91U7Zs\nN0knAX2Atwp1IGmKpEZJjRs3bizBkMzMzKpPKe7ypQLLdt9MVdIQ4H8BF0VES6EOImIWMAugoaEh\nOuPOKdVwZxmonjrAtVSqaqmlWuoA11KpUr3LVxNwRN7rWmAtgKRDgceA/xYRL5Tgs8zMzLqtUoT2\no8DXs7PITwa2RsQ6SX2Ah8kd7/6XEnyOmZlZt9bu7nFJ9wMTgBpJTcCNQG+AiJgJPA6cBawkd8b4\nJdlb/xNwCjBY0sXZsosj4pUSjt/MzKzbKObs8QvaWR/AlQWW3wfcd+BDMzMzs3y+IpqZmVkiHNpm\nZmaJcGibmZklwqFtZmaWCIe2mZlZIhzaZmZmiXBom5mZJcKhbWZmlgiHtpmZWSIc2mZmZolwaJuZ\nmSXCoW1mZpYIh7aZmVkiHNpmZmaJcGibmZklwqFtZmaWCIe2mZlZIhzaZmZmiXBom5mZJcKhbWZm\nloh2Q1vSbEnrJS1tY70kzZC0UtKrkk7IW3eRpBXZ46JSDrxbmTuX/qNH89lBg6CuDubOLfeIDDwv\nlchzUpk8LyVTzJb2HODMfayfBBydPaYAvwCQNAi4Efi3wEnAjZIGdmSwB+LFF3vy05/2ZdGirv7k\nEpk7F6ZMoWdTE4qAd9+FKVOS/qZPfk7A81KJPCeVyfNSUoqI9htJdcD/iYhjC6z7R+DZiLg/e/0m\nMGHXIyIuK9SuLWPGjImFCxfuVxFtefHFnpx11mdoboYePWDUqGYOPbQkXXeZni+9hD75y17Lo09f\nmseOLcOIOubDD2HZsp60tKQ7J+B5qUSek8pUrfMSAf36wW9+s52TTmrucL8DBw5cHBEN7bUrxTHt\nocDqvNdN2bK2lu9F0hRJjZIaN27cWIIh5fz+971oaQEQLS3w4YcqWd9dpdA3+76WV7oPP1TycwKe\nl0rkOalM1TovEeKTT3I505VK8WmFvpNiH8v3XhgxC5gF0NDQEAMGDCjBsGDSJPjJT4JPPgn69BH3\n39+TceNK0nXXqbs4tzuptWHD4PfvdPVoOmzRIjj99MTnBDwvlchzUpmqfF4mTTqIAQMO6rLPL8WW\ndhNwRN7rWmDtPpZ3mXHjcrsubrjhY55+mjS/4adPh4MP3nPZwQfnlieoKuYEPC+VyHNSmTwvpRUR\n7T6AOmBpG+vOBp4gt2V9MvBitnwQ8DYwMHu8DQxq77NOPPHEKKUPPvggPvjgg5L22eXuuy921tZG\nixQxbFjEffeVe0QdUhVzEuF5qUSek8rkeWkX0BhF5HG7J6JJup/cSWU1wPvkzgjvnQX+TEkC7iB3\nhvlHwCUR0Zi991Lghqyr6RFxd3u/RDQ0NERjY2MRv24UZ8uWLQCUapd7uVRLHeBaKlW11FItdYBr\nqVSdUYukok5Ea/eYdkRc0M76AK5sY91sYHZ7n2FmZmbt8xXRzMzMEuHQNjMzS4RD28zMLBEObTMz\ns0Q4tM3MzBLh0DYzM0uEQ9vMzCwRDm0zM7NEOLTNzMwS4dA2MzNLhEPbzMwsEQ5tMzOzRDi0zczM\nEuHQNjMzS4RD28zMLBEObTMzs0Q4tM3MzBLh0DYzM0uEQ9vMzCwRDm0zM7NEOLTNzMwSUVRoSzpT\n0puSVkq6vsD6YZKelvSqpGcl1eatu03SMknLJc2QpFIWYGZm1l20G9qSegJ3ApOAkcAFkka2avYT\n4N6IqAduAm7J3vu3wHigHjgWGAucWrLRm5mZdSO9imhzErAyIlYBSPoVcC7wel6bkcA12fOFwG+y\n5wH0A/oAAnoD7+/rw5qbm9myZUux42/Xtm3bStZXOVVLHeBaKlW11FItdYBrqVTlrKWY3eNDgdV5\nr5uyZfmWAOdnz78M9Jc0OCIWkQvxddljXkQsb/0BkqZIapTUuHHjxv2twczMrFsoZku70DHoaPX6\nWuAOSRcDzwFrgJ2S/hoYAew6xv2UpFMi4rk9OouYBcwCaGhoiAEDBhRfQZE6o89yqJY6wLVUqmqp\npVrqANdSqcpRSzGh3QQckfe6Flib3yAi1gJ/ByDpM8D5EbFV0hTghYjYnq17AjiZXLCbmZnZfihm\n9/hLwNGShkvqA3wVeDS/gaQaSbv6mgrMzp6/B5wqqZek3uROQttr97iZmZm1r93QjoidwFXAPHKB\n+2BELJN0k6RzsmYTgDcl/Qn4HDA9W/4Q8BbwGrnj3ksi4relLcHMzKx7KGb3OBHxOPB4q2U/yHv+\nELmAbv2+ZuCyDo7RzMzM8BXRzMzMkuHQNjMzS4RD28zMLBEObTMzs0Q4tM3MzBLh0DYzM0uEQ9vM\nzCwRDm0zM7NEOLTNzMwS4dA2MzNLhEPbzMwsEQ5tMzOzRDi0zczMEuHQNjMzS4RD28zMLBEObTMz\ns0Q4tM3MzBLh0DYzM0uEQ9vMzCwRDm0zM7NEOLTNzMwSUVRoSzpT0puSVkq6vsD6YZKelvSqpGcl\n1eatO1LSfEnLJb0uqa50wzczM+s+2g1tST2BO4FJwEjgAkkjWzX7CXBvRNQDNwG35K27F7g9IkYA\nJwHrSzFwMzOz7qZXEW1OAlZGxCoASb8CzgVez2szErgme74Q+E3WdiTQKyKeAoiI7e19WHNzM1u2\nbCm6gPZs27atZH2VU7XUAa6lUlVLLdVSB7iWSlXOWorZPT4UWJ33uilblm8JcH72/MtAf0mDgb8B\ntkj635JelnR7tuW+B0lTJDVKaty4ceP+V2FmZtYNFLOlrQLLotXra4E7JF0MPAesAXZm/X8BGAO8\nBzwAXAz8zz06i5gFzAJoaGiIAQMGFF1AsTqjz3KoljrAtVSqaqmlWuoA11KpylFLMVvaTcARea9r\ngbX5DSJibUT8XUSMAaZly7Zm7305IlZFxE5yu81PKMnIzczMupliQvsl4GhJwyX1Ab4KPJrfQFKN\npF19TQVm5713oKTDstf/jj2PhZuZmVmR2g3tbAv5KmAesBx4MCKWSbpJ0jlZswnAm5L+BHwOmJ69\nt5ncrvOnJb1Gblf7P5W8CjMzs26gmGPaRMTjwOOtlv0g7/lDwENtvPcpoL4DYzQzMzN8RTQzM7Nk\nOLTNzMwS4dA2MzNLhEPbzMwsEQ5tMzOzRDi0zczMEuHQNjMzS4RD28zMLBEObTMzs0Q4tM3MzBLh\n0DYzM0uEQ9vMzCwRDm0zM7NEOLTNzMwS4dA2MzNLhEPbzMwsEYqIco9hD5I2AO+WuNsaYGOJ+yyH\naqkDXEulqpZaqqUOcC2VqtS1DIuIw9prVHGh3RkkNUZEQ7nH0VHVUge4lkpVLbVUSx3gWipVuWrx\n7nEzM7NEOLTNzMwS0V1Ce1a5B1Ai1VIHuJZKVS21VEsd4FoqVVlq6RbHtM3MzKpBd9nSNjMzS55D\n28zMLBFVE9qSzpT0pqSVkq4vsL6vpAey9f9XUl3Xj3L/SBok6SlJK7KvA9to1yzplezxaFePsxiS\n/qOkZZJaJLX5ZxLtzWMl2I9a3pH0WjYvjV05xmJJul3SG5JelfSwpAFttKvoedmPOlKYk5uzOl6R\nNF/SX7XR7qLsZ8MKSRd19TiLsR+1VPzPsF0kXSspJNW0sb5z5yUikn8APYG3gKOAPsASYGSrNlcA\nM7PnXwUeKPe4i6jrNuD67Pn1wI/baLe93GMtopYRwDHAs0DDgc5jJTyKqSVr9w5QU+7xtlPLRKBX\n9vzHhb7HUpiXYupIaE4OzXv+X3f93GrVZhCwKvs6MHs+sNxjP5BasnUV/zMsG+cRwDxyFwDb6/uo\nK+alWra0TwJWRsSqiPgE+BVwbqs25wL3ZM8fAk6XpC4c44HIH/M9wHllHEuHRMTyiHiznWbFzGPZ\nFVlLEiJifkTszF6+ANQWaFbx81JkHUmIiA/zXh4CFDpb+AzgqYjYHBEfAE8BZ3bF+PZHkbWk5GfA\n92i7jk6fl2oJ7aHA6rzXTdmygm2y/9xbgcFdMroD97mIWAeQfT28jXb9JDVKekFSssFOcfOYkgDm\nS1osaUq5B1OES4EnCixPbV7aqgMSmRNJ0yWtBiYDPyjQJJk5KaIWSOBnmKRzgDURsWQfzTp9XnqV\nsrMyKrTF3Po3oWLadDlJC4B/U2DVtP3o5siIWCvpKOAZSa9FxFulGWHx9lVLRDxSTBcFlpVljkpQ\nC8D4bF4OB56S9EZEPFe6URanmFokTQN2AnMLdVFgWZfPSwnqgETmJCKmAdMkTQWuAm5s3UWB91bk\n/5UiaoEEfoYBN5A7DLPPLgosK+m8VEtoN5E71rBLLbC2jTZNknoBnwU2d83w2hYR/76tdZLelzQk\nItZJGgKsb6OPtdnXVZKeBcaQOwbZpfZVS5GKmccuUYJa8udlvaSHye1m7vKAaK+W7GSZ/wCcHtmB\nuVYqYl5KUEcyc5Lnl8Bj7B10TcCEvNe15M6x6HIlqKXif4ZJGg0MB5ZkR1VrgT9KOiki/l9e006f\nl2rZPf4ScLSk4ZL6kDvRrPUZiI8Cu87k+wrwTFv/sStI/pgvAvbawpM0UFLf7HkNMB54vctGWFrF\nzGMSJB0iqf+u5+R+Q19a3lHtTdKZwHXAORHxURvNKn5eiqkjoTk5Ou/lOcAbBZrNAyZm//8Hkqtl\nXleMb38UU0sKP8Mi4rWIODwi6iKijlw4n9AqsKEr5qXcZ+OV6gGcBfyJ3G9n07JlN5H7TwzQD/gX\nYCXwInBUucdcRE2DgaeBFdnXQdnyBuCfs+d/C7xG7oze14D/Uu5xt1HLl7Nv9L8A7wPzsuV/BTy+\nr3mstEcxtZA703pJ9lhWwbWsJHcM7pXssesvLJKal2LqSGhOfk3ul4lXgd8CQ7Plu//fZ68vzepe\nCVxS7nEfaC2p/AxrVdc7ZGePd/W8+DKmZmZmiaiW3eNmZmZVz6FtZmaWCIe2mZlZIhzaZmZmiXBo\nm5mZJcKhbdbNSXoouxJVse1HS5rTiUMyszY4tM26MUmjgJ4RsarY90TEa0CtpCM7b2RmVohD26xK\nSbo87x7Fb0taWKDZZPKutCdpu6QfZzfUWCDpJEnPSlqV3TBhl9+SuzKamXUhh7ZZlYqImRFxPDCW\n3BXcflqg2Xhgcd7rQ4BnI+JEYBvwD8AXyV0F7qa8do3AFzpj3GbWtmq5YYiZte1/kLvW/m8LrBsC\nbMh7/QnwZPb8NeAvEbFD0mtAXV679eQuEWpmXcihbVbFJF0MDCN3S8RC/kzuuvy77IhPr23cQu76\n6kRES3Z3vF36Ze81sy7k3eNmVUrSicC1wIUR0dJGs+XAXx9A939DBd4hy6zaObTNqtdVwCBgYXYy\n2j8XaPMYe97/t1inZe81sy7ku3yZdWOSDgIWAuMjornI9/QF/hX4fETs7MzxmdmeHNpm3ZykM4Dl\nEfFeke2PJndf5Gc7dWBmtheHtpmZWSJ8TNvMzCwRDm0zM7NEOLTNzMwS4dA2MzNLhEPbzMwsEf8f\nnAOD2L/y/1MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1156cff10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,3))\n",
    "\n",
    "mesh.plotGrid(centers=True, faces=True, ax=ax)\n",
    "ax.invert_xaxis()  # put the surface of the earth on the left. \n",
    "ax.set_xlabel('z (m)')\n",
    "ax.grid(which=\"both\", linewidth=0.4)\n",
    "plt.legend((\"centers\", \"faces\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To count, we will use \n",
    "$$\n",
    "i = 0, 1, 2, ..., N\n",
    "$$\n",
    "to denote cell centers, so faces are at $i \\pm 1/2$.\n",
    "\n",
    "To discretize our system of equations, we put the physical properties, $\\sigma$, $\\mu$, $\\epsilon$ at cell centers and stagger the electric and magnetic fields so that $E_x$ is on cell centers and $H_y$ is at cell faces.\n",
    "\n",
    "Our physical properties are described by the discrete vectors\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\sigma} = [\\sigma_0, \\sigma_1, \\sigma_2, ..., \\sigma_N]^\\top\n",
    "$$\n",
    "$$\n",
    "\\boldsymbol{\\mu} = [\\mu_0, \\mu_1, \\mu_2, ..., \\mu_N]^\\top\n",
    "$$\n",
    "$$\n",
    "\\boldsymbol{\\epsilon} = [\\epsilon_0, \\epsilon_1, \\epsilon_2, ..., \\epsilon_N]^\\top\n",
    "$$\n",
    "\n",
    "and \n",
    "$$\n",
    "\\boldsymbol{\\hat{\\sigma}} = \\boldsymbol{\\sigma} + \\imath \\omega \\boldsymbol{\\epsilon}\n",
    "$$\n",
    " \n",
    "is also defined at cell centers. \n",
    "\n",
    "Our fields are described by the discrete vectors\n",
    "$$\n",
    "\\mathbf{e_x} = [e_0, e_1, e_2, ..., e_N]^\\top\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{h_y} = [h_{-1/2}, h_{1/2}, h_{1+1/2}, h_{2+1/2}, ..., h_{N+1/2}]^\\top\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretizing Ampere's law\n",
    "\n",
    "Lets start by exmining Ampere's law (the second equation), \n",
    "\n",
    "$$-\\frac{\\partial H_y}{\\partial z} +  \\hat{\\sigma} E_x = 0$$\n",
    "\n",
    "To approximate the derivative of $H_y$ (which is defined on faces) with respect to $z$, we use centered differences, so \n",
    "\n",
    "$$\n",
    "\\frac{\\partial H_y}{\\partial z} \\bigg\\rvert_i \\simeq \\frac{h_{i+1/2} - h_{i-1/2}}{\\Delta z_i}\n",
    "$$\n",
    "\n",
    "where $\\Delta z_i$ is the width of the cell, and the approximation of the derivative lands on the cell center. We repeat this operation for each cell in our mesh. You could do this in a for loop, but it is often benificial to work with the matrix form, so we will do that here. \n",
    "\n",
    "The differential operator matrix that takes the derivative of a variable defined on faces is the face divergence operator.\n",
    "\n",
    "$$\n",
    "\\frac{\\partial H_y}{\\partial z} \\simeq \\mathbf{Div} ~\\mathbf{h_y}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.  1.  0.  0.  0.]\n",
      " [ 0. -1.  1.  0.  0.]\n",
      " [ 0.  0. -1.  1.  0.]\n",
      " [ 0.  0.  0. -1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "Div = mesh.faceDiv\n",
    "print(Div.todense())  # operators are stored as sparse matrices in SimPEG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the physical properties $\\boldsymbol{\\hat{\\sigma}}$ is defined at cell centers which is in the same location as $E_x$, we can simply multiply them. In matrix form, we use a diagonal matrix, \n",
    "\n",
    "$$\\mathbf{M^{cc}_{\\boldsymbol{\\hat{\\sigma}}}} = \\mathbf{diag}(\\boldsymbol{\\hat{\\sigma}})$$ \n",
    "\n",
    "so the product is given by\n",
    "\n",
    "$$\n",
    "\\hat{\\sigma} E_x \\simeq \\mathbf{M^{cc}_{\\boldsymbol{\\hat{\\sigma}}}} ~\\mathbf{e_x}\n",
    "$$\n",
    "\n",
    "in the example that follows, we will assume that $\\sigma \\gg \\imath \\omega \\epsilon$, so $\\hat{\\sigma} \\simeq \\sigma$. In the more general implementation later on in the tutorial, we will use the full definition of $\\hat{\\sigma}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.1  0.   0.   0. ]\n",
      " [ 0.   0.1  0.   0. ]\n",
      " [ 0.   0.   0.1  0. ]\n",
      " [ 0.   0.   0.   0.1]]\n"
     ]
    }
   ],
   "source": [
    "sigma = 1e-1 * np.ones(mesh.nC)\n",
    "Mcc_sigma = Utils.sdiag(sigma)\n",
    "print(Mcc_sigma.todense())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have taken\n",
    "\n",
    "$$-\\frac{\\partial H_y}{\\partial z} +  \\hat{\\sigma} E_x = 0$$\n",
    "\n",
    "and discretized to \n",
    "\n",
    "$$\n",
    "- \\mathbf{Div} ~ \\mathbf{h_y} + \\mathbf{M^{cc}_{\\boldsymbol{\\hat{\\sigma}}}} ~ \\mathbf{e_x} = \\mathbf{0}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretizing Faraday's Law\n",
    "\n",
    "Next, we examine Faraday's law:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}  + \\imath \\omega \\mu H_y = 0\n",
    "$$\n",
    "\n",
    "Over one cell, the discrete approximation of the derivative of $E_x$ with respect to $z$ is\n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{i+1/2} = \\frac{e_{i+1} - e_{i}}{\\Delta z_{i+1/2}}\n",
    "$$\n",
    "\n",
    "where $\\Delta z_{1+1/2}$ is the distance (m) from the cell center $z_{i}$ to $z_{i+1}$. Notice here that we are going from cell centers to cell faces. So in this case we need to handle the boundary conditions, what do we do at \n",
    "the top and the bottom:\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{-1/2}, \\quad\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{nC+1/2} \n",
    "$$\n",
    "\n",
    "we somehow need to define \"ghost points\" $e_{-1}$ and $e_{N+1}$ so that we can solve  \n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{-1/2} = \\frac{e_{0} - e_{-1}}{\\Delta z_{-1/2}}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{N+1/2} = \\frac{e_{N+1} - e_{N}}{\\Delta z_{N+1/2}}\n",
    "$$\n",
    "\n",
    "#### Boundary Conditions\n",
    "\n",
    "Lets start with the bottom boundary - we know that MT fields and fluxes are diffusive and decay as they travel through conductive media, so if our boundary is sufficiently far away\n",
    "\n",
    "$$E_x (z=-\\infty) = 0$$\n",
    "\n",
    "Clearly we can't discretize to infinity... but we know approximately how quickly the fields decay, this is captured by the skin depth\n",
    "\n",
    "$$\n",
    "\\delta \\simeq \\frac{500}{\\sqrt{\\sigma f}}\n",
    "$$\n",
    "\n",
    "So as long as we define our mesh such that we are a few skin depths from the surface, then we can safely assume that the fields will have decayed to zero (dirichlet boundary condition). In our discrete world, this means that we want to enforce\n",
    "\n",
    "$$\n",
    "E_x \\big|_{nC+{1/2}} = 0\n",
    "$$\n",
    "\n",
    "The elements of $e$ are defined on cell centers and our boundary is a face, so we choose our ghost point $e_N$ such that when we average across the boundary, the average is 0, eg\n",
    "\n",
    "$$\n",
    "\\frac{1}{2} (e_{N-1} + e_{N}) = 0\n",
    "$$\n",
    "\n",
    "which means\n",
    "\n",
    "$$\n",
    "e_N = - e_{N-1}\n",
    "$$\n",
    "\n",
    "and our discrete approximation of the derivative at this boundary is\n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{N+1/2} = \\frac{e_{N+1} - e_{N}}{\\Delta z_{N+1/2}} = \\frac{-2 e_{N}}{\\Delta z_{N+1/2}}\n",
    "$$\n",
    "\n",
    "At the top boundary is where our incoming plane wave is, so we specify an electric field at the surface of \n",
    "\n",
    "$$E_x (z=0) = 1$$\n",
    "\n",
    "So this means we want to define our ghost point $e_{-1}$ such that\n",
    "\n",
    "$$\n",
    "\\frac{1}{2}(e_{-1} + e_0) = 1\n",
    "$$\n",
    "\n",
    "or \n",
    "\n",
    "$$\n",
    "(e_{-1}) = 2 - e_0\n",
    "$$\n",
    "\n",
    "and the derivative is \n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}\\bigg\\rvert_{-1/2} = \\frac{e_{0} - e_{-1}}{\\Delta z_{-1/2}} = \\underbrace{\\frac{2 e_{0}}{\\Delta z_{-1/2}}}_{\\text{due to dirichlet BC}} - \\underbrace{\\frac{2}{\\Delta z_{-1/2}}}_{\\text{due to non-homogeneous BC}}\n",
    "$$\n",
    "\n",
    "For conveienence, when we discretize, we first employ dirichlet boundary conditions on each boundary, and add the second term, due to a non-homogeneous boundary condition. \n",
    "\n",
    "The differential operator matrix that takes the derivative of a variable defined on faces is the cell gradient operator, so the discrete derivative of $E_x$ is given by\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z} \\simeq \\mathbf{Grad} ~ \\mathbf{e_x} + \\mathbf{B} ~ \\mathbf{e_x}^{BC}\n",
    "$$ \n",
    "\n",
    "where $\\mathbf{Grad}$ includes dirichlet boundary conditions, and $\\mathbf{B}$ is a $\\text{nC}\\times2$ matrix that accounts for the non-homogeneous boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.  0.  0.  0.]\n",
      " [-1.  1.  0.  0.]\n",
      " [ 0. -1.  1.  0.]\n",
      " [ 0.  0. -1.  1.]\n",
      " [ 0.  0.  0. -2.]]\n"
     ]
    }
   ],
   "source": [
    "# Grad matrix with dirichlet boundary conditions\n",
    "mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # set the boundary conditions\n",
    "Grad = mesh.cellGrad\n",
    "print(Grad.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.  0.]\n",
      " [ 0.  0.]\n",
      " [ 0.  0.]\n",
      " [ 0.  0.]\n",
      " [ 0.  2.]]\n"
     ]
    }
   ],
   "source": [
    "# deal with the boundary conditions\n",
    "ex_bc = np.r_[0., 1.] # bottom boundary, fields decay to zero, top is source\n",
    "B = mesh.cellGradBC  # a matrix for boundary conditions\n",
    "print(B.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.  0.  0.  0.  2.]\n"
     ]
    }
   ],
   "source": [
    "# B * e_BC describes what we need to add to Grad e in order to addount for \n",
    "# the boundary conditions\n",
    "print(B*ex_bc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last piece we need to define is how to take the product $\\imath \\omega \\mu H_y$. $\\imath$ and $\\omega$ are scalars, so they are easy. The tricky part is $\\mu H_y$ since $\\mathbf{\\mu}$ is defined at cell centers (there are $\\text{nC}$ of them) and $\\mathbf{h}$ is at faces (there are $\\text{nC+1}$ of them). So to take this product, we will average the magnetic permeability to faces, and again stick it in a diagonal matrix \n",
    "\n",
    "$$\\mathbf{M^{f}_{\\mu}} = \\mathbf{diag}(\\mathbf{Av^{cc2f} \\mathbf{\\mu}})$$\n",
    "\n",
    "so the product is then\n",
    "\n",
    "$$\n",
    "\\imath\\omega\\mu H_y \\simeq \\imath\\omega\\mathbf{M^{f}_{\\mu}} ~\\mathbf{h_y}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.   0.   0.   0. ]\n",
      " [ 0.5  0.5  0.   0. ]\n",
      " [ 0.   0.5  0.5  0. ]\n",
      " [ 0.   0.   0.5  0.5]\n",
      " [ 0.   0.   0.   1. ]]\n"
     ]
    }
   ],
   "source": [
    "# Averaging matrix\n",
    "AvCC2F = mesh.aveCC2F\n",
    "print(AvCC2F.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.25663706e-06   0.00000000e+00   0.00000000e+00   0.00000000e+00\n",
      "    0.00000000e+00]\n",
      " [  0.00000000e+00   1.25663706e-06   0.00000000e+00   0.00000000e+00\n",
      "    0.00000000e+00]\n",
      " [  0.00000000e+00   0.00000000e+00   1.25663706e-06   0.00000000e+00\n",
      "    0.00000000e+00]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.25663706e-06\n",
      "    0.00000000e+00]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00\n",
      "    1.25663706e-06]]\n"
     ]
    }
   ],
   "source": [
    "mu = mu_0*np.ones(mesh.nC)\n",
    "Mfmu = Utils.sdiag(AvCC2F * mu)\n",
    "print(Mfmu.todense())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have taken Faraday's law\n",
    "\n",
    "$$\n",
    "\\frac{\\partial E_x}{\\partial z}  + \\imath \\omega \\mu H_y = 0\n",
    "$$\n",
    "\n",
    "and arrived at the discrete system\n",
    "\n",
    "$$\n",
    "\\mathbf{Grad} ~ \\mathbf{e_x} + \\mathbf{B} ~ \\mathbf{e_x}^{BC} + \\imath\\omega \\mathbf{M^f_\\mu} \\mathbf{h_y} = 0\n",
    "$$\n",
    "\n",
    "since the boundary conditions are known, we can move them to the right hand side\n",
    "\n",
    "$$\n",
    "\\mathbf{Grad} ~ \\mathbf{e_x}  + \\imath\\omega \\mathbf{M^f_\\mu} \\mathbf{h_y} = - \\mathbf{B} ~ \\mathbf{e_x}^{BC}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two equations, Two unknowns\n",
    "\n",
    "Our discrete Maxwell system is \n",
    "\n",
    "$$\n",
    "\\mathbf{Grad} ~ \\mathbf{e_x}  + \\imath\\omega \\mathbf{M^f_\\mu} \\mathbf{h_y} = - \\mathbf{B} ~ \\mathbf{e_x}^{BC}\n",
    "$$\n",
    "\n",
    "$$\n",
    "- \\mathbf{Div} ~ \\mathbf{h_y} + \\mathbf{M^{cc}_{\\boldsymbol{\\hat{\\sigma}}}} ~ \\mathbf{e_x} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "For convienence, lets re-arrage... \n",
    "$$\n",
    "\\mathbf{Grad} ~ \\mathbf{e_x}  + \\imath\\omega \\mathbf{M^f_\\mu} \\mathbf{h_y} = - \\mathbf{B} ~ \\mathbf{e_x}^{BC}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{M^{cc}_{\\boldsymbol{\\hat{\\sigma}}}} ~ \\mathbf{e_x} - \\mathbf{Div} ~ \\mathbf{h_y} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "and assemble into a single matrix system\n",
    "\n",
    "$$\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "        \\mathbf{Grad} & \\imath \\omega \\mathbf{M}^{f2cc}_{\\mu} \\\\[0.3em]\n",
    "        \\mathbf{M}^{cc}_{\\hat{\\sigma}} & \\mathbf{Div}           \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{A}}\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "       \\mathbf{e}_x \\\\[0.3em]\n",
    "       \\mathbf{h}_y \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{u}}\n",
    "=\n",
    "\\underbrace{\n",
    "    \\begin{bmatrix}\n",
    "       - \\mathbf{B}\\mathbf{E}_x^{BC} \\\\[0.3em]\n",
    "       \\boldsymbol{0} \\\\[0.3em]\n",
    "    \\end{bmatrix}\n",
    "}_{\\mathbf{rhs}}\n",
    "$$\n",
    "\n",
    "with \n",
    "\n",
    "- $\\mathbf{e}_x$: Discrete $E_x$ $[\\text{nC} \\times 1]$\n",
    "\n",
    "- $\\mathbf{e}_y$: Dicrete $H_x$ $[(\\text{nC}+1) \\times 1]$\n",
    "\n",
    "- $ \\mathbf{Grad}$: Discrete gradient operator with dirichlet boundary conditions $[\\text{nC} \\times (\\text{nC}+1)]$\n",
    "\n",
    "- $ \\mathbf{Div}$: Discrete divergence operator $[(\\text{nC}+1) \\times \\text{nC}]$\n",
    "\n",
    "- $\\mathbf{M}^{f}_{\\boldsymbol{\\mu}} = \\mathbf{diag}(\\mathbf{Av^{cc2f}}  \\boldsymbol{\\mu})$ $[(\\text{nC}+1) \\times (\\text{nC}+1)]$\n",
    "\n",
    "- $\\mathbf{M}^{cc}_{\\boldsymbol{\\hat{\\sigma}}} = \\mathbf{diag}(\\boldsymbol{\\hat{\\sigma}})$ $[\\text{nC} \\times \\text{nC}]$\n",
    "\n",
    "- $\\mathbf{B} \\mathbf{e_x}^{BC}$ handles the boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now we have all of the pieces\n",
    "\n",
    "Here, lets create a larger mesh, and we can go ahead and asseble the system of equations to solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sigma_halfspace = 1e-2\n",
    "freq = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mesh extends 1.2e+05m, is that far enough? (should be at least 5.0e+03m away)\n"
     ]
    }
   ],
   "source": [
    "cs = 39.  # core cell size\n",
    "npad = 25  # number of padding cells\n",
    "ncz = 100   # number of core cells\n",
    "\n",
    "# define a tensor mesh\n",
    "hz = [(cs, npad, -1.3), (cs, ncz)] \n",
    "mesh = Mesh.TensorMesh([hz], x0='N')  # put the origin at the surface \n",
    "\n",
    "print(\n",
    "    \"The mesh extends {:1.1e}m, is that far enough? (should be at least {:1.1e}m away)\".format(\n",
    "        mesh.hx.sum(),\n",
    "        skin_depth(sigma_halfspace, freq)\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# physical properties\n",
    "sigma = np.ones(mesh.nC)*sigma_halfspace # conductivity values for all cells\n",
    "mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n",
    "epsilon = np.ones(mesh.nC)*epsilon_0 # dielectric constant values for all cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Grad \n",
    "mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n",
    "Grad = mesh.cellGrad # Gradient matrix\n",
    "\n",
    "# MfMu\n",
    "Mmu = Utils.sdiag(mesh.aveCC2F * mu) \n",
    "\n",
    "# Mccsigma\n",
    "Msighat = Utils.sdiag(sigmahat) \n",
    "\n",
    "# Div\n",
    "Div = mesh.faceDiv # Divergence matrix\n",
    "\n",
    "# Right Hand Side\n",
    "B = mesh.cellGradBC  # a matrix for boundary conditions\n",
    "Exbc = np.r_[0., 1.] # boundary values for Ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Assemble the matrix\n",
    "\n",
    "# A-matrix\n",
    "A = sp.vstack([\n",
    "    sp.hstack([Grad, 1j*omega*Mmu]), # Top row of A matrix\n",
    "    sp.hstack((Msighat, Div)) # Bottom row of A matrix\n",
    "])\n",
    "\n",
    "# Right-hand side\n",
    "rhs = np.r_[\n",
    "    -B*Exbc, \n",
    "    np.zeros(mesh.nC)\n",
    "] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.65 ms, sys: 990 µs, total: 2.64 ms\n",
      "Wall time: 1.71 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Ainv = Solver(A) # Factorize A matrix\n",
    "sol = Ainv*rhs   # Solve A^-1 rhs = sol\n",
    "Ex = sol[:mesh.nC] # Extract Ex from solution vector u\n",
    "Hy = sol[mesh.nC:mesh.nC+mesh.nN] # Extract Hy from solution vector u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Impedance: 6.2e-01 + 6.4e-01i\n",
      "or in terms of Amplidude: 8.9e-01 and phase: 45.9 degrees\n"
     ]
    }
   ],
   "source": [
    "Zxy = - 1./Hy[-1] # Impedance at the surface\n",
    "\n",
    "print(\"Impedance: {:1.1e} + {:1.1e}i\".format(Zxy.real, Zxy.imag))\n",
    "print(\n",
    "    \"or in terms of Amplidude: {:1.1e} and phase: {:1.1f} degrees\".format(\n",
    "    np.absolute(Zxy), np.rad2deg(np.arctan(Zxy.imag / Zxy.real))\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}