{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Constable and Weiss, 2006\n\nReproducing Figure 3 of Constable and Weiss, 2006, Geophysics. This is a marine\nCSEM example.\n\n**Reference**\n\n- **Constable, S., and C. J. Weiss,  2006**, Mapping thin resistors and\n  hydrocarbons with marine EM methods: Insights from 1D modeling: Geophysics,\n  71, G43-G51; DOI: [10.1190/1.2187748](http://dx.doi.org/10.1190/1.2187748).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import empymod\nimport numpy as np\nfrom copy import deepcopy as dc\nimport matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Computation\n\nNote: Exact reproduction is not possible, as source and receiver depths are\nnot explicitly specified in the publication. I made a few checks, and it\nlooks like a source-depth of 900 meter gives good accordance. Receivers are\non the sea-floor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Offsets\nx = np.linspace(0, 20000, 101)\n\n# TG model\ninp3 = {'src': [0, 0, 900],\n        'rec': [x, np.zeros(x.shape), 1000],\n        'depth': [0, 1000, 2000, 2100],\n        'res': [2e14, 0.3, 1, 100, 1],\n        'freqtime': 1,\n        'verb': 1}\n\n# HS model\ninp4 = dc(inp3)\ninp4['depth'] = inp3['depth'][:2]\ninp4['res'] = inp3['res'][:3]\n\n# Compute radial responses\nrhs = empymod.dipole(**inp4)  # Step, HS\nrhs = empymod.utils.EMArray(np.nan_to_num(rhs))\nrtg = empymod.dipole(**inp3)  # \" \"   Target\nrtg = empymod.utils.EMArray(np.nan_to_num(rtg))\n\n# Compute azimuthal response\nahs = empymod.dipole(**inp4, ab=22)  # Step, HS\nahs = empymod.utils.EMArray(np.nan_to_num(ahs))\natg = empymod.dipole(**inp3, ab=22)  # \" \"   Target\natg = empymod.utils.EMArray(np.nan_to_num(atg))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(9, 13))\nplt.subplots_adjust(wspace=.3, hspace=.3)\noldsettings = np.geterr()\n_ = np.seterr(all='ignore')\n\n# Radial amplitude\nplt.subplot(321)\nplt.title('(a) Radial mode fields')\nplt.plot(x/1000, np.log10(rtg.amp()), 'k', label='Model')\nplt.plot(x/1000, np.log10(rhs.amp()), 'k-.', label='Half-space response')\nplt.axis([0, 20, -18, -8])\nplt.xlabel('Range (km)')\nplt.ylabel(r'Log$_{10}$(E-field magnitude, V/Am$^2$)')\nplt.legend()\n\n# Radial phase\nplt.subplot(323)\nplt.title('(b) Radial mode phase')\nplt.plot(x/1000, rtg.pha(deg=True), 'k')\nplt.plot(x/1000, rhs.pha(deg=True), 'k-.')\nplt.axis([0, 20, -500, 0])\nplt.xlabel('Range (km)')\nplt.ylabel('Phase (degrees)')\n\n# Azimuthal amplitude\nplt.subplot(325)\nplt.title('(c) Azimuthal mode fields')\nplt.plot(x/1000, np.log10(atg.amp()), 'k', label='Model')\nplt.plot(x/1000, np.log10(ahs.amp()), 'k-.', label='Half-space response')\nplt.axis([0, 20, -18, -8])\nplt.xlabel('Range (km)')\nplt.ylabel(r'Log$_{10}$(E-field magnitude, V/Am$^2$)')\nplt.legend()\n\n# Azimuthal phase\nplt.subplot(322)\nplt.title('(d) Azimuthal mode phase')\nplt.plot(x/1000, atg.pha(deg=True)+180, 'k')\nplt.plot(x/1000, ahs.pha(deg=True)+180, 'k-.')\nplt.axis([0, 20, -500, 0])\nplt.xlabel('Range (km)')\nplt.ylabel('Phase (degrees)')\n\n# Normalized\nplt.subplot(324)\nplt.title('(e) Normalized E-field magnitude')\nplt.plot(x/1000, np.abs(rtg/rhs), 'k', label='Radial')\nplt.plot(x/1000, np.abs(atg/ahs), 'k--', label='Azimuthal')\nplt.axis([0, 20, 0, 70])\nplt.xlabel('Range (km)')\nplt.legend()\n\nplt.show()\n_ = np.seterr(**oldsettings)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Original Figure\n\nFigure 3 of Constable and Weiss, 2006, Geophysics:\n\n<img src=\"file://../../_static/figures/Constable2006.jpg\">\n\n\n"
      ]
    },
    {
      "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
}