seidart.routines.definitions module¶

seidart.routines.definitions.agc(ts: ndarray, k: int, agctype: str) → ndarray¶

Applies auto-gain control (AGC) normalization to a time series using a specified window length and type. This function normalizes the input time series using a running window approach, with the normalization type determined by agctype. It supports standard deviation (“std”), mean (“mean”), or median (“median”) based normalization. The function modifies the series to have a uniform amplitude characteristic over time.

Parameters:

ts (np.ndarray) – The input time series as a numpy array.
k (int) – The length of the window used for normalization.
agctype (str) – The type of normalization to apply, choices are “std”, “mean”, or “median”.

Returns:

The time series after applying AGC normalization.

Return type:

np.ndarray

seidart.routines.definitions.airsurf(material, domain, N: int = 2) → ndarray¶

Generates a gradient matrix to simulate the air surface in a domain, based on material properties.

This function calculates a gradient matrix representing the transition between air and non-air materials within a domain. It is particularly useful for creating a surface with variable density or other properties near the air-material interface.

Parameters:

material (Material) – The material object containing material lists and properties.
domain (Domain) – The domain object containing the geometry of the simulation area.
N (int) – The number of gradational steps in the air surface simulation, defaults to 2.

Returns:

A gradient matrix representing the air surface within the domain.

Return type:

np.ndarray

The function uses the domain’s geometry to determine air interfaces and applies a gradational approach to simulate the air surface effect. The gradient matrix can be used to adjust physical properties like density near the surface.

seidart.routines.definitions.coherdt(alpha: float, v: float, n: int, dx: float, loc: str = 'left') → ndarray¶

Calculates a vector of time delays based on the angle of incidence, velocity, and spatial discretization.

This function determines the time delay for coherent addition of traces based on the propagation angle (alpha), velocity (v), and the discretization in space (dx) along the domain length.

Parameters:

alpha (float) – The angle of incidence in degrees, counter clockwise from the x-plane. Range: (-90, 90).
v (float) – Velocity in meters/second.
n (int) – The number of spatial points in the domain.
dx (float) – The spatial discretization in meters.
loc (str) – Specifies the reference point’s location for time delay calculation, defaults to ‘left’.

Returns:

A numpy array containing the time delays for each spatial point in the domain.

Return type:

np.ndarray

The function adjusts angles outside the (-90, 90) range and calculates time delays accordingly.

seidart.routines.definitions.coherstf(timediff: ndarray, source_function: ndarray, dt: float, m: int, n: int, cpml: int, bottom: bool = False) → None¶

Applies time shifts to a source time function along the domain and saves it as an m-by-n-by-p matrix.

Parameters:

timediff (np.ndarray) – A numpy array containing time differences for time shifting the source function.
sf (np.ndarray) – The source time function as a numpy array.
dt (float) – The time interval between successive samples in the source function.
m (int) – Number of indices in the y-direction (height of the domain).
n (int) – Number of indices in the x-direction (width of the domain).
cpml (int) – The size of the Convolutional Perfectly Matched Layer (CPML) padding.
bottom (bool) – Specifies whether to place the source at the bottom of the domain. Defaults to False.

Returns:

None

This function time shifts the source time function for each node in the x-direction, considering CPML adjustments.

‘topleft’ origin is (0,0), ‘topright’ origin is (x_n, 0) ‘bottomleft’ origin is (0, y_n), ‘bottomright’ origin is (x_n, y_n)

seidart.routines.definitions.complex2str(complex_vector: ndarray) → str¶

Convert a numpy array of complex numbers into a string representation, with the first element representing an identifier and the subsequent elements representing complex numbers.

This function extracts the real and imaginary parts of each complex number in the input vector, converts them to strings, and concatenates them with ‘j’ as the delimiter. The first real number of the input vector is treated as an identifier (ID) and is converted separately.

Parameters:: complex_vector (np.ndarray) – A numpy array of complex numbers where the first element is an ID.
Returns:: A string representation of the ID followed by the complex numbers in the vector.
Return type:: str

The returned string format is ‘ID,real1jimag1,real2jimag2,…’ where ID is the integer ID, and realN and imagN are the real and imaginary parts of the Nth complex number, respectively.

seidart.routines.definitions.compute_dispersion(seismograms: ndarray[tuple[Any, ...], dtype[floating]], dx: float, dt: float, fmin: float, fmax: float, nfreq: int = 50) → Tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]¶

Estimate phase-velocity dispersion from a 2D seismogram gather.

A linear phase-vs-distance trend is fitted for each selected frequency bin, with \(2\pi\) branch correction applied to unwrap the overall phase slope. The resulting wavenumber is converted to phase velocity, which is then clamped to a physically plausible range.

Parameters:

seismograms (numpy.typing.NDArray[numpy.floating]) – Input seismogram gather with shape (nt, nx), where nt is the number of time samples and nx is the number of spatial traces.
dx (float) – Spatial sampling interval in meters.
dt (float) – Time sampling interval in seconds.
fmin (float) – Minimum frequency in Hz for dispersion estimation.
fmax (float) – Maximum frequency in Hz for dispersion estimation.
nfreq (int) – Number of frequency samples between fmin and fmax to evaluate. Default is 50.

Returns:

A tuple (freqs, vs) where: freqs is a 1D array of frequencies in Hz, and vs is a 1D array of phase velocities in m/s for each frequency, with invalid or out-of-range estimates set to NaN.

