Module utils

Utility functions to interact with telewavesim modules.

telewavesim.utils.set_iso_tensor(a, b)

Function to generate tensor for isotropic material.

Parameters

• a (float) – P-wave velocity (km/s)

• b (float) – S-wave velocity (km/s)

Returns

cc: Elastic tensor (GPa /density) (shape (3, 3, 3, 3))

Return type

(np.ndarray)

telewavesim.utils.set_tri_tensor(a, b, tr, pl, ani)

Function to generate tensor for transverse isotropy. The tensor is rotated using the trend and plunge of the symmetry axis.

Parameters

• a (float) – P-wave velocity (km/s)

• b (float) – S-wave velocity (km/s)

• tr (float) – Trend angle of symmetry axis (degree)

• pl (float) – Plunge angle of symmetry axis (degree)

• ani (float) – Percent anisotropy

Returns

cc: Elastic tensor (GPa / kg/m^3) (shape (3, 3, 3, 3))

Return type

(np.ndarray)

telewavesim.utils.set_aniso_tensor(tr, pl, typ='atg')

Function to generate tensor for anisotropic minerals. The tensor is rotated using the trend and plunge of the symmetry axis.

Parameters

• tr (float) – Trend angle of symmetry axis (degree)

• pl (float) – Plunge angle of symmetry axis (degree)

• type (str, optional) – Type of elastic material

Returns

Tuple containing:

• cc (np.ndarray): Elastic tensor (GPa / kg/m^3) (shape (3, 3, 3, 3))

• rho (float): Density (kg/m^3)

Return type

(tuple)

telewavesim.utils.full_3x3_to_Voigt_6_index(i, j)

Conversion of tensor to Voigt notation for indices

telewavesim.utils.voigt2cc(C)

Convert the Voigt representation of the stiffness matrix to the full 3x3x3x3 tensor representation.

Parameters

C (np.ndarray) – Stiffness matrix (shape (6, 6))

Returns

cc: Elastic tensor (shape (3, 3, 3, 3))

Return type

(np.ndarray)

telewavesim.utils.cc2voigt(cc)

Convert from the full 3x3x3x3 tensor representation to the Voigt notation of the stiffness matrix.

Parameters

cc (np.ndarray) – Elastic tensor (shape (3, 3, 3, 3))

Returns

C: Stiffness matrix (shape (6, 6))

Return type

(np.ndarray)

telewavesim.utils.VRH_average(C)

Performs a Voigt-Reuss-Hill average of the anisotropic stifness matrix to the bulk modulus K and the shear modulus G.

Parameters

C (np.ndarray) – Stiffness matrix (shape (6, 6))

Returns

Tuple containing:

• Kvoigt (float): Voigt average bulk modulus (GPa)

• Gvoigt (float): Voigt average shear modulus (GPa)

• Kreuss (float): Reuss average bulk modulus (GPa)

• Greuss (float): Reuss average shear modulus (GPa)

• Kvrh (float): Voigt-Reuss-Hill average bulk modulus (GPa)

• Gvrh (float): Voigt-Reuss-Hill average shear modulus (GPa)

Return type

(tuple)

Example

>>> from telewavesim import utils
>>> cc, rho = utils.set_aniso_tensor(0., 0., typ='atg')
>>> C = utils.cc2voigt(cc)
>>> utils.VRH_average(C*rho)
(75655555555.555557, 48113333333.333336, 61245706544.967415,
28835098086.844658, 68450631050.26149, 38474215710.088997)

telewavesim.utils.mod2vel(K, G, rho)

Calculates the isotropic P and S wave velocities from given bulk (K) and shear (G) moduli and density (rho) in kg/m^3

Parameters

• K (float) – Bulk modulus (GPa)

• G (float) – Shear modulus (GPa)

• rho (float) – Density (kg/m^3)

Returns

tuple containing:

• Vp (float): P-wave velocity (m/s)

• Vs (float): S-wave velocity (m/s)

Return type

(tuple)

Example

>>> from telewavesim import utils
>>> cc, rho = utils.set_aniso_tensor(0., 0., typ='atg')
>>> C = utils.cc2voigt(cc)
>>> K, G = utils.VRH_average(C*rho)[4:6]
>>> utils.mod2vel(K, G, rho)
(6760.617471753726, 3832.0771334254896)

telewavesim.utils.rot_tensor(a, alpha, beta, gam)

Performs a rotation of the tensor cc (c_ijkl) about three angles (alpha, beta, gamma)

Parameters

• a (np.ndarray) – Elastic tensor with shape (3, 3, 3, 3)

• alpha (float) – Angle in radians

• beta (float) – Angle in radians

• gam (float) – Angle in radians

Returns

