Main modelling routines¶

Electromagnetic modeller to model electric or magnetic responses due to a three-dimensional electric or magnetic source in a layered-earth model with vertical transverse isotropic (VTI) resistivity, VTI electric permittivity, and VTI magnetic permeability, from very low frequencies (DC) to very high frequencies (GPR). The calculation is carried out in the wavenumber-frequency domain, and various Hankel- and Fourier-transform methods are included to transform the responses into the space-frequency and space-time domains.

`empymod.model` – Model EM-responses¶

EM-modelling routines. The implemented routines might not be the fastest solution to your specific problem. Use these routines as template to create your own, problem-specific modelling routine!

Principal routines:

bipole()
dipole()
loop()

The main routine is bipole(), which can model bipole source(s) and bipole receiver(s) of arbitrary direction, for electric or magnetic sources and receivers, both in frequency and in time. A subset of bipole() is dipole(), which models infinitesimal small dipoles along the principal axes x, y, and z. The third routine, loop(), can be used if the source or the receivers are loops instead of dipoles.

Further routines are:

analytical(): Calculate analytical fullspace and halfspace solutions.
dipole_k(): Calculate the electromagnetic wavenumber-domain solution.
gpr(): Calculate the Ground-Penetrating Radar (GPR) response.

The dipole_k() routine can be used if you are interested in the wavenumber-domain result, without Hankel nor Fourier transform. It calls straight the empymod.kernel. The gpr()-routine convolves the frequency-domain result with a wavelet, and applies a gain to the time-domain result. This function is still experimental.

The modelling routines make use of the following two core routines:

fem(): Calculate wavenumber-domain electromagnetic field and carry out the Hankel transform to the frequency domain.
tem(): Carry out the Fourier transform to time domain after fem().

empymod.model.bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, msrc=False, srcpts=1, mrec=False, recpts=1, strength=0, **kwargs)[source]¶

Return EM fields due to arbitrary rotated, finite length EM dipoles.

Calculate the electromagnetic frequency- or time-domain field due to arbitrary rotated, finite electric or magnetic bipole sources, measured by arbitrary rotated, finite electric or magnetic bipole receivers. By default, the electromagnetic response is normalized to source and receiver of 1 m length, and source strength of 1 A.

Parameters:

src, rec : list of floats or arrays

Source and receiver coordinates (m):

[x0, x1, y0, y1, z0, z1] (bipole of finite length)
[x, y, z, azimuth, dip] (dipole, infinitesimal small)

Dimensions:

The coordinates x, y, and z (dipole) or x0, x1, y0, y1, z0, and z1 (bipole) can be single values or arrays.
The variables x and y (dipole) or x0, x1, y0, and y1 (bipole) must have the same dimensions.
The variables z, azimuth, and dip (dipole) or z0 and z1 (bipole) must either be single values or having the same dimension as the other coordinates.

Angles (coordinate system is either left-handed with positive z down or right-handed with positive z up; East-North-Depth):

azimuth (°): horizontal deviation from x-axis, anti-clockwise.
+/-dip (°): vertical deviation from xy-plane down/up-wards.

Sources or receivers placed on a layer interface are considered in the upper layer.

depth : list

Absolute layer interfaces z (m); #depth = #res - 1 (excluding +/- infinity).

res : array_like

Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.

Alternatively, res can be a dictionary. See the main manual of empymod too see how to exploit this hook to re-calculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the Cole-Cole model for IP.

freqtime : array_like

Frequencies f (Hz) if signal==None, else times t (s); (f, t > 0).

signal : {None, 0, 1, -1}, optional

Source signal, default is None:

None: Frequency-domain response
-1 : Switch-off time-domain response
0 : Impulse time-domain response
+1 : Switch-on time-domain response

aniso : array_like, optional

Anisotropies lambda = sqrt(rho_v/rho_h) (-); #aniso = #res. Defaults to ones.

epermH, epermV : array_like, optional

Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (-); #epermH = #epermV = #res. If epermH is provided but not epermV, isotropic behaviour is assumed. Default is ones.

mpermH, mpermV : array_like, optional

Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (-); #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed. Default is ones.

msrc, mrec : bool, optional

If True, source/receiver (msrc/mrec) is magnetic, else electric. Default is False.

srcpts, recpts : int, optional

Number of integration points for bipole source/receiver, default is 1:

srcpts/recpts < 3 : bipole, but calculated as dipole at centre
srcpts/recpts >= 3 : bipole

strength : float, optional

Source strength (A):