Return type:

tuple

Note

For each target frequency, the complex spectrum is extracted across all traces, the phase is unwrapped, and a linear fit of phase versus distance is computed. A \(2\pi\) branch correction ensures that the fitted slope is consistent over the aperture.

The maximum physically resolvable velocity is taken as dx / dt based on the spatial and temporal sampling. Phase-velocity estimates outside (0, max_vel] are rejected and stored as NaN in the output.

# Example: compute a dispersion curve from a synthetic gather
nt, nx = 2048, 48
dt, dx = 0.004, 5.0
seismograms = np.random.randn(nt, nx)

freqs, vs = compute_dispersion(
    seismograms,
    dx=dx,
    dt=dt,
    fmin=2.0,
    fmax=40.0,
    nfreq=60,
)

seidart.routines.definitions.compute_dispersion_image(seismograms: ndarray[tuple[Any, ...], dtype[floating]], dx: float, dt: float, fmin: float, fmax: float, nv: int = 100, nfreq: int = 100, vmin: float | None = None, vmax: float | None = None) → Tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]¶

Compute a MASW-style dispersion image and phase-velocity picks.

A multi-channel analysis of surface waves (MASW) slant-stack is formed over a range of trial phase velocities and frequencies. The resulting f–v image is normalized, smoothed, and used to perform sub-bin parabolic peak picking of the dominant dispersion curve.

Parameters:

seismograms (numpy.typing.NDArray[numpy.floating]) – Input seismogram gather with shape (nt, nx), where nt is the number of time samples and nx is the number of spatial traces.
dx (float) – Spatial sampling interval in meters.
dt (float) – Time sampling interval in seconds.
fmin (float) – Minimum frequency in Hz for the dispersion image.
fmax (float) – Maximum frequency in Hz for the dispersion image.
nv (int) – Number of phase-velocity samples between vmin and vmax. Default is 100.
nfreq (int) – Number of frequency samples between fmin and fmax for the image and picks. Default is 100.
vmin (float or None) – Minimum phase velocity in m/s for the image. If None, it is estimated from initial phase picks.
vmax (float or None) – Maximum phase velocity in m/s for the image. If None, it is estimated from initial phase picks.

Returns:

A tuple (freqs, vs, image, picks) where: freqs is a 1D array of frequencies in Hz, vs is a 1D array of trial phase velocities in m/s, image is a 2D normalized dispersion image with shape (len(freqs), len(vs)), and picks is a 1D array of sub-bin phase-velocity picks (m/s) at each frequency, with out-of-bounds values set to NaN.

Return type:

tuple

Raises:

ValueError – If the automatically or manually determined velocity bounds satisfy vmin >= vmax.

Note

Initial phase-velocity estimates are obtained via compute_dispersion() and used to set automatic bounds when vmin and vmax are not provided. If these picks are unavailable, fallback bounds are based on the sampling ratio dx / dt.

For each frequency, complex amplitudes are slant-stacked over trial velocities using an exponential steering term. The image is normalized per frequency, smoothed with a Gaussian filter, and a parabolic fit around the maximum value in each row is used to obtain sub-bin phase-velocity picks.

# Example: build a dispersion image and picks from a gather
nt, nx = 2048, 48
dt, dx = 0.004, 5.0
seismograms = np.random.randn(nt, nx)

freqs, vs, image, picks = compute_dispersion_image(
    seismograms,
    dx=dx,
    dt=dt,
    fmin=2.0,
    fmax=40.0,
    nv=120,
    nfreq=80,
)

Compute the amplitude envelope of a 1D signal using the Hilbert transform.

The signal is first padded by reflection at both ends to reduce edge effects in the Hilbert transform. The envelope is then computed on the padded signal and the padding is removed before returning.

Parameters:

signal (ArrayLike) – Real-valued input signal.
pad_width (int) – Number of samples to pad on each side using reflection. The input length must be at least 2 * pad_width.

Returns:

Amplitude envelope of the input signal with the same length as the original (unpadded) signal.

Return type:

numpy.typing.NDArray[numpy.floating]

Raises:

ValueError – If the input signal is shorter than 2 * pad_width.

Note

Reflection padding helps mitigate boundary artifacts introduced by the Hilbert transform. If the signal is very short or strongly non-stationary near the edges, you may need to adjust pad_width to obtain a stable envelope.

seidart.routines.definitions.compute_fk_spectrum(data: ndarray[tuple[Any, ...], dtype[floating]], dt: float, dx: float, taper: bool = True, ntfft: int | None = None, nxfft: int | None = None, fmin: float = 0.0, fmax: float | None = None, kmin: float | None = None, kmax: float | None = None, wavenumber_units: str = 'cycles') → Tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]¶

Compute the slant-stack f–k power spectrum of a 2D shot gather.

The input gather is assumed to have time along the first axis and space (offset) along the second axis. A 2D FFT is computed, and the power spectrum is returned as a function of frequency and spatial wavenumber with optional band-limiting in both dimensions.

Parameters:

