{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Magnetotelluric data\n\nThe magnetotelluric (MT) method is a passive method using as a source\nvariations in Earth's magnetic field, which create telluric currents in the\nEarth. The variation of Earth's magnetic field has many origins, e.g.,\nlightning or the interaction between the Earth's magnetic field and solar wind.\n\nUbiquitous in MT is the plane-wave approximation, hence, that the source signal\nis a plane wave hitting the Earth's surface. Having a 1D source (vertical plane\nwave) for a layered, 1D model reduces the computation of the impedances and\nfrom there the apparent resistivity and apparent phase a simple recursion\nalgorithm. As such it does not make sense to use a full EM wavefield algorithm\nwith three-dimensional sources such as empymod to compute MT responses.\nHowever, it is still interesting to see if we can compute MT impedances with a\nthree-dimensional point source.\n\nAs background theory we reproduce here Equations (11) to (17) from Pedersen and\nHermance (1986), with surrounding text. For a more in-depth read we refer to\nChave and Jones (2012).\n\n--------\n\nIf we define the impedance as\n\n\\begin{align}:label: ph-11\n\n    Z = E_x / H_y \\, ,\\end{align}\n\n[...]  we can develop a recursive relationship for the impedance at the top of\nthe j-th layer looking down\n\n\\begin{align}:label: ph-12\n\n    Z_j = z_{oj} \\frac{1-R_j \\exp(-2\\gamma_j t_j)}{1+R_j \\exp(-2\\gamma_j t_j)}\n        \\, , \\quad\n    j = N-1, \\dots, 1 \\, ,\\end{align}\n\nwhere\n\n\\begin{align}:label: ph-13\n\n    z_{oj} \\equiv \\text{intrinsic impedance}\n        \\equiv \\sqrt{\\rm{i}\\omega \\mu \\rho_j} \\, ,\\end{align}\n\n\\begin{align}:label: ph-14\n\n    R_j \\equiv \\text{reflection coefficient}\n        \\equiv \\frac{z_{oj} - Z_{j+1}}{z_{oj} + Z_{j+1}} \\, ,\\end{align}\n\n\\begin{align}t_j = \\text{thickness of layer} j \\, ,\\end{align}\n\nand the impedance at the surface of the deepest layer ($j=N$) is given by\n\n\\begin{align}:label: ph-15\n\n    Z_N = z_{oN}\\, .\\end{align}\n\nThe surface impedance, $Z_j$, is found by applying Equation :eq:`ph-12`\nrecursively from the top of the bottom half-space, $j = N$, and\npropagating upwards. From the surface impedance, $Z_1$, we can then\ncalculate the apparent resistivity, $\\rho_a$, and phase,\n$\\theta_a$, as\n\n\\begin{align}:label: ph-16\n\n    \\rho_a = \\frac{|Z_1|^2}{\\omega \\mu} \\, ,\\end{align}\n\n\\begin{align}:label: ph-17\n\n    \\theta_a = \\tan^{-1}\\frac{\\Im(Z1)}{\\Re(Z_1)} \\, .\\end{align}\n\nThis calculation is repeated for a range of periods and is used to model the\nmagnetotelluric response of the layered structure.\n\n--------\n\n**Reference**:\n\n- Chave, A., and Jones, A. (Eds.). (2012). The Magnetotelluric Method: Theory\n  and Practice. Cambridge: Cambridge University Press;\n  https://doi.org/10.1017/CBO9781139020138.\n- Pedersen, J., Hermance, J.F. Least squares inversion of one-dimensional\n  magnetotelluric data: An assessment of procedures employed by Brown\n  University. Surv Geophys 8, 187\u2013231 (1986);\n  https://doi.org/10.1007/BF01902413.\n- This example was strongly motivated by Andrew Pethicks blog post\n  https://www.digitalearthlab.com/tutorial/tutorial-1d-mt-forward.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import empymod\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.constants import mu_0\nplt.style.use('ggplot')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define model parameter and frequencies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "resistivities = np.array([2e14, 300, 2500, 0.8, 3000, 2500])\ndepths = np.array([0, 200, 600, 640, 1140])\nfrequencies = np.logspace(-4, 5, 101)\nomega = 2 * np.pi * frequencies"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 1D-MT recursion following Pedersen & Hermance\n\nUsing the variable names as in the paper.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Initiate recursive formula with impedance of the deepest layer.\nZ_j = np.sqrt(2j * np.pi * frequencies * mu_0 * resistivities[-1])\n\n# Move up the stack of layers till the top (without air).\nfor j in range(depths.size-1, 0, -1):\n\n    # Thickness\n    t_j = depths[j] - depths[j-1]\n\n    # Intrinsic impedance\n    z_oj = np.sqrt(1j * omega * mu_0 * resistivities[j])\n\n    # Reflection coefficient\n    R_j = (z_oj - Z_j) / (z_oj + Z_j)\n\n    # Exponential factor\n    gamma_j = np.sqrt(1j * omega * mu_0 / resistivities[j])\n    exp_j = np.exp(-2 * gamma_j * t_j)\n\n    # Impedance at this layer\n    Z_j = z_oj * (1 - R_j * exp_j) / (1 + R_j * exp_j)\n\n# Step 3. Compute apparent resistivity last impedance\napres_mt1d = abs(Z_j)**2/(omega * mu_0)\nphase_mt1d = np.arctan2(Z_j.imag, Z_j.real)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1D MT using empymod\n\nThe above derivation and code assume plane waves. We can \"simulate\" plane\nwaves by putting the source _very_ far away. In this example, we set it one\nmillion km away in all directions. As the impedance is the ratio of the Ex\nand Hy fields the source type (electric or magnetic) and the source\norientation do not matter; hence, it can be an arbitrarily rotated electric\nor magnetic source, the apparent resistivity and apparent phase will always\nbe the same.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dist = 1_000_000_000  # 1 million kilometer (!)\ninp = {\n    'src': (-dist, -dist, -dist),\n    'rec': (0, 0, 0.1),\n    'res': resistivities,\n    'depth': depths,\n    'freqtime': frequencies,\n    'verb': 1,\n}\n\n# Get Ex, Hy.\nex = empymod.dipole(ab=11, **inp)\nhy = empymod.dipole(ab=51, **inp)\n\n# Impedance.\nZ = ex/hy\n\n# Apparent resistivity and apparent phase.\napres_empy = abs(Z)**2 / (omega * mu_0)\nphase_empy = np.arctan2(Z.imag, Z.real)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot results\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\n\nax1.set_title('Apparent resistivity')\nax1.plot(frequencies, apres_mt1d, label='MT-1D')\nax1.plot(frequencies, apres_empy, '--', label='empymod')\nax1.set_xscale('log')\nax1.set_yscale('log')\nax1.set_xlabel('Frequency (Hz)')\nax1.set_ylabel(r'Apparent resistivity ($\\Omega\\,$m)')\nax1.legend()\n\nax2.set_title('Phase')\nax2.plot(frequencies, phase_mt1d*180/np.pi)\nax2.plot(frequencies, phase_empy*180/np.pi, '--')\nax2.set_xscale('log')\nax2.yaxis.tick_right()\nax2.set_xlabel('Frequency (Hz)')\nax2.set_ylabel('Phase (degree)')\nax2.yaxis.set_label_position(\"right\")\n\nfig.tight_layout()\nfig.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "empymod.Report()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}