aa: Rotated tensor with shape (3, 3, 3, 3)

Return type

(np.ndarray)

Note

The three angles (alpha, beta, gam) correspond to rotation about the x_2, x_3, x_1 axes. Note that the sequence of the rotation is important: (AB ~= BA). In this case we rotate about x_2 first, x_3 second and x_1 third.

For trend and plunge of symmetry axis (e.g., tri_tensor):

alpha = plunge

beta = trend

telewavesim.utils.rotate_zrt_pvh(trZ, trR, trT, slow, vp=None, vs=None)

Rotates traces from Z-R-T orientation to P-SV-SH wave mode.

Parameters

• trZ (obspy.trace) – Vertical component

• trR (obspy.trace) – Radial component

• trT (obspy.trace) – Transverse component

• slow (float) – slowness of wave

• vp (float, optional) – P-wave velocity used for rotation

• vs (float, optional) – S-wave velocity used for rotation

Returns

tuple containing:

• trP (obspy.trace): Compressional (P) wave mode

• trV (obspy.trace): Vertically polarized shear (SV) wave mode

• trH (obspy.trace): Horizontally polarized shear (SH) wave mode

Return type

(tuple)

telewavesim.utils.stack_all(st1, st2, pws=False)

Stacks all traces in two Stream objects.

Parameters

• st1 (obspy.stream) – Stream 1

• st2 (obspy.stream,) – Stream 2

• pws (bool, optional) – Enables Phase-Weighted Stacking

Returns

tuple containing:

• stack1 (obspy.trace): Stacked trace for Stream 1

• stack2 (obspy.trace): Stacked trace for Stream 2

Return type

(tuple)

telewavesim.utils.calc_ttime(model, slow, wvtype='P')

Calculates total propagation time through model given the corresponding 'P' or 'S' wave type. The bottom layer is irrelevant in this calculation. All 'S' wave types will return the same predicted time. This function is useful mainly for plotting purposes. For example, to show the first phase arrival at a time of 0, the traces can be shifted by the total propagation time through the model.

Parameters

• model (Model) – Model object

• slow (float) – Slowness value (s/km)

• wvtype (str) – Incident wavetype ('P', 'SV', 'SH', 'Si')

Returns

t1: Time in seconds

Return type

(float)

Example

>>> from telewavesim import utils
>>> # Define two-layer model with identical material
>>> model = utils.Model([10, 0], None, 0, 0, 'atg', 0, 0, 0)
>>> # Only topmost layer is useful for travel time calculation
>>> wvtype = 'P'
>>> slow = 0.06     # s/km
>>> utils.calc_ttime(model, slow, wvtype)
1.3519981570791182

class telewavesim.utils.Model(thickn, rho, vp, vs, isoflg='iso', ani=None, tr=None, pl=None)

model parameters:

• thickn (np.ndarray): Thickness of layers (km) (shape (nlay))

• rho (np.ndarray): Density (kg/m^3) (shape (nlay))

• vp (np.ndarray): P-wave velocity (km/s) (shape (nlay))

• vs (np.ndarray): S-wave velocity (km/s) (shape (nlay))

• isoflg (list of str, optional, defaut: 'iso'):

  Flags for type of layer material (dimension nlay)

• ani (np.ndarray, optional): Anisotropy (percent) (shape (nlay))

• tr (np.ndarray, optional):

  Trend of symmetry axis (degree) (shape (nlay))

• pl (np.ndarray, optional):

  Plunge of symmetry axis (degree) (shape (nlay))

• nlay (int): Number of layers

• a (np.ndarray): Elastic thickness (shape (3, 3, 3, 3, nlay))

update_tensor()

Update the elastic thickness tensor a.

Needs to be called when model parameters change.

telewavesim.utils.read_model(modfile, encoding=None)

Reads model parameters from file and returns a Model object.

Returns

Model object

telewavesim.utils.model2for(model)

Passes global model variables to Fortran conf module.

Returns

None

Variables to pass are a, rho, thickn, isoflg

telewavesim.utils.wave2for(dt, slow, baz)

Passes global wavefield variables to Fortran conf module.

Returns

None

Variables to pass are dt, slow, baz

telewavesim.utils.obs2for(dp, c, rhof)

Passes global OBS-related variables to Fortran conf module.

Returns

None

Variables to pass are dp, c, rhof

telewavesim.utils.run_plane(model, slow, npts, dt, baz=0, wvtype='P', obs=False, dp=100.0, c=1.5, rhof=1027)

Function to run the plane module and return 3-component seismograms as an obspy.Stream object. This function builds the seismic response spectrum by frequency using the matrix propagation approach. Required arguments are the seismic model, slowness value, and the sampling properties (maximum number of samples and the sampling distance in seconds).

