{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Ward and Hohmann, 1988, SEG\n\nFrequency and time-domain modelling of magnetic loop sources and magnetic\nreceivers.\n\nReproducing Figures 2.2-2.5, 4.2-4.5, and 4.7-4.8 of Ward and Hohmann (1988):\nFrequency- and time-domain isotropic solutions for a full-space (2.2-2.5) and a\nhalf-space (4.2-4.5, 4.7-4.8), where source and receiver are at the interface.\nSource is a loop, receiver is a magnetic dipole.\n\n**Reference**\n\n- **Ward, S. H., and G. W. Hohmann, 1988**, Electromagnetic theory for\n geophysical applications, Chapter 4 of Electromagnetic Methods in Applied\n Geophysics: SEG, Investigations in Geophysics No. 3, 130--311; DOI:\n [10.1190/1.9781560802631.ch4](https://doi.org/10.1190/1.9781560802631.ch4).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import empymod\nimport numpy as np\nfrom scipy.special import erf\nimport matplotlib.pyplot as plt\nfrom scipy.constants import mu_0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ward and Hohmann, 1988, Fig 4.4\n\nWard and Hohmann (1988), Equations 4.69a and 4.70:\n\n\\begin{align}h_z = \\frac{m}{4\\pi r^3} \\left[\n \\frac{9}{2\\theta^2 r^2} \\rm{erf}(\\theta r) - \\rm{erf}(\\theta r) -\n \\frac{1}{\\pi^{1/2}} \\left(\\frac{9}{\\theta r} + 4\\theta r\\right)\n \\exp(-\\theta^2 r^2) \\right] \\, , \\qquad (4.69\\rm{a})\\end{align}\n\nand\n\n\\begin{align}\\frac{\\partial h_z}{\\partial t} = -\\frac{m\\rho}{2\\pi\\mu_0 r^5} \\left[\n 9\\rm{erf}(\\theta r) -\n \\frac{2\\theta r}{\\pi^{1/2}} \\left(9 + 6\\theta^2 r^2 +\n 4\\theta^4 r^4\\right) \\exp(-\\theta^2 r^2) \\right] \\, , \\qquad (4.70)\\end{align}\n\nwhere\n\n\\begin{align}\\theta = \\sqrt{\\frac{\\mu_0}{4t\\rho}} \\, ,\\end{align}\n\n$t$ is time in s, $\\rho$ is resistivity in\n$\\Omega\\,\\text{m}$, $r$ is offset in m, and $m$ the\nmagnetic moment in $\\text{A}\\,\\text{m}^2$ .\n\n### Analytical solutions\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def hz(t, res, r, m=1.):\n r\"\"\"Return equation 4.69a, Ward and Hohmann, 1988.\n\n Switch-off response (i.e., Hz(t)) of a homogeneous isotropic half-space,\n where the vertical magnetic source and receiver are at the interface.\n\n Parameters\n ----------\n t : array\n Times (t)\n res : float\n Halfspace resistivity (Ohm.m)\n r : float\n Offset (m)\n m : float, optional\n Magnetic moment, default is 1.\n\n Returns\n -------\n hz : array\n Vertical magnetic field (A/m)\n\n \"\"\"\n theta = np.sqrt(mu_0/(4*res*t))\n theta_r = theta*r\n\n s = 9/(2*theta_r**2)*erf(theta_r) - erf(theta_r)\n s -= (9/theta_r+4*theta_r)*np.exp(-theta_r**2)/np.sqrt(np.pi)\n s *= m/(4*np.pi*r**3)\n\n return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def dhzdt(t, res, r, m=1.):\n r\"\"\"Return equation 4.70, Ward and Hohmann, 1988.\n\n Impulse response (i.e., dHz(t)/dt) of a homogeneous isotropic half-space,\n where the vertical magnetic source and receiver are at the interface.\n\n Parameters\n ----------\n t : array\n Times (t)\n res : float\n Halfspace resistivity (Ohm.m)\n r : float\n Offset (m)\n m : float, optional\n Magnetic moment, default is 1.\n\n Returns\n -------\n dhz : array\n Time-derivative of the vertical magnetic field (A/m/s)\n\n \"\"\"\n theta = np.sqrt(mu_0/(4*res*t))\n theta_r = theta*r\n\n s = (9 + 6 * theta_r**2 + 4 * theta_r**4) * np.exp(-theta_r**2)\n s *= -2 * theta_r / np.sqrt(np.pi)\n s += 9 * erf(theta_r)\n s *= -(m*res)/(2*np.pi*mu_0*r**5)\n\n return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def pos(data):\n \"\"\"Return positive data; set negative data to NaN.\"\"\"\n return np.array([x if x > 0 else np.nan for x in data])\n\n\ndef neg(data):\n \"\"\"Return -negative data; set positive data to NaN.\"\"\"\n return np.array([-x if x < 0 else np.nan for x in data])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Survey parameters\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "time = np.logspace(-8, 0, 301)\n\nsrc = [0, 0, 0, 0, 90]\nrec = [100, 0, 0, 0, 90]\ndepth = 0\nres = [2e14, 100]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analytical result\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hz_ana = hz(time, res[1], rec[0])\ndhz_ana = dhzdt(time, res[1], rec[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical result\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "eperm = [0, 0] # Reduce early time numerical noise (diffusive approx for air)\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': time, 'verb': 1, 'xdirect': True, 'epermH': eperm}\n\nhz_num = empymod.loop(signal=-1, **inp)\ndhz_num = empymod.loop(signal=0, **inp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the result\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(5, 6))\n\nplt.plot(time*1e3, pos(dhz_ana), 'k-', lw=2, label='Ward & Hohmann')\nplt.plot(time*1e3, neg(dhz_ana), 'k--', lw=2)\nplt.plot(time*1e3, pos(dhz_num), 'C1-', label='empymod; dHz/dt')\nplt.plot(time*1e3, neg(dhz_num), 'C1--')\n\nplt.plot(time*1e3, pos(hz_ana), 'k-', lw=2)\nplt.plot(time*1e3, neg(hz_ana), 'k--', lw=2)\nplt.plot(time*1e3, pos(hz_num), 'C0-', label='empymod; Hz')\nplt.plot(time*1e3, neg(hz_num), 'C0--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-5, 1e3])\nplt.yticks(10**np.arange(-11., 0))\nplt.ylim([1e-11, 1e-1])\nplt.xlabel('time (ms)')\nplt.legend()\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\nNote that $h_z$ has the opposite sign in the original figure, which is\nprobably a typo (it is not what their equation yields).\n\n\nThe following examples are just compared to the figures, without the provided\nanalytical solutions.\n\n## Ward and Hohmann, 1988, Fig 4.2\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nfreq = np.logspace(-1, 5, 301)\nsrc = [0, 0, 0, 0, 90]\nrec = [100, 0, 0, 0, 90]\ndepth = 0\nres = [2e14, 100]\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': freq, 'verb': 1}\nfhz_num = empymod.loop(**inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nplt.plot(freq, pos(fhz_num.real), 'C0-', label='Real')\nplt.plot(freq, neg(fhz_num.real), 'C0--')\n\nplt.plot(freq, pos(fhz_num.imag), 'C1-', label='Imaginary')\nplt.plot(freq, neg(fhz_num.imag), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-1, 1e5])\nplt.ylim([1e-12, 1e-6])\nplt.xlabel('FREQUENCY (Hz)')\nplt.ylabel('$H_z$ (A/m)')\nplt.legend()\n\nplt.tight_layout()\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 4.3\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nfreq = np.logspace(-1, 5, 301)\nsrc = [0, 0, 0, 0, 90]\nrec = [100, 0, 0, 0, 0]\ndepth = 0\nres = [2e14, 100]\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': freq, 'verb': 1}\nfhz_num = empymod.loop(**inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nplt.plot(freq, pos(fhz_num.real), 'C0-', label='Real')\nplt.plot(freq, neg(fhz_num.real), 'C0--')\n\nplt.plot(freq, pos(fhz_num.imag), 'C1-', label='Imaginary')\nplt.plot(freq, neg(fhz_num.imag), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-1, 1e5])\nplt.ylim([1e-12, 1e-6])\nplt.xlabel('FREQUENCY (Hz)')\nplt.ylabel(r'$H_{\\rho}$ (A/m)')\nplt.legend()\n\nplt.tight_layout()\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n## Ward and Hohmann, 1988, Fig 4.5\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\ntime = np.logspace(-6, 0.5, 301)\nsrc = [0, 0, 0, 0, 90]\nrec = [100, 0, 0, 0, 0]\ndepth = 0\nres = [2e14, 100]\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'epermH': eperm, 'freqtime': time, 'verb': 1}\nfhz_num = empymod.loop(signal=1, **inp)\nfdhz_num = empymod.loop(signal=0, **inp)\n\n# Figure\nplt.figure(figsize=(5, 6))\n\nax1 = plt.subplot(111)\nplt.plot(time*1e3, pos(fdhz_num), 'C0-', label='dHz/dt')\nplt.plot(time*1e3, neg(fdhz_num), 'C0--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-3, 2e3])\nplt.yticks(10**np.arange(-11., -1))\nplt.ylim([1e-11, 1e-1])\nplt.xlabel('time (ms)')\nplt.ylabel(r'$\\frac{\\partial h_{\\rho}}{\\partial t}$ (A/m-s)')\nplt.legend(loc=8)\n\nax2 = ax1.twinx()\n\nplt.plot(time*1e3, pos(fhz_num), 'C1-', label='Hz')\nplt.plot(time*1e3, neg(fhz_num), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-3, 2e3])\nplt.yticks(10**np.arange(-16., -7))\nplt.ylim([1e-17, 1e-7])\nplt.ylabel(r'$h_{\\rho}$ (A/m)')\nplt.legend(loc=5)\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 4.7\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nradius = 50\narea = radius**2*np.pi\nfreq = np.logspace(-1, np.log10(250000), 301)\nsrc = [radius, 0, 0, 90, 0]\nrec = [0, 0, 0, 0, 90]\ndepth = 0\nres = [2e14, 100]\nstrength = area/(radius/2)\nmrec = True\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': freq, 'strength': strength, 'mrec': mrec,\n 'verb': 1}\nfhz_num = empymod.bipole(**inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nplt.plot(freq, pos(fhz_num.real), 'C0-', label='Real')\nplt.plot(freq, neg(fhz_num.real), 'C0--')\n\nplt.plot(freq, pos(fhz_num.imag), 'C1-', label='Imaginary')\nplt.plot(freq, neg(fhz_num.imag), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-1, 1e6])\nplt.ylim([1e-8, 1e-1])\nplt.xlabel('frequency (Hz)')\nplt.ylabel('$H_z$ (A/m)')\nplt.legend()\n\nplt.tight_layout()\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 4.8\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nradius = 50\narea = radius**2*np.pi\ntime = np.logspace(-7, -1, 301)\nsrc = [radius, 0, 0, 90, 0]\nrec = [0, 0, 0, 0, 90]\ndepth = 0\nres = [2e14, 100]\nstrength = area/(radius/2)\nmrec = True\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': time, 'strength': strength, 'mrec': mrec,\n 'epermH': eperm, 'verb': 1}\n\nfhz_num = empymod.bipole(signal=-1, **inp)\nfdhz_num = empymod.bipole(signal=0, **inp)\n\n# Figure\nplt.figure(figsize=(4, 6))\n\nax1 = plt.subplot(111)\nplt.plot(time*1e3, pos(fdhz_num), 'C0-', label=r'dhz/dt (A/m-s)')\nplt.plot(time*1e3, neg(fdhz_num), 'C0--')\n\nplt.plot(time*1e3, pos(fhz_num), 'C1-', label='hz (A/m)')\nplt.plot(time*1e3, neg(fhz_num), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-4, 1e2])\nplt.yticks(10**np.arange(-7., 4))\nplt.ylim([1e-7, 5e3])\n\nplt.xlabel('time (ms)')\nplt.legend()\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 2.2\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nfreq = np.logspace(-2, 5, 301)\nsrc = [0, 0, 0, 0, 0]\nrec = [0, 100, 0, 0, 0]\ndepth = []\nres = 100\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': freq, 'verb': 1}\nfhz_num = empymod.loop(**inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nplt.plot(freq, pos(fhz_num.real), 'C0-', label='Real')\nplt.plot(freq, neg(fhz_num.real), 'C0--')\n\nplt.plot(freq, pos(fhz_num.imag), 'C1-', label='Imaginary')\nplt.plot(freq, neg(fhz_num.imag), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-2, 1e5])\nplt.ylim([1e-13, 1e-6])\nplt.xlabel('frequency (Hz)')\nplt.ylabel('$H_z$ (A/m)')\nplt.legend()\n\nplt.tight_layout()\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 2.3\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nfreq = np.logspace(-2, 5, 301)\nsrc = [0, 0, 0, 0, 0]\nrec = [100, 0, 0, 0, 0]\ndepth = []\nres = 100\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'freqtime': freq, 'verb': 1}\nfhz_num = empymod.loop(**inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nplt.plot(freq, pos(fhz_num.real), 'C0-', label='Real')\nplt.plot(freq, neg(fhz_num.real), 'C0--')\n\nplt.plot(freq, pos(fhz_num.imag), 'C1-', label='Imaginary')\nplt.plot(freq, neg(fhz_num.imag), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-2, 1e5])\nplt.ylim([1e-13, 1e-6])\nplt.xlabel('Frequency (Hz)')\nplt.ylabel(r'$H_{\\rho}$ (A/m)')\nplt.legend()\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 2.4\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\ntime = np.logspace(-7, 0, 301)\nsrc = [0, 0, 0, 0, 0]\nrec = [0, 100, 0, 0, 0]\ndepth = []\nres = 100\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'xdirect': True, 'freqtime': time, 'verb': 1}\nfhz_num = empymod.loop(signal=1, **inp)\nfdhz_num = empymod.loop(signal=0, **inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nax1 = plt.subplot(111)\n\nplt.plot(time*1e3, pos(fdhz_num), 'C0-', label='dHz/dt')\nplt.plot(time*1e3, neg(fdhz_num), 'C0--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-4, 1e3])\nplt.yticks(10**np.arange(-12., -1))\nplt.ylim([1e-12, 1e-2])\nplt.xlabel('time (ms)')\nplt.ylabel(r'$\\frac{\\partial h_{\\rho}}{\\partial t}$ (A/m-s)')\nplt.legend(loc=8)\n\nax2 = ax1.twinx()\n\nplt.plot(time*1e3, pos(fhz_num), 'C1-', label='Hz')\nplt.plot(time*1e3, neg(fhz_num), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-4, 1e3])\nplt.yticks(10**np.arange(-14., -3))\nplt.ylim([1e-14, 1e-4])\nplt.ylabel(r'$h_{\\rho}$ (A/m)')\nplt.legend(loc=5)\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\n\n\n## Ward and Hohmann, 1988, Fig 2.5\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\ntime = np.logspace(-7, 0, 301)\nsrc = [0, 0, 0, 0, 0]\nrec = [100, 0, 0, 0, 0]\ndepth = []\nres = 100\n\n# Computation\ninp = {'src': src, 'rec': rec, 'depth': depth, 'res': res,\n 'xdirect': True, 'freqtime': time, 'verb': 1}\nfhz_num = empymod.loop(signal=1, **inp)\nfdhz_num = empymod.loop(signal=0, **inp)\n\n# Figure\nplt.figure(figsize=(5, 5))\n\nax1 = plt.subplot(111)\n\nplt.plot(time*1e3, pos(fdhz_num), 'C0-', label='dHz/dt')\nplt.plot(time*1e3, neg(fdhz_num), 'C0--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-4, 1e3])\nplt.yticks(10**np.arange(-12., -1))\nplt.ylim([1e-12, 1e-2])\nplt.xlabel('time (ms)')\nplt.ylabel(r'$\\frac{\\partial h_{\\rho}}{\\partial t}$ (A/m-s)')\nplt.legend(loc=8)\n\nax2 = ax1.twinx()\n\nplt.plot(time*1e3, pos(fhz_num), 'C1-', label='Hz')\nplt.plot(time*1e3, neg(fhz_num), 'C1--')\n\nplt.xscale('log')\nplt.yscale('log')\nplt.xlim([1e-4, 1e3])\nplt.yticks(10**np.arange(-16., -5))\nplt.ylim([1e-16, 1e-6])\nplt.ylabel(r'$h_{\\rho}$ (A/m)')\nplt.legend(loc=5)\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Original Figure\n\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 }