data (numpy.typing.NDArray[numpy.floating]) – Input shot gather with shape (nt, nx) where nt is the number of time samples and nx is the number of spatial traces.
dt (float) – Time sampling interval in seconds.
dx (float) – Spatial sampling interval in meters.
taper (bool) – If True, apply a separable Hann taper in time and space to reduce spectral leakage. Default is True.
ntfft (int or None) – Number of time samples for the FFT. If None, use nt (no zero-padding).
nxfft (int or None) – Number of spatial samples for the FFT. If None, use nx (no zero-padding).
fmin (float) – Minimum frequency in Hz to include in the output.
fmax (float or None) – Maximum frequency in Hz to include in the output. If None, include up to the Nyquist frequency.
kmin (float or None) – Minimum absolute wavenumber to include in the output, expressed in the units specified by wavenumber_units. If None, defaults to 0.0.
kmax (float or None) – Maximum absolute wavenumber to include in the output, expressed in the units specified by wavenumber_units. If None, defaults to the maximum available wavenumber.
wavenumber_units (str) – Units of the output wavenumber axis: "cycles" for cycles per meter (default) or "radians" for radians per meter.

Returns:

A tuple (freqs, ks, P) where: freqs is a 1D array of frequencies in Hz, ks is a 1D array of spatial wavenumbers in the requested units, and P is a 2D power spectrum array with shape (len(freqs), len(ks)).

Return type:

tuple

Note

The temporal frequency axis is constructed using numpy.fft.rfftfreq and restricted to the range [fmin, fmax]. The spatial wavenumber axis is built from numpy.fft.fftfreq and fftshift, and optionally converted to radians per meter.

Only wavenumbers whose absolute value lies within [kmin, kmax] are retained, and the corresponding columns are selected from the shifted power spectrum.

# Example: compute an f–k spectrum from synthetic data
nt, nx = 1024, 64
dt, dx = 0.004, 10.0  # 4 ms, 10 m
data = np.random.randn(nt, nx)

freqs, ks, P = compute_fk_spectrum(
    data,
    dt=dt,
    dx=dx,
    fmin=2.0,
    fmax=60.0,
    kmin=0.0,
    kmax=0.2,
    wavenumber_units="cycles",
)

Compute phase shift and gain between two time series at a specific frequency.

The function computes the real FFT of both input time series, locates the frequency bin closest to f0, and evaluates the complex transfer function between the spectra. From this transfer function, it returns the phase shift (in degrees) and gain (in decibels).

Parameters:

ts_ref (ArrayLike) – Reference time series (1D array).
ts (ArrayLike) – Comparison time series, with the same shape as ts_ref.
f0 (float) – Center frequency in Hz at which to compute phase and gain.
dt (float) – Sampling interval in seconds.

Returns:

A tuple (phase_deg, gain_db) where: phase_deg is the phase shift of ts relative to ts_ref at f0, in degrees, and gain_db is the gain in decibels, computed as 20 * log10(|H|) with H the complex spectral ratio.

Return type:

tuple[float, float]

Raises:

ValueError – If the input time series do not have the same shape or are empty.

Note

The frequency bin used for the calculation is chosen as the one whose center frequency is closest to f0 in the real FFT frequency grid (as returned by numpy.fft.rfftfreq). For more accurate estimates at arbitrary frequencies, a longer time series or additional spectral interpolation may be required.

If the reference spectrum is (near) zero at the selected bin, the transfer function may be ill-conditioned and the resulting gain or phase can be unstable. It is recommended to check that the reference signal has sufficient energy near f0.

Compute the integrated energy of the horizontal field projection.

The horizontal field is projected onto multiple polarization angles for a 3-component time series (using the x and y components), and the energy is integrated over time for each projection angle.

Parameters:

E_x (ArrayLike) – 1D array of the x-component time series.
E_y (ArrayLike) – 1D array of the y-component time series.
dt (float) – Time-step interval in seconds.
angles (ArrayLike) – 1D array of polarization angles in radians at which to project the horizontal field.
to_db (bool) – If True, convert energies to dB scale using 10 * log10(E / E.max()). Default is False.
corrected (bool) – If True, apply a cosine-squared baseline correction based on the projection geometry and return the corrected energies in dB. Default is False.

Returns:

Array of shape (len(angles),) containing the energy at each projection angle. If to_db is True and corrected is False, the values are in dB relative to the maximum energy. If corrected is True, the returned values are baseline-corrected and expressed in dB.

Return type:

numpy.typing.NDArray[numpy.floating]

Note

For each angle \(\theta\), the horizontal field is projected as \(E_\theta = E_x \cos(\theta) + E_y \sin(\theta)\), and the energy is computed as the discrete integral \(\sum |E_\theta|^2 \Delta t\).

When corrected is True, the raw energy pattern is divided by a cosine-squared baseline (cos(angles)**2) to account for purely geometric projection effects before converting to dB. This can be useful when comparing polarization strength across angles.

seidart.routines.definitions.correct_geometric_spreading(time_series: ndarray, distances: ndarray, exponent: float = 1.0) → ndarray¶

Corrects a seismic time series for geometric spreading.

Parameters:

time_series (np.ndarray) – The seismic time series data. Should be a numpy array of shape (n_samples, n_traces).
distances (np.ndarray) – The distances from the source to each trace. Should be a numpy array of shape (n_traces,).
exponent (float, optional) – The exponent used in the geometric spreading correction. Defaults to 1.0.

Raises:

ValueError – If the length of distances does not match the number of traces.

Returns:

The corrected seismic time series. Will be a numpy array of shape (n_samples, n_traces).

Return type:

np.ndarray

Example usage:

