comet.pyhed.hed package

Submodules

comet.pyhed.hed.hed_bib module

Part of comet/pyhed/hed

comet.pyhed.hed.hed_bib.btp(u, model, rho, d, f, mode)[source]

Airspace only, internal function, for imput see downout.

Do not call directly.

comet.pyhed.hed.hed_bib.calcField(polar, rho, d, f, Ids, ftype, mode)[source]

Calculates field for a given dipole on given polar coords. Internally.

Internally used by makeField. Please be referred to the docstrings of makeField. And please use makeField directly!!!

comet.pyhed.hed.hed_bib.downout(u, model, rho, d, f, mode)[source]

Overall call function for recursive calculation.

Parameters:
  • u (np.ndarray) – Horizontal wavenumbers based on Hankel factors and horizontal tx-rx distance. Shape: (Hankel, n_rx)
  • model (np.ndarray) – Polar coords of the receiver pos (3, n).
  • rho (np.ndarray) – Resistivities for each layer (Ohm*m).
  • d (np.ndarray) – Thicknesses of each layer (m).
  • f (float) – Frequency (Hz).
  • mode (str) – Calculation for ‘te’, ‘tm’ or ‘tetm’ possible. Mode ‘te’ for closed loops and ‘tetm’ for grounded wires. Single ‘tm’ is only for debug.
Returns:

  • aa (np.ndarray) – Ratio of the partial wave amplitude A(z,u)/A(0,u)
  • aap (np.andarray) – Ratio of the partial wave amplitude A’(z,u)/A’(0,u)
  • bt (np.ndarray) – Admittance at the surface of the layerd halfspace
comet.pyhed.hed.hed_bib.downward(u, model, rho, d, f, mode)[source]

Downward attenuation.

Parameters:
  • u (np.ndarray) – Horizontal wavenumbers based on Hankel factors and horizontal tx-rx distance. Shape: (Hankel, n_rx)
  • model (np.ndarray) – Polar coords of the receiver pos (3, n).
  • rho (np.ndarray) – Resistivities for each layer (Ohm*m).
  • d (np.ndarray) – Thicknesses of each layer (m).
  • f (float) – Frequency (Hz).
Returns:

  • aa (np.ndarray) – Ratio of the partial wave amplitude A(z,u)/A(0,u)
  • aap (np.andarray) – Ratio of the partial wave amplitude A’(z,u)/A’(0,u)
  • bt (np.ndarray) – Admittance at the surface of the layerd halfspace
comet.pyhed.hed.hed_bib.efield_3D_hed_te(polar, u, aa, aap, bt, f, Ids)[source]

Calculation of electric field for transversal electric mode.

Computes the transversal electric induced electric field of a x-directed dipole at (0, 0, 0). Field shape (3, n) with x, y, z components for each reciever point in polar.

Internal function. Called by makeField if ftype == ‘E’ and mode in (‘te’, ‘tetm’).

Parameters:
  • polar (np.ndarray) – Polar coords of the receiver pos (3, n).
  • u (np.ndarray) – Horizontal wavenumbers based on Hankel factors and horizontal tx-rx distance. Shape: (Hankel, n_rx)
  • aa (np.ndarray) – Ratio of the partial wave amplitude A(z,u)/A(0,u)
  • aap (np.andarray) – Ratio of the partial wave amplitude A’(z,u)/A’(0,u)
  • bt (np.ndarray) – Admittance at the surface of the layerd halfspace
  • f (float) – Frequency (Hz).
  • Ids (float) – Dipole current * dipole length.
Returns:

field – Transversal electric component of the electric field of a x-directed dipole at (0, 0, 0). field.shape = polar.shape.
Return type:

np.ndarray
comet.pyhed.hed.hed_bib.hankelfc(order)[source]

Filter coefficients for hankel transformation by Anderson (1980)

comet.pyhed.hed.hed_bib.hfield_3D_hed_te(polar, u, aa, aap, bt, f, Ids)[source]

Calculation of magnetic field for transversal electric mode.

Computes the transversal electric induced magnetic field of a x-directed dipole at (0, 0, 0). Field shape (3, n) with x, y, z components for each reciever point in polar.

Internal function. Called by makeField if ftype == ‘H’ and mode in (‘te’, ‘tetm’).

Parameters:
  • polar (np.ndarray) – Polar coords of the receiver pos (3, n).
  • u (np.ndarray) – Horizontal wavenumbers based on Hankel factors and horizontal tx-rx distance. Shape: (Hankel, n_rx)
  • aa (np.ndarray) – Ratio of the partial wave amplitude A(z,u)/A(0,u)
  • aap (np.andarray) – Ratio of the partial wave amplitude A’(z,u)/A’(0,u)
  • bt (np.ndarray) – Admittance at the surface of the layerd halfspace
  • f (float) – Frequency (Hz).
  • Ids (float) – Dipole current * dipole length.
