{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nDifference between magnetic dipole and loop sources\n===================================================\n\nIn this example we look at the differences between an electric loop loop, which\nresults in a magnetic source, and a magnetic dipole source.\n\nThe derivation of the electromagnetic field in Hunziker et al. (2015) is for\nelectric and magnetic point-dipole sources and receivers. The magnetic field\ndue to a magnetic source ($mm$) is obtain from the electric field due to\nan electric source ($ee$) using the duality principle, given in their\nEquation (11),\n\n\\begin{align}\\hat{G}^{mm}_{pq}(\\mathbf{x}, \\mathbf{x'}, s, \\eta_{kr}, \\zeta_{ij}) =\n -\\hat{G}^{ee}_{pq}(\\mathbf{x}, \\mathbf{x'}, s, -\\zeta_{kr}, -\\eta_{ij}) \\,\n . \\qquad (1)\\end{align}\n\nWithout going into the details of the different parameters, we can focus on the\ndifference between the $mm$ and $ee$ fields for a homogeneous,\nisotropic fullspace by simplifying this further to\n\n\\begin{align}\\mathbf{G}^{mm}_\\text{dip-dip} = \\frac{\\eta}{\\zeta}\\mathbf{G}^{ee} \\quad\n \\xrightarrow{\\text{diff. approx}} \\quad \\frac{\\sigma}{\\mathrm{i}\\omega\n \\mu}\\mathbf{G}^{ee}_\\text{dip-dip} \\, . \\qquad (2)\\end{align}\n\nHere, $\\sigma$ is conductivity (S/m), $\\omega=2\\pi f$ is angular\nfrequency (Hz), and $\\mu$ is the magnetic permeability (H/m). So from\nEquation (2) we see that the $mm$ field differs from the $ee$\nfield by a factor $\\sigma/(\\mathrm{i}\\omega\\mu)$.\n\nA magnetic dipole source has a moment of $I^mds$; however, a magnetic\ndipole source is basically never used in geophysics. Instead a loop of an\nelectric wire is used, which generates a magnetic field. The moment generated\nby this loop is given by $I^m = \\mathrm{i}\\omega\\mu N A I^e$, where\n$A$ is the area of the loop (m$^2$), and $N$ the number of\nturns of the loop. So the difference between a unit magnetic dipole and a unit\nloop ($A=1, N=1$) is the factor $\\mathrm{i}\\omega\\mu$, hence\nEquation (2) becomes\n\n\\begin{align}\\mathbf{G}^{mm}_\\text{loop-dip} =\n \\mathrm{i}\\omega\\mu\\mathbf{G}^{mm}_\\text{dip-dip} =\n \\sigma\\,\\mathbf{G}^{ee}_\\text{dip-dip} \\, . \\qquad (3)\\end{align}\n\nThis notebook shows this relation in the frequency domain, as well as for\nimpulse, step-on, and step-off responses in the time domain.\n\nWe can actually model an **electric loop** instead of adjusting the magnetic\ndipole solution to correspond to a loop source. This is shown in the second\npart of the notebook.\n\n**References**\n\n- Hunziker, J., J. Thorbecke, and E. Slob, 2015, The electromagnetic response\n in a layered vertical transverse isotropic medium: A new look at an old\n problem: Geophysics, 80(1), F1\u2013F18; DOI: `10.1190/geo2013-0411.1\n `_.\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": [ "1. Using the magnetic dipole solution\n-------------------------------------\n\nSurvey parameters\n~~~~~~~~~~~~~~~~~\n\n- Homogenous fullspace of $\\sigma$ = 0.01 S/m.\n- Source at the origin, x-directed.\n- Inline receiver with offset of 100 m, x-directed.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "freq = np.logspace(-1, 5, 301) # Frequencies (Hz)\ntime = np.logspace(-6, 0, 301) # Times (s)\nsrc = [0, 0, 0, 0, 0] # x-dir. source at the origin [x, y, z, azimuth, dip]\nrec = [100, 0, 0, 0, 0] # x-dir. receiver 100m away from source, inline\ncond = 0.01 # Conductivity (S/m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Computation using ``empymod``\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Collect common parameters\ninp = {'src': src, 'rec': rec, 'depth': [], 'res': 1/cond, 'verb': 1}\n\n# Frequency domain\ninp['freqtime'] = freq\nfee_dip_dip = empymod.bipole(**inp)\nfmm_dip_dip = empymod.bipole(msrc=True, mrec=True, **inp)\nf_loo_dip = empymod.loop(**inp)\n\n# Time domain\ninp['freqtime'] = time\n\n# ee\nee_dip_dip_of = empymod.bipole(signal=-1, **inp)\nee_dip_dip_im = empymod.bipole(signal=0, **inp)\nee_dip_dip_on = empymod.bipole(signal=1, **inp)\n\n# mm dip-dip\ndip_dip_of = empymod.bipole(signal=-1, msrc=True, mrec=True, **inp)\ndip_dip_im = empymod.bipole(signal=0, msrc=True, mrec=True, **inp)\ndip_dip_on = empymod.bipole(signal=1, msrc=True, mrec=True, **inp)\n\n# mm loop-dip\nloo_dip_of = empymod.loop(signal=-1, **inp)\nloo_dip_im = empymod.loop(signal=0, **inp)\nloo_dip_on = empymod.loop(signal=1, **inp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the result\n~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fs = 16 # Fontsize\n\n# Figure\nfig = plt.figure(figsize=(12, 8))\n\n# Frequency Domain\nplt.subplot(231)\nplt.title(r'$G^{ee}_{\\rm{dip-dip}}$', fontsize=fs)\nplt.plot(freq, fee_dip_dip.real, 'C0-', label='Real')\nplt.plot(freq, -fee_dip_dip.real, 'C0--')\nplt.plot(freq, fee_dip_dip.imag, 'C1-', label='Imag')\nplt.plot(freq, -fee_dip_dip.imag, 'C1--')\nplt.xscale('log')\nplt.yscale('log')\nplt.ylim([5e-8, 2e-5])\n\nax1 = plt.subplot(232)\nplt.title(r'$G^{mm}_{\\rm{dip-dip}}$', fontsize=fs)\nplt.plot(freq, fmm_dip_dip.real, 'C0-', label='Real')\nplt.plot(freq, -fmm_dip_dip.real, 'C0--')\nplt.plot(freq, fmm_dip_dip.imag, 'C1-', label='Imag')\nplt.plot(freq, -fmm_dip_dip.imag, 'C1--')\nplt.xscale('log')\nplt.yscale('log')\nplt.xlabel('Frequency (Hz)', fontsize=fs-2)\nplt.legend()\n\nplt.subplot(233)\nplt.title(r'$G^{mm}_{\\rm{loop-dip}}$', fontsize=fs)\nplt.plot(freq, f_loo_dip.real, 'C0-', label='Real')\nplt.plot(freq, -f_loo_dip.real, 'C0--')\nplt.plot(freq, f_loo_dip.imag, 'C1-', label='Imag')\nplt.plot(freq, -f_loo_dip.imag, 'C1--')\nplt.xscale('log')\nplt.yscale('log')\nplt.ylim([5e-10, 2e-7])\n\nplt.text(1.05, 0.5, \"Frequency Domain\", {'fontsize': fs},\n horizontalalignment='left', verticalalignment='center',\n rotation=-90, clip_on=False, transform=plt.gca().transAxes)\n\n# Time Domain\nplt.subplot(234)\nplt.plot(time, ee_dip_dip_of, 'C0-', label='Step-Off')\nplt.plot(time, -ee_dip_dip_of, 'C0--')\nplt.plot(time, ee_dip_dip_im, 'C1-', label='Impulse')\nplt.plot(time, -ee_dip_dip_im, 'C1--')\nplt.plot(time, ee_dip_dip_on, 'C2-', label='Step-On')\nplt.plot(time, -ee_dip_dip_on, 'C2--')\nplt.xscale('log')\nplt.yscale('log')\n\nplt.subplot(235)\nplt.plot(time, dip_dip_of, 'C0-', label='Step-Off')\nplt.plot(time, -dip_dip_of, 'C0--')\nplt.plot(time, dip_dip_im, 'C1-', label='Impulse')\nplt.plot(time, -dip_dip_im, 'C1--')\nplt.plot(time, dip_dip_on, 'C2-', label='Step-On')\nplt.plot(time, -dip_dip_on, 'C2--')\nplt.xscale('log')\nplt.yscale('log')\nplt.xlabel('Time (s)', fontsize=fs-2)\nplt.legend()\n\nplt.subplot(236)\nplt.plot(time, loo_dip_of, 'C0-', label='Step-Off')\nplt.plot(time, -loo_dip_of, 'C0--')\nplt.plot(time, loo_dip_im, 'C1-', label='Impulse')\nplt.plot(time, -loo_dip_im, 'C1--')\nplt.plot(time, loo_dip_on, 'C2-', label='Step-On')\nplt.plot(time, -loo_dip_on, 'C2--')\nplt.xscale('log')\nplt.yscale('log')\n\nplt.text(1.05, 0.5, \"Time Domain\", {'fontsize': fs},\n horizontalalignment='left', verticalalignment='center',\n rotation=-90, clip_on=False, transform=plt.gca().transAxes)\n\nfig.text(-0.01, 0.5, 'Amplitude; e-rec (V/m); m-rec (A/m)',\n va='center', rotation='vertical', fontsize=fs, color='.4')\n\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure shows the main points of Equations (2) and (3):\n\n- The magnetic dipole-dipole response differs by a factor\n $\\sigma/(\\mathrm{i}\\omega\\mu)$ from the electric dipole-dipole\n response. That means for the time-domain that the magnetic response looks\n more like the time derivative of the electric response (e.g., the magnetic\n impulse responses resembles the electric step-on response).\n- The magnetic loop-dipole response differs only by $\\sigma$ from the\n electric dipole-dipole response, hence a factor of 0.01.\n\nThe units of the response only depend on the receiver, what the receiver\nactually measures. So if we change the source from a dipole to a loop it does\nnot change the units of the received responses.\n\n2. Using an electric loop\n-------------------------\n\nWe can use ``empymod`` to model arbitrary shaped sources by simply adding\npoint dipole sources together. This is what ``empymod`` does internally to\nmodel a finite length dipole (``empymod.bipole``), where it uses a Gaussian\nquadrature with a few points.\n\nHere, we are going to compare the result from ``loop``, as presented above,\nwith two different simulations of an electric loop source, assuming a square\nloop which sides are 1 m long, so the area correspond to one square meter.\n\nPlotting routines\n~~~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def discrete_cmap(N, base_cmap=None):\n \"\"\"Create an N-bin discrete colormap from the specified input map\n https://gist.github.com/jakevdp/91077b0cae40f8f8244a\n \"\"\"\n base = plt.cm.get_cmap(base_cmap)\n color_list = base(np.linspace(0, 1, N))\n cmap_name = base.name + str(N)\n return base.from_list(cmap_name, color_list, N)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_result(data1, data2, x, title, vmin=-15., vmax=-7., rx=0):\n \"\"\"Plot result.\"\"\"\n fig = plt.figure(figsize=(18, 10))\n\n def setplot(name):\n \"\"\"Plot settings\"\"\"\n plt.title(name)\n plt.xlim(rx.min(), rx.max())\n plt.ylim(rx.min(), rx.max())\n plt.axis(\"equal\")\n\n # Plot Re(data)\n ax1 = plt.subplot(231)\n setplot(r\"(a) |Re(magn.dip*iwu)|\")\n cf0 = plt.pcolormesh(rx, rx, np.log10(np.abs(data1.real)), linewidth=0,\n rasterized=True, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n\n ax2 = plt.subplot(232)\n setplot(r\"(b) |Re(el. square)|\")\n plt.pcolormesh(rx, rx, np.log10(np.abs(data2.real)), linewidth=0,\n rasterized=True, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n\n ax3 = plt.subplot(233)\n setplot(r\"(c) Error real part\")\n error_r = np.abs((data1.real-data2.real)/data1.real)*100\n cf2 = plt.pcolormesh(rx, rx, np.log10(error_r), vmin=-2, vmax=2,\n linewidth=0, rasterized=True,\n cmap=discrete_cmap(8, \"RdBu_r\"))\n\n # Plot Im(data)\n ax4 = plt.subplot(234)\n setplot(r\"(d) |Im(magn.dip*iwu)|\")\n plt.pcolormesh(rx, rx, np.log10(np.abs(data1.imag)), linewidth=0,\n rasterized=True, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n\n ax5 = plt.subplot(235)\n setplot(r\"(e) |Im(el. square)|\")\n plt.pcolormesh(rx, rx, np.log10(np.abs(data2.imag)), linewidth=0,\n rasterized=True, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n\n ax6 = plt.subplot(236)\n setplot(r\"(f) Error imag part\")\n error_i = np.abs((data1.imag-data2.imag)/data1.imag)*100\n plt.pcolormesh(rx, rx, np.log10(error_i), vmin=-2, vmax=2,\n linewidth=0, rasterized=True,\n cmap=discrete_cmap(8, \"RdBu_r\"))\n\n # Colorbars\n fig.colorbar(cf0, ax=[ax1, ax2, ax3], label=r\"$\\log_{10}$ Amplitude (A/m)\")\n cbar = fig.colorbar(cf2, ax=[ax4, ax5, ax6], label=r\"Relative Error\")\n cbar.set_ticks([-2, -1, 0, 1, 2])\n cbar.ax.set_yticklabels([r\"$0.01\\,\\%$\", r\"$0.1\\,\\%$\", r\"$1\\,\\%$\",\n r\"$10\\,\\%$\", r\"$100\\,\\%$\"])\n\n # Axis label\n fig.text(0.4, 0.05, \"Inline Offset (m)\", fontsize=14)\n fig.text(0.08, 0.5, 'Crossline Offset (m)', rotation=90, fontsize=14)\n\n # Title\n fig.suptitle(title, y=.95, fontsize=20)\n plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Model parameters\n~~~~~~~~~~~~~~~~\n\n- Resistivity: $1 \\Omega$ m fullspace\n\nSurvey\n~~~~~~\n\n- Source at [0, 0, 0]\n- Receivers at [x, y, 10]\n- frequencies: 100 Hz.\n- Offsets: -250 m - 250 m\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Survey parameters\nx = ((np.arange(502))-250.5)\nrx = np.repeat([x, ], np.size(x), axis=0)\nry = rx.transpose()\nrxx = rx.ravel()\nryy = ry.ravel()\n\n# Model\nmodel = {\n 'depth': [], # Fullspace\n 'res': 1., # 1 Ohm.m\n 'freqtime': 100, # 100 Hz\n 'htarg': {'pts_per_dec': -1},\n 'verb': 1,\n}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute ``empymod.loop`` result\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "epm_loop = empymod.loop(src=[0, 0, 0, 0, 90], rec=[rxx, ryy, 10, 0, 0],\n **model).reshape(np.shape(rx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2.1 Point dipoles at (x, y) using ``empymod.dipole``\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n - (0.5, 0), ab=42\n - (0, 0.5), ab=41\n - (-0.5, 0), ab=-42\n - (0, -0.5), ab=-41\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rec_dip = [rxx, ryy, 10]\n\nsquare_pts = +empymod.dipole(src=[+0.5, +0.0, 0], rec=rec_dip, ab=42,\n **model).reshape(np.shape(rx))\nsquare_pts += empymod.dipole(src=[+0.0, +0.5, 0], rec=rec_dip, ab=41,\n **model).reshape(np.shape(rx))\nsquare_pts -= empymod.dipole(src=[-0.5, +0.0, 0], rec=rec_dip, ab=42,\n **model).reshape(np.shape(rx))\nsquare_pts -= empymod.dipole(src=[+0.0, -0.5, 0], rec=rec_dip, ab=41,\n **model).reshape(np.shape(rx))\n\nplot_result(epm_loop, square_pts, x, 'Loop made of four points',\n vmin=-13, vmax=-5, rx=x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2.2 Finite length dipoles using ``empymod.bipole``\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\nEach simulated with a 5pt Gaussian quadrature. The dipoles are:\n\n - (-0.5, -0.5) to (+0.5, -0.5)\n - (+0.5, -0.5) to (+0.5, +0.5)\n - (+0.5, +0.5) to (-0.5, +0.5)\n - (-0.5, +0.5) to (-0.5, -0.5)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inp_dip = {\n 'rec': [rxx, ryy, 10, 0, 0],\n 'mrec': True,\n 'srcpts': 5 # Gaussian quadr. with 5 pts to simulate a finite length dip.\n}\n\nsquare_dip = +empymod.bipole(src=[+0.5, +0.5, -0.5, +0.5, 0, 0],\n **inp_dip, **model)\nsquare_dip += empymod.bipole(src=[+0.5, -0.5, +0.5, +0.5, 0, 0],\n **inp_dip, **model)\nsquare_dip += empymod.bipole(src=[-0.5, -0.5, +0.5, -0.5, 0, 0],\n **inp_dip, **model)\nsquare_dip += empymod.bipole(src=[-0.5, +0.5, -0.5, -0.5, 0, 0],\n **inp_dip, **model)\nsquare_dip = square_dip.reshape(np.shape(rx))\n\nplot_result(epm_loop, square_dip, x, 'Loop made of four dipoles',\n vmin=-13, vmax=-5, rx=x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Close to the source the results between\n\n- (1) a magnetic dipole,\n- (2) an electric loop conisting of four point sources, and\n- (3) an electric loop consisting of four finite length dipoles,\n\ndiffer, as expected. However, for the vast majority they are identical. Skin\ndepth for our example with $\\rho=1\\Omega$ m and $f=100\\,$ Hz is\nroughly 50 m, so the results are basically identical for 4-5 skin depths,\nafter which the signal is very low.\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.7.3" } }, "nbformat": 4, "nbformat_minor": 0 }