# Example seismic time series (n_samples x n_traces)
time_series = np.random.rand(1000, 10)  # Replace with your actual data

# Example distances (n_traces,)
distances = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])  # Replace with your actual distances

# Apply geometric spreading correction
corrected_time_series = correct_geometric_spreading(time_series, distances)

print(corrected_time_series)

seidart.routines.definitions.cpmlcompute(modelclass, domain, half: bool = False, velocity_scaling_factor=1.0, pml_smoothing: bool = 'global', smoothing_window: int = 7, smoothing_passes: int = 1) → None¶

Computes CPML parameters for a given direction and updates model/domain.

For 2 / 2.5-D domains the output arrays are 2-D (nx+2c, nz+2c) and are written as Fortran binary records (row-major transpose).

For true 3-D domains (domain.dim == 3) separate damping profiles are built along x, y, and z independently and the output arrays are 3-D (nx+2c, ny+2c, nz+2c), written in Fortran-array order.

Parameters:

modelclass (Model) – The model class instance to update.
domain (Domain) – The domain class instance to update.
half (bool) – Flag to compute half-step CPML parameters. Defaults to False.
pml_smoothing – Smoothing strategy for the PML boundary values. Options are None, 'global' (default), or 'smoothed'.

seidart.routines.definitions.exponential_gain(time_series: ndarray, dt: float, alpha: float) → ndarray¶

Applies an exponential gain to a time series to correct for attenuation.

Parameters:

time_series (np.ndarray) – The original time series data. Should be a numpy array of shape (n_samples,).
alpha (float) – The attenuation coefficient.

Returns:

The corrected time series.

Return type:

np.ndarray

seidart.routines.definitions.indvar(model_object, domain) → tuple[ndarray, ndarray | None, ndarray, ndarray]¶

Generates spatial and temporal grids based on domain and model object attributes.

Parameters:

model_object – An object containing model-specific parameters, including time steps and time interval (dt).
domain – An object containing domain-specific parameters such as spatial dimensions (nx, ny, nz) and discretizations (dx, dy, dz).

Returns:

A tuple of numpy arrays representing the grids in x, y (None if ny is not set), z, and t dimensions.

Return type:

tuple

The function calculates the spatial grid points in the x, z, and optionally y dimensions, and temporal grid points based on the provided domain and model parameters. The y-dimension grid is generated only if the domain parameters include ‘ny’.

seidart.routines.definitions.loadproject(project_file: str, domain, material, seismic, electromag)¶

Loads project settings from a project file and assigns values to domain, material, seismic, and electromagnetic objects.

This function parses a specified project file, reading configurations and parameters for the domain, materials, seismic modeling, and electromagnetic properties. These parameters are then assigned to the provided objects.

Parameters:

project_file (str) – The path to the project configuration file.
domain – An object to hold domain-specific parameters and configurations.
material – An object to hold material-specific parameters and configurations.
seismic – An object to hold seismic modeling parameters and configurations.
electromag – An object to hold electromagnetic modeling parameters and configurations.

Returns:

A tuple containing the updated domain, material, seismic, and electromagnetic objects.

Return type:

tuple

The function assumes the project file contains specific prefixes for different types of configurations (e.g., ‘I’ for images, ‘D’ for domain configurations, ‘S’ for seismic modeling parameters, etc.).

seidart.routines.definitions.movingsrc(sf: ndarray, txlocs: ndarray) → None¶

Simulates a moving source given a stationary source time function.

Parameters:

sf (np.ndarray) – The stationary source time function as a numpy array.
txlocs (np.ndarray) – A numpy array containing the transmission locations over time.

Returns:

None

This function is intended as a placeholder to simulate a moving source by manipulating the stationary source time function according to the provided transmission locations.

seidart.routines.definitions.plot_3c(time, data1, data2=None, data3=None, data_label1='Data1', data_label2='Data2', data_label3='Data3', source_wavelet=None, arrivals1=None, arrivals2=None, arrivals3=None, component_labels=None, arrival_labels1=None, arrival_labels2=None, arrival_labels3=None, arrival_colors1=None, arrival_colors2=None, arrival_colors3=None, color_data1='black', color_data2='blue', color_data3='green', envelope_color1='lightgray', envelope_color2='lightgray', envelope_color3='lightgray', plot_envelope=True, figsize=(10, 8), show_legend=True, xlims=None, yaxis_rotation=90, envelope_pad_width=50)¶

Plot 3-component seismograms with vertical dashed lines for expected arrival times.

For each component channel, the function will:

Optionally cross-correlate the channel with a provided source wavelet.
Compute a power envelope (via Hilbert transform with padding to reduce edge effects).
Plot the envelope as a shaded region and overlay the cross-correlated seismogram.

The function supports up to three datasets:

data1 is required.
data2 and data3 are optional.

Separate arrival times, labels, and line colors can be specified for each dataset.

Parameters:

time (numpy.ndarray) – 1D array of time values.
data1 (numpy.ndarray) – 2D array (n_samples x 3) for the first set of seismogram components.
data2 (numpy.ndarray, optional) – 2D array (n_samples x 3) for the second set of seismogram components.
data3 (numpy.ndarray, optional) – 2D array (n_samples x 3) for the third set of seismogram components.
source_wavelet (numpy.ndarray, optional) – 1D array containing the source wavelet. If provided, cross correlation is performed.
arrivals1 (dict or list/tuple, optional) – Arrival times for data1. If a list/tuple, must have 3 elements in order [P, S1, S2].
arrivals2 (dict or list/tuple, optional) – Arrival times for data2. If a list/tuple, must have 3 elements in order [P, S1, S2].
arrivals3 (dict or list/tuple, optional) – Arrival times for data3. If a list/tuple, must have 3 elements in order [P, S1, S2].
component_labels (list of str, optional) – Labels for the 3 components. Default is [‘Vertical’, ‘Radial’, ‘Transverse’].
arrival_labels1 (dict, optional) – Arrival labels for data1. Default is {‘P’: ‘P1’, ‘S1’: ‘S1_1’, ‘S2’: ‘S2_1’}.
arrival_labels2 (dict, optional) – Arrival labels for data2. Default is {‘P’: ‘P2’, ‘S1’: ‘S1_2’, ‘S2’: ‘S2_2’}.
arrival_labels3 (dict, optional) – Arrival labels for data3. Default is {‘P’: ‘P3’, ‘S1’: ‘S1_3’, ‘S2’: ‘S2_3’}.
arrival_colors1 (dict, optional) – Colors for vertical dashed lines for data1 arrivals. Default is {‘P’: ‘red’, ‘S1’: ‘red’, ‘S2’: ‘red’}.
arrival_colors2 (dict, optional) – Colors for vertical dashed lines for data2 arrivals. Default is {‘P’: ‘blue’, ‘S1’: ‘blue’, ‘S2’: ‘blue’}.
arrival_colors3 (dict, optional) – Colors for vertical dashed lines for data3 arrivals. Default is {‘P’: ‘green’, ‘S1’: ‘green’, ‘S2’: ‘green’}.
color_data1 (str, optional) – Color for the data1 trace (default ‘black’).
color_data2 (str, optional) – Color for the data2 trace (default ‘blue’).
color_data3 (str, optional) – Color for the data3 trace (default ‘green’).
envelope_color1 (str, optional) – Fill color for the envelope of data1 (default ‘lightgray’).
envelope_color2 (str, optional) – Fill color for the envelope of data2 (default ‘lightgray’).
envelope_color3 (str, optional) – Fill color for the envelope of data3 (default ‘lightgray’).
figsize (tuple, optional) – Figure size.
envelope_pad_width (int, optional) – Padding width (in samples) for the Hilbert transform envelope computation.

Returns:

fig – The figure handle.

Return type:

matplotlib.figure.Figure

seidart.routines.definitions.plot_dispersion(freqs: ndarray[tuple[Any, ...], dtype[floating]], velocities: ndarray[tuple[Any, ...], dtype[floating]]) → None¶

Plot a dispersion curve of frequency versus phase velocity.

A simple line plot is created showing phase velocity as a function of frequency, suitable for visualizing dispersion picks or curves.

Parameters:

freqs (numpy.typing.NDArray[numpy.floating]) – 1D array of frequencies in Hz.
velocities (numpy.typing.NDArray[numpy.floating]) – 1D array of phase velocities in m/s, with the same shape as freqs.

Returns:

This function does not return anything; it creates a Matplotlib figure and displays it.

Return type:

None

Note

This function calls matplotlib.pyplot.show(), which will block execution in some environments. If you need more control over the figure (e.g., saving to file or embedding in a GUI), consider modifying the code to return the figure and axes instead of calling plt.show() directly.

# Example: plot dispersion picks
freqs = np.linspace(2.0, 40.0, 50)
velocities = 300.0 + 5.0 * freqs  # synthetic trend

plot_dispersion(freqs, velocities)

seidart.routines.definitions.plot_dispersion_image(freqs: ndarray[tuple[Any, ...], dtype[floating]], vs: ndarray[tuple[Any, ...], dtype[floating]], image: ndarray[tuple[Any, ...], dtype[floating]], picks_freqs: ndarray[tuple[Any, ...], dtype[floating]] | None = None, picks_vs: ndarray[tuple[Any, ...], dtype[floating]] | None = None, colormap: str = 'inferno') → Tuple[matplotlib.figure.Figure, matplotlib.axes.Axes]¶

Plot a MASW-style dispersion image with an optional picked curve.

The function displays a frequency–velocity (f–v) image of normalized energy and, if provided, overlays a picked dispersion curve on top of the image.

Parameters:

freqs (numpy.typing.NDArray[numpy.floating]) – 1D array of frequencies in Hz with shape (Nf,).
vs (numpy.typing.NDArray[numpy.floating]) – 1D array of phase velocities in m/s with shape (Nv,).
image (numpy.typing.NDArray[numpy.floating]) – 2D normalized energy map with shape (Nf, Nv), typically the output of compute_dispersion_image().
picks_freqs (numpy.typing.NDArray[numpy.floating] or None) – Optional 1D array of frequencies in Hz for a picked dispersion curve.
picks_vs (numpy.typing.NDArray[numpy.floating] or None) – Optional 1D array of phase velocities in m/s for a picked dispersion curve, matching picks_freqs.
colormap (str) – Matplotlib colormap name used for the image (default "inferno").

Returns:

The figure and axes containing the dispersion image plot.

Return type:

tuple[matplotlib.figure.Figure, matplotlib.axes.Axes]

Note

The energy map is plotted using matplotlib.axes.Axes.pcolormesh() with transposed data so that frequency appears on the horizontal axis and phase velocity on the vertical axis.

