Code¶
model
– Model EMresponses¶
EMmodelling routines. The implemented routines might not be the fastest solution to your specific problem. Use these routines as template to create your own, problemspecific modelling routine!
 Principal routines:
bipole
dipole
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.
Further routines are:
analytical
: Calculate analytical fullspace and halfspace solutions.dipole_k
: Calculate the electromagnetic wavenumberdomain solution.gpr
: Calculate the GroundPenetrating Radar (GPR) response.
The dipole_k
routine can be used if you are interested in the
wavenumberdomain result, without Hankel nor Fourier transform. It calls
straight the kernel
. The gpr
routine convolves the frequencydomain
result with a wavelet, and applies a gain to the timedomain result. This
function is still experimental.
 The modelling routines make use of the following two core routines:
fem
: Calculate wavenumberdomain electromagnetic field and carry out the Hankel transform to the frequency domain.
tem
: Carry out the Fourier transform to time domain afterfem
.

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, xdirect=False, ht='fht', htarg=None, ft='sin', ftarg=None, opt=None, loop=None, verb=2)[source]¶ Return the electromagnetic field due to an electromagnetic source.
Calculate the electromagnetic frequency or timedomain field due to arbitrary finite electric or magnetic bipole sources, measured by arbitrary finite electric or magnetic bipole receivers. By default, the electromagnetic response is normalized to 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 variable z (dipole) or z0 and z1 (bipole) must either be single values or having the same dimension as the other coordinates.
 The variables azimuth and dip must be single values. If they have different angles, you have to use the bipolemethod (with srcpts/recpts = 1, so it is calculated as dipoles).
Angles (coordinate system is lefthanded, positive z down (EastNorthDepth):
 azimuth (°): horizontal deviation from xaxis, anticlockwise.
 dip (°): vertical deviation from xyplane downwards.
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 recalculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the ColeCole 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: Frequencydomain response
 1 : Switchoff timedomain response
 0 : Impulse timedomain response
 +1 : Switchon timedomain 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. Default is ones.
 mpermH, mpermV : array_like, optional
Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (); #mpermH = #mpermV = #res. Default is ones.
 msrc, mrec : boolean, 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.
 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.
 ht : {‘fht’, ‘qwe’, ‘quad’}, optional
Flag to choose either the Digital Linear Filter method (FHT, Fast Hankel Transform), the QuadratureWithExtrapolation (QWE), or a simple Quadrature (QUAD) for the Hankel transform. Defaults to ‘fht’.
 htarg : dict or list, optional
 Depends on the value for
ht
: If
ht
= ‘fht’: [fhtfilt, pts_per_dec]: fhtfilt: string of filter name in
empymod.filters
or  the filter method itself.
(default:
empymod.filters.key_201_2009()
)
 fhtfilt: string of filter name in
 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, atol, nquad, maxint, pts_per_dec, diff_quad, a, b, limit]:
 rtol: relative tolerance (default: 1e12)
 atol: absolute tolerance (default: 1e30)
 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 opt=’spline’) (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
If
ht
= ‘quad’: [atol, rtol, limit, lmin, lmax, pts_per_dec]: rtol: relative tolerance (default: 1e12)
 atol: absolute tolerance (default: 1e20)
 limit: An upper bound on the number of subintervals used in the adaptive algorithm (default: 500)
 lmin: Minimum wavenumber (default 1e6)
 lmax: Maximum wavenumber (default 0.1)
 pts_per_dec: points per decade (default: 40)
The values can be provided as dict with the keywords, or as list. However, if provided as list, you have to follow the order given above. A few examples, assuming
ht
=qwe
: Only changing rtol:
 {‘rtol’: 1e4} or [1e4] or 1e4
 Changing rtol and nquad:
 {‘rtol’: 1e4, ‘nquad’: 101} or [1e4, ‘’, 101]
 Only changing diff_quad:
 {‘diffquad’: 10} or [‘’, ‘’, ‘’, ‘’, ‘’, 10]
 Depends on the value for
 ft : {‘sin’, ‘cos’, ‘qwe’, ‘fftlog’, ‘fft’}, optional
Only used if
signal
!= None. Flag to choose either the Digital Linear Filter method (Sine or CosineFilter), the QuadratureWithExtrapolation (QWE), the FFTLog, or the FFT for the Fourier transform. Defaults to ‘sin’. ftarg : dict or list, optional
 Only used if
signal
!=None. Depends on the value forft
: If
ft
= ‘sin’ or ‘cos’: [fftfilt, pts_per_dec]: fftfilt: string of filter name in
empymod.filters
or  the filter method itself.
(Default:
empymod.filters.key_201_CosSin_2012()
)
 fftfilt: string of filter name in
 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, atol, nquad, maxint, pts_per_dec]: rtol: relative tolerance (default: 1e8)
 atol: absolute tolerance (default: 1e20)
 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, add_dec, q]: 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, nfreq, ntot]: dfreq: Linear stepsize 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.
The values can be provided as dict with the keywords, or as list. However, if provided as list, you have to follow the order given above. See
htarg
for a few examples. Only used if
 opt : {None, ‘parallel’}, optional
 Optimization flag. Defaults to None:
 None: Normal case, no parallelization nor interpolation is used.
 If ‘parallel’, the package