By default, the function uses a back-azimuth value of 0 degree, which is suitable for events coming from the North pole or isotropic seismic velocity models (i.e., those that do not vary with direction of incoming waves). For anisotropic velocity models, users need to specify the back-azimuth value in degrees. Furthermore, the default type of the incoming teleseismic body wave is 'P' for compressional wave. Other options are 'SV', 'SH', or 'Si' for vertically-polarized shear wave, horizontally-polarized shear wave or isotropic shear wave, respectively. Wave modes cannot be mixed.

Finally, it is possible to simulate the seismic response for ocean-bottom seismic (OBS) stations using the flag obs=True. If the flag is set to True, the user can specify the water depth below sea level (in meters, positive value) as well as the properties of the sea water (defaults are acoustic wavespeed of 1.5 km/s and density of 1027 kg/m^3).

Parameters

• model (Model) – Instance of the Model class that contains the physical properties of subsurface layers.

• slow (float) – Slowness (s/km)

• baz (float, optional) – Back-azimuth (degree)

• npts (int) – Number of samples in time series

• dt (float) – Sampling distance (s)

• baz – Back-azimuth (degree)

• (str, optional, default (wvtype) – 'P'): Incident wavetype ('P', 'SV', 'SH', 'Si')

• obs (bool, optional) – Whether or not the analysis is done for an OBS stations

• dp (float, optional) – Deployment depth below sea level (m)

• c (float, optional) – P-wave velocity of salt water (default = 1.5 km/s)

• rhof (float, optional) – Density of salt water (default = 1027.0 kg/m^3)

Returns

trxyz: Stream containing 3-component displacement seismograms

Return type

(obspy.stream)

Example

Basic example:

>>> from telewavesim import utils
>>> # Define three-layer model with isotropic crust and antigorite upper mantle layer over isotropic half-space
>>> model = utils.Model([20, 10, 0], [2800., None, 3300.], [4.6, 0, 6.0], [2.6, 0, 3.6], ['iso', 'atg', 'iso'], [0, 0, 0], [0, 0, 0], [0, 0, 0])
>>> slow = 0.06     # s/km
>>> npts = 1500
>>> dt = 0.025      # s
>>> st = utils.run_plane(model, slow, npts, dt)

>>> type(st)
<class 'obspy.core.stream.Stream'>
>>> print(st)
3 Trace(s) in Stream:
...N | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
...E | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
...Z | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
>>> st.plot(size=(600, 450))

OBS station:

>>> from telewavesim import utils
>>> # Define two-layer model with foliated eclogitic crust over isotropic half-space
>>> model = utils.Model([20, 0], [None, 3300.], [0, 6.0], [0, 3.6], ['EC_f', 'iso'], [0, 0], [0, 0], [0, 0])
>>> slow = 0.06     # s/km
>>> npts = 3000
>>> dt = 0.01      # s
>>> wvtype = 'SV'
>>> baz = 45.
>>> dp = 1000.
>>> st = utils.run_plane(model, slow, npts, dt, baz=baz, wvtype=wvtype, obs=True, dp=dp)
>>> st.plot(size=(600, 450))

telewavesim.utils.get_trxyz(yx, yy, yz, npts, dt, slow, baz, wvtype)

Function to store displacement seismograms into obspy.Trace objects and then an obspy Stream object.

Parameters

• ux (np.ndarray) – x-component displacement seismogram

• uy (np.ndarray) – y-component displacement seismogram

• uz (np.ndarray) – z-component displacement seismogram

Returns

trxyz: Stream containing 3-component displacement

seismograms

Return type

(obspy.stream)

telewavesim.utils.tf_from_xyz(trxyz, pvh=False, vp=None, vs=None)

Function to generate transfer functions from displacement traces.

Parameters

• trxyz (obspy.stream) – Obspy Stream object in cartesian coordinate system

• pvh (bool, optional) – Whether to rotate from Z-R-T coordinate system to P-SV-SH wave mode

• vp (float, optional) – Vp velocity at surface for rotation to P-SV-SH system

• vs (float, optional) – Vs velocity at surface for rotation to P-SV-SH system

Returns

tfs: Stream containing Radial and Transverse transfer functions

Return type

(obspy.stream)

telewavesim.utils.update_stats(tr, dt, slow, baz, wvtype, cha)

Updates the stats doctionary from an obspy Trace object.

Parameters

• tr (obspy.trace) – Trace object to update

• nt (int) – Number of samples

• dt (float) – Sampling rate

• slow (float) – Slowness value (s/km)

• baz (float) – Back-azimuth value (degree)

Returns

tr: Trace with updated stats

Return type

(obspy.trace)