When both picks_freqs and picks_vs are provided, the picked curve is overlaid as red scatter points and a legend is added in the upper-right corner.

# Example: visualize a dispersion image and picked curve
freqs, vs, image, picks = compute_dispersion_image(
    seismograms,
    dx=5.0,
    dt=0.004,
    fmin=2.0,
    fmax=40.0,
)

fig, ax = plot_dispersion_image(
    freqs,
    vs,
    image,
    picks_freqs=freqs,
    picks_vs=picks,
    colormap="inferno",
)

plt.show()

seidart.routines.definitions.plot_fk_spectrum(freqs: ndarray[tuple[Any, ...], dtype[floating]], ks: ndarray[tuple[Any, ...], dtype[floating]], P: ndarray[tuple[Any, ...], dtype[floating]], fig: matplotlib.figure.Figure | None = None, ax: matplotlib.axes.Axes | None = None, vmin: float | None = None, vmax: float | None = None, log_scale: bool = False, smooth: bool = False, smooth_sigma: Tuple[float, float] = (1.0, 1.0), mode_lines: Dict[str, float | Tuple[float, str]] | None = None, wavenumber_units: str = 'cycles') → Tuple[matplotlib.figure.Figure, matplotlib.axes.Axes]¶

Display an f–k power spectrum.

The f–k power spectrum is shown as a 2D color mesh of power versus frequency and wavenumber, with options for logarithmic scaling, Gaussian smoothing, and overlaying theoretical mode-velocity lines in either cycles/m or rad/m.

Parameters:

freqs (numpy.typing.NDArray[numpy.floating]) – 1D array of frequencies in Hz, typically the output from compute_fk_spectrum().
ks (numpy.typing.NDArray[numpy.floating]) – 1D array of spatial wavenumbers in either cycles/m or rad/m, matching the wavenumber_units setting.
P (numpy.typing.NDArray[numpy.floating]) – 2D power spectrum array with shape (len(freqs), len(ks)).
fig (matplotlib.figure.Figure or None) – Existing Matplotlib matplotlib.figure.Figure to draw into. If None, a new figure is created.
ax (matplotlib.axes.Axes or None) – Existing Matplotlib matplotlib.axes.Axes to draw into. If None, a new axes is created.
vmin (float or None) – Minimum value for the color scale. If None, it is determined automatically by Matplotlib.
vmax (float or None) – Maximum value for the color scale. If None, it is determined automatically by Matplotlib.
log_scale (bool) – If True, display power in decibels using 10 * log10(P). Default is False.
smooth (bool) – If True, apply a Gaussian filter to P before plotting. Default is False.
smooth_sigma (tuple[float, float]) – Standard deviations of the Gaussian smoothing kernel in the (frequency, wavenumber) directions.
mode_lines (dict[str, float or tuple[float, str]] or None) – Optional dictionary defining theoretical mode lines to overlay. Keys are labels; values are either velocities (float) or (velocity, color) tuples. Velocities are in m/s.
wavenumber_units (str) – Units of ks and the horizontal axis: "cycles" for cycles/m (default) or "radians" for rad/m.

Returns:

The figure and axes used for the plot.

Return type:

tuple[matplotlib.figure.Figure, matplotlib.axes.Axes]

Note

When log_scale is True, a small epsilon is added to the power spectrum before converting to dB to avoid taking the logarithm of zero. Smoothing can help visualize coherent energy trends but may blur fine spectral features.

Mode-velocity lines are drawn as pairs of symmetric curves \(\pm k(f)\) for each velocity, computed using \(k = 2 \pi f / v\) for wavenumbers in radians/m or \(k = f / v\) for cycles/m, depending on wavenumber_units.

# Example: plot an f–k spectrum with a fundamental mode line
freqs, ks, P = compute_fk_spectrum(data, dt=0.004, dx=10.0)

fig, ax = plot_fk_spectrum(
    freqs,
    ks,
    P,
    log_scale=True,
    smooth=True,
    smooth_sigma=(1.0, 2.0),
    mode_lines={"fundamental": 300.0},  # velocity in m/s
)

plt.show()

seidart.routines.definitions.plot_hodogram(data, dt, window=None, components=('L', 'Q', 'T'), title='Hodogram')¶

Plots hodograms (particle motion plots) from given seismic data and returns the figure and axes objects.

Parameters:

data (ndarray (N x 3)) – Seismic data array with columns corresponding to [L, Q, T] components.
dt (float) – Time interval between samples (seconds).
window (tuple (start, end), optional) – Time window (in seconds) to select data. If None, uses full data.
components (tuple of str) – Labels for components, default (‘L’, ‘Q’, ‘T’).
title (str) – Plot title.

Returns:

fig (matplotlib.figure.Figure) – The figure object.
axs (ndarray of matplotlib.axes.Axes) – Array containing the axes objects.

seidart.routines.definitions.polarization_analysis(data: ndarray[tuple[Any, ...], dtype[floating]], dt: float, M: int, alpha: float = 0.1) → tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[bool]]]¶

Compute rectilinearity, backazimuth, and incidence angle using a sliding Tukey window.

The analysis uses a sliding window over 3-component motion data to estimate polarization attributes via covariance eigen-decomposition. Near the signal edges, shorter windows are used, and a cone-of-influence mask is returned to indicate where the full window length is available.

Parameters:

data (numpy.typing.NDArray[numpy.floating]) – 3-component motion data with shape (N, 3) (e.g., LQT, ZRT, or XYZ).
dt (float) – Time step in seconds.
M (int) – Window length in samples.
alpha (float) – Tukey window shape parameter in [0, 1]. Default is 0.1.

Returns:

A tuple (times, rectilinearity, backazimuth, incidence, cone_mask) where: times is the time vector in seconds, rectilinearity contains rectilinearity values, backazimuth contains backazimuth angles in degrees from the y-axis, incidence contains incidence angles in degrees from the vertical, cone_mask is a boolean array that is True where a full-length window is used.

Return type:

tuple

Note

The polarization attributes are computed from the eigenvectors and eigenvalues of the windowed covariance matrix. The principal eigenvector is forced to point in a consistent (upward) direction by flipping it when its vertical component is negative.

Near the series boundaries, the effective window is shorter than M samples. In these regions, the estimates may be less stable; the cone_mask output can be used to restrict interpretation to times where the full window length is available.

seidart.routines.definitions.rcxgen(rgb: list, domain, material, filename: str = 'receivers.xyz')¶

Creates a receiver list from a specified RGB value found in the model image, outputting to a CSV file.

This function searches the domain’s geometry for occurrences of a given RGB value (specified as a list of integers representing red, green, and blue components). It then generates a list of receiver coordinates based on the locations where the RGB value is found. Currently, the function is set up for 2D receivers only, with Y-coordinates set to zero.

Parameters:

rgb (list) – A list of three integers representing the RGB value to search for in the domain’s geometry.
domain – An object representing the domain, including its geometry and spatial discretization parameters (dz, dx).
material – An object containing the material list and associated RGB values.
filename (str) – The name of the output CSV file containing the receiver coordinates, defaults to ‘receivers.xyz’.

Returns:

A numpy array of receiver coordinates where each row represents [Z, Y, X] coordinates.

Return type:

np.ndarray

The function converts the RGB list into a string format, looks up the: corresponding integer value from the
material list, finds all occurrences in the domain’s geometry, and outputs: the coordinates to a specified CSV file.

seidart.routines.definitions.read_dat(fn: str, channel: str, domain, single: bool = False) → ndarray¶

Read data from an unformatted Fortran binary file and return it as a numpy array. This function supports reading data in single or double precision, and can handle complex data.

The domain parameter is expected to be an object with dim, nx, ny, and nz attributes. These values determine the shape of the returned array.

Parameters:

fn (str) – The filename of the data file to read.
channel (str) – Specifies the data channel (‘Ex’, ‘Ey’, ‘Ez’, etc.) to read.
domain (Domain) – The domain object that contains dimensions (dim, nx, ny, nz) of the data.
single (bool) – A flag indicating if the data should be read in single precision. Defaults to False (double precision).

Returns:

The data read from the file, reshaped according to the domain dimensions.

Return type:

np.ndarray

seidart.routines.definitions.rotate_sdr(strike: float, dip: float, rake: float, is_degree: bool = True) → ndarray[tuple[Any, ...], dtype[floating]]¶

Create a 3D rotation matrix from strike, dip, and rake angles.

The rotation matrix is constructed from three sequential rotations about the global axes corresponding to strike, dip, and rake angles. This is commonly used to transform between fault-plane and global coordinate systems in seismological applications.

Parameters:

strike (float) – Strike angle. Interpreted as degrees if is_degree is True, otherwise radians.
dip (float) – Dip angle. Interpreted as degrees if is_degree is True, otherwise radians.
rake (float) – Rake (slip) angle. Interpreted as degrees if is_degree is True, otherwise radians.
is_degree (bool) – If True, input angles are treated as degrees and converted to radians internally. If False, angles are assumed to be in radians.

Returns:

Rotation matrix of shape (3, 3) formed as Rz(strike) @ Ry(dip) @ Rx(rake).

Return type:

numpy.typing.NDArray[numpy.floating]

Note

The rotation is applied as a product of three basic rotations: about the z-axis by strike, about the y-axis by dip, and about the x-axis by rake. The precise physical interpretation depends on your coordinate and fault-plane conventions.

When using this matrix to rotate vectors v in the original coordinate system, the transformed vectors can be obtained as v_rot = R @ v.

# Example: build a rotation matrix for a given focal mechanism
strike, dip, rake = 210.0, 30.0, 90.0  # degrees
R = rotate_sdr(strike, dip, rake)

v = np.array([1.0, 0.0, 0.0])
v_rot = R @ v  # rotated vector

Rotate 3-component data from Cartesian (x, y, z) to (L, Q, T) components.

The input 3-component data are rotated into the longitudinal (L), quasi-vertical (Q), and transverse (T) directions using geometry defined by the source and receiver locations and, optionally, given backazimuth and incidence angles.

Parameters:

data (numpy.typing.NDArray[numpy.floating]) – Input data with shape (N, 3) and columns [vx, vy, vz].
source_location (array_like) – Source coordinates [x, y, z].
receiver_location (array_like) – Receiver coordinates [x, y, z].
backazimuth (float or None) – Backazimuth angle in radians (clockwise from the positive y-axis). If None, it is computed from the source–receiver geometry.
incidence (float or None) – Incidence (inclination) angle in radians, positive downward from horizontal. If None, it is computed from the source–receiver geometry.

Returns:

A tuple (ba, inc, lqt) where: ba is the backazimuth angle in radians, inc is the inclination angle in radians, lqt is the rotated data with shape (N, 3) and columns [L, Q, T].

Return type:

tuple[float, float, numpy.typing.NDArray[numpy.floating]]

Note

When backazimuth or incidence is not provided, both are derived from the displacement vector receiver_location - source_location. The inclination is taken as the arctangent of the vertical offset over the horizontal distance, and the backazimuth is measured clockwise from the positive y-axis.

The rotation matrix is constructed so that the L direction is aligned with the propagation direction, Q lies in the vertical plane containing the propagation direction, and T is the horizontal transverse (SH) direction.

# Example: rotate from vx, vy, vz to L, Q, T using geometry
data = np.random.randn(1000, 3)  # [vx, vy, vz]
src = np.array([0.0, 0.0, 0.0])
rcv = np.array([10.0, 5.0, -2.0])

ba, inc, lqt = rotate_to_qlt(data, src, rcv)

Rotate seismogram data from Cartesian components to Z, R, T components.

The input 3-component seismogram data (vx, vy, vz) are rotated into a local coordinate system defined by the vertical, radial, and transverse directions. The radial direction is determined either from a propagation direction vector or from the source–receiver geometry.

Parameters:

data (numpy.typing.NDArray[numpy.floating]) – Seismogram data in Cartesian components, with shape (n_samples, 3) and components ordered as (vx, vy, vz).
source (array_like or None) – Source coordinates as a 3-element array [x, y, z]. Used together with receiver when direction is not provided.
receiver (array_like or None) – Receiver coordinates as a 3-element array [x, y, z]. Used together with source when direction is not provided.
direction (array_like or None) – Propagation direction vector (3 elements). If provided, its horizontal projection defines the radial direction.

Returns:

Rotated seismogram data with shape (n_samples, 3) and components ordered as (Z, R, T).

Return type:

numpy.typing.NDArray[numpy.floating]

Raises:

ValueError – If neither direction nor a valid source and receiver pair is provided, if the direction vector has zero length, or if the source–receiver horizontal separation is effectively zero.

Note

If direction is provided, it is normalized and its horizontal \((x, y)\) components are used to define the radial direction.

If direction is not provided, both source and receiver must be given. The horizontal difference receiver - source is used to determine the azimuth and radial direction.

The vertical unit vector is assumed to be \([0, 0, 1]\), and the transverse direction is computed as the cross product of the vertical and radial directions, ensuring an approximately right-handed coordinate system.

# Example: rotate from vx, vy, vz to Z, R, T using source/receiver
data = np.random.randn(1000, 3)  # (vx, vy, vz)
source = [0.0, 0.0, 0.0]
receiver = [10.0, 5.0, 0.0]

zrt = rotate_to_zrt(data, source=source, receiver=receiver)

seidart.routines.definitions.sponge_boundary(domain, eta_max: float, eta: ndarray[tuple[Any, ...], dtype[floating]]) → ndarray[tuple[Any, ...], dtype[floating]]¶

Apply a polynomial sponge boundary to a 2D damping field.

A polynomial sponge profile is added near the domain boundaries to gradually damp outgoing waves. The strength and thickness of the sponge are controlled by attributes of domain.

Parameters:

domain (Any) – Simulation domain object providing sponge parameters. Must define cpml (sponge thickness in grid points) and NPS (polynomial order for the sponge profile).
eta_max (float) – Maximum sponge damping coefficient applied at the boundaries.
eta (numpy.typing.NDArray[numpy.floating]) – Existing 2D damping field to be modified in place by adding the sponge boundary contribution.

Returns:

The updated damping field with sponge boundary applied.

Return type:

numpy.typing.NDArray[numpy.floating]

seidart.routines.definitions.stfvec2mat(source_function: ndarray, xind: int, zind: int, m: int, n: int, cpml: int, yind: int = None) → None¶

Converts a source time function vector to a matrix for 2D or 3D simulations.

Parameters:

sf (np.ndarray) – The source time function vector.
xind (int) – The index in the x-direction where the source is located.
zind (int) – The index in the z-direction where the source is located.
m (int) – The number of indices in the y-direction (height of the domain).
n (int) – The number of indices in the x-direction (width of the domain).
cpml (int) – The size of the Convolutional Perfectly Matched Layer (CPML) padding.
yind (Optional[int]) – The index in the y-direction where the source is located, for 3D simulations. Defaults to None for 2D simulations.

Returns:

None

This function arranges the source time function in the spatial domain and applies CPML padding.

seidart.routines.definitions.str2complex(strsplit: list) → ndarray¶

Convert a list of strings into a numpy array of complex numbers, where the first string is treated as an identifier (ID) and the subsequent strings represent complex numbers.

The function splits each string by ‘j’ to separate the real and imaginary parts of each complex number, then constructs a numpy array of complex numbers from these parts. The first element of the array is the ID.

Parameters:: strsplit (list) – A list of strings, where the first string is an ID and the remaining strings represent complex numbers in ‘realjimag’ format.
Returns:: A numpy array of complex numbers with the ID as the first element.
Return type:: np.ndarray

Note that the input list should follow the format [‘ID’, ‘real1jimag1’, ‘real2jimag2’, …], where ‘ID’ is an integer identifier, and ‘realNjimagN’ are the real and imaginary parts of the Nth complex number.