numexpr
is used to evaluate the most expensive statements. Always check if it actually improves performance for a specific problem. It can speed up the calculation for big arrays, but will most likely be slower for small arrays. It will use all available cores for these specific statements, which all containGamma
in one way or another, which has dimensions (#frequencies, #offsets, #layers, #lambdas), therefore can grow pretty big. The modulenumexpr
uses by default all available cores up to a maximum of 8. You can change this behaviour to your desired number of threadsnthreads
withnumexpr.set_num_threads(nthreads)
.  The value ‘spline’ is deprecated and will be removed. See
htarg
instead for the interpolated versions.
The option ‘parallel’ only affects speed and memory usage, whereas ‘spline’ also affects precision! Please read the note in the README documentation for more information.
 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 ifopt = 'spline'
. 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! 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:  EM : ndarray, (nfreq, nrec, nsrc)
 Frequency or timedomain EM field (depending on
signal
):  If rec is electric, returns E [V/m].
 If rec is magnetic, returns H [A/m].
However, source and receiver are normalised (unless strength != 0). So for instance in the electric case the source strength is 1 A and its length is 1 m. So the electric field could also be written as [V/(A.m2)].
In the magnetic case the source strength is given by \(i\omega\mu_0 A I^e\), where A is the loop area (m2), and \(I^e\) the electric source strength. For the normalized magnetic source \(A=1m^2\) and \(I^e=1 Ampere\). A magnetic source is therefore frequency dependent.
The shape of EM is (nfreq, nrec, nsrc). However, single dimensions are removed.
 Frequency or timedomain EM field (depending on
Examples
>>> import numpy as np >>> from empymod import bipole >>> # xdirected bipole source: x0, x1, y0, y1, z0, z1 >>> src = [50, 50, 0, 0, 100, 100] >>> # xdirected dipole sourcearray: 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 = True (default)] >>> EMfield = bipole(src, rec, depth, res, freq, verb=4) :: empymod START :: ~ 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 frequency [Hz] : 1 Hankel : DLF (Fast Hankel Transform) > Filter : Key 201 (2009) > DLF type : Standard Kernel Opt. : None Loop over : None (all vectorized) Source(s) : 1 bipole(s) > intpts : 1 (as dipole) > length [m] : 100 > 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 [minmax; #] : 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 > y [m] : 0  0 : 10 [minmax; #] : 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.68809346e10 3.08303130e10j 8.77189179e12 3.76920235e11j 3.46654704e12 4.87133683e12j 3.60159726e13 1.12434417e12j 1.87807271e13 6.21669759e13j 1.97200208e13 4.38210489e13j 1.44134842e13 3.17505260e13j 9.92770406e14 2.33950871e13j 6.75287598e14 1.74922886e13j 4.62724887e14 1.32266600e13j]

empymod.model.
dipole
(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, xdirect=False, ht='fht', htarg=None, ft='sin', ftarg=None, opt=None, loop=None, verb=2)[source]¶ Return the electromagnetic field due to a dipole source.
Calculate the electromagnetic frequency or timedomain 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 bybipole
(all there is to do is translateab
intomsrc
,mrec
,azimuth
’s anddip
’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 (m): [x, y, z]. The x and ycoordinates can be arrays, z is a single value. The x and ycoordinates 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 recalculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the ColeCole 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: Frequencydomain response
 1 : Switchoff timedomain response
 0 : Impulse timedomain response
 +1 : Switchon timedomain response
 ab : int, optional
Sourcereceiver 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. Default is ones.
 mpermH, mpermV : array_like, optional
Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (); #mpermH = #mpermV = #res. Default is ones.
 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.
 ht : {‘fht’, ‘qwe’, ‘quad’}, optional
Flag to choose either the Digital Linear Filter method (FHT, Fast Hankel Transform), the QuadratureWithExtrapolation (QWE), or a simple Quadrature (QUAD) for the Hankel transform. Defaults to ‘fht’.
 htarg : dict or list, optional
 Depends on the value for
ht
: If
ht
= ‘fht’: [fhtfilt, pts_per_dec]: fhtfilt: string of filter name in
empymod.filters
or  the filter method itself.
(default:
empymod.filters.key_201_2009()
)
 fhtfilt: string of filter name in
 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, atol, nquad, maxint, pts_per_dec, diff_quad, a, b, limit]:
 rtol: relative tolerance (default: 1e12)
 atol: absolute tolerance (default: 1e30)
 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 opt=’spline’) (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
If
ht
= ‘quad’: [atol, rtol, limit, lmin, lmax, pts_per_dec]: rtol: relative tolerance (default: 1e12)
 atol: absolute tolerance (default: 1e20)
 limit: An upper bound on the number of subintervals used in the adaptive algorithm (default: 500)
 lmin: Minimum wavenumber (default 1e6)
 lmax: Maximum wavenumber (default 0.1)
 pts_per_dec: points per decade (default: 40)
The values can be provided as dict with the keywords, or as list. However, if provided as list, you have to follow the order given above. A few examples, assuming
ht
=qwe
: Only changing rtol:
 {‘rtol’: 1e4} or [1e4] or 1e4
 Changing rtol and nquad:
 {‘rtol’: 1e4, ‘nquad’: 101} or [1e4, ‘’, 101]
 Only changing diff_quad:
 {‘diffquad’: 10} or [‘’, ‘’, ‘’, ‘’, ‘’, 10]
 Depends on the value for
 ft : {‘sin’, ‘cos’, ‘qwe’, ‘fftlog’, ‘fft’}, optional
Only used if
signal
!= None. Flag to choose either the Digital Linear Filter method (Sine or CosineFilter), the QuadratureWithExtrapolation (QWE), the FFTLog, or the FFT for the Fourier transform. Defaults to ‘sin’. ftarg : dict or list, optional
 Only used if
signal
!=None. Depends on the value forft
: If
ft
= ‘sin’ or ‘cos’: [fftfilt, pts_per_dec]: fftfilt: string of filter name in
empymod.filters
or  the filter method itself.
(Default:
empymod.filters.key_201_CosSin_2012()
)
 fftfilt: string of filter name in
 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, atol, nquad, maxint, pts_per_dec]: rtol: relative tolerance (default: 1e8)
 atol: absolute tolerance (default: 1e20)
 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, add_dec, q]: 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, nfreq, ntot]: dfreq: Linear stepsize 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.
The values can be provided as dict with the keywords, or as list. However, if provided as list, you have to follow the order given above. See
htarg
for a few examples. Only used if
 opt : {None, ‘parallel’}, optional
 Optimization flag. Defaults to None:
 None: Normal case, no parallelization nor interpolation is used.
 If ‘parallel’, the package