If 0, output is normalized to source and receiver of 1 m length, and source strength of 1 A.
If != 0, output is returned for given source and receiver length, and source strength.

Default is 0.

verb : {0, 1, 2, 3, 4}, optional

Level of verbosity, default is 2:

0: Print nothing.
1: Print warnings.
2: Print additional runtime and kernel calls
3: Print additional start/stop, condensed parameter information.
4: Print additional full parameter information

ht : {‘dlf’, ‘qwe’, ‘quad’}, optional

Flag to choose either the Digital Linear Filter (DLF) method, the Quadrature-With-Extrapolation (QWE), or a simple Quadrature (QUAD) for the Hankel transform. Defaults to ‘dlf’.

htarg : dict, optional

Possible parameters depends on the value for ht:

If ht=’dlf’:
- dlf: string of filter name in empymod.filters or the filter method itself. (default: empymod.filters.key_201_2009())
- pts_per_dec: points per decade; (default: 0):
  - If 0: Standard DLF.
  - If < 0: Lagged Convolution DLF.
  - If > 0: Splined DLF
If ht=’qwe’:
- rtol: relative tolerance (default: 1e-12)
- atol: absolute tolerance (default: 1e-30)
- nquad: order of Gaussian quadrature (default: 51)
- maxint: maximum number of partial integral intervals (default: 40)
- pts_per_dec: points per decade; (default: 0)
  - If 0, no interpolation is used.
  - If > 0, interpolation is used.
- diff_quad: criteria when to swap to QUAD (only relevant if pts_per_dec=-1) (default: 100)
- a: lower limit for QUAD (default: first interval from QWE)
- b: upper limit for QUAD (default: last interval from QWE)
- limit: limit for quad (default: maxint)
If ht=’quad’:
- rtol: relative tolerance (default: 1e-12)
- atol: absolute tolerance (default: 1e-20)
- limit: An upper bound on the number of subintervals used in the adaptive algorithm (default: 500)
- lmin: Minimum wavenumber (default 1e-6)
- lmax: Maximum wavenumber (default 0.1)
- pts_per_dec: points per decade (default: 40)

ft : {‘dlf’, ‘sin’, ‘cos’, ‘qwe’, ‘fftlog’, ‘fft’}, optional

Only used if signal!=None. Flag to choose either the Digital Linear Filter method (Sine- or Cosine-Filter), the Quadrature-With-Extrapolation (QWE), the FFTLog, or the FFT for the Fourier transform. Defaults to ‘dlf’ (which is ‘sin’ if signal>=0, else ‘cos’).

ftarg : dict, optional

Only used if signal!=None. Possible parameters depends on the value for ft:

If ft=’dlf’, ‘sin’, or ‘cos’:
- dlf: string of filter name in empymod.filters or the filter method itself. (Default: empymod.filters.key_201_CosSin_2012())
- pts_per_dec: points per decade; (default: -1)
  - If 0: Standard DLF.
  - If < 0: Lagged Convolution DLF.
  - If > 0: Splined DLF
If ft=’qwe’:
- rtol: relative tolerance (default: 1e-8)
- atol: absolute tolerance (default: 1e-20)
- nquad: order of Gaussian quadrature (default: 21)
- maxint: maximum number of partial integral intervals (default: 200)
- pts_per_dec: points per decade (default: 20)
- diff_quad: criteria when to swap to QUAD (default: 100)
- a: lower limit for QUAD (default: first interval from QWE)
- b: upper limit for QUAD (default: last interval from QWE)
- limit: limit for quad (default: maxint)
If ft=’fftlog’:
- pts_per_dec: sampels per decade (default: 10)
- add_dec: additional decades [left, right] (default: [-2, 1])
- q: exponent of power law bias (default: 0); -1 <= q <= 1
If ft=’fft’:
- dfreq: Linear step-size of frequencies (default: 0.002)
- nfreq: Number of frequencies (default: 2048)
- ntot: Total number for FFT; difference between nfreq and ntot is padded with zeroes. This number is ideally a power of 2, e.g. 2048 or 4096 (default: nfreq).
- pts_per_dec: points per decade (default: None)
Padding can sometimes improve the result, not always. The default samples from 0.002 Hz - 4.096 Hz. If pts_per_dec is set to an integer, calculated frequencies are logarithmically spaced with the given number per decade, and then interpolated to yield the required frequencies for the FFT.

xdirect : bool or None, optional

Direct field calculation (only if src and rec are in the same layer):

If True, direct field is calculated analytically in the frequency domain.
If False, direct field is calculated in the wavenumber domain.
If None, direct field is excluded from the calculation, and only reflected fields are returned (secondary field).