Returns:

field – Transversal electric induced magnetic field of a x-directed dipole at (0, 0, 0). field.shape = polar.shape.
Return type:

np.ndarray
comet.pyhed.hed.hed_bib.hfield_3D_hed_tm(polar, u, aa, aap, bt, f, Ids)[source]

Calculation of magnetic field for transversal magnetic mode.

Computes the transversal magnetic component of the magnetic field of a x-directed dipole at (0, 0, 0). Field shape (3, n) with x, y, z components for each reciever point in polar.

Internal function. Called by makeField if ftype == ‘H’ and mode in (‘tm’, ‘tetm’).

Parameters:
  • polar (np.ndarray) – Polar coords of the receiver pos (3, n).
  • u (np.ndarray) – Horizontal wavenumbers based on Hankel factors and horizontal tx-rx distance. Shape: (Hankel, n_rx)
  • aa (np.ndarray) – Ratio of the partial wave amplitude A(z,u)/A(0,u)
  • aap (np.andarray) – Ratio of the partial wave amplitude A’(z,u)/A’(0,u)
  • bt (np.ndarray) – Admittance at the surface of the layerd halfspace
  • f (float) – Frequency (Hz).
  • Ids (float) – Dipole current * dipole length.
Returns:

field – Transversal magnetic component of the magnetic field of a x-directed dipole at (0, 0, 0). field.shape = polar.shape.
Return type:

np.ndarray
comet.pyhed.hed.hed_bib.makeField(coords, rho_in, d_in, f=2000, Ids=1, pos=(0, 0, 0), angle=0, mode='te', inputType='M', ftype='B', cell_center=False, drop_tol=0.01, src_z=-0.01)[source]

Calculation of the electric or magnetic field of a horizontal electric dipole at position pos, pointing in a direction defined by angle on given cartesian coordinates.

Parameters:
  • coords (np.ndarray or string) – Reciever coords. Possible input types are numpy ndarrays for direct cartesian coordinates, ranges for (x, y, z) or pygimli Meshes.
  • rho_in (float or np.ndarray) – Float or Array of resistivity values for the 1d layered earth model. Airspace is at the level of the source dipole.
  • d_in (float or np.ndarray) – Layer thicknesses in m. As the lower halfspace is considered to have an infinite thickness, d_in is always one value short of rho_in (an empty list ar array or a 0 for homogenous halfspace.)
  • f (float [ 2000 ]) – Frequency (Hz).
  • Ids (float [ 1 ]) – Dipole current * dipole length. Used for simple scaling of the calcualted field.
  • pos (tuple of length 3 [ (0 ,0, 0) ]) – Absolute position of the source dipole in cartesian coordinates. Values for z are used for a shift of the airspace. Currently only sources at the upper halfpsce boundary are permitted.
  • angle (float [ 0 ]) – Rotation of the dipole with respect of an x-directed dipole counting positive clockwise.
  • mode (string [ ‘te’ ]) – For a closed loop consisting of a finite number of dipoles the total field can be seen as superposition of the transversal electric components of the single dipole fields (‘te’). For grounded dipoles ‘tetm’ is needed.
  • inputType (string [ ‘M’ ]) – Specifier for input coordinates. Possible choices are ‘M’ if coords is a pygimli mesh or file path to a pygimli mesh, ‘C’ if coords is a np.ndarray with ranges to build a meshgrid, or ‘V’ to indicate that coords is a vector of cartesian coordinates.
  • ftype (string [ ‘B’ ]) – Flag to control calculated field type. Possible choices are ‘E’, ‘B’ or ‘H’ (asuming B = 4e-7 * pi * H).
  • verbose (boolean [ False ]) – Turn on verbose mode.
  • cell_center (boolean [ False ]) – If coords is a pygimli mesh, there is the additional possibility to calulate the fields in the cell Centers, instead of the node positions.
  • drop_tol (float [ 1e-2 ]) – Singularity fix. All horizontal distances between drop_tol and the transmitter dipole are placed between the first reciever outside the tolerance and the tolerance, maintaining the correct order and angle. This has been very useful for later usage of the fields in FEM approaches.
  • src_z (float [-0.001]) – This is only used if grounded terms for an electric field are used. In this case the source has to be buried in order to get the correct results. Default is 1 cm (remember: z defined positive upwards). So in most cases this value should be negative.

comet.pyhed.hed.hed_para module

Part of comet/pyhed/hed

comet.pyhed.hed.hed_para.InterpolationWorker(num, pos_queue, out_queue, data, srcmeshName, outmeshName, outtype, verbose)[source]

MPI Worker used to interpolate fields to target source location.

comet.pyhed.hed.hed_para.SummationWorker(queueIn, queueSum, queueEnd, verbose)[source]

MPI Worker used to sum up single fields.