numexpr
is used to evaluate the most expensive statements. Always check if it actually improves performance for a specific problem. It can speed up the calculation for big arrays, but will most likely be slower for small arrays. It will use all available cores for these specific statements, which all containGamma
in one way or another, which has dimensions (#frequencies, #offsets, #layers, #lambdas), therefore can grow pretty big. The modulenumexpr
uses by default all available cores up to a maximum of 8. You can change this behaviour to your desired number of threadsnthreads
withnumexpr.set_num_threads(nthreads)
.  The value ‘spline’ is deprecated and will be removed. See
htarg
instead for the interpolated versions.
The option ‘parallel’ only affects speed and memory usage, whereas ‘spline’ also affects precision! Please read the note in the README documentation for more information.
 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 ifopt = 'spline'
. 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! 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:  EM : ndarray, (nfreq, nrec, nsrc)
 Frequency or timedomain EM field (depending on
signal
):  If rec is electric, returns E [V/m].
 If rec is magnetic, returns H [A/m].
However, source and receiver are normalised. So for instance in the electric case the source strength is 1 A and its length is 1 m. So the electric field could also be written as [V/(A.m2)].
The shape of EM is (nfreq, nrec, nsrc). However, single dimensions are removed.
 Frequency or timedomain EM field (depending on
See also
Examples
>>> import numpy as np >>> from empymod import dipole >>> 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 = dipole(src, rec, depth, res, freqtime=1, verb=0) >>> print(EMfield) [ 1.68809346e10 3.08303130e10j 8.77189179e12 3.76920235e11j 3.46654704e12 4.87133683e12j 3.60159726e13 1.12434417e12j 1.87807271e13 6.21669759e13j 1.97200208e13 4.38210489e13j 1.44134842e13 3.17505260e13j 9.92770406e14 2.33950871e13j 6.75287598e14 1.74922886e13j 4.62724887e14 1.32266600e13j]

empymod.model.
analytical
(src, rec, res, freqtime, solution='fs', signal=None, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, verb=2)[source]¶ Return the analytical full or halfspace solution.
Calculate the electromagnetic frequency or timedomain 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 airinterface is located at z = 0 m.
You can call the functions
fullspace
andhalfspace
inkernel.py
directly. This interface is just to provide a consistent interface with the same input parameters as for instance fordipole
.This function yields the same result if
solution='fs'
asdipole
, if the model is a fullspace. Included are:
 Full fullspace solution (
solution='fs'
) for ee, me, em, mmfields, only frequency domain, [HuTS15].  Diffusive fullspace solution (
solution='dfs'
) for eefields, [SlHM10].  Diffusive halfspace solution (
solution='dhs'
) for eefields, [SlHM10].  Diffusive direct and reflected field and airwave
(
solution='dsplit'
) for eefields, [SlHM10].  Diffusive direct and reflected field and airwave
(
solution='dtetm'
) for eefields, split into TE and TM mode [SlHM10].
 Full fullspace solution (
Parameters:  src, rec : list of floats or arrays
Source and receiver coordinates (m): [x, y, z]. The x and ycoordinates can be arrays, z is a single value. The x and ycoordinates 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 recalculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the ColeCole 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, mmfields); fdomain.
 ‘dfs’ : Diffusive fullspace solution (eefields only).
 ‘dhs’ : Diffusive halfspace solution (eefields only).
 ‘dsplit’ : Diffusive direct and reflected field and airwave
 (eefields only).
 ‘dtetm’ : as dsplit, but direct fielt TE, TM; reflected field TE, TM,
 and airwave (eefields only).
 signal : {None, 0, 1, 1}, optional
 Source signal, default is None:
 None: Frequencydomain response
 1 : Switchoff timedomain response
 0 : Impulse timedomain response
 +1 : Switchon timedomain response
 ab : int, optional
Sourcereceiver 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 (); default is one. Ignored for the diffusive solution.
 mpermH, mpermV : float, optional
Relative horizontal/vertical magnetic permeability mu_h/mu_v (); 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 : ndarray, (nfreq, nrec, nsrc)
 Frequency or timedomain EM field (depending on
signal
):  If rec is electric, returns E [V/m].
 If rec is magnetic, returns H [A/m].
However, source and receiver are normalised. So for instance in the electric case the source strength is 1 A and its length is 1 m. So the electric field could also be written as [V/(A.m2)].
The shape of EM is (nfreq, 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. Frequency or timedomain EM field (depending on
Examples
>>> import numpy as np >>> from empymod import analytical >>> src = [0, 0, 0] >>> rec = [np.arange(1, 11)*500, np.zeros(10), 200] >>> res = 50 >>> EMfield = analytical(src, rec, res, freqtime=1, verb=0) >>> print(EMfield) [ 4.03091405e08 9.69163818e10j 6.97630362e09 4.88342150e10j 2.15205979e09 2.97489809e10j 8.90394459e10 1.99313433e10j 4.32915802e10 1.40741644e10j 2.31674165e10 1.02579391e10j 1.31469130e10 7.62770461e11j 7.72342470e11 5.74534125e11j 4.61480481e11 4.36275540e11j 2.76174038e11 3.32860932e11j]

empymod.model.
gpr
(src, rec, depth, res, freqtime, cf, gain=None, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, xdirect=False, ht='quad', htarg=None, ft='fft', ftarg=None, opt=None, loop=None, verb=2)[source]¶ Return the GroundPenetrating 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 frequencydomain calculation. It subsequently convolves the response with a Ricker wavelet with central frequencycf
. If signal!=None, it carries out the Fourier transform and applies a gain to the response.For input parameters see the function
dipole
, except for:Parameters:  cf : float
Centre frequency of GPRsignal, 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.
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, verb=2)[source]¶ Return the electromagnetic wavenumberdomain field.
Calculate the electromagnetic wavenumberdomain 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 (m): [x, y, z]. The x and ycoordinates can be arrays, z is a single value. The x and ycoordinates must have the same dimension. The x and ycoordinates only matter for the angledependent 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
Sourcereceiver 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. Default is ones.
 mpermH, mpermV : array_like, optional
Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (); #mpermH = #mpermV = #res. 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
 Wavenumberdomain EM responses:
 PJ0: Wavenumberdomain solution for the kernel with a Bessel function of the first kind of order zero.
 PJ1: Wavenumberdomain solution for the kernel with a Bessel function of the first kind of order one.
See also
Examples
>>> import numpy as np >>> from empymod.model import dipole_k >>> src = [0, 0, 100] >>> rec = [5000, 0, 200] >>> depth = [0, 300, 1000, 1050] >>> res = [1e20, .3, 1, 50, 1] >>> freq = 1 >>> wavenrs = np.logspace(3.7, 3.6, 10) >>> PJ0, PJ1 = dipole_k(src, rec, depth, res, freq, wavenrs, verb=0) >>> print(PJ0) [ 1.02638329e08 +4.91531529e09j 1.05289724e08 +5.04222413e09j 1.08009148e08 +5.17238608e09j 1.10798310e08 +5.30588284e09j 1.13658957e08 +5.44279805e09j 1.16592877e08 +5.58321732e09j 1.19601897e08 +5.72722830e09j 1.22687889e08 +5.87492067e09j 1.25852765e08 +6.02638626e09j 1.29098481e08 +6.18171904e09j] >>> print(PJ1) [ 1.79483705e10 6.59235332e10j 1.88672497e10 6.93749344e10j 1.98325814e10 7.30068377e10j 2.08466693e10 7.68286748e10j 2.19119282e10 8.08503709e10j 2.30308887e10 8.50823701e10j 2.42062030e10 8.95356636e10j 2.54406501e10 9.42218177e10j 2.67371420e10 9.91530051e10j 2.80987292e10 1.04342036e09j]

empymod.model.
fem
(ab, off, angle, zsrc, zrec, lsrc, lrec, depth, freq, etaH, etaV, zetaH, zetaV, xdirect, isfullspace, ht, htarg, use_ne_eval, msrc, mrec, loop_freq, loop_off, conv=True)[source]¶ Return the electromagnetic frequencydomain response.
This function is called from one of the above modelling routines. No inputcheck is carried out here. See the main description of
model
for information regarding input and output parameters.This function can be directly used if you are sure the provided input is in the correct format. This is useful for inversion routines and similar, as it can speedup the calculation by omitting inputchecks.

empymod.model.
tem
(fEM, off, freq, time, signal, ft, ftarg, conv=True)[source]¶ Return the timedomain response of the frequencydomain response fEM.
This function is called from one of the above modelling routines. No inputcheck is carried out here. See the main description of
model
for information regarding input and output parameters.This function can be directly used if you are sure the provided input is in the correct format. This is useful for inversion routines and similar, as it can speedup the calculation by omitting inputchecks.
kernel
– Kernel calculation¶
Kernel of empymod
, calculates the wavenumberdomain electromagnetic
response. Plus analytical full and halfspace solutions.
The functions wavenumber
, angle_factor
, fullspace
, greenfct
,
reflections
, and fields
are based on source files (specified in each
function) from the source code distributed with [HuTS15], which can be found
at software.seg.org/2015/0001. These
functions are (c) 2015 by Hunziker et al. and the Society of Exploration
Geophysicists, http://software.seg.org/disclaimer.txt. Please read the
NOTICEfile in the root directory for more information regarding the involved
licenses.

empymod.kernel.
wavenumber
(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd, ab, xdirect, msrc, mrec, use_ne_eval)[source]¶ Calculate wavenumber domain solution.
Return the wavenumber domain solutions
PJ0
,PJ1
, andPJ0b
, which have to be transformed with a Hankel transform to the frequency domain.PJ0
/PJ0b
andPJ1
have to be transformed with Bessel functions of order 0 (\(J_0\)) and 1 (\(J_1\)), respectively.This function corresponds loosely to equations 105–107, 111–116, 119–121, and 123–128 in [HuTS15], and equally loosely to the file
kxwmod.c
.[HuTS15] uses Bessel functions of orders 0, 1, and 2 (\(J_0, J_1, J_2\)). The implementations of the Fast Hankel Transform and the QuadraturewithExtrapolation in
transform
are setup with Bessel functions of order 0 and 1 only. This is achieved by applying the recurrence formula\[J_2(kr) = \frac{2}{kr} J_1(kr)  J_0(kr) \ .\]Note
PJ0
andPJ0b
could theoretically be added here into one, and then be transformed in one go. However,PJ0b
has to be multiplied byfactAng
later. This has to be done after the Hankel transform for methods which make use of spline interpolation, in order to work for offsets that are not in line with each other.This function is called from one of the Hankel functions in
transform
. Consult the modelling routines inmodel
for a description of the input and output parameters.If you are solely interested in the wavenumberdomain solution you can call this function directly. However, you have to make sure all input arguments are correct, as no checks are carried out here.

empymod.kernel.
angle_factor
(angle, ab, msrc, mrec)[source]¶ Return the angledependent factor.
The whole calculation in the wavenumber domain is only a function of the distance between the source and the receiver, it is independent of the angel. The angledependency is this factor, which can be applied to the corresponding parts in the wavenumber or in the frequency domain.
The
angle_factor
corresponds to the sine and cosinefunctions in Eqs 105107, 111116, 119121, 123128.This function is called from one of the Hankel functions in
transform
. Consult the modelling routines inmodel
for a description of the input and output parameters.

empymod.kernel.
fullspace
(off, angle, zsrc, zrec, etaH, etaV, zetaH, zetaV, ab, msrc, mrec)[source]¶ Analytical fullspace solutions in the frequency domain.
\[\hat{G}^{ee}_{\alpha\beta}, \hat{G}^{ee}_{3\alpha}, \hat{G}^{ee}_{33}, \hat{G}^{em}_{\alpha\beta}, \hat{G}^{em}_{\alpha 3}\]This function corresponds to equations 45–50 in [HuTS15], and loosely to the corresponding files
Gin11.F90
,Gin12.F90
,Gin13.F90
,Gin22.F90
,Gin23.F90
,Gin31.F90
,Gin32.F90
,Gin33.F90
,Gin41.F90
,Gin42.F90
,Gin43.F90
,Gin51.F90
,Gin52.F90
,Gin53.F90
,Gin61.F90
, andGin62.F90
.This function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.

empymod.kernel.
greenfct
(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd, ab, xdirect, msrc, mrec, use_ne_eval)[source]¶ Calculate Green’s function for TM and TE.
\[\tilde{g}^{tm}_{hh}, \tilde{g}^{tm}_{hz}, \tilde{g}^{tm}_{zh}, \tilde{g}^{tm}_{zz}, \tilde{g}^{te}_{hh}, \tilde{g}^{te}_{zz}\]This function corresponds to equations 108–110, 117/118, 122; 89–94, A18–A23, B13–B15; 97–102 A26–A31, and B16–B18 in [HuTS15], and loosely to the corresponding files
Gamma.F90
,Wprop.F90
,Ptotalx.F90
,Ptotalxm.F90
,Ptotaly.F90
,Ptotalym.F90
,Ptotalz.F90
, andPtotalzm.F90
.The Green’s functions are multiplied according to Eqs 105107, 111116, 119121, 123128; with the factors inside the integrals.
This function is called from the function
kernel.wavenumber
.

empymod.kernel.
reflections
(depth, e_zH, Gam, lrec, lsrc, use_ne_eval)[source]¶ Calculate Rp, Rm.
\[R^\pm_n, \bar{R}^\pm_n\]This function corresponds to equations 64/65 and A11/A12 in [HuTS15], and loosely to the corresponding files
Rmin.F90
andRplus.F90
.This function is called from the function
kernel.greenfct
.

empymod.kernel.
fields
(depth, Rp, Rm, Gam, lrec, lsrc, zsrc, ab, TM, use_ne_eval)[source]¶ Calculate Pu+, Pu, Pd+, Pd.
\[P^{u\pm}_s, P^{d\pm}_s, \bar{P}^{u\pm}_s, \bar{P}^{d\pm}_s; P^{u\pm}_{s1}, P^{u\pm}_n, \bar{P}^{u\pm}_{s1}, \bar{P}^{u\pm}_n; P^{d\pm}_{s+1}, P^{d\pm}_n, \bar{P}^{d\pm}_{s+1}, \bar{P}^{d\pm}_n\]This function corresponds to equations 81/82, 95/96, 103/104, A8/A9, A24/A25, and A32/A33 in [HuTS15], and loosely to the corresponding files
Pdownmin.F90
,Pdownplus.F90
,Pupmin.F90
, andPdownmin.F90
.This function is called from the function
kernel.greenfct
.

empymod.kernel.
halfspace
(off, angle, zsrc, zrec, etaH, etaV, freqtime, ab, signal, solution='dhs')[source]¶ Return frequency or timespace domain VTI halfspace solution.
Calculates the frequency or timespace domain electromagnetic response for a halfspace below air using the diffusive approximation, as given in [SlHM10], where the electric source is located at [0, 0, zsrc], and the electric receiver at [xco, yco, zrec].
It can also be used to calculate the fullspace solution or the separate fields: direct field, reflected field, and airwave; always using the diffusive approximation. See
solution
parameter.This function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and solution parameters.
transform
– Hankel and Fourier Transforms¶
Methods to carry out the required Hankel transform from wavenumber to frequency domain and Fourier transform from frequency to time domain.
The functions for the QWE and DLF Hankel and Fourier transforms are based on source files (specified in each function) from the source code distributed with [Key12], which can be found at software.seg.org/2012/0003. These functions are (c) 2012 by Kerry Key and the Society of Exploration Geophysicists, http://software.seg.org/disclaimer.txt. Please read the NOTICEfile in the root directory for more information regarding the involved licenses.

empymod.transform.
fht
(zsrc, zrec, lsrc, lrec, off, factAng, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, fhtarg, use_ne_eval, msrc, mrec)[source]¶ Hankel Transform using the Digital Linear Filter method.
The Digital Linear Filter method was introduced to geophysics by [Ghos70], and made popular and widespread by [Ande75], [Ande79], [Ande82]. The DLF is sometimes referred to as the Fast Hankel Transform FHT, from which this routine has its name.
This implementation of the DLF follows [Key12], equation 6. Without going into the mathematical details (which can be found in any of the above papers) and following [Key12], the DLF method rewrites the Hankel transform of the form
\[F(r) = \int^\infty_0 f(\lambda)J_v(\lambda r)\ \mathrm{d}\lambda\]as
\[F(r) = \sum^n_{i=1} f(b_i/r)h_i/r \ ,\]where \(h\) is the digital filter.The Filter abscissae b is given by
\[b_i = \lambda_ir = e^{ai}, \qquad i = l, l+1, \cdots, l \ ,\]with \(l=(n1)/2\), and \(a\) is the spacing coefficient.
This function is loosely based on
get_CSEM1D_FD_FHT.m
from the source code distributed with [Key12].The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.Returns:  fEM : array
Returns frequencydomain EM response.
 kcount : int
Kernel count. For DLF, this is 1.
 conv : bool
Only relevant for QWE/QUAD.

empymod.transform.
hqwe
(zsrc, zrec, lsrc, lrec, off, factAng, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, qweargs, use_ne_eval, msrc, mrec)[source]¶ Hankel Transform using QuadratureWithExtrapolation.
QuadratureWithExtrapolation was introduced to geophysics by [Key12]. It is one of many socalled ISE methods to solve Hankel Transforms, where ISE stands for Integration, Summation, and Extrapolation.
Following [Key12], but without going into the mathematical details here, the QWE method rewrites the Hankel transform of the form
\[F(r) = \int^\infty_0 f(\lambda)J_v(\lambda r)\ \mathrm{d}\lambda\]as a quadrature sum which form is similar to the DLF (equation 15),
\[F_i \approx \sum^m_{j=1} f(x_j/r)w_j g(x_j) = \sum^m_{j=1} f(x_j/r)\hat{g}(x_j) \ ,\]but with various bells and whistles applied (using the socalled Shanks transformation in the form of a routine called \(\epsilon\)algorithm ([Shan55], [Wynn56]; implemented with algorithms from [Tref00] and [Weni89]).
This function is based on
get_CSEM1D_FD_QWE.m
,qwe.m
, andgetBesselWeights.m
from the source code distributed with [Key12].In the splineversion,
hqwe
checks how steep the decay of the wavenumberdomain result is, and calls QUAD for the very steep interval, for which QWE is not suited.The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.Returns:  fEM : array
Returns frequencydomain EM response.
 kcount : int
Kernel count.
 conv : bool
If true, QWE/QUAD converged. If not, <htarg> might have to be adjusted.

empymod.transform.
hquad
(zsrc, zrec, lsrc, lrec, off, factAng, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, quadargs, use_ne_eval, msrc, mrec)[source]¶ Hankel Transform using the
QUADPACK
library.This routine uses the
scipy.integrate.quad
module, which in turn makes use of the Fortran libraryQUADPACK
(qagse
).It is massively (orders of magnitudes) slower than either
fht
orhqwe
, and is mainly here for completeness and comparison purposes. It always uses interpolation in the wavenumber domain, hence it generally will not be as precise as the other methods. However, it might work in some areas where the others fail.The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.Returns:  fEM : array
Returns frequencydomain EM response.
 kcount : int
Kernel count. For HQUAD, this is 1.
 conv : bool
If true, QUAD converged. If not, <htarg> might have to be adjusted.

empymod.transform.
ffht
(fEM, time, freq, ftarg)[source]¶ Fourier Transform using the Digital Linear Filter method.
It follows the Filter methodology [Ande75], using Cosine and Sinefilters; see
fht
for more information.The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.This function is based on
get_CSEM1D_TD_FHT.m
from the source code distributed with [Key12].Returns:  tEM : array
Returns timedomain EM response of
fEM
for giventime
. conv : bool
Only relevant for QWE/QUAD.

empymod.transform.
fqwe
(fEM, time, freq, qweargs)[source]¶ Fourier Transform using QuadratureWithExtrapolation.
It follows the QWE methodology [Key12] for the Hankel transform, see
hqwe
for more information.The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.This function is based on
get_CSEM1D_TD_QWE.m
from the source code distributed with [Key12].fqwe
checks how steep the decay of the frequencydomain result is, and calls QUAD for the very steep interval, for which QWE is not suited.Returns:  tEM : array
Returns timedomain EM response of
fEM
for giventime
. conv : bool
If true, QWE/QUAD converged. If not, <ftarg> might have to be adjusted.

empymod.transform.
fftlog
(fEM, time, freq, ftarg)[source]¶ Fourier Transform using FFTLog.
FFTLog is the logarithmic analogue to the Fast Fourier Transform FFT. FFTLog was presented in Appendix B of [Hami00] and published at <http://casa.colorado.edu/~ajsh/FFTLog>.
This function uses a simplified version of
pyfftlog
, which is a pythonversion ofFFTLog
. For more details regardingpyfftlog
see <https://github.com/prisae/pyfftlog>.Not the full flexibility of
FFTLog
is available here: Only the logarithmic FFT (fftl
inFFTLog
), not the Hankel transform (fht
inFFTLog
). Furthermore, the following parameters are fixed:kr
= 1 (initial value)kropt
= 1 (silently adjustskr
)dir
= 1 (forward)
Furthermore,
q
is restricted to 1 <= q <= 1.The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.Returns:  tEM : array
Returns timedomain EM response of
fEM
for giventime
. conv : bool
Only relevant for QWE/QUAD.

empymod.transform.
fft
(fEM, time, freq, ftarg)[source]¶ Fourier Transform using the Fast Fourier Transform.
The function is called from one of the modelling routines in
model
. Consult these modelling routines for a description of the input and output parameters.Returns:  tEM : array
Returns timedomain EM response of
fEM
for giventime
. conv : bool
Only relevant for QWE/QUAD.

empymod.transform.
dlf
(signal, points, out_pts, filt, pts_per_dec, kind=None, factAng=None, ab=None, int_pts=None)[source]¶ Digital Linear Filter method.
This is the kernel of the DLF method, used for the Hankel (
fht
) and the Fourier (ffht
) Transforms. Seefht
for an extensive description.For the Hankel transform, signal contains 3 complex wavenumberdomain signals: (PJ0, PJ1, PJ0b), as returned from kernel.wavenumber. The Hankel DLF has two additional, optional parameters: factAng, as returned from kernel.angle_factor, and ab. The PJ0kernel is the part of the wavenumberdomain calculation which contains a zerothorder Bessel function and does NOT depend on the angle between source and receiver, only on offset. PJ0b and PJ1 are the parts of the wavenumberdomain calculation which contain a zeroth and firstorder Bessel function, respectively, and can depend on the angle between source and receiver. PJ0, PJ1, or PJ0b can also be None, if they are not used.
For the Fourier transform, signal is a complex frequencydomain signal. The Fourier DLF requires one additional parameter, kind, which will be ‘cos’ or ‘sin’.

empymod.transform.
qwe
(rtol, atol, maxint, inp, intervals, lambd=None, off=None, factAng=None)[source]¶ QuadratureWithExtrapolation.
This is the kernel of the QWE method, used for the Hankel (
hqwe
) and the Fourier (fqwe
) Transforms. Seehqwe
for an extensive description.This function is based on
qwe.m
from the source code distributed with [Key12].
filters
– Digital Linear Filters¶
Filters for the Digital Linear Filter (DLF) method for the Hankel [Ghos70]) and the Fourier ([Ande75]) transforms.
To calculate the dlf.factor
I used
np.around(np.average(dlf.base[1:]/dlf.base[:1]), 15)
The filters kong_61_2007
and kong_241_2007
from [Kong07], and
key_101_2009
, key_201_2009
, key_401_2009
, key_81_CosSin_2009
,
key_241_CosSin_2009
, and key_601_CosSin_2009
from [Key09] are taken
from DIPOLE1D, [Key09], which can be downloaded at
http://marineemlab.ucsd.edu/Projects/Occam/1DCSEM (1DCSEM). DIPOLE1D is
distributed under the license GNU GPL version 3 or later. Kerry Key gave his
written permission to redistribute the filters under the Apache License,
Version 2.0 (email from Kerry Key to Dieter Werthmüller, 21 November 2016).
The filters anderson_801_1982
from [Ande82] and key_51_2012
,
key_101_2012
, key_201_2012
, key_101_CosSin_2012
, and
key_201_CosSin_2012
, all from [Key12], are taken from the software
distributed with [Key12] and available at http://software.seg.org/2012/0003
(SEG2012003). These filters are distributed under the SEG license.
The filter wer_201_2018
was designed with the addon fdesign
, see
https://github.com/empymod/articlefdesign.

class
empymod.filters.
DigitalFilter
(name, savename=None)[source]¶ Simple Class for Digital Linear Filters.

fromfile
(self, path='filters')[source]¶ Load filter values from asciifiles.
Load filter base and filter coefficients from ascii files in the directory path; path can be a relative or absolute path.
Examples
>>> import empymod >>> # Create an empty filter; >>> # Name has to be the base of the text files >>> filt = empymod.filters.DigitalFilter('myfilter') >>> # Load the asciifiles >>> filt.fromfile() >>> # This will load the following three files: >>> # ./filters/myfilter_base.txt >>> # ./filters/myfilter_j0.txt >>> # ./filters/myfilter_j1.txt >>> # and store them in filt.base, filt.j0, and filt.j1.

tofile
(self, path='filters')[source]¶ Save filter values to asciifiles.
Store the filter base and the filter coefficients in separate files in the directory path; path can be a relative or absolute path.
Examples
>>> import empymod >>> # Load a filter >>> filt = empymod.filters.wer_201_2018() >>> # Save it to pure asciifiles >>> filt.tofile() >>> # This will save the following three files: >>> # ./filters/wer_201_2018_base.txt >>> # ./filters/wer_201_2018_j0.txt >>> # ./filters/wer_201_2018_j1.txt


empymod.filters.
kong_61_2007
()[source]¶ Kong 61 pt Hankel filter, as published in [Kong07].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
kong_241_2007
()[source]¶ Kong 241 pt Hankel filter, as published in [Kong07].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
key_101_2009
()[source]¶ Key 101 pt Hankel filter, as published in [Key09].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
key_201_2009
()[source]¶ Key 201 pt Hankel filter, as published in [Key09].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
key_401_2009
()[source]¶ Key 401 pt Hankel filter, as published in [Key09].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
anderson_801_1982
()[source]¶ Anderson 801 pt Hankel filter, as published in [Ande82].
Taken from file
wa801Hankel.txt
provided with SEG2012003.License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_51_2012
()[source]¶ Key 51 pt Hankel filter, as published in [Key12].
Taken from file
kk51Hankel.txt
provided with SEG2012003.License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_101_2012
()[source]¶ Key 101 pt Hankel filter, as published in [Key12].
Taken from file
kk101Hankel.txt
provided with SEG2012003.License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_201_2012
()[source]¶ Key 201 pt Hankel filter, as published in [Key12].
Taken from file
kk201Hankel.txt
provided with SEG2012003.License: http://software.seg.org/disclaimer.txt.

empymod.filters.
wer_201_2018
()[source]¶ Werthmüller 201 pt Hankel filter, 2018.
Designed with the empymod addon
fdesign
, see https://github.com/empymod/articlefdesign.License: Apache License, Version 2.0,.

empymod.filters.
key_81_CosSin_2009
()[source]¶ Key 81 pt CosSin filter, as published in [Key09].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
key_241_CosSin_2009
()[source]¶ Key 241 pt CosSin filter, as published in [Key09].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
key_601_CosSin_2009
()[source]¶ Key 601 pt CosSin filter, as published in [Key09].
Taken from file
FilterModules.f90
provided with 1DCSEM.License: Apache License, Version 2.0,.

empymod.filters.
key_101_CosSin_2012
()[source]¶ Key 101 pt CosSin filter, as published in [Key12].
Taken from file
kk101CosSin.txt
provided with SEG2012003.License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_201_CosSin_2012
()[source]¶ Key 201 pt CosSin filter, as published in [Key12].
Taken from file
kk201CosSin.txt
provided with SEG2012003.License: http://software.seg.org/disclaimer.txt.
utils
– Utilites¶
Utilities for model
such as checking input parameters.
 This module consists of four groups of functions:
 General settings
 Class EMArray
 Input parameter checks for modelling
 Internal utilities

class
empymod.utils.
EMArray
[source]¶ Subclassing an ndarray: add amplitude <amp> and phase <pha>.
Parameters:  realpart : array
 Real part of input, if input is real or complex.
 Imaginary part of input, if input is pure imaginary.
 Complex input.
In cases 2 and 3,
imagpart
must be None. imagpart: array, optional
Imaginary part of input. Defaults to None.
Examples
>>> import numpy as np >>> from empymod.utils import EMArray >>> emvalues = EMArray(np.array([1,2,3]), np.array([1, 0, 1])) >>> print('Amplitude : ', emvalues.amp) Amplitude : [ 1.41421356 2. 3.16227766] >>> print('Phase : ', emvalues.pha) Phase : [ 45. 0. 18.43494882]
Attributes:  amp : ndarray
Amplitude of the input data.
 pha : ndarray
Phase of the input data, in degrees, lagdefined (increasing with increasing offset.) To get leaddefined phases, multiply
imagpart
by 1 before passing through this function.

empymod.utils.
check_time_only
(time, signal, verb)[source]¶ Check time and signal parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  time : array_like
Times t (s).
 signal : {None, 0, 1, 1}
 Source signal:
 None: Frequencydomain response
 1 : Switchoff timedomain response
 0 : Impulse timedomain response
 +1 : Switchon timedomain response
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  time : float
Time, checked for size and assured min_time.

empymod.utils.
check_time
(time, signal, ft, ftarg, verb)[source]¶ Check time domain specific input parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  time : array_like
Times t (s).
 signal : {None, 0, 1, 1}
 Source signal:
 None: Frequencydomain response
 1 : Switchoff timedomain response
 0 : Impulse timedomain response
 +1 : Switchon timedomain response
 ft : {‘sin’, ‘cos’, ‘qwe’, ‘fftlog’, ‘fft’}
Flag for Fourier transform.
 ftarg : str or filter from empymod.filters or array_like,
Only used if
signal
!=None. Depends on the value forft
: verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  time : float
Time, checked for size and assured min_time.
 freq : float
Frequencies required for given times and ftsettings.
 ft, ftarg
Checked if valid and set to defaults if not provided, checked with signal.

empymod.utils.
check_model
(depth, res, aniso, epermH, epermV, mpermH, mpermV, xdirect, verb)[source]¶ Check the model: depth and corresponding layer parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  depth : list
Absolute layer interfaces z (m); #depth = #res  1 (excluding +/ infinity).
 res : array_like
Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.
 aniso : array_like
Anisotropies lambda = sqrt(rho_v/rho_h) (); #aniso = #res.
 epermH, epermV : array_like
Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (); #epermH = #epermV = #res.
 mpermH, mpermV : array_like
Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (); #mpermH = #mpermV = #res.
 xdirect : bool, optional
If True and source and receiver are in the same layer, the direct field is calculated analytically in the frequency domain, if False it is calculated in the wavenumber domain.
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  depth : array
Depths of layer interfaces, adds infty at beginning if not present.
 res : array
As input, checked for size.
 aniso : array
As input, checked for size. If None, defaults to an array of ones.
 epermH, epermV : array_like
As input, checked for size. If None, defaults to an array of ones.
 mpermH, mpermV : array_like
As input, checked for size. If None, defaults to an array of ones.
 isfullspace : bool
If True, the model is a fullspace (res, aniso, epermH, epermV, mpermM, and mpermV are in all layers the same).

empymod.utils.
check_frequency
(freq, res, aniso, epermH, epermV, mpermH, mpermV, verb)[source]¶ Calculate frequencydependent parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  freq : array_like
Frequencies f (Hz).
 res : array_like
Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.
 aniso : array_like
Anisotropies lambda = sqrt(rho_v/rho_h) (); #aniso = #res.
 epermH, epermV : array_like
Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (); #epermH = #epermV = #res.
 mpermH, mpermV : array_like
Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (); #mpermH = #mpermV = #res.
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  freq : float
Frequency, checked for size and assured min_freq.
 etaH, etaV : array
Parameters etaH/etaV, same size as provided resistivity.
 zetaH, zetaV : array
Parameters zetaH/zetaV, same size as provided resistivity.

empymod.utils.
check_hankel
(ht, htarg, verb)[source]¶ Check Hankel transform parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  ht : {‘fht’, ‘qwe’, ‘quad’}
Flag to choose the Hankel transform.
 htarg : str or filter from empymod.filters or array_like,
Depends on the value for
ht
. verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  ht, htarg
Checked if valid and set to defaults if not provided.

empymod.utils.
check_opt
(opt, loop, ht, htarg, verb)[source]¶ Check optimization parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  opt : {None, ‘parallel’}
Optimization flag; use
numexpr
or not. loop : {None, ‘freq’, ‘off’}
Loop flag.
 ht : str
Flag to choose the Hankel transform.
 htarg : array_like,
Depends on the value for
ht
. verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  use_ne_eval : bool
Boolean if to use
numexpr
. loop_freq : bool
Boolean if to loop over frequencies.
 loop_off : bool
Boolean if to loop over offsets.

empymod.utils.
check_dipole
(inp, name, verb)[source]¶ Check dipole parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  inp : list of floats or arrays
Pole coordinates (m): [polex, poley, polez].
 name : str, {‘src’, ‘rec’}
Poletype.
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  inp : list
List of pole coordinates [x, y, z].
 ninp : int
Number of inpelements

empymod.utils.
check_bipole
(inp, name)[source]¶ Check di/bipole parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  inp : list of floats or arrays
Coordinates of inp (m): [dipolex, dipoley, dipolez, azimuth, dip] or. [bipolex0, bipolex1, bipoley0, bipoley1, bipolez0, bipolez1].
 name : str, {‘src’, ‘rec’}
Poletype.
Returns:  inp : list
As input, checked for type and length.
 ninp : int
Number of inp.
 ninpz : int
Number of inp depths (ninpz is either 1 or ninp).
 isdipole : bool
True if inp is a dipole.

empymod.utils.
check_ab
(ab, verb)[source]¶ Check sourcereceiver configuration.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  ab : int
Sourcereceiver configuration.
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  ab_calc : int
Adjusted sourcereceiver configuration using reciprocity.
 msrc, mrec : bool
If True, src/rec is magnetic; if False, src/rec is electric.

empymod.utils.
check_solution
(solution, signal, ab, msrc, mrec)[source]¶ Check required solution with parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  solution : str
String to define analytical solution.
 signal : {None, 0, 1, 1}
 Source signal:
 None: Frequencydomain response
 1 : Switchoff timedomain response
 0 : Impulse timedomain response
 +1 : Switchon timedomain response
 msrc, mrec : bool
True if src/rec is magnetic, else False.

empymod.utils.
get_abs
(msrc, mrec, srcazm, srcdip, recazm, recdip, verb)[source]¶ Get required ab’s for given angles.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  msrc, mrec : bool
True if src/rec is magnetic, else False.
 srcazm, recazm : float
Horizontal source/receiver angle (azimuth).
 srcdip, recdip : float
Vertical source/receiver angle (dip).
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  ab_calc : array of int
ab’s to calculate for this bipole.

empymod.utils.
get_geo_fact
(ab, srcazm, srcdip, recazm, recdip, msrc, mrec)[source]¶ Get required geometrical scaling factor for given angles.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  ab : int
Sourcereceiver configuration.
 srcazm, recazm : float
Horizontal source/receiver angle.
 srcdip, recdip : float
Vertical source/receiver angle.
Returns:  fact : float
Geometrical scaling factor.

empymod.utils.
get_azm_dip
(inp, iz, ninpz, intpts, isdipole, strength, name, verb)[source]¶ Get angles, interpolation weights and normalization weights.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  inp : list of floats or arrays
 Input coordinates (m):
 [x0, x1, y0, y1, z0, z1] (bipole of finite length)
 [x, y, z, azimuth, dip] (dipole, infinitesimal small)
 iz : int
Index of current di/bipole depth ().
 ninpz : int
Total number of di/bipole depths (ninpz = 1 or npinz = nsrc) ().
 intpts : int
Number of integration points for bipole ().
 isdipole : bool
Boolean if inp is a dipole.
 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.
 name : str, {‘src’, ‘rec’}
Poletype.
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  tout : list of floats or arrays
Dipole coordinates x, y, and z (m).
 azm : float or array of floats
Horizontal angle (azimuth).
 dip : float or array of floats
Vertical angle (dip).
 g_w : float or array of floats
Factors from Gaussian interpolation.
 intpts : int
As input, checked.
 inp_w : float or array of floats
Factors from source/receiver length and source strength.

empymod.utils.
get_off_ang
(src, rec, nsrc, nrec, verb)[source]¶ Get depths, offsets, angles, hence spatial input parameters.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  src, rec : list of floats or arrays
Source/receiver dipole coordinates x, y, and z (m).
 nsrc, nrec : int
Number of sources/receivers ().
 verb : {0, 1, 2, 3, 4}
Level of verbosity.
Returns:  off : array of floats
Offsets
 angle : array of floats
Angles

empymod.utils.
get_layer_nr
(inp, depth)[source]¶ Get number of layer in which inp resides.
Note: If zinp is on a layer interface, the layer above the interface is chosen.
This checkfunction is called from one of the modelling routines in
model
. Consult these modelling routines for a detailed description of the input parameters.Parameters:  inp : list of floats or arrays
Dipole coordinates (m)
 depth : array
Depths of layer interfaces.
Returns:  linp : int or array_like of int
Layer number(s) in which inp resides (plural only if bipole).
 zinp : float or array
inp[2] (depths).

empymod.utils.
printstartfinish
(verb, inp=None, kcount=None)[source]¶ Print start and finish with time measure and kernel count.

empymod.utils.
conv_warning
(conv, targ, name, verb)[source]¶ Print error if QWE/QUAD did not converge at least once.

empymod.utils.
set_minimum
(min_freq=None, min_time=None, min_off=None, min_res=None, min_angle=None)[source]¶ Set minimum values of parameters.
The given parameters are set to its minimum value if they are smaller.
Parameters:  min_freq : float, optional
Minimum frequency [Hz] (default 1e20 Hz).
 min_time : float, optional
Minimum time [s] (default 1e20 s).
 min_off : float, optional
Minimum offset [m] (default 1e3 m). Also used to round src & reccoordinates.
 min_res : float, optional
Minimum horizontal and vertical resistivity [Ohm.m] (default 1e20).
 min_angle : float, optional
Minimum angle [] (default 1e10).

empymod.utils.
get_minimum
()[source]¶ Return the current minimum values.
Returns:  min_vals : dict
Dictionary of current minimum values with keys
 min_freq : float
 min_time : float
 min_off : float
 min_res : float
 min_angle : float
For a full description of these options, see set_minimum.

empymod.utils.
spline_backwards_hankel
(ht, htarg, opt)[source]¶ Check opt if deprecated ‘spline’ is used.
Returns corrected htarg, opt. r

class
empymod.utils.
Report
(add_pckg=None, ncol=3, text_width=80, sort=False)[source]¶ Print date, time, and version information.
Use
scooby
to print date, time, and package version information in any environment (Jupyter notebook, IPython console, Python console, QT console), either as htmltable (notebook) or as plain text (anywhere).Always shown are the OS, number of CPU(s),
numpy
,scipy
,empymod
,sys.version
, and time/date.Additionally shown are, if they can be imported,
numexpr
,IPython
, andmatplotlib
. It also shows MKL information, if available.All modules provided in
add_pckg
are also shown.Parameters:  add_pckg : packages, optional
Package or list of packages to add to output information (must be imported beforehand).
 ncol : int, optional
Number of packagecolumns in html table (no effect in textversion); Defaults to 3.
 text_width : int, optional
The text width for nonHTML display modes
 sort : bool, optional
Sort the packages when the report is shown
Examples
>>> import pytest >>> import dateutil >>> from emg3d import Report >>> Report() # Default values >>> Report(pytest) # Provide additional package >>> Report([pytest, dateutil], ncol=5) # Set nr of columns