{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Improve land CSEM computation\n\n## The problem\n\nThere exists a numerical singularity in the wavenumber-frequency domain. This\nsingularity is always present, but it is usually neglectible except when the\nresistivity is very high (like air; hence conductivity goes to zero and\ntherefore the real part of $\\eta_H$, $\\eta_V$ goes to zero) and\nsource and receiver are close to this boundary. So unfortunately exactly in the\ncase of a land CSEM survey.\n\nThis singularity leads to noise at very high frequencies and therefore at very\nearly times because the Hankel transform cannot handle the singularity\ncorrectly (or, if you would choose a sufficiently precise quadrature, it would\ntake literally forever to compute it).\n\n## The \"solution\"\n\nElectric permittivity ($\\varepsilon_H$, $\\varepsilon_V$) are set to\n1 by default. They are not important for the frequency range of CSEM. By\nsetting the electric permittivity of the air-layer to 0, the singularity\ndisapears, which subsquently improves a lot the time-domain result for land\nCSEM. It therefore uses the diffusive approximation for the air layer, but\nagain, that doesn't matter for the frequencies required for CSEM.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import empymod\nimport numpy as np\nimport matplotlib.pyplot as plt\nplt.style.use('ggplot')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Times (s)\nt = np.logspace(-2, 1, 301)\n\n# Model\nmodel = {\n    'src': [0, 0, 0.001],           # Src at origin, slightly below interface\n    'rec': [6000, 0, 0.0001],       # 6 km off., in-line, slightly bel. interf.\n    'depth': [0, 2000, 2100],       # Target of 100 m at 2 km depth\n    'res': [2e14, 10, 100, 10],     # Res: [air, overb., target, half-space]\n    'epermH': [1, 1, 1, 1],         # El. permittivity: default values\n    'freqtime': t,                  # Times\n    'signal': 0,                    # 0: Impulse response\n    'ftarg': {'dlf': 'key_81_CosSin_2009'},  # Shorter filter then the default\n    'verb': 1,                      # Verbosity; set to 3 to see all parameters\n}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute with default eperm_air = 1\nres_1 = empymod.dipole(**model)\n\n# Set horizontal and vertical electric permittivity of air to 0\nmodel['epermH'][0] = 0\n# Note that for empymod < v2.0.0 you have to set `epermH` AND `epermV`. From\n# v2.0.0 onwards `eperm` is assumed isotropic if `epermV` is not provided, and\n# `epermV` is therefore internally a copy of `epermH`.\n\n# Compute with eperm_air = 0\nres_0 = empymod.dipole(**model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot result\n\nAs it can be seen, setting $\\varepsilon_{air} =$ 0 improves the land\nCSEM result significantly for earlier times, where the signal should be zero.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure()\n\nplt.title('Time domain impulse response')\nplt.semilogx(t, res_1, label=r'$\\varepsilon_{air}=1$')\nplt.semilogx(t, res_0, label=r'$\\varepsilon_{air}=0$')\nplt.xlabel(r'Time $[s]$')\nplt.ylabel(r'Amplitude $[V/(m\\,s)]$')\nplt.legend()\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Version 1.7.0 and older\n\nIf you use a version of `empymod` that is smaller than `1.7.1` and set\n$\\varepsilon_H$, $\\varepsilon_V$ to zero you will see the\nfollowing warning,\n\n::\n\n    * WARNING :: Parameter epermH < 1e-20  are set to 1e-20 !\n    * WARNING :: Parameter epermV < 1e-20  are set to 1e-20 !\n\nand the values will be re-set to the defined minimum value, which is by\ndefault 1e-20. Using a value of 1e-20 for `epermH`/`epermV`  is also working\njust fine for land CSEM.\n\nIt is possible to change the minimum to zero for these old versions of\n`empymod`. However, there is a caveat: In `empymod v1.7.0` and older, the\n`min_param`-parameter also checks frequency and anisotropy. If you set\n`min_param=0`, and provide `empymod` with resistivities or anisotropies\nequals to zero, `empymod` will crash due to division by zero errors (avoiding\ndivision by zero is the purpose behind the `min_param`-parameter).\n\nTo change the `min_param`-parameter do:\n\n::\n\n    import empymod\n    empymod.set_minimum(min_param=0)\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
}