{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nCole-Cole\n=========\n\nThere are various different definitions of a Cole-Cole model, see for instance\nTarasov and Titov (2013). We try a few different ones here, but you can supply\nyour preferred version.\n\nThe original Cole-Cole (1940) model was formulated for the complex dielectric\npermittivity. It is reformulated to conductivity to use it for IP,\n\n\\begin{align}\\sigma(\\omega) = \\sigma_\\infty + \\frac{\\sigma_0 - \\sigma_\\infty}{1 +\n    (i\\omega\\tau)^C}\\ . \\qquad\\qquad\\qquad (1)\\end{align}\n\nAnother, similar model is given by Pelton et al. (1978),\n\n\\begin{align}\\rho(\\omega) = \\rho_\\infty + \\frac{\\rho_0 - \\rho_\\infty}{1 +\n    (i\\omega\\tau)^C}\\ . \\qquad\\qquad\\qquad (2)\\end{align}\n\nEquation (2) is just like equation (1), but replaces $\\sigma$ by\n$\\rho$. However, mathematically they are not the same. Substituting\n$\\rho = 1/\\sigma$ in the latter and resolving it for $\\sigma$ will\nnot yield the former. Equation (2) is usually written in the following form,\nusing the chargeability $m = (\\rho_0-\\rho_\\infty)/\\rho_0$,\n\n\\begin{align}\\rho(\\omega) = \\rho_0 \\left[1 - m \\left(1- \\frac{1}{1 + (i\\omega\\tau)^C}\n    \\right)\\right]\\ . \\quad (3)\\end{align}\n\nIn all cases we add the part coming from the dielectric permittivity\n(displacement currents), even tough it usually doesn't matter in the frequency\nrange of IP.\n\n**References**\n\n- **Cole, K.S., and R.H. Cole, 1941**, Dispersion and adsorption in\n  dielectrics. I. Alternating current characteristics; *Journal of Chemical\n  Physics*, Volume 9, Pages 341-351, doi:\n  `10.1063/1.1750906 <https://doi.org/10.1063/1.1750906>`_.\n- **Pelton, W.H., S.H. Ward, P.G. Hallof, W.R. Sill, and P.H. Nelson, 1978**,\n  Mineral discrimination and removal of inductive coupling with multifrequency\n  IP, *Geophysics*, Volume 43, Pages 588-609, doi:\n  `10.1190/1.1440839 <https://doi.org/10.1190/1.1440839>`_.\n- **Tarasov, A., and K. Titov, 2013**, On the use of the Cole\u2013Cole equations in\n  spectral induced polarization; *Geophysical Journal International*, Volume\n  195, Issue 1, Pages 352-356, doi:\n  `10.1093/gji/ggt251 <https://doi.org/10.1093/gji/ggt251>`_.\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": [
        "Use empymod with user-def. func. to adjust $\\eta$ and $\\zeta$\n-------------------------------------------------------------------------\n\nIn principal it is always best to write your own modelling routine if you\nwant to adjust something. Just copy ``empymod.dipole`` or ``empymod.bipole``\nas a template, and modify it to your needs. Since ``empymod v1.7.4``,\nhowever, there is a hook which allows you to modify $\\eta_h, \\eta_v,\n\\zeta_h$, and $\\zeta_v$ quite easily.\n\nThe trick is to provide a dictionary (we name it ``inp`` here) instead of the\nresistivity vector in ``res``. This dictionary, ``inp``, has two mandatory\nplus optional entries: - ``res``: the resistivity vector you would have\nprovided normally (mandatory).\n\n- A function name, which has to be either or both of (mandatory):\n\n  - ``func_eta``: To adjust ``etaH`` and ``etaV``, or\n  - ``func_zeta``: to adjust ``zetaH`` and ``zetaV``.\n\n- In addition, you have to provide all parameters you use in\n  ``func_eta``/``func_zeta`` and are not already provided to ``empymod``.\n  All additional parameters must have #layers elements.\n\nThe functions ``func_eta`` and ``func_zeta`` must have the following\ncharacteristics:\n\n- The signature is ``func(inp, p_dict)``, where\n\n  - ``inp`` is the dictionary you provide, and\n  - ``p_dict`` is a dictionary that contains all parameters so far computed\n    in empymod [``locals()``].\n\n- It must return ``etaH, etaV`` if ``func_eta``, or ``zetaH, zetaV`` if\n  ``func_zeta``.\n\nDummy example\n~~~~~~~~~~~~~\n\n::\n\n    def my_new_eta(inp, p_dict):\n        # Your computations, using the parameters you provided\n        # in ``inp`` and the parameters from empymod in ``p_dict``.\n        # In the example below, we provide, e.g., inp['tau']\n        return etaH, etaV\n\nAnd then you call ``empymod`` with ``res = {'res': res-array, 'tau': tau,\n'func_eta': my_new_eta}``.\n\nDefine the Cole-Cole model\n--------------------------\n\nIn this notebook we exploit this hook in empymod to compute $\\eta_h$\nand $\\eta_v$ with the Cole-Cole model. By default, $\\eta_h$ and\n$\\eta_v$ are computed like this:\n\n\\begin{align}\\eta_h = \\frac{1}{\\rho} + j\\omega \\varepsilon_{r;h}\\varepsilon_0 \\ ,\n    \\qquad (4)\\\\ \\eta_v = \\frac{1}{\\rho \\lambda^2} +\n    j\\omega\\varepsilon_{r;v}\\varepsilon_0 \\ . \\qquad (5)\\end{align}\n\n\nWith this function we recompute it. We replace the real part, the resistivity\n$\\rho$, in equations (4) and (5) by the complex, frequency-dependent\nCole-Cole resistivity [$\\rho(\\omega)$], as given, for instance, in\nequations (1)-(3). Then we add back the imaginary part coming from thet\ndielectric permittivity (basically zero for low frequencies).\n\nNote that in this notebook we use this hook to model relaxation in the low\nfrequency spectrum for IP measurements, replacing $\\rho$ by a\nfrequency-dependent model $\\rho(\\omega)$. However, this could also be\nused to model dielectric phenomena in the high frequency spectrum, replacing\n$\\varepsilon_r$ by a frequency-dependent formula\n$\\varepsilon_r(\\omega)$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def cole_cole(inp, p_dict):\n    \"\"\"Cole and Cole (1941).\"\"\"\n\n    # Compute complex conductivity from Cole-Cole\n    iotc = np.outer(2j*np.pi*p_dict['freq'], inp['tau'])**inp['c']\n    condH = inp['cond_8'] + (inp['cond_0']-inp['cond_8'])/(1+iotc)\n    condV = condH/p_dict['aniso']**2\n\n    # Add electric permittivity contribution\n    etaH = condH + 1j*p_dict['etaH'].imag\n    etaV = condV + 1j*p_dict['etaV'].imag\n\n    return etaH, etaV\n\n\ndef pelton_et_al(inp, p_dict):\n    \"\"\" Pelton et al. (1978).\"\"\"\n\n    # Compute complex resistivity from Pelton et al.\n    iotc = np.outer(2j*np.pi*p_dict['freq'], inp['tau'])**inp['c']\n    rhoH = inp['rho_0']*(1 - inp['m']*(1 - 1/(1 + iotc)))\n    rhoV = rhoH*p_dict['aniso']**2\n\n    # Add electric permittivity contribution\n    etaH = 1/rhoH + 1j*p_dict['etaH'].imag\n    etaV = 1/rhoV + 1j*p_dict['etaV'].imag\n\n    return etaH, etaV"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Example\n-------\n\nTwo half-space model, air above earth:\n\n- x-directed sourcer at the surface\n- x-directed receiver, also at the surface, inline at an offset of 500 m.\n- Switch-on time-domain response\n- Isotropic\n- Model [air, subsurface]\n\n    - $\\rho_\\infty = 1/\\sigma_\\infty =$ [2e14, 10]\n    - $\\rho_0 = 1/\\sigma_0 =$ [2e14, 5]\n    - $\\tau =$ [0, 1]\n    - $c =$ [0, 0.5]\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Times\ntimes = np.logspace(-2, 2, 101)\n\n# Model parameter which apply for all\nmodel = {\n    'src': [0, 0, 1e-5, 0, 0],\n    'rec': [500, 0, 1e-5, 0, 0],\n    'depth': 0,\n    'freqtime': times,\n    'signal': 1,\n    'verb': 1\n}\n\n# Collect Cole-Cole models\nres_0 = np.array([2e14, 10])\nres_8 = np.array([2e14, 5])\ntau = [0, 1]\nc = [0, 0.5]\nm = (res_0-res_8)/res_0\n\ncole_model = {'res': res_0, 'cond_0': 1/res_0, 'cond_8': 1/res_8,\n              'tau': tau, 'c': c, 'func_eta': cole_cole}\npelton_model = {'res': res_0, 'rho_0': res_0, 'm': m,\n                'tau': tau, 'c': c, 'func_eta': pelton_et_al}\n\n# Compute\nout_bipole = empymod.bipole(res=res_0, **model)\nout_cole = empymod.bipole(res=cole_model, **model)\nout_pelton = empymod.bipole(res=pelton_model, **model)\n\n# Plot\nplt.figure()\nplt.title('Switch-off')\nplt.plot(times, out_bipole, label='Regular Bipole')\nplt.plot(times, out_cole, '--', label='Cole and Cole (1941)')\nplt.plot(times, out_pelton, '-.', label='Pelton et al. (1978)')\nplt.legend()\nplt.yscale('log')\nplt.xscale('log')\nplt.xlabel('time (s)')\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.7.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}