Info¶
The latest version of this documentation can be found at https://empymod.readthedocs.io.
Installation & requirements¶
The easiest way to install the latest stable version of empymod is via conda:
> conda install c prisae empymod
or via pip:
> pip install empymod
Alternatively, you can download the latest version from GitHub and either add the path to empymod to your pythonpath variable, or install it in your python distribution via:
> python setup.py install
Required are python version 3.4 or higher and the modules NumPy and SciPy. If you want to run parts of the kernel in parallel, the module numexpr is required additionally.
If you are new to Python I recommend using a Python distribution, which will ensure that all dependencies are met, specifically properly compiled versions of NumPy and SciPy; I recommend using Anaconda (version 3.x; continuum.io/downloads). If you install Anaconda you can simply start the Anaconda Navigator, add the channel prisae and empymod will appear in the package list and can be installed with a click.
Usage¶
The main modelling routines is bipole, which can 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. The model is defined by horizontal resistivity and anisotropy, horizontal and vertical electric permittivities and horizontal and vertical magnetic permeabilities. By default, the electromagnetic response is normalized to to source and receiver of 1 m length, and source strength of 1 A.
A simple frequencydomain example, with most of the parameters left at the default value:
>>> 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)
>>> 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
Hankel : Fast Hankel Transform
> Filter : Key 201 (2009)
Hankel 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  0 : 10 [minmax; #]
: 0 0 0 0 0 0 0 0 0 0
> dip [°] : 0  0 : 10 [minmax; #]
: 0 0 0 0 0 0 0 0 0 0
Required ab's : 11
~
:: empymod END; runtime = 0:00:00.022349 :: 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]
Have a look at the notebooks for more extensive examples including figures.
Structure¶
 model.py: EM modelling routines.
 utils.py: Utilities for model such as checking input parameters.
 kernel.py: Kernel of empymod, calculates the wavenumberdomain electromagnetic response. Plus analytical, frequencydomain full and halfspace solutions.
 transform.py: Methods to carry out the required Hankel transform from wavenumber to frequency domain and Fourier transform from frequency to time domain.
 filters.py: Filters for the Fast Hankel Transform (FHT, [Anderson_1982]) and the Fourier Sine and Cosine Transforms [Anderson_1975].
Missing features¶
A list of things that should or could be added and improved:
 Kernel
 Include scipy.integrate.quad as an additional Hankel transform. There are cases when both QWE and FHT struggle, e.g. at very short offsets with very high frequencies (GPR).
 A cython or numba (pure C?) implementation of the kernel and the transform modules. Maybe not worth it, as it may improve speed, but decrease accessibility. Both at the same time would be nice. A fast Cversion for calculations (inversions), and a Pythonversion to tinker with for interested folks. (Probably combined with default parallelisation, removing the numexpr variant.)
 More modelling routines:
 Convolution with a wavelet for GPR (proper version of model.gpr, needs another HT/FT).
 Various sourcereceiver arrangements (loops etc).
 Load and save functions to store and load model, together with all information.
 Module to create Hankel filters (nice to have addition, mainly for educational purposes).
 GUI.
 Add a benchmark suite, e.g. http://asv.readthedocs.io, in addition to the testing suite.
Testing¶
The modeller comes with a test suite using pytest. If you want to run the tests, just install pytest and run it within the empymodtopdirectory.
> conda install pytest
> # or
> pip install pytest
> # and then
> cd to/the/empymod/folder
> ls
CHANGELOG.md empymod notebooks paper requirements.txt tests
docs LICENSE NOTICE README.rst setup.py THANKS.md
> # pytest will find the tests, which are located in the testsfolder.
> # simply run
> pytest
It should run all tests successfully. Please let me know if not!
Note that the installations via conda/pip do not have the testsuite included. To run the testsuite you must download empymod from GitHub.
Citation¶
I am in the process of publishing an article regarding empymod, and I will put the info here once it is reality. If you publish results for which you used empymod, please consider citing this article. Also consider citing [Hunziker_et_al_2015] and [Key_2012], without which empymod would not exist.
License¶
Copyright 20162017 Dieter Werthmüller
Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at
Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “AS IS” BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
See the LICENSEfile in the root directory for a full reprint of the Apache License.
Notice¶
This product includes software that was initially (till 01/2017) developed at The Mexican Institute of Petroleum IMP (Instituto Mexicano del Petróleo, http://www.imp.mx). The project was funded through The Mexican National Council of Science and Technology (Consejo Nacional de Ciencia y Tecnología, http://www.conacyt.mx).
This product is a derivative work of [Hunziker_et_al_2015] and [Key_2012], and their publicly available software:
 Hunziker, J., J. Thorbecke, and E. Slob, 2015, The electromagnetic response in a layered vertical transverse isotropic medium: A new look at an old problem: Geophysics, 80, F1F18; DOI: 10.1190/geo20130411.1; Software: software.seg.org/2015/0001.
 Key, K., 2012, Is the fast Hankel transform faster than quadrature?: Geophysics, 77, F21F30; DOI: 10.1190/GEO20110237.1; Software: software.seg.org/2012/0003.
Both pieces of software are published under the SEG disclaimer. Parts of the modeller emmod from Hunziker et al, 2015, is furthermore released under the Common Public License Version 1.0 (CPL). See the NOTICEfile in the root directory for more information and a reprint of the SEG disclaimer and the CPL.
Note on speed, memory, and accuracy¶
There is the usual tradeoff between speed, memory, and accuracy. Very generally speaking we can say that the FHT is faster than QWE, but QWE is much easier on memory usage. I doubt you will ever run into memory issues with QWE, whereas for FHT you might for ten thousands of offsets or hundreds of layers. Furthermore, QWE allows you to control the accuracy.
There are two optimisation possibilities included via the opt
flag:
parallelisation (opt='parallel'
) and spline interpolation
(opt='spline'
). They are switched off by default. The optimization
opt='parallel'
only affects speed and memory usage, whereas
opt='spline'
also affects precision!
I am sure empymod could be made much faster with cleverer coding style or with the likes of cython or numba. Suggestions and contributions are welcomed!
Depths, Rotation, and Bipole¶
Depths: Calculation of many source and receiver positions is fastest if they remain at the same depth, as they can be calculated in one kernelcall. If depths do change, one has to loop over them.
Rotation: Sources and receivers aligned along the principal axes x, y, and z can be calculated in one kernel call. For arbitrary aligned di or bipoles, 3 kernel calls are required. If source and receiver are arbitrary aligned, 3x3 hence 9 kernel calls are required.
Bipole: Bipole increase the calculation time by the amount of integration points used. For a source and a receiver bipole with each 5 integration points you need 5x5 hence 25 kernel calls. You can calculate it in 1 kernel call if you set both integration points to 1, and hence calculate the bipole as if they were dipoles at their centre.
Example: For 1 source and 10 receivers, all at the same depth, 1 kernel call is required. If all receivers are at different depths, 10 kernel calls are required. If you make source and receivers bipoles with 5 integration points, 250 kernel calls are required. If you rotate the source arbitrary horizontally, 500 kernel calls are required. If you rotate the receivers too, in the horizontal plane, 1‘000 kernel calls are required. If you rotate the receivers also vertically, 1‘500 kernel calls are required. If you rotate the source vertically too, 2‘250 kernel calls are required. So your calculation will take 2‘250 times longer! No matter how fast the kernel is, this will take a long time. Therefore carefully plan how precise you want to define your source and receiver bipoles.
source bipole  receiver bipole  

kernel calls  intpts  azimuth  dip  intpts  azimuth  dip  diff. z 
1  1  0/90  0/90  1  0/90  0/90  1 
10  1  0/90  0/90  1  0/90  0/90  10 
250  5  0/90  0/90  5  0/90  0/90  10 
500  5  arb.  0/90  5  0/90  0/90  10 
1000  5  arb.  0/90  5  arb.  0/90  10 
1500  5  arb.  0/90  5  arb.  arb.  10 
2250  5  arb.  arb.  5  arb.  arb.  10 
Parallelisation¶
If opt = 'parallel'
, a good dozen of the most timeconsuming statements are
calculated by using the numexpr package
(https://github.com/pydata/numexpr/wiki/NumexprUsersGuide). These statements
are all in the kernelfunctions greenfct, reflections, and fields, and
all involve \(\Gamma\) in one way or another, often calculating square
roots or exponentials. As \(\Gamma\) has dimensions (#frequencies,
#offsets, #layers, #lambdas), it can become fairly big.
The module numexpr uses by default all available cores up to a maximum of 8. You can change this behaviour to a lower or a higher value with the following command (in the example it is changed to 4):
>>> import numexpr
>>> numexpr.set_num_threads(4)
This parallelisation will make empymod faster if you calculate a lot of offsets/frequencies at once, but slower for few offsets/frequencies. Best practice is to check first which one is faster. (You can use the included jupyter notebookbenchmark.)
Spline interpolation¶
If opt = 'spline'
, the socalled lagged convolution or splined variant
of the FHT (depending on htarg
) or the splined version of the QWE are
applied. The spline option should be used with caution, as it is an
interpolation and therefore less precise than the nonspline version. However,
it significantly speeds up QWE, and massively speeds up FHT. (The
numexprversion of the spline option is slower than the pure spline one, and
therefore it is only possible to have either 'parallel'
or 'spline'
on.)
Setting opt = 'spline'
is generally faster. Good speedup is achieved for
QWE by setting maxint
as low as possible. Also, the higher nquad
is,
the higher the speedup will be. The variable pts_per_dec
has also some
influence. For FHT, big improvements are achieved for long FHTfilters and
for many offsets/frequencies (thousands). Additionally, spline minimizes
memory requirements a lot. Speedup is greater if all sourcereceiver angles
are identical.
FHT: Default for pts_per_dec = None
, which is the original lagged
convolution, where the spacing is defined by the filterbase, the transform is
carried out first followed by splineinterpolation. You can set this parameter
to an integer, which defines the number of points to evaluate per decade. In
this case the splineinterpolation is carried out first, followed by the
transformation. The original lagged convolution is generally the fastest for
a very good precision. However, by setting pts_per_dec
appropriately one
can achieve higher precision, normally at the cost of speed.
Warning
Keep in mind that it uses interpolation, and is therefore not as accurate as the nonspline version. Use with caution and always compare with the nonspline version if you can apply the splineversion to your problem at hand!
Be aware that the QWE and the FHTVersions for the frequencytotime transformation always use the splined version and always loop over offsets.
Looping¶
By default, you can calculate many offsets and many frequencies all in one go,
vectorized (for the FHT), which is the default. The loop
parameter gives
you the possibility to force looping over frequencies or offsets. This
parameter can have severe effects on both runtime and memory usage. Play around
with this factor to find the fastest version for your problem at hand. It
ALWAYS loops over frequencies if ht = 'QWE'
or if opt = 'spline'
. All
vectorized is very fast if there are few offsets or few frequencies. If there
are many offsets and many frequencies, looping over the smaller of the two will
be faster. Choosing the right looping together with opt = 'parallel'
can
have a huge influence.
Vertical components¶
It is advised to use xdirect = True
(the default) if source and receiver
are in the same layer to calculate
 the vertical electric field due to a vertical electric source,
 configurations that involve vertical magnetic components (source or receiver),
 all configurations when source and receiver depth are exactly the same.
The Hankel transforms methods are having sometimes difficulties transforming these functions.
FFTLog¶
FFTLog is the logarithmic analogue to the Fast Fourier Transform FFT originally proposed by [Talman_1978]. The code used by empymod was published in Appendix B of [Hamilton_2000] and is publicly available at casa.colorado.edu/~ajsh/FFTLog. From the FFTLogwebsite:
FFTLog is a set of fortran subroutines that compute the fast Fourier or Hankel (= FourierBessel) transform of a periodic sequence of logarithmically spaced points.
FFTlog can be used for the Hankel as well as for the Fourier Transform, but currently empymod uses it only for the Fourier transform. It uses a simplified version of the python implementation of FFTLog, pyfftlog (github.com/prisae/pyfftlog).
References ¶
[Anderson_1975]  (1, 2, 3, 4) Anderson, W.L., 1975, Improved digital filters for evaluating Fourier and Hankel transform integrals: USGS Unnumbered Series; http://pubs.usgs.gov/unnumbered/70045426/report.pdf. 
[Anderson_1979]  Anderson, W. L., 1979, Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering: Geophysics, 44, 1287–1305; DOI: 10.1190/1.1441007. 
[Anderson_1982]  (1, 2, 3, 4, 5, 6) Anderson, W. L., 1982, Fast Hankel transforms using related and lagged convolutions: ACM Trans. on Math. Softw. (TOMS), 8, 344–368; DOI: 10.1145/356012.356014. 
[Gosh_1971]  Ghosh, D. P., 1971, The application of linear filter theory to the direct interpretation of geoelectrical resistivity sounding measurements: Geophysical Prospecting, 19, 192–217; DOI: 10.1111/j.13652478.1971.tb00593.x. 
[Hamilton_2000]  (1, 2) Hamilton, A. J. S., 2000, Uncorrelated modes of the nonlinear power spectrum: Monthly Notices of the Royal Astronomical Society, 312, pages 257284; DOI: 10.1046/j.13658711.2000.03071.x; Website of FFTLog: casa.colorado.edu/~ajsh/FFTLog. 
[Hunziker_et_al_2015]  (1, 2, 3, 4, 5, 6, 7, 8, 9) Hunziker, J., J. Thorbecke, and E. Slob, 2015, The electromagnetic response in a layered vertical transverse isotropic medium: A new look at an old problem: Geophysics, 80, F1–F18; DOI: 10.1190/geo20130411.1; Software: software.seg.org/2015/0001. 
[Key_2009]  (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22) Key, K., 2009, 1D inversion of multicomponent, multifrequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers: Geophysics, 74, F9–F20; DOI: 10.1190/1.3058434. Software: marineemlab.ucsd.edu/Projects/Occam/1DCSEM. 
[Key_2012]  (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26) Key, K., 2012, Is the fast Hankel transform faster than quadrature?: Geophysics, 77, F21–F30; DOI: 10.1190/GEO20110237.1; Software: software.seg.org/2012/0003. 
[Kong_2007]  (1, 2, 3, 4, 5) Kong, F. N., 2007, Hankel transform filters for dipole antenna radiation in a conductive medium: Geophysical Prospecting, 55, 83–89; DOI: 10.1111/j.13652478.2006.00585.x. 
[Shanks_1955]  Shanks, D., 1955, Nonlinear transformations of divergent and slowly convergent sequences: Journal of Mathematics and Physics, 34, 1–42; DOI: 10.1002/sapm19553411. 
[Slob_et_al_2010]  (1, 2) Slob, E., J. Hunziker, and W. A. Mulder, 2010, Green’s tensors for the diffusive electric field in a VTI halfspace: PIER, 107, 1–20: DOI: 10.2528/PIER10052807. 
[Talman_1978]  Talman, J. D., 1978, Numerical Fourier and Bessel transforms in logarithmic variables: Journal of Computational Physics, 29, pages 3548; DOI: 10.1016/00219991(78)901079. 
[Trefethen_2000]  Trefethen, L. N., 2000, Spectral methods in MATLAB: Society for Industrial and Applied Mathematics (SIAM), volume 10 of Software, Environments, and Tools, chapter 12, page 129; DOI: 10.1137/1.9780898719598.ch12. 
[Weniger_1989]  Weniger, E. J., 1989, Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series: Computer Physics Reports, 10, 189–371; arXiv: abs/math/0306302. 
[Wynn_1956]  Wynn, P., 1956, On a device for computing the \(e_m(S_n)\) tranformation: Math. Comput., 10, 91–96; DOI: 10.1090/S00255718195600840566. 
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.
 These principal 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 after fem.
Two further routines are shortcuts for frequency and timedomain dipoles, respectively, and mainly in for legacy reasons:
 frequency: Shortcut of dipole for frequencydomain calculation.
 time: Shortcut of dipole for timedomain calculation.
Two more routines are more kind of examples and cannot be regarded stable; they can serve as template to create your own routines:
 gpr: Calculate the GroundPenetrating Radar (GPR) response.
 wavenumber: Calculate the electromagnetic wavenumberdomain solution.

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=True, ht='fht', htarg=None, ft='sin', ftarg=None, opt=None, loop=None, verb=2)¶ 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.
depth : list
Absolute layer interfaces z (m); #depth = #res  1 (excluding +/ infinity).
res : array_like
Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.
freqtime : array_like
Frequencies f (Hz) if signal == None, else times t (s).
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, 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. Defaults to True.
ht : {‘fht’, ‘qwe’}, optional
Flag to choose either the Fast Hankel Transform (FHT) or the QuadratureWithExtrapolation (QWE) for the Hankel transform. Defaults to ‘fht’.
htarg : str or filter from empymod.filters or array_like, optional
 Depends on the value for ht:
If ht = ‘fht’: array containing: [filter, pts_per_dec]:
 filter: string of filter name in empymod.filters or
 the filter method itself. (default: empymod.filters.key_201_2009())
 pts_per_dec: points per decade (only relevant if spline=True)
 If none, standard lagged convolution is used.
 (default: None)
If ht = ‘qwe’: array containing: [rtol, atol, nquad, maxint, pts_per_dec]:
 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 (only relevant if opt=’spline’) (default: 80)
All are optional, you only have to maintain the order. To only change nquad to 11 and use the defaults otherwise, you can provide htarg=[‘’, ‘’, 11].
ft : {‘sin’, ‘cos’, ‘qwe’, ‘fftlog’}, optional
Only used if signal != None. Flag to choose either the Sine or CosineFilter, the QuadratureWithExtrapolation (QWE), or FFTLog for the Fourier transform. Defaults to ‘sin’.
ftarg : str or filter from empymod.filters or array_like, optional
 Only used if signal !=None. Depends on the value for ft:
If ft = ‘sin’ or ‘cos’: array containing: [filter, pts_per_dec]:
 filter: 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. If none, standard lagged
 convolution is used. (Default: None)
If ft = ‘qwe’: array containing: [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 (only relevant if spline=True) (default: 20)
All are optional, you only have to maintain the order. To only change nquad to 11 and use the defaults otherwise, you can provide ftarg=[‘’, ‘’, 11].
If ft = ‘fftlog’: array containing: [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
All are optional, you only have to maintain the order. To only change add_dec to [1, 1] and use the defaults otherwise, you can provide ftarg=[‘’, [1, 1]].
opt : {None, ‘parallel’, ‘spline’}, 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 contain Gamma in one way or another, which has dimensions (#frequencies, #offsets, #layers, #lambdas), therefore can grow pretty big. The module numexpr uses by default all available cores up to a maximum of 8. You can change this behaviour to your desired number of threads nthreads with numexpr.set_num_threads(nthreads).
 If ‘spline’, the lagged convolution or splined variant of the FHT or the splined version of the QWE are used. Use with caution and check with the nonspline version for a specific problem. (Can be faster, slower, or plainly wrong, as it uses interpolation.) If spline is set it will make use of the parameter pts_per_dec that can be defined in htarg. If pts_per_dec is not set for FHT, then the lagged version is used, else the splined.
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 B [T] (not H [A/m]!).
In the case of the impulse timedomain response, the unit is further divided by seconds [1/s].
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)].
The shape of EM is (nfreq, nrec, nsrc). However, single dimensions are removed.
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 Hankel : Fast Hankel Transform > Filter : Key 201 (2009) Hankel 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  0 : 10 [minmax; #] : 0 0 0 0 0 0 0 0 0 0 > dip [°] : 0  0 : 10 [minmax; #] : 0 0 0 0 0 0 0 0 0 0 Required ab's : 11 ~ :: empymod END; runtime = 0:00:00.022349 :: 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=True, ht='fht', htarg=None, ft='sin', ftarg=None, opt=None, loop=None, verb=2)¶ 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 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 (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.
depth : list
Absolute layer interfaces z (m); #depth = #res  1 (excluding +/ infinity).
res : array_like
Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.
freqtime : array_like
Frequencies f (Hz) if signal == None, else times t (s).
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, 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. Defaults to True.
ht : {‘fht’, ‘qwe’}, optional
Flag to choose either the Fast Hankel Transform (FHT) or the QuadratureWithExtrapolation (QWE) for the Hankel transform. Defaults to ‘fht’.
htarg : str or filter from empymod.filters or array_like, optional
 Depends on the value for ht:
If ht = ‘fht’: array containing: [filter, pts_per_dec]:
 filter: string of filter name in empymod.filters or
 the filter method itself. (default: empymod.filters.key_201_2009())
 pts_per_dec: points per decade (only relevant if spline=True)
 If none, standard lagged convolution is used.
 (default: None)
If ht = ‘qwe’: array containing: [rtol, atol, nquad, maxint, pts_per_dec]:
 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 (only relevant if opt=’spline’) (default: 80)
All are optional, you only have to maintain the order. To only change nquad to 11 and use the defaults otherwise, you can provide htarg=[‘’, ‘’, 11].
ft : {‘sin’, ‘cos’, ‘qwe’, ‘fftlog’}, optional
Only used if signal != None. Flag to choose either the Sine or CosineFilter, the QuadratureWithExtrapolation (QWE), or FFTLog for the Fourier transform. Defaults to ‘sin’.
ftarg : str or filter from empymod.filters or array_like, optional
 Only used if signal !=None. Depends on the value for ft:
If ft = ‘sin’ or ‘cos’: array containing: [filter, pts_per_dec]:
 filter: 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. If none, standard lagged
 convolution is used. (Default: None)
If ft = ‘qwe’: array containing: [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 (only relevant if spline=True) (default: 20)
All are optional, you only have to maintain the order. To only change nquad to 11 and use the defaults otherwise, you can provide ftarg=[‘’, ‘’, 11].
If ft = ‘fftlog’: array containing: [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
All are optional, you only have to maintain the order. To only change add_dec to [1, 1] and use the defaults otherwise, you can provide ftarg=[‘’, [1, 1]].
opt : {None, ‘parallel’, ‘spline’}, 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 contain Gamma in one way or another, which has dimensions (#frequencies, #offsets, #layers, #lambdas), therefore can grow pretty big. The module numexpr uses by default all available cores up to a maximum of 8. You can change this behaviour to your desired number of threads nthreads with numexpr.set_num_threads(nthreads).
 If ‘spline’, the lagged convolution or splined variant of the FHT or the splined version of the QWE are used. Use with caution and check with the nonspline version for a specific problem. (Can be faster, slower, or plainly wrong, as it uses interpolation.) If spline is set it will make use of the parameter pts_per_dec that can be defined in htarg. If pts_per_dec is not set for FHT, then the lagged version is used, else the splined.
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 B [T] (not H [A/m]!).
In the case of the impulse timedomain response, the unit is further divided by seconds [1/s].
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.
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.
frequency
(src, rec, depth, res, freq, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, xdirect=True, ht='fht', htarg=None, opt=None, loop=None, verb=2)¶ Return the frequencydomain EM field due to a dipole source.
This is a shortcut for frequencydomain modelling using dipole (mainly for legacy reasons).
See dipole for info and a description of input and output parameters. Only difference is that frequency here corresponds to freqtime in dipole.
See also
Examples
>>> import numpy as np >>> from empymod import frequency >>> 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 = frequency(src, rec, depth, res, freq=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.
time
(src, rec, depth, res, time, ab=11, signal=0, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, xdirect=True, ht='fht', htarg=None, ft='sin', ftarg=None, opt=None, loop='off', verb=2)¶ Return the timedomain EM field due to a dipole source.
This is a shortcut for timedomain modelling using dipole (mainly for legacy reasons).
See dipole for info and a description of input and output parameters. Only difference is that time here corresponds to freqtime in dipole.
See also
Examples
>>> import numpy as np >>> from empymod import time >>> 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 = time(src, rec, depth, res, time=1, verb=0) >>> print(EMfield) [ 4.23754930e11 3.13805193e11 1.98884433e11 1.14387827e11 6.34605628e12 3.54905259e12 2.03906739e12 1.20569287e12 7.31746271e13 4.55825907e13]

empymod.model.
gpr
(src, rec, depth, res, fc=250, ab=11, gain=None, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, xdirect=True, ht='fht', htarg=None, opt=None, loop='off', verb=2)¶ Return the GroundPenetrating Radar signal.
THIS FUNCTION IS IN DEVELOPMENT, USE WITH CAUTION.
Or in other words it is merely an example how one could calculate the GPRresponse. However, the currently included FHT and QWE struggle for these high frequencies, and another Hankel transform has to be included to make GPR work properly (e.g. scipy.integrate.quad).
 QWE is slow, but does a pretty good job except for very short offsets: only direct wave for offset < 0.1 m, trianglelike noise at later times.
 FHT is fast. Airwave, direct wave and first reflection are well visible, but afterwards it is very noisy.
A lot is still hardcoded in this routine, for instance the frequencyrange used to calculate the response.
For input parameters see frequency, except for:
Parameters: fc : float
Centre frequency of GPRsignal (MHz). Sensible values are between 10 MHz and 3000 MHz.
gain : float
Power of gain function. If None, no gain is applied.
Returns: t : array
Times (s)
gprEM : ndarray
GPR response

empymod.model.
wavenumber
(src, rec, depth, res, freq, wavenumber, ab=11, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, xdirect=True, verb=2)¶ 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.
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.
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. Defaults to True.
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 wavenumber >>> 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 = wavenumber(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_spline, use_ne_eval, msrc, mrec, loop_freq, loop_off, conv=True)¶ 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)¶ 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, frequencydomain 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 [Hunziker_et_al_2015], 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)¶ Calculate wavenumber domain solution.
Return the wavenumber domain solutions PJ0, PJ1, and PJ0b, which have to be transformed with a Hankel transform to the frequency domain. PJ0/PJ0b and PJ1 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 [Hunziker_et_al_2015], and equally loosely to the file kxwmod.c.
[Hunziker_et_al_2015] 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 and PJ0b could theoretically be added here into one, and then be transformed in one go. However, PJ0b has to be multiplied by factAng 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)¶ 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)¶ 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 [Hunziker_et_al_2015], 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, and Gin62.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)¶ 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 [Hunziker_et_al_2015], and loosely to the corresponding files Gamma.F90, Wprop.F90, Ptotalx.F90, Ptotalxm.F90, Ptotaly.F90, Ptotalym.F90, Ptotalz.F90, and Ptotalzm.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)¶ Calculate Rp, Rm.
\[R^\pm_n, \bar{R}^\pm_n\]This function corresponds to equations 64/65 and A11/A12 in [Hunziker_et_al_2015], and loosely to the corresponding files Rmin.F90 and Rplus.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)¶ 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 [Hunziker_et_al_2015], and loosely to the corresponding files Pdownmin.F90, Pdownplus.F90, Pupmin.F90, and Pdownmin.F90.
This function is called from the function
kernel.greenfct
.

empymod.kernel.
halfspace
(xco, yco, zsrc, zrec, res, freq, aniso=1, ab=11)¶ Return frequencyspace domain VTI halfspace solution.
Calculates the frequencyspace domain electromagnetic response for a halfspace below air using the diffusive approximation, as given in [Slob_et_al_2010].
This routine is not strictly part of empymod and not used by it. However, it can be useful to compare the code to the analytical solution.
There are a few known typos in the equations of [Slob_et_al_2010]. Write the authors to receive an updated version!
This could be integrated into empymod by checking if the toplayer is a very resistive layer, hence air, and the rest is a halfspace, and then calling this function instead of wavenumber. (Similar to the way fullspace is incorporated if all layer parameters are identical.) The timespace domain solution could be implemented as well.
Parameters: xco, yco : array
Inline and crossline coordinates (m)
zsrc, zrec : float
Source and receiver depth (m)
res : float or array
Halfspace resistivity (Ohm.m)
freq : float
Frequency (Hz)
aniso : float, optional
Anisotropy (), default = 1
ab : int, optional
SrcRec config, default = 11; {11, 12, 13, 21, 22, 23, 31, 32, 33}
Returns: EM halfspace solution
Examples
>>> from empymod.kernel import halfspace >>> EM = halfspace(1000, 0, 10, 1, 10, 1) >>> print('HS response : ', EM) HS response : (3.02186073352e093.87322421836e10j)
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 FHT Hankel and Fourier transforms are based on source files (specified in each function) from the source code distributed with [Key_2012], 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, angle, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, fhtarg, use_spline, use_ne_eval, msrc, mrec)¶ Hankel Transform using the Fast Hankel Transform.
The Fast Hankel Transform is a Digital Filter Method, introduced to geophysics by [Gosh_1971], and made popular and widespread by [Anderson_1975], [Anderson_1979], [Anderson_1982].
This implementation of the FHT follows [Key_2012], equation 6. Without going into the mathematical details (which can be found in any of the above papers) and following [Key_2012], the FHT 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 [Key_2012].
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 FHT, this is 1.
conv : bool
Only relevant for QWE, not for FHT.

empymod.transform.
hqwe
(zsrc, zrec, lsrc, lrec, off, angle, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, qweargs, use_spline, use_ne_eval, msrc, mrec)¶ Hankel Transform using QuadratureWithExtrapolation.
QuadratureWithExtrapolation was introduced to geophysics by [Key_2012]. It is one of many socalled ISE methods to solve Hankel Transforms, where ISE stands for Integration, Summation, and Extrapolation.
Following [Key_2012], 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 FHT (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 ([Shanks_1955], [Wynn_1956]; implemented with algorithms from [Trefethen_2000] and [Weniger_1989]).
This function is based on get_CSEM1D_FD_QWE.m, qwe.m, and getBesselWeights.m from the source code distributed with [Key_2012].
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 converged. If not, maxint might have to be set higher.

empymod.transform.
fft
(fEM, time, freq, ftarg)¶ Fourier Transform using a Cosine or a Sinefilter.
It follows the Filter methodology [Anderson_1975], 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 [Key_2012].
Returns: tEM : array
Returns timedomain EM response of fEM for given time.
conv : bool
Only relevant for QWE, not for FFT.

empymod.transform.
fqwe
(fEM, time, freq, qweargs)¶ Fourier Transform using QuadratureWithExtrapolation.
It follows the QWE methodology [Key_2012] 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 [Key_2012].
Returns: tEM : array
Returns timedomain EM response of fEM for given time.
conv : bool
If true, QWE converged. If not, maxint might have to be set higher.

empymod.transform.
fftlog
(fEM, time, freq, ftarg)¶ Fourier Transform using FFTLog.
FFTLog is the logarithmic analogue to the Fast Fourier Transform FFT. FFTLog was presented in Appendix B of [Hamilton_2000] and published at <http://casa.colorado.edu/~ajsh/FFTLog>.
This function uses a simplified version of pyfftlog, which is a pythonversion of FFTLog. For more details regarding pyfftlog see <https://github.com/prisae/pyfftlog>.
Not the full flexibility of FFTLog is available here: Only the logarithmic FFT (fftl in FFTLog), not the Hankel transform (fht in FFTLog). Furthermore, the following parameters are fixed:
 mu = 0.5 (sinetransform)
 kr = 1 (initial value)
 kropt = 1 (silently adjusts kr)
 dir = 1 (forward)
Furthermore, q is restricted to 1 <= q <= 1.
I am trying to get FFTLog into scipy. If this happens the current implementation will be replaced by the scipy.fftpack.fftlogversion.
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 given time.
conv : bool
Only relevant for QWE, not for FFTLog.

empymod.transform.
qwe
(rtol, atol, maxint, inp, intervals, lambd=None, off=None, factAng=None)¶ QuadratureWithExtrapolation.
This is the kernel of the QWE method, used for the Hankel (hqwe) and the Fourier (fqwe) Transforms. See hqwe for an extensive description.
This function is based on qwe.m from the source code distributed with [Key_2012].

empymod.transform.
get_spline_values
(filt, inp, nr_per_dec=None)¶ Return required calculation points.

empymod.transform.
fhti
(rmin, rmax, n, q)¶ Return parameters required for FFTLog.
filters
– Digital Filters for FHT¶
Filters for the Fast Hankel Transform (FHT, [Anderson_1982]) and the Fourier Sine and Cosine Transforms [Anderson_1975].
To calculate the fhtfilter.factor I used
np.around(np.average(fhtfilter.base[1:]/fhtfilter.base[:1]), 15)
The filters kong_61_2007 and kong_241_2007 from [Kong_2007], and key_101_2009, key_201_2009, key_401_2009, key_81_CosSin_2009, key_241_CosSin_2009, and key_601_CosSin_2009 from [Key_2009] are taken from DIPOLE1D, [Key_2009], which can be downloaded at marineemlab.ucsd.edu/Projects/Occam/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 [Anderson_1982] and key_51_2012, key_101_2012, key_201_2012, key_101_CosSin_2012, and key_201_CosSin_2012, all from [Key_2012], are taken from the software distributed with [Key_2012] and available at software.seg.org/2012/0003. These filters are distributed under the SEG license.

class
empymod.filters.
DigitalFilter
(name)¶ Simple Class for Digital Filters.

empymod.filters.
anderson_801_1982
()¶ Anderson 801: [Anderson_1982].
Anderson 801 pt filter, as published in [Anderson_1982]; taken from file wa801Hankel.txt from [Key_2012], published by the Society of Exploration Geophysicists; software.seg.org/2012/0003. License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_101_2009
()¶ Key 101 2009: [Key_2009].
Key 101 pt filter, as published in [Key_2009]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
key_101_2012
()¶ Key 101 2012: [Key_2012].
Key 101 pt filter, taken from file kk101Hankel.txt from [Key_2012], published by the Society of Exploration Geophysicists; software.seg.org/2012/0003. License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_101_CosSin_2012
()¶ Key 101 CosSin 2012: [Key_2012].
Key 101 pt filter, taken from file kk101CosSin.txt from [Key_2012], published by the Society of Exploration Geophysicists; software.seg.org/2012/0003. License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_201_2009
()¶ Key 201 2009: [Key_2009].
Key 201 pt filter, as published in [Key_2009]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
key_201_2012
()¶ Key 201 2012: [Key_2012].
Key 201 pt filter, taken from file kk201Hankel.txt from [Key_2012], published by the Society of Exploration Geophysicists; software.seg.org/2012/0003. License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_201_CosSin_2012
()¶ Key 201 CosSin 2012: [Key_2012].
Key 201 pt filter, taken from file kk201CosSin.txt from [Key_2012], published by the Society of Exploration Geophysicists; software.seg.org/2012/0003. License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_241_CosSin_2009
()¶ Key 241 CosSin 2009: [Key_2009].
Key 241 pt filter, as published in [Key_2009]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
key_401_2009
()¶ Key 401 2009: [Key_2009].
Key 401 pt filter, as published in [Key_2009]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
key_51_2012
()¶ Key 51 2012: [Key_2012].
Key 51 pt filter, taken from file kk51Hankel.txt from [Key_2012], published by the Society of Exploration Geophysicists; software.seg.org/2012/0003. License: http://software.seg.org/disclaimer.txt.

empymod.filters.
key_601_CosSin_2009
()¶ Key 601 CosSin 2009: [Key_2009].
Key 601 pt filter, as published in [Key_2009]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
key_81_CosSin_2009
()¶ Key 81 CosSin 2009: [Key_2009].
Key 81 pt filter, as published in [Key_2009]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
kong_241_2007
()¶ Kong 241: [Kong_2007].
Kong 241 pt filter, as published in [Kong_2007]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.

empymod.filters.
kong_61_2007
()¶ Kong 61: [Kong_2007].
Kong 61 pt filter, as published in [Kong_2007]; taken from file FilterModules.f90 from [Key_2009], available on marineemlab.ucsd.edu/Projects/Occam/1DCSEM. License: Apache License, Version 2.0, http://www.apache.org/licenses/LICENSE2.0.
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
¶ 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
(time, signal, ft, ftarg, verb)¶ 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’}
Flag for Fourier transform, only used if signal != None.
ftarg : str or filter from empymod.filters or array_like,
Only used if signal !=None. Depends on the value for ft:
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, verb)¶ 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.
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)¶ 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)¶ 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’}
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)¶ 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’, ‘spline’}
Optimization flag.
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_spline : bool
Boolean if to use spline interpolation.
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)¶ 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)¶ 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)¶ 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.
get_abs
(msrc, mrec, srcazm, srcdip, recazm, recdip, verb)¶ 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)¶ 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)¶ 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)¶ 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)¶ 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)¶ Print start and finish with time measure and kernel count.

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