{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Definition of the coordinate system in empymod\n\n## Short version\n\nThe used coordinate system is either a\n\n- Left-Handed System (LHS), where Easting is the $x$-direction, Northing\n  the $y$-direction, and positive $z$ is pointing downwards;\n- Right-Handed System (RHS), where Easting is the $x$-direction, Northing\n  the $y$-direction, and positive $z$ is pointing upwards.\n\n\n## In more detail\n\nThe derivation on which ``empymod`` is based ([HuTS15]_) uses a right-handed\nsystem with $x$ to the East, $y$ to the South, and $z$\ndownwards (ESD). In the actual original implementation of ``empymod`` this was\nchanged to a left-handed system with $x$ to the East, $y$ to the\nNorth, and $z$ downwards (END). However, ``empymod`` can equally well be\nused for a coordinate system where positive $z$ is pointing up by just\nflipping $z$, resulting in $x$ to the East, $y$ to the North,\nand $z$ upwards (ENU).\n\n\n+----------------+--------------------+---------------------+\n|                | Left-handed system | Right-handed system |\n+================+====================+=====================+\n| $x$      | Easting            | Easting             |\n+----------------+--------------------+---------------------+\n| $y$      | Northing           | Northing            |\n+----------------+--------------------+---------------------+\n| $z$      | Down               | Up                  |\n+----------------+--------------------+---------------------+\n| $\\theta$ | Angle E-N          | Angle E-N           |\n+----------------+--------------------+---------------------+\n| $\\varphi$| Angle down         | Angle up            |\n+----------------+--------------------+---------------------+\n\nThere are a few other important points to keep in mind when switching between\ncoordinate systems:\n\n- The interfaces (``depth``) have to be defined continuously increasing or\n  decreasing. E.g., think of a simple five-layer model with the sea-surface at\n  0 m, a 100 m water column, and a target of 50 m 900 m below the seafloor,\n  with the following resistivities:\n\n  - ``res = [1e12, 0.3, 1, 50, 1]``: air, water, background, target,\n    background.\n\n  We can now define the depths in two different ways:\n\n  - ``depth = [0, 100, 1000, 1050]``: **LHS** (+z down)\n  - ``depth = [0, -100, -1000, -1050]``: **RHS** (+z up)\n\n- A source or a receiver *exactly on* a boundary is taken as being in the lower\n  layer. Hence, if $z_\\rm{rec} = z_0$, where $z_0$ is the surface,\n  then the receiver is taken as in the air in the LHS, but as in the subsurface\n  in the RHS. Similarly, if $z_\\rm{rec} = z_\\rm{seafloor}$, then the\n  receiver is taken as in the sea in the LHS, but as in the subsurface in the\n  RHS. This can be avoided by never placing it exactly on a boundary, but\n  slightly (e.g., 1 mm) in the layer where you want to have it.\n\n- Sign switches:\n\n  - In ``bipole``, the ``dip`` switches sign.\n  - In ``dipole``, the ``ab``'s containing vertical directions switch the sign\n    for each vertical component.\n  - Sign switches also occur for magnetic sources or receivers.\n\n  These sign switches multiply; e.g., ``ab=34`` has no sign switch, as the\n  switch due to the vertical source and the switch due to the magnetic receiver\n  cancel each other. Here an overview of the sign switches (--) for ``dipole``\n  when changing from the default LHS to a RHS:\n\n  +---------------+-------+------+------+------+------+------+------+\n  |                       | electric  source   | magnetic source    |\n  +===============+=======+======+======+======+======+======+======+\n  |                       | **x**| **y**| **z**| **x**| **y**| **z**|\n  +---------------+-------+------+------+------+------+------+------+\n  |               | **x** |      |      |  --  |  --  |  --  |      |\n  + **electric**  +-------+------+------+------+------+------+------+\n  |               | **y** |      |      |  --  |  --  |  --  |      |\n  + **receiver**  +-------+------+------+------+------+------+------+\n  |               | **z** |  --  |  --  |      |      |      |  --  |\n  +---------------+-------+------+------+------+------+------+------+\n  |               | **x** |  --  |  --  |      |      |      |  --  |\n  + **magnetic**  +-------+------+------+------+------+------+------+\n  |               | **y** |  --  |  --  |      |      |      |  --  |\n  + **receiver**  +-------+------+------+------+------+------+------+\n  |               | **z** |      |      |  --  |  --  |  --  |      |\n  +---------------+-------+------+------+------+------+------+------+\n\n- The depths can also be defined the other way around, starting at +/-1050\n  going to 0. You have to change the order of all model parameters (``res``,\n  ``aniso``, ``{e;m}perm{H;V}``) accordingly.\n\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>In a two-layer scenario with only one interface ``depth = z`` it always\n  assumes a **LHS**, as it is not possible to detect the direction from only\n  one interface. To force a **RHS** you have to add ``-np.infty`` for the\n  down-most interface:\n\n  - **LHS**: ``depth = z`` (+z down); default, corresponds to ``[z, np.infty]``\n  - **RHS**: ``depth = [z, -np.infty]`` (+z up)</p></div>\n\n\nIn this example we first create a sketch of the LHS and RHS for visualization,\nfollowed by a few examples using ``dipole`` and ``bipole`` to demonstrate the\ntwo possibilities.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import empymod\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom mpl_toolkits.mplot3d import proj3d\nfrom matplotlib.patches import FancyArrowPatch\nplt.style.use('ggplot')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## RHS vs LHS\n\nComparison of the right-handed system with positive $z$ downwards and\nthe left-handed system with positive $z$ upwards. Easting is always\n$x$, and Northing is $y$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class Arrow3D(FancyArrowPatch):\n    \"\"\"https://github.com/matplotlib/matplotlib/issues/21688\"\"\"\n\n    def __init__(self, xs, ys, zs):\n        super().__init__(\n                (0, 0), (0, 0), mutation_scale=20, lw=1.5,\n                arrowstyle='-|>', color='.2', zorder=100)\n        self._verts3d = xs, ys, zs\n\n    def do_3d_projection(self, renderer=None):\n        xs3d, ys3d, zs3d = self._verts3d\n        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.M)\n        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))\n        return np.min(zs)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def repeated(ax, pm):\n    \"\"\"These are all repeated for the two subplots.\"\"\"\n\n    # Coordinate system\n    # The first three are not visible, but for the aspect ratio of the plot.\n    ax.plot([-2, 12], [0, 0], [0, 0], c='w')\n    ax.plot([0, 0], [-2, 12], [0, 0], c='w')\n    ax.plot([0, 0], [0, 0], [-pm*2, pm*12], c='w')\n    ax.add_artist(Arrow3D([-2, 14], [0, 0], [0, 0]))\n    ax.add_artist(Arrow3D([0, 0], [-2, 14], [0, 0]))\n    ax.add_artist(Arrow3D([0, 0], [0, 0], [-pm*2, pm*14]))\n\n    # Annotate it\n    ax.text(12, 2, 0, r'$x$')\n    ax.text(0, 12, 2, r'$y$')\n    ax.text(-2, 0, pm*12, r'$z$')\n\n    # Helper lines\n    ax.plot([0, 10], [0, 10], [0, 0], '--', c='.6')\n    ax.plot([0, 10], [0, 0], [0, -10], '--', c='.6')\n    ax.plot([10, 10], [0, 10], [0, 0], ':', c='.6')\n    ax.plot([10, 10], [10, 10], [0, -10], ':', c='.6')\n    ax.plot([10, 10], [0, 0], [0, -10], ':', c='.6')\n    ax.plot([10, 10], [0, 10], [-10, -10], ':', c='.6')\n\n    # Resulting trajectory\n    ax.plot([0, 10], [0, 10], [0, -10], 'C0')\n\n    # Theta\n    azimuth = np.linspace(np.pi/4, np.pi/2, 31)\n    ax.plot(np.sin(azimuth)*5, np.cos(azimuth)*5, 0, c='C5')\n    ax.text(3, 5, 0, r\"$\\theta$\", color='C5', fontsize=14)\n\n    # Phi\n    ax.plot(np.sin(azimuth)*7, azimuth*0, -np.cos(azimuth)*7, c='C1')\n\n    ax.view_init(azim=-60, elev=20)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create figure\nfig = plt.figure(figsize=(8, 3.5))\n\n# Left-handed system\nax1 = fig.add_subplot(121, projection=Axes3D.name, facecolor='w')\nax1.axis('off')\nplt.title('Left-handed system (LHS)\\nfor positive $z$ downwards', fontsize=12)\nax1.text(7, 0, -5, r\"$\\varphi$\", color='C1', fontsize=14)\n\nrepeated(ax1, -1)\n\n# Right-handed  system\nax2 = fig.add_subplot(122, projection='3d', facecolor='w', sharez=ax1)\nax2.axis('off')\nplt.title('Right-handed system (RHS)\\nfor positive $z$ upwards', fontsize=12)\nax2.text(7, 0, -5, r\"$-\\varphi$\", color='C1', fontsize=14)\n\nrepeated(ax2, 1)\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Dipole\n\nA simple example using ``dipole``. It is a marine case with 300 meter water\ndepth and a 50 m thick target 700 m below the seafloor.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "off = np.linspace(500, 10000, 301)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### LHS\n\nIn the left-handed system positive $z$ is downwards. So we have to\ndefine our model by beginning with the air layer, followed by water,\nbackground, target, and background again. This means that all our\ndepth-values are positive, as the air-interface $z_0$ is at 0 m.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lhs = empymod.dipole(\n        src=[0, 0, 100],\n        rec=[off, np.zeros(off.size), 200],\n        depth=[0, 300, 1000, 1050],\n        res=[1e20, 0.3, 1, 50, 1],\n        # depth=[1050, 1000, 300, 0],  # Alternative way, LHS high to low.\n        # res=[1, 50, 1, 0.3, 1e20],   # \" \" \"\n        freqtime=1,\n        verb=0\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### RHS\n\nIn the right-handed system positive $z$ is upwards. So we have to\ndefine our model by beginning with the background, followed by the target,\nbackground again, water, and air. This means that all our depth-values are\nnegative.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rhs = empymod.dipole(\n        src=[0, 0, -100],\n        rec=[off, np.zeros(off.size), -200],\n        depth=[0, -300, -1000, -1050],\n        res=[1e20, 0.3, 1, 50, 1],\n        # depth=[-1050, -1000, -300, 0],  # Alternative way, RHS low to high.\n        # res=[1, 50, 1, 0.3, 1e20],      # \" \" \"\n        freqtime=1,\n        verb=0\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Compare\n\nPlotting the two confirms that the results agree, no matter if we use the LHS\nor the RHS definition.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(9, 4))\n\nax1 = plt.subplot(121)\nplt.title('Real')\nplt.plot(off/1e3, lhs.real, 'C0', label='LHS +')\nplt.plot(off/1e3, -lhs.real, 'C4', label='LHS -')\nplt.plot(off/1e3, rhs.real, 'C1--', label='RHS +')\nplt.plot(off/1e3, -rhs.real, 'C2--', label='RHS -')\nplt.yscale('log')\nplt.xlabel('Offset (km)')\nplt.ylabel('$E_x$ (V/m)')\nplt.legend()\n\nax2 = plt.subplot(122, sharey=ax1)\nplt.title('Imaginary')\nplt.plot(off/1e3, lhs.imag, 'C0', label='LHS +')\nplt.plot(off/1e3, -lhs.imag, 'C4', label='LHS -')\nplt.plot(off/1e3, rhs.imag, 'C1-.', label='RHS +')\nplt.plot(off/1e3, -rhs.imag, 'C2:', label='RHS -')\nplt.yscale('log')\nplt.xlabel('Offset (km)')\nplt.ylabel('$E_x$ (V/m)')\nplt.legend()\nax2.yaxis.set_label_position(\"right\")\nax2.yaxis.tick_right()\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Bipole [x1, x2, y1, y2, z1, z2]\n\nA time-domain example using rotated bipoles, where we define them as\n$[x_1, x_2, y_1, y_2, z_1, z_2]$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = np.linspace(0.1, 10, 301)\ninp = {'freqtime': times, 'signal': 0, 'verb': 0}\n\nlhs = empymod.bipole(\n        src=[-50, 50, -10, 10, 100, 110],\n        rec=[6000, 6100, 20, -20, 220, 200],\n        depth=[0, 300, 1000, 1050],\n        res=[1e20, 0.3, 1, 50, 1],\n        **inp\n)\n\nrhs = empymod.bipole(\n        src=[-50, 50, -10, 10, -100, -110],\n        rec=[6000, 6100, 20, -20, -220, -200],\n        depth=[0, -300, -1000, -1050],\n        res=[1e20, 0.3, 1, 50, 1],\n        **inp\n)\n\nplt.figure()\n\nplt.plot(times, lhs, 'C0', label='LHS')\nplt.plot(times, rhs, 'C1--', label='RHS')\nplt.xlabel('Time (s)')\nplt.ylabel('$E_x$ (V/m)')\nplt.legend()\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Bipole [x, y, z, azimuth, dip]\n\nA very similar time-domain example using rotated bipoles, but this time\ndefining them as $[x, y, z, \\theta, \\varphi]$. Note that\n$\\varphi$ has to change the sign, while $\\theta$ does not.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lhs = empymod.bipole(\n        src=[0, 0, 100, 10, 20],\n        rec=[6000, 0, 200, -5, 15],\n        depth=[0, 300, 1000, 1050],\n        res=[1e20, 0.3, 1, 50, 1],\n        **inp\n)\n\nrhs = empymod.bipole(\n        src=[0, 0, -100, 10, -20],\n        rec=[6000, 0, -200, -5, -15],\n        depth=[0, -300, -1000, -1050],\n        res=[1e20, 0.3, 1, 50, 1],\n        **inp\n)\n\nplt.figure()\n\nplt.plot(times, lhs, 'C0', label='LHS')\nplt.plot(times, rhs, 'C1--', label='RHS')\nplt.xlabel('Time (s)')\nplt.ylabel('$E_x$ (V/m)')\nplt.legend()\n\nplt.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
}