comet.pyhed.hed.hed_para.multiInterpolation(DipoleDataName, SrcMeshName, OutMesh, DipolePos=None, verbose=False)[source]

Call function for multiprocessing interpoaltion of dipole fields.

comet.pyhed.hed.libHED module

Part of comet/pyhed/hed

Earth class for calculation of dipole (HED) fields for 1d layered earth.

The algorithms in method calcFieldForLayer of HED class is partly taken from Kerry Key Dipole1D.f90 after the algorithms published in [Key2009G] (Appendix A).

Hankel factors of Hankelfc are based on the original values of Anderson (1990).

References:

[Key2009G](1, 2, 3) Key, K., 2009, 1D inversion of multicomponent, multifrequency marine (CSEM) data: Methodology and synthetic studies for resolving thin resistive layers: Geophysics.
class comet.pyhed.hed.libHED.HED(src_z=-0.01, src_theta=0.0, src_ids=1.0, config=None, timer=None, debug=False)[source]

Bases: object

calcFieldForLayer(rx_lay, lay_indices, lam, lam_2, lam_c, lam_c2, R_p, R_m, S_p, S_m, exp)[source]
calculate()[source]

Calculates the 1d layered earth recursive formula.

Calculates the recursive attenuation and reflection coefficients for each layer on basis of the given set of cylindrical coordinates.

Fills the variables R_p, R_m, r_p, r_s, S_p, S_m, s_p, s_m, hem_a, hem_b, hem_c, and hem_d. The used formulas correspond to equations A-6 to A-13 in [Key2009G] (Appendix A).

reflectionCoefficients(rx_lay, lam, lam_2, lam_c, lam_c2, exp)[source]

Calculation of the general reflection coefficients R+, R-, S+, and S- as stated in [Key2009G] (Appendix A, equations A-06 to A-09).

Computed from the air and halfspace, respectively, inward to the source layer.

setCoords(cartesian, nodes=True, drop_tol=0.01)[source]

Sets coordinates of the receiver for calculation.

All calcualtions will be performed in cylindrically coordinates.

Parameters:
  • cartesian (np.ndarray) – Cartesian coordinates (N points) of the receiver points with shape (3, N). Z is defined positive upwards.
  • drop_tol (float) – Tolerance in meter, where the horizontal src distance is capped to ensure a safe division (singularity fix). Distances smaller than drop_tol are distributed between droptol and 20% of the first value outside the droptol. Raises Exception if all points within drop_tol. Default value of 1cm.
setTheta(theta)[source]
class comet.pyhed.hed.libHED.World1D(rho=1000.0, thk=None, airspace_interface=0.0, f=2000.0)[source]

Bases: object

evalSrcIdx(src_depth)[source]

Evaluates in which layer the source is considered.

setFrequency(freq)[source]

Simple setter for frequency + implicit omega/sigma calculation

setRes(rho=1000.0, thk=None, air_resistivty=10000000000000.0)[source]

Sets the resisitvity model for the dipoles + calc sigma complex.

Parameters:
  • rho (float or array_like) – Resistivity distribution in Ohm*m. Airspace is considered to have 0 Ohm*m. The first entry of rho correspond to the first layer of the subsurface. The airspace interface is considered to be at z = 0 m which simplyfies the calculations. For offsets in z, a cordinate transformation has to be performed externally.
  • thk (float or array_like) – Layer thicknesses of each subsurface layer except the last.
class comet.pyhed.hed.libHED.hankelfc[source]

Bases: object

getFactors(string)[source]

Returns the requested set of Hankel factors.

Parameters:string ([str]) – Evaluates which Hankel factors the wavenumbers is calculated. Possible choices are sin, cos, j0, or j1.
Returns:factors – Hankel factors.
Return type:[np.ndarray]
getWavenumbes(string)[source]

Calcualtes the wavenumbers for the requested set of Hankel factors.

Parameters:string ([str]) – Evaluates for which Hankel factors the wavenumbers is calculated. Possible choices are sin, cos, j0, or j1.
Returns:wavenumbers – Normed wavenumber factors for evaluation of the Hankel integral. Divide by horizontal distance of the receiver to get horizontal wavenumber lambda = sqrt(k_x² + k_y²).
Return type:[np.ndarray]
class comet.pyhed.hed.libHED.wer_201_2018[source]

Bases: object

Hankel factors after Werthmüller 2018 implemented from the empymod package after consultation with Dieter Werthmüller. Thank you very much!

The filter coefficient are published in:

Werthmüller, D., K. Key, and E. Slob, 2019, A tool for designing digital filters for the Hankel and Fourier transforms in potential, diffusive, and wavefield modeling: 84(2), F47-F56; DOI: 10.1190/geo2018-0069.1

under the Apache 2.0 license.

Module contents

module comet/pyhed/hed

Init file for HED calculation routines inside pyhed.