Defaults to False.

loop : {None, ‘freq’, ‘off’}, optional

Define if to calculate everything vectorized or if to loop over frequencies (‘freq’) or over offsets (‘off’), default is None. It always loops over frequencies if ht=’qwe’ or if pts_per_dec=-1. Calculating everything vectorized is fast for few offsets OR for few frequencies. However, if you calculate many frequencies for many offsets, it might be faster to loop over frequencies. Only comparing the different versions will yield the answer for your specific problem at hand!

Returns:

EM : EMArray, (nfreqtime, nrec, nsrc)

Frequency- or time-domain EM field (depending on signal):

If rec is electric, returns E [V/m].
If rec is magnetic, returns H [A/m].

EMArray is a subclassed ndarray with .pha and .amp attributes (only relevant for frequency-domain data).

The shape of EM is (nfreqtime, nrec, nsrc). However, single dimensions are removed.

See also

dipole(): EM fields due to infinitesimal small EM dipoles.
loop(): EM fields due to a magnetic source loop.

Examples

>>> import empymod
>>> import numpy as np
>>> # x-directed bipole source: x0, x1, y0, y1, z0, z1
>>> src = [-50, 50, 0, 0, 100, 100]
>>> # x-directed dipole receiver-array: x, y, z, azimuth, dip
>>> rec = [np.arange(1, 11)*500, np.zeros(10), 200, 0, 0]
>>> # layer boundaries
>>> depth = [0, 300, 1000, 1050]
>>> # layer resistivities
>>> res = [1e20, .3, 1, 50, 1]
>>> # Frequency
>>> freq = 1
>>> # Calculate electric field due to an electric source at 1 Hz.
>>> # [msrc = mrec = False (default)]
>>> EMfield = empymod.bipole(src, rec, depth, res, freq, verb=4)
~
:: empymod START  ::  v2.0.0
~
   depth       [m] :  0 300 1000 1050
   res     [Ohm.m] :  1E+20 0.3 1 50 1
   aniso       [-] :  1 1 1 1 1
   epermH      [-] :  1 1 1 1 1
   epermV      [-] :  1 1 1 1 1
   mpermH      [-] :  1 1 1 1 1
   mpermV      [-] :  1 1 1 1 1
   direct field    :  Comp. in wavenumber domain
   frequency  [Hz] :  1
   Hankel          :  DLF (Fast Hankel Transform)
     > Filter      :  Key 201 (2009)
     > DLF type    :  Standard
   Loop over       :  None (all vectorized)
   Source(s)       :  1 bipole(s)
     > intpts      :  1 (as dipole)
     > length  [m] :  100
     > strength[A] :  0
     > x_c     [m] :  0
     > y_c     [m] :  0
     > z_c     [m] :  100
     > azimuth [°] :  0
     > dip     [°] :  0
   Receiver(s)     :  10 dipole(s)
     > x       [m] :  500 - 5000 : 10  [min-max; #]
                   :  500 1000 1500 2000 2500 3000 3500 4000 4500 5000
     > y       [m] :  0 - 0 : 10  [min-max; #]
                   :  0 0 0 0 0 0 0 0 0 0
     > z       [m] :  200
     > azimuth [°] :  0
     > dip     [°] :  0
   Required ab's   :  11
~
:: empymod END; runtime = 0:00:00.005536 :: 1 kernel call(s)
~
>>> print(EMfield)
[  1.68809346e-10 -3.08303130e-10j  -8.77189179e-12 -3.76920235e-11j
  -3.46654704e-12 -4.87133683e-12j  -3.60159726e-13 -1.12434417e-12j
   1.87807271e-13 -6.21669759e-13j   1.97200208e-13 -4.38210489e-13j
   1.44134842e-13 -3.17505260e-13j   9.92770406e-14 -2.33950871e-13j
   6.75287598e-14 -1.74922886e-13j   4.62724887e-14 -1.32266600e-13j]

empymod.model.dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, **kwargs)[source]¶

Return EM fields due to infinitesimal small EM dipoles.

Calculate the electromagnetic frequency- or time-domain field due to infinitesimal small electric or magnetic dipole source(s), measured by infinitesimal small electric or magnetic dipole receiver(s); sources and receivers are directed along the principal directions x, y, or z, and all sources are at the same depth, as well as all receivers are at the same depth.

Use the functions bipole() to calculate dipoles with arbitrary angles or bipoles of finite length and arbitrary angle.

The function dipole() could be replaced by bipole() (all there is to do is translate ab into msrc, mrec, azimuth’s and dip’s). However, dipole() is kept separately to serve as an example of a simple modelling routine that can serve as a template.

Parameters:

src, rec : list of floats or arrays

Source and receiver coordinates [x, y, z] (m):

The x- and y-coordinates can be arrays, z is a single value.
The x- and y-coordinates must have the same dimension.

Sources or receivers placed on a layer interface are considered in the upper layer.

depth : list

Absolute layer interfaces z (m); #depth = #res - 1 (excluding +/- infinity).

res : array_like

Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.

Alternatively, res can be a dictionary. See the main manual of empymod too see how to exploit this hook to re-calculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the Cole-Cole model for IP.

freqtime : array_like

Frequencies f (Hz) if signal==None, else times t (s); (f, t > 0).

signal : {None, 0, 1, -1}, optional

Source signal, default is None:

None: Frequency-domain response
-1 : Switch-off time-domain response
0 : Impulse time-domain response
+1 : Switch-on time-domain response

ab : int, optional

Source-receiver configuration, defaults to 11.

		electric source			magnetic source
		x	y	z	x	y	z
electric receiver	x	11	12	13	14	15	16
	y	21	22	23	24	25	26
	z	31	32	33	34	35	36
magnetic receiver	x	41	42	43	44	45	46
	y	51	52	53	54	55	56
	z	61	62	63	64	65	66

aniso : array_like, optional

Anisotropies lambda = sqrt(rho_v/rho_h) (-); #aniso = #res. Defaults to ones.

epermH, epermV : array_like, optional

Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (-); #epermH = #epermV = #res. If epermH is provided but not epermV, isotropic behaviour is assumed. Default is ones.

mpermH, mpermV : array_like, optional

Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (-); #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed. Default is ones.

verb : {0, 1, 2, 3, 4}, optional

Level of verbosity, default is 2:

0: Print nothing.
1: Print warnings.
2: Print additional runtime and kernel calls
3: Print additional start/stop, condensed parameter information.
4: Print additional full parameter information

ht, htarg, ft, ftarg, xdirect, loop : settings, optinal

See docstring of bipole() for a description.

Returns:

EM : EMArray, (nfreqtime, nrec, nsrc)

Frequency- or time-domain EM field (depending on signal):

If rec is electric, returns E [V/m].
If rec is magnetic, returns H [A/m].

EMArray is a subclassed ndarray with .pha and .amp attributes (only relevant for frequency-domain data).

The shape of EM is (nfreqtime, nrec, nsrc). However, single dimensions are removed.

See also

bipole(): EM fields due to arbitrary rotated, finite length EM dipoles.
loop(): EM fields due to a magnetic source loop.

Examples

>>> import empymod
>>> import numpy as np
>>> src = [0, 0, 100]
>>> rec = [np.arange(1, 11)*500, np.zeros(10), 200]
>>> depth = [0, 300, 1000, 1050]
>>> res = [1e20, .3, 1, 50, 1]
>>> EMfield = empymod.dipole(src, rec, depth, res, freqtime=1, verb=0)
>>> print(EMfield)
[  1.68809346e-10 -3.08303130e-10j  -8.77189179e-12 -3.76920235e-11j
  -3.46654704e-12 -4.87133683e-12j  -3.60159726e-13 -1.12434417e-12j
   1.87807271e-13 -6.21669759e-13j   1.97200208e-13 -4.38210489e-13j
   1.44134842e-13 -3.17505260e-13j   9.92770406e-14 -2.33950871e-13j
   6.75287598e-14 -1.74922886e-13j   4.62724887e-14 -1.32266600e-13j]

empymod.model.loop(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, mrec=True, recpts=1, strength=0, **kwargs)[source]¶

Return EM fields due to a magnetic source loop.

Calculate the electromagnetic frequency- or time-domain field due to an arbitrary rotated, magnetic source consisting of an electric loop, measured by arbitrary rotated, finite electric or magnetic bipole receivers or arbitrary rotated magnetic receivers consisting of electric loops. By default, the electromagnetic response is normalized to source loop area of 1 m2 and receiver length or area of 1 m or 1 m2, respectively, and source strength of 1 A.

A magnetic dipole, as used in dipole() and bipole(), has a moment of \(I^m ds\). However, if the magnetic dipole is generated by an electric-wire loop, this changes to \(I^m = i\omega\mu A I^e\), where A is the area of the loop. The same factor \(i\omega\mu A\), applies to the receiver, if it consists of an electric-wire loop.

The current implementation only handles loop sources and receivers in layers where \(\mu_r^h=\mu_r^v\); the horizontal magnetic permeability is used, and a warning is thrown if the vertical differs from the horizontal one.

Note that the kernel internally still calculates dipole sources and receivers, the moment is a factor that is multiplied in the frequency domain. The logs will therefore still inform about bipoles and dipoles.

Parameters:

src, rec : list of floats or arrays

Source and receiver coordinates (m):

[x0, x1, y0, y1, z0, z1] (bipole of finite length)
[x, y, z, azimuth, dip] (dipole, infinitesimal small)

Dimensions:

The coordinates x, y, and z (dipole) or x0, x1, y0, y1, z0, and z1 (bipole) can be single values or arrays.
The variables x and y (dipole) or x0, x1, y0, and y1 (bipole) must have the same dimensions.
The variables z, azimuth, and dip (dipole) or z0 and z1 (bipole) must either be single values or having the same dimension as the other coordinates.

Angles (coordinate system is either left-handed with positive z down or right-handed with positive z up; East-North-Depth):

azimuth (°): horizontal deviation from x-axis, anti-clockwise.
+/-dip (°): vertical deviation from xy-plane down/up-wards.

Sources or receivers placed on a layer interface are considered in the upper layer.

depth : list

Absolute layer interfaces z (m); #depth = #res - 1 (excluding +/- infinity).

res : array_like

Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.

Alternatively, res can be a dictionary. See the main manual of empymod too see how to exploit this hook to re-calculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the Cole-Cole model for IP.

freqtime : array_like

Frequencies f (Hz) if signal==None, else times t (s); (f, t > 0).

signal : {None, 0, 1, -1}, optional

Source signal, default is None:

None: Frequency-domain response
-1 : Switch-off time-domain response
0 : Impulse time-domain response
+1 : Switch-on time-domain response

aniso : array_like, optional

Anisotropies lambda = sqrt(rho_v/rho_h) (-); #aniso = #res. Defaults to ones.

epermH, epermV : array_like, optional

Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (-); #epermH = #epermV = #res. If epermH is provided but not epermV, isotropic behaviour is assumed. Default is ones.

mpermH, mpermV : array_like, optional

Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (-); #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed. Default is ones.

Note that the relative horizontal and vertical magnetic permeabilities in layers with loop sources or receivers will be set to 1.

mrec : bool or string, optional

Receiver options; default is True:

True: Magnetic dipole receiver;
False: Electric dipole receiver;
‘loop’: Magnetic receiver consisting of an electric-wire loop.

recpts : int, optional

Number of integration points for bipole receiver, default is 1:

recpts < 3 : bipole, but calculated as dipole at centre
recpts >= 3 : bipole

Note that if mrec=’loop’, recpts will be set to 1.

strength : float, optional

Source strength (A):

If 0, output is normalized to source of 1 m2 area and receiver of 1 m length or 1 m2 area, and source strength of 1 A.
If != 0, output is returned for given source strength and receiver length (if mrec!=’loop’).

The strength is simply a multiplication factor. It can also be used to provide the source and receiver loop area, or also to multiply by :math:mu_0`, if you want the B-field instead of the H-field.

Default is 0.

verb : {0, 1, 2, 3, 4}, optional

Level of verbosity, default is 2:

0: Print nothing.
1: Print warnings.
2: Print additional runtime and kernel calls
3: Print additional start/stop, condensed parameter information.
4: Print additional full parameter information

ht, htarg, ft, ftarg, xdirect, loop : settings, optinal

See docstring of bipole() for a description.

Returns:

EM : EMArray, (nfreqtime, nrec, nsrc)

Frequency- or time-domain EM field (depending on signal):

If rec is electric, returns E [V/m].
If rec is magnetic, returns H [A/m].

EMArray is a subclassed ndarray with .pha and .amp attributes (only relevant for frequency-domain data).

See also

dipole(): EM fields due to infinitesimal small EM dipoles.
bipole(): EM fields due to arbitrary rotated, finite length EM dipoles.

Examples

>>> import empymod
>>> import numpy as np
>>> # z-directed loop source: x, y, z, azimuth, dip
>>> src = [0, 0, 0, 0, 90]
>>> # z-directed magnetic dipole receiver-array: x, y, z, azimuth, dip
>>> rec = [np.arange(1, 11)*500, np.zeros(10), 200, 0, 90]
>>> # layer boundaries
>>> depth = [0, 300, 500]
>>> # layer resistivities
>>> res = [2e14, 10, 500, 10]
>>> # Frequency
>>> freq = 1
>>> # Calculate magnetic field due to a loop source at 1 Hz.
>>> # [mrec = True (default)]
>>> EMfield = empymod.loop(src, rec, depth, res, freq, verb=4)
~
:: empymod START  ::  w2.0.0
~
   depth       [m] :  0 300 500
   res     [Ohm.m] :  2E+14 10 500 10
   aniso       [-] :  1 1 1 1
   epermH      [-] :  1 1 1 1
   epermV      [-] :  1 1 1 1
   mpermH      [-] :  1 1 1 1
   mpermV      [-] :  1 1 1 1
   direct field    :  Comp. in wavenumber domain
   frequency  [Hz] :  1
   Hankel          :  DLF (Fast Hankel Transform)
     > Filter      :  Key 201 (2009)
     > DLF type    :  Standard
   Loop over       :  None (all vectorized)
   Source(s)       :  1 dipole(s)
     > x       [m] :  0
     > y       [m] :  0
     > z       [m] :  0
     > azimuth [°] :  0
     > dip     [°] :  90
   Receiver(s)     :  10 dipole(s)
     > x       [m] :  500 - 5000 : 10  [min-max; #]
                   :  500 1000 1500 2000 2500 3000 3500 4000 4500 5000
     > y       [m] :  0 - 0 : 10  [min-max; #]
                   :  0 0 0 0 0 0 0 0 0 0
     > z       [m] :  200
     > azimuth [°] :  0
     > dip     [°] :  90
   Required ab's   :  33
~
:: empymod END; runtime = 0:00:00.005025 :: 1 kernel call(s)

>>> print(EMfield)
[ -3.05449848e-10 -2.00374185e-11j  -7.12528991e-11 -5.37083268e-12j
  -2.52076501e-11 -1.62732412e-12j  -1.18412295e-11 -8.99570998e-14j
  -6.44054097e-12 +5.61150066e-13j  -3.77109625e-12 +7.89022722e-13j
  -2.28484774e-12 +8.08897623e-13j  -1.40021365e-12 +7.32151174e-13j
  -8.55487532e-13 +6.18402706e-13j  -5.15642408e-13 +4.99091919e-13j]

empymod.model.analytical(src, rec, res, freqtime, solution='fs', signal=None, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, **kwargs)[source]¶

Return analytical full- or half-space solution.

Calculate the electromagnetic frequency- or time-domain field due to infinitesimal small electric or magnetic dipole source(s), measured by infinitesimal small electric or magnetic dipole receiver(s); sources and receivers are directed along the principal directions x, y, or z, and all sources are at the same depth, as well as all receivers are at the same depth.

In the case of a halfspace the air-interface is located at z = 0 m.

You can call the functions empymod.kernel.fullspace() and empymod.kernel.halfspace() in empymod.kernel directly. This interface is just to provide a consistent interface with the same input parameters as for instance for dipole().

This function yields the same result if solution=’fs’ as dipole(), if the model is a fullspace.

Included are:

Full fullspace solution (solution=’fs’) for ee-, me-, em-, mm-fields, only frequency domain, [HuTS15].
Diffusive fullspace solution (solution=’dfs’) for ee-fields, [SlHM10].
Diffusive halfspace solution (solution=’dhs’) for ee-fields, [SlHM10].
Diffusive direct- and reflected field and airwave (solution=’dsplit’) for ee-fields, [SlHM10].
Diffusive direct- and reflected field and airwave (solution=’dtetm’) for ee-fields, split into TE and TM mode [SlHM10].

Parameters:

src, rec : list of floats or arrays

Source and receiver coordinates [x, y, z] (m):

The x- and y-coordinates can be arrays, z is a single value.
The x- and y-coordinates must have the same dimension.

res : float

Horizontal resistivity rho_h (Ohm.m).

Alternatively, res can be a dictionary. See the main manual of empymod too see how to exploit this hook to re-calculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the Cole-Cole model for IP.

freqtime : array_like

Frequencies f (Hz) if signal==None, else times t (s); (f, t > 0).

solution : str, optional

Defines which solution is returned:

‘fs’ : Full fullspace solution (ee-, me-, em-, mm-fields); f-domain.
‘dfs’ : Diffusive fullspace solution (ee-fields only).
‘dhs’ : Diffusive halfspace solution (ee-fields only).
‘dsplit’ : Diffusive direct- and reflected field and airwave (ee-fields only).
‘dtetm’ : as dsplit, but direct fielt TE, TM; reflected field TE, TM, and airwave (ee-fields only).

signal : {None, 0, 1, -1}, optional

Source signal, default is None:

None: Frequency-domain response
-1 : Switch-off time-domain response
0 : Impulse time-domain response
+1 : Switch-on time-domain response

ab : int, optional

Source-receiver configuration, defaults to 11.

		electric source			magnetic source
		x	y	z	x	y	z
electric receiver	x	11	12	13	14	15	16
	y	21	22	23	24	25	26
	z	31	32	33	34	35	36
magnetic receiver	x	41	42	43	44	45	46
	y	51	52	53	54	55	56
	z	61	62	63	64	65	66

aniso : float, optional

Anisotropy lambda = sqrt(rho_v/rho_h) (-); defaults to one.

epermH, epermV : float, optional

Relative horizontal/vertical electric permittivity epsilon_h/epsilon_v (-). If epermH is provided but not epermV, isotropic behaviour is assumed. Default is one; ignored for the diffusive solution.

mpermH, mpermV : float, optional

Relative horizontal/vertical magnetic permeability mu_h/mu_v (-); #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed. Default is one; ignored for the diffusive solution.

verb : {0, 1, 2, 3, 4}, optional

Level of verbosity, default is 2:

0: Print nothing.
1: Print warnings.
2: Print additional runtime
3: Print additional start/stop, condensed parameter information.
4: Print additional full parameter information

Returns:

EM : EMArray, (nfreqtime, nrec, nsrc)

Frequency- or time-domain EM field (depending on signal):

If rec is electric, returns E [V/m].
If rec is magnetic, returns H [A/m].

EMArray is a subclassed ndarray with .pha and .amp attributes (only relevant for frequency-domain data).

The shape of EM is (nfreqtime, nrec, nsrc). However, single dimensions are removed.

If solution=’dsplit’, three ndarrays are returned: direct, reflect, air.

If solution=’dtetm’, five ndarrays are returned: direct_TE, direct_TM, reflect_TE, reflect_TM, air.

Examples

>>> import empymod
>>> import numpy as np
>>> src = [0, 0, 0]
>>> rec = [np.arange(1, 11)*500, np.zeros(10), 200]
>>> res = 50
>>> EMfield = empymod.analytical(src, rec, res, freqtime=1, verb=0)
>>> print(EMfield)
[  4.03091405e-08 -9.69163818e-10j   6.97630362e-09 -4.88342150e-10j
   2.15205979e-09 -2.97489809e-10j   8.90394459e-10 -1.99313433e-10j
   4.32915802e-10 -1.40741644e-10j   2.31674165e-10 -1.02579391e-10j
   1.31469130e-10 -7.62770461e-11j   7.72342470e-11 -5.74534125e-11j
   4.61480481e-11 -4.36275540e-11j   2.76174038e-11 -3.32860932e-11j]

empymod.model.gpr(src, rec, depth, res, freqtime, cf, gain=None, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, **kwargs)[source]¶

Return Ground-Penetrating Radar signal.

THIS FUNCTION IS EXPERIMENTAL, USE WITH CAUTION.

It is rather an example how you can calculate GPR responses; however, DO NOT RELY ON IT! It works only well with QUAD or QWE (quad, qwe) for the Hankel transform, and with FFT (fft) for the Fourier transform.

It calls internally dipole() for the frequency-domain calculation. It subsequently convolves the response with a Ricker wavelet with central frequency cf. If signal!=None, it carries out the Fourier transform and applies a gain to the response.

Parameters:

src, rec, freqtime : survey parameters: See docstring of dipole() for a description.
depth, res, aniso, epermH, epermV, mpermH, mpermV : model parameters: See docstring of dipole() for a description.
cf : float: Centre frequency of GPR-signal, in Hz. Sensible values are between 10 MHz and 3000 MHz.
gain : float: Power of gain function. If None, no gain is applied. Only used if signal!=None.
ht, htarg, ft, ftarg, xdirect, loop : settings, optinal: See docstring of bipole() for a description.

Returns:

EM : ndarray: GPR response

empymod.model.dipole_k(src, rec, depth, res, freq, wavenumber, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, **kwargs)[source]¶

Return electromagnetic wavenumber-domain field.

Calculate the electromagnetic wavenumber-domain field due to infinitesimal small electric or magnetic dipole source(s), measured by infinitesimal small electric or magnetic dipole receiver(s); sources and receivers are directed along the principal directions x, y, or z, and all sources are at the same depth, as well as all receivers are at the same depth.

Parameters:

src, rec : list of floats or arrays

Source and receiver coordinates [x, y, z] (m):

The x- and y-coordinates can be arrays, z is a single value.
The x- and y-coordinates must have the same dimension.
The x- and y-coordinates only matter for the angle-dependent factor.

Sources or receivers placed on a layer interface are considered in the upper layer.

depth : list

Absolute layer interfaces z (m); #depth = #res - 1 (excluding +/- infinity).

res : array_like

Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.

freq : array_like

Frequencies f (Hz), used to calculate etaH/V and zetaH/V.

wavenumber : array

Wavenumbers lambda (1/m)

ab : int, optional

Source-receiver configuration, defaults to 11.

		electric source			magnetic source
		x	y	z	x	y	z
electric receiver	x	11	12	13	14	15	16
	y	21	22	23	24	25	26
	z	31	32	33	34	35	36
magnetic receiver	x	41	42	43	44	45	46
	y	51	52	53	54	55	56
	z	61	62	63	64	65	66

aniso : array_like, optional

Anisotropies lambda = sqrt(rho_v/rho_h) (-); #aniso = #res. Defaults to ones.

epermH, epermV : array_like, optional

Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (-); #epermH = #epermV = #res. If epermH is provided but not epermV, isotropic behaviour is assumed. Default is ones.

mpermH, mpermV : array_like, optional

Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (-); #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed. Default is ones.

verb : {0, 1, 2, 3, 4}, optional

Level of verbosity, default is 2:

0: Print nothing.
1: Print warnings.
2: Print additional runtime and kernel calls
3: Print additional start/stop, condensed parameter information.
4: Print additional full parameter information

Returns:

PJ0, PJ1 : array

Wavenumber-domain EM responses:

PJ0: Wavenumber-domain solution for the kernel with a Bessel function of the first kind of order zero.
PJ1: Wavenumber-domain solution for the kernel with a Bessel function of the first kind of order one.

See also

dipole(): EM fields due to infinitesimal small EM dipoles.
bipole(): EM fields due to arbitrary rotated, finite length EM dipoles.
loop(): EM fields due to a magnetic source loop.

Examples

>>> import empymod
>>> import numpy as np
>>> src = [0, 0, 100]
>>> rec = [5000, 0, 200]
>>> depth = [0, 300, 1000, 1050]
>>> res = [1e20, .3, 1, 50, 1]
>>> freq = 1
>>> wavenr = np.logspace(-3.7, -3.6, 10)
>>> PJ0, PJ1 = empymod.dipole_k(src, rec, depth, res, freq, wavenr, verb=0)
>>> print(PJ0)
[ -1.02638329e-08 +4.91531529e-09j  -1.05289724e-08 +5.04222413e-09j
  -1.08009148e-08 +5.17238608e-09j  -1.10798310e-08 +5.30588284e-09j
  -1.13658957e-08 +5.44279805e-09j  -1.16592877e-08 +5.58321732e-09j
  -1.19601897e-08 +5.72722830e-09j  -1.22687889e-08 +5.87492067e-09j
  -1.25852765e-08 +6.02638626e-09j  -1.29098481e-08 +6.18171904e-09j]
>>> print(PJ1)
[  1.79483705e-10 -6.59235332e-10j   1.88672497e-10 -6.93749344e-10j
   1.98325814e-10 -7.30068377e-10j   2.08466693e-10 -7.68286748e-10j
   2.19119282e-10 -8.08503709e-10j   2.30308887e-10 -8.50823701e-10j
   2.42062030e-10 -8.95356636e-10j   2.54406501e-10 -9.42218177e-10j
   2.67371420e-10 -9.91530051e-10j   2.80987292e-10 -1.04342036e-09j]

empymod.model.fem(ab, off, angle, zsrc, zrec, lsrc, lrec, depth, freq, etaH, etaV, zetaH, zetaV, xdirect, isfullspace, ht, htarg, msrc, mrec, loop_freq, loop_off, conv=True)[source]¶

Return electromagnetic frequency-domain response.

This function is called from one of the modelling routines empymod.model. Consult those for more details regarding the input and output parameters.

This function can be used directly if you are sure the provided input is in the correct format. This is useful for inversion routines and similar, as it can speed-up the calculation by omitting input-checks.

empymod.model.tem(fEM, off, freq, time, signal, ft, ftarg, conv=True)[source]¶

Return time-domain response of the frequency-domain response fEM.

This function is called from one of the modelling routines empymod.model. Consult those for more details regarding the input and output parameters.

This function can be used directly if you are sure the provided input is in the correct format. This is useful for inversion routines and similar, as it can speed-up the calculation by omitting input-checks.

Main modelling routines¶

empymod.model – Model EM-responses¶

`empymod.model` – Model EM-responses¶