.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/geophysical_em.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_geophysical_em.py: Geophysical Electromagnetic modelling ===================================== In this example we use `pyfftlog` to obtain time-domain EM data from frequency-domain data and vice versa. We do this by using analytical halfspace solution in both domains, and comparing the transformed responses to the true result. The analytical halfspace solutions are computed using `empymod` (see https://empymod.github.io). .. GENERATED FROM PYTHON SOURCE LINES 11-18 .. code-block:: Python import empymod import pyfftlog import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import InterpolatedUnivariateSpline as iuSpline .. GENERATED FROM PYTHON SOURCE LINES 19-21 Model and Survey parameters --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 21-63 .. code-block:: Python # Impulse response (in the time domain) signal = 0 # x-directed electric source and receiver point-dipoles ab = 11 # We use the same range of times (s) and frequencies (Hz) ftpts = np.logspace(-4, 4, 301) # Source and receiver src = [0, 0, 100] # At the origin, 100 m below surface rec = [6000, 0, 200] # At an inline offset of 6 km, 200 m below surface # Resistivity depth = [0] # Interface at z = 0, default for empymod.analytical res = [2e14, 1] # Horizontal resistivity [air, subsurface] aniso = [1, 2] # Anisotropy [air, subsurface] # Collect parameters analytical = { 'src': src, 'rec': rec, 'res': res[1], 'aniso': aniso[1], 'solution': 'dhs', # Diffusive half-space solution 'verb': 2, 'ab': ab, } dipole = { 'src': src, 'rec': rec, 'depth': depth, 'res': res, 'aniso': aniso, 'ht': 'dlf', 'verb': 2, 'ab': ab, } .. GENERATED FROM PYTHON SOURCE LINES 64-66 Analytical solutions -------------------- .. GENERATED FROM PYTHON SOURCE LINES 66-74 .. code-block:: Python # Frequency Domain f_ana = empymod.analytical(**analytical, freqtime=ftpts) # Time Domain t_ana = empymod.analytical(**analytical, freqtime=ftpts, signal=signal) .. rst-class:: sphx-glr-script-out .. code-block:: none :: empymod END; runtime = 0:00:00.002036 :: :: empymod END; runtime = 0:00:00.001036 :: .. GENERATED FROM PYTHON SOURCE LINES 75-77 FFTLog ------ .. GENERATED FROM PYTHON SOURCE LINES 77-141 .. code-block:: Python # FFTLog parameters pts_per_dec = 5 # Increase if not precise enough add_dec = [-2, 2] # e.g. [-2, 2] to add 2 decades on each side q = 0 # -1 - +1; can improve results # Compute minimum and maximum required inputs rmin = np.log10(1/ftpts.max()) + add_dec[0] rmax = np.log10(1/ftpts.min()) + add_dec[1] n = int(rmax - rmin)*pts_per_dec # Pre-allocate output f_resp = np.zeros(ftpts.shape, dtype=complex) # Loop over Sine, Cosine transform. for mu in [0.5, -0.5]: # Central point log10(r_c) of periodic interval logrc = (rmin + rmax)/2 # Central index (1/2 integral if n is even) nc = (n + 1)/2. # Log spacing of points dlogr = (rmax - rmin)/n dlnr = dlogr*np.log(10.) # Compute required input x-values pts_req = 10**(logrc + (np.arange(1, n+1) - nc)*dlogr)/2/np.pi # Initialize FFTLog kr, xsave = pyfftlog.fhti(n, mu, dlnr, q, kr=1, kropt=1) # Compute pts_out with adjusted kr logkc = np.log10(kr) - logrc pts_out = 10**(logkc + (np.arange(1, n+1) - nc)*dlogr) # rk = r_c/k_r; adjust for Fourier transform scaling rk = 10**(logrc - logkc)*np.pi/2 # Compute required times/frequencies with the analytical solution t2f_t_resp = empymod.analytical(**analytical, freqtime=pts_req, signal=signal) f2t_f_resp = empymod.analytical(**analytical, freqtime=pts_req) # Carry out FFTLog t2f_f_coarse = pyfftlog.fftl(t2f_t_resp, xsave.copy(), rk, 1) if mu > 0: f2t_t_coarse = pyfftlog.fftl(f2t_f_resp.imag, xsave.copy(), rk, 1) else: f2t_t_coarse = pyfftlog.fftl(f2t_f_resp.real, xsave.copy(), rk, 1) # Interpolate for required frequencies/times t2f_f_spline = iuSpline(np.log(pts_out), t2f_f_coarse) f2t_t_spline = iuSpline(np.log(pts_out), f2t_t_coarse) if mu > 0: f_resp += -1j*t2f_f_spline(np.log(ftpts))/np.pi/2 t_resp_sin = -f2t_t_spline(np.log(ftpts))/np.pi*2 else: f_resp += t2f_f_spline(np.log(ftpts))/np.pi/2 t_resp_cos = f2t_t_spline(np.log(ftpts))/np.pi*2 .. rst-class:: sphx-glr-script-out .. code-block:: none :: empymod END; runtime = 0:00:00.000764 :: :: empymod END; runtime = 0:00:00.000989 :: :: empymod END; runtime = 0:00:00.000751 :: :: empymod END; runtime = 0:00:00.000958 :: .. GENERATED FROM PYTHON SOURCE LINES 142-144 Comparison ---------- .. GENERATED FROM PYTHON SOURCE LINES 144-172 .. code-block:: Python fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(9, 4)) # TIME DOMAIN ax0.set_title(r'Frequency domain') ax0.set_xlabel('Frequency (Hz)') ax0.set_ylabel('Amplitude (V/m)') ax0.semilogx(ftpts, f_ana.real, 'k-', label='Analytical') ax0.semilogx(ftpts, f_ana.imag, 'k-') ax0.semilogx(ftpts, f_resp.real, 'C3--', label=r'FFTLog, $\mu=-0.5$') ax0.semilogx(ftpts, f_resp.imag, 'C2--', label=r'FFTLog, $\mu=+0.5$') ax0.legend(loc='best') ax0.grid(which='both', c='.95') # TIME DOMAIN ax1.set_title(r'Time domain') ax1.set_xlabel('Time (s)') ax1.set_ylabel('Amplitude (V/m)') ax1.semilogx(ftpts, t_ana, 'k', label='Analytical') ax1.semilogx(ftpts, t_resp_cos, 'C3--', label=r'FFTLog, $\mu=-0.5$') ax1.semilogx(ftpts, t_resp_sin, 'C2-.', label=r'FFTLog, $\mu=+0.5$') ax1.legend(loc='best') ax1.yaxis.set_label_position("right") ax1.yaxis.tick_right() ax1.grid(which='both', c='.95') fig.tight_layout() fig.show() .. image-sg:: /examples/images/sphx_glr_geophysical_em_001.png :alt: Frequency domain, Time domain :srcset: /examples/images/sphx_glr_geophysical_em_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.956 seconds) .. _sphx_glr_download_examples_geophysical_em.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: geophysical_em.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: geophysical_em.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: geophysical_em.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_