API Reference#

LayTracer — fast two-point ray tracing in 1-D layered media.

Public API#

Model#

`LayerStack`(h, vp, vs[, rho, qp, qs])	Layers traversed by a ray between source and receiver depths.
`ModelArrays`(depths, vp, vs[, rho, qp, qs])	Pre-extracted numpy arrays from a velocity DataFrame.
`build_layer_stack`(vel_model, z_src, z_rcv)	Extract the layer stack between source and receiver depths.

Solver#

`solve`(h, v, segments, interactions, ...[, ...])	Solve the two-point ray tracing problem for an arbitrary path.
`RayResult`(travel_time, ray_path, ray_parameter)	Result of a single two-point ray trace.

Multi-ray#

`trace_rays`(sources, receivers, velocity_df)	Trace rays for all source-receiver pairs.
`TraceResult`(travel_times[, rays, ...])	Container for multi-ray tracing results.

Amplitude#

`transmission_normal`(v1, rho1, v2, rho2)	Normal-incidence displacement transmission coefficient.
`psv_rt_coefficients`(p, vp1, vs1, rho1, vp2, ...)	Compute all eight P-SV reflection/transmission coefficients.
`find_brewster_angles`(rt_coefficients, angles)	Find Brewster-like angles (deep minima) in R/T coefficient curves.
`normalize_rt_coefficient`(R_bar, p, v_in, ...)	Apply Červený (2001) normalization to a displacement R/T coefficient.

Visualisation#

plot

Standalone visualisation helpers for LayTracer.

Model#

class laytracer.LayerStack(h, vp, vs, rho=None, qp=None, qs=None)[source]#

Layers traversed by a ray between source and receiver depths.

Parameters:

hnumpy.ndarray: Layer thicknesses (m), shape (N,). Ordered from the shallowest traversed depth to the deepest.
vpnumpy.ndarray: P-wave velocities (m/s), shape (N,).
vsnumpy.ndarray: S-wave velocities (m/s), shape (N,).
rhonumpy.ndarray or None: Densities (kg/m³), shape (N,), or None.
qpnumpy.ndarray or None: P-wave quality factors, shape (N,), or None.
qsnumpy.ndarray or None: S-wave quality factors, shape (N,), or None.

property n_layers#: Number of layers in the stack.

v(vel_type='Vp')[source]#

Return the velocity array for the requested wave type.

Parameters:

vel_typestr: 'Vp' or 'Vs'.

q_factor(vel_type='Vp')[source]#: Return the Q-factor array for the requested wave type.

class laytracer.ModelArrays(depths, vp, vs, rho=None, qp=None, qs=None)[source]#

Pre-extracted numpy arrays from a velocity DataFrame.

Use from_dataframe to construct once and pass to build_layer_stack to avoid repeated DataFrame column access when tracing many rays.

classmethod from_dataframe(vel_df)[source]#: Extract arrays once from vel_df.

laytracer.build_layer_stack(vel_model, z_src, z_rcv)[source]#

Extract the layer stack between source and receiver depths.

The returned LayerStack contains the layers traversed by a ray connecting z_src and z_rcv, with partial thicknesses at the source and receiver layers. Layers are always ordered from the shallowest point to the deepest, regardless of which endpoint is the source.

Parameters:

vel_modelpandas.DataFrame or ModelArrays: Velocity model. A DataFrame must have columns Depth, Vp, Vs (optional: Rho, Qp, Qs). A ModelArrays instance avoids repeated column extraction when tracing many rays.
z_srcfloat: Source depth (positive downward, m).
z_rcvfloat: Receiver depth (positive downward, m).

Returns:

LayerStack

Solver#

laytracer.solve(h, v, segments, interactions, epicentral_dist, z_src, z_rcv, requested=None, return_ray_path=True, need_ray_parameter=True, need_tstar=False, need_spreading=False, need_trans_product=False, transcoef_method='standard', tol=0.0001, max_iter=10)[source]#

Solve the two-point ray tracing problem for an arbitrary path.

Parameters:

hnumpy.ndarray: Concatenated layer thicknesses for the entire path (m).
vnumpy.ndarray: Concatenated phase velocities (Vp or Vs) for the entire path (m/s).
segmentslist of dict: Metadata for each logical monotonic segment (used for reconstruction). Dict keys: ‘h’, ‘v’, ‘phase’, ‘start_z’, ‘end_z’.
interactionslist of dict: Metadata for interactions (reflections/refractions) impacting amplitude. Dict keys: ‘type’, ‘depth’, ‘in_phase’, ‘out_phase’, ‘seg_idx’.
interactions are assumed to occur at the END of the defined segment.

class laytracer.RayResult(travel_time, ray_path, ray_parameter, tstar=None, spreading=None, trans_product=None)[source]#

Result of a single two-point ray trace.

Attributes:

travel_timefloat: Total travel time (s).
ray_pathnumpy.ndarray or None: Ray coordinates in the 2-D ray plane, shape (M, 2) with columns [x, z]. None if not requested.
ray_parameterfloat or None: Horizontal slowness p (s/m).
tstarfloat or None: Attenuation operator \(t^*\) (s), if requested.
spreadingfloat or None: Relative geometrical spreading factor \(\mathcal{L}\), if requested.
trans_productfloat or None: Product of transmission-coefficient magnitudes along the ray, if requested.

laytracer.offset(q, h, lmd)[source]#

Total horizontal range \(X(q)\).

\[X(q) = \sum_{k=1}^{N} \frac{q\,\lambda_k\,h_k} {\sqrt{1 + (1 - \lambda_k^2)\,q^2}}\]

Parameters:

qfloat: Dimensionless ray parameter.
hnumpy.ndarray: Layer thicknesses (m), shape (N,).
lmdnumpy.ndarray: Velocity ratios \(\lambda_k = v_k / v_{\max}\), shape (N,).

Returns:

float

laytracer.offset_dq(q, h, lmd)[source]#: First derivative \(\mathrm{d}X/\mathrm{d}q\).

\[\frac{\mathrm{d}X}{\mathrm{d}q} = \sum_{k=1}^{N} \frac{\lambda_k\,h_k} {\bigl[1 + (1 - \lambda_k^2)\,q^2\bigr]^{3/2}}\]

laytracer.offset_dq2(q, h, lmd)[source]#: Second derivative \(\mathrm{d}^2X/\mathrm{d}q^2\).

\[\frac{\mathrm{d}^2X}{\mathrm{d}q^2} = -3\,q\,\sum_{k=1}^{N} \frac{(1 - \lambda_k^2)\,\lambda_k\,h_k} {\bigl[1 + (1 - \lambda_k^2)\,q^2\bigr]^{5/2}}\]

laytracer.q_from_p(p, vmax)[source]#: Convert horizontal slowness p to dimensionless parameter q.

\[q = \frac{p\,v_{\max}}{\sqrt{1 - p^2\,v_{\max}^2}}\]

laytracer.p_from_q(q, vmax)[source]#: Convert dimensionless parameter q to horizontal slowness p.

\[p = \frac{q}{v_{\max}\,\sqrt{1 + q^2}}\]

laytracer.initial_q(X_target, h, lmd)[source]#

Asymptotic initial estimate \(q_0\) for the Newton iteration.

Two linear asymptotes of \(X(q)\):

Near-field (\(q\to 0\)): \(X \approx m_0\,q\), where \(m_0 = \sum \lambda_k h_k\).
Far-field (\(q\to\infty\)): \(X \approx m_\infty q + b_\infty\), where \(m_\infty = \sum_{k:\,\lambda_k=1} h_k\) and \(b_\infty = \sum_{k:\,\lambda_k<1} \frac{\lambda_k h_k}{\sqrt{1-\lambda_k^2}}\).

The initial estimate is chosen as:

\(q_0 = X / m_0\) when \(X \le X^*\) (near-field),
\(q_0 = (X - b_\infty) / m_\infty\) otherwise.

where \(X^* = m_0\,q^*\) and \(q^* = b_\infty/(m_0-m_\infty)\).

laytracer.newton_step(q_i, X_target, h, lmd)[source]#

One second-order (quadratic) Newton step.

Solves the local quadratic approximation

\[\tfrac{1}{2}\,X''(q_i)\,\Delta q^2 + X'(q_i)\,\Delta q + \bigl[X(q_i) - X_R\bigr] = 0\]

and selects the root that minimises \(|X(q_i+\Delta q)-X_R|\).

Parameters:

q_ifloat: Current iterate.
X_targetfloat: Desired horizontal range.
h, lmdnumpy.ndarray: Layer thicknesses and velocity ratios.

Returns:

q_newfloat: Updated iterate.
X_newfloat: Offset at q_new.

Multi-ray interface#

laytracer.trace_rays(sources, receivers, velocity_df, source_phase='P', reflection=None, refraction=None, requested=('travel_times', 'rays', 'ray_parameters'), transcoef_method='standard', n_jobs=-1, backend='loky', sequential_limit=10000, rays_per_chunk=None, tol=0.0001, max_iter=10, verbose=True)[source]#

Trace rays for all source-receiver pairs.

Every source is paired with every receiver, producing n_src x n_rcv rays (each source traced to all receivers).

Parameters:

sourcesnumpy.ndarray: Source coordinates, shape (n_src, 3) or (3,).
receiversnumpy.ndarray: Receiver coordinates, shape (n_rcv, 3) or (3,).
velocity_dfpandas.DataFrame: Velocity model with columns Depth, Vp, Vs and optionally Rho, Qp, Qs.
source_phasestr: Initial wave phase at source: 'P' or 'S'.
reflectionlist of (depth, phase), optional: Reflection points as (depth, out_phase) tuples.
refractionlist of (depth, phase), optional: Refraction / mode-conversion points as (depth, out_phase) tuples.
requestedsequence of str, optional: Explicit set of requested outputs. Valid names are travel_times, rays, ray_parameters, tstar, spreading, and trans_product. The set must include travel_times.
transcoef_methodstr: 'standard' (Zoeppritz) or 'normalized'.
n_jobsint: Number of parallel jobs (-1 = all physical cores).
backendstr: Joblib parallel backend (default 'loky').
sequential_limitint: If the total number of rays is below this threshold, run sequentially to avoid parallel overhead.
rays_per_chunkint or None: Maximum number of rays to process per memory-bounded chunk.
tolfloat: Newton convergence tolerance (m).
max_iterint: Maximum Newton iterations.
verbosebool: If True, print progress information for chunked processing.

Returns:

TraceResult

class laytracer.TraceResult(travel_times, rays=None, ray_parameters=None, tstar=None, spreading=None, trans_product=None)[source]#

Container for multi-ray tracing results.

Attributes:

travel_timesnumpy.ndarray: Travel times (s), shape (n_rays,).
rayslist of numpy.ndarray or None: Ray paths; each element is shape (M_i, 3) in the original 3-D coordinate system. None if not requested.
ray_parametersnumpy.ndarray or None: Horizontal slowness p for each ray, shape (n_rays,).
tstarnumpy.ndarray or None: Attenuation operator \(t^*\) for each ray, shape (n_rays,).
spreadingnumpy.ndarray or None: Relative geometrical spreading factor for each ray, shape (n_rays,).
trans_productnumpy.ndarray or None: Product of transmission coefficients along each ray.

Amplitude#

laytracer.psv_rt_coefficients(p, vp1, vs1, rho1, vp2, vs2, rho2)[source]#

Compute all eight P-SV reflection/transmission coefficients.

Direct port of the Ammon MATLAB PSVRTmatrix function (Lay & Wallace / Aki & Richards formulation).

For an incident P-wave the system unknowns are \([R_{PP},\; R_{PS},\; T_{PP},\; T_{PS}]\). For an incident SV-wave the unknowns are \([R_{SP},\; R_{SS},\; T_{SP},\; T_{SS}]\).

Parameters:

pfloat or array_like: Ray parameter (horizontal slowness, s/m). Scalar or 1-D array.
vp1, vs1, rho1float: P-velocity, S-velocity, density of the incident medium.
vp2, vs2, rho2float: Same for the transmitted medium.

Returns:

dict: Keys 'Rpp', 'Rps', 'Rss', 'Rsp', 'Tpp', 'Tps', 'Tss', 'Tsp'. Each value has the same shape as p (complex).

laytracer.normalize_rt_coefficient(R_bar, p, v_in, rho_in, v_out, rho_out)[source]#

Apply Červený (2001) normalization to a displacement R/T coefficient.

Eq. 5.3.10:

\[R_{mn} = \bar{R}_{mn} \left( \frac{V(\tilde Q)\,\rho(\tilde Q)\,P(\tilde Q)} {V(Q)\,\rho(Q)\,P(Q)} \right)^{1/2}\]

where \(P(Q) = (1 - V^2 p^2)^{1/2}\).

Parameters:

R_barcomplex or float: Standard (unnormalized) displacement coefficient.
pfloat: Ray parameter (horizontal slowness, s/m).
v_infloat: Phase velocity of the incident wave (m/s).
rho_infloat: Density of the incident medium (kg/m³).
v_outfloat: Phase velocity of the outgoing (reflected/transmitted) wave (m/s).
rho_outfloat: Density of the outgoing medium (kg/m³).

Returns:

complex or float: Normalized displacement coefficient.

laytracer.find_brewster_angles(rt_coefficients, angles, keys=None, threshold=0.05, order=20)[source]#

Find Brewster-like angles (deep minima) in R/T coefficient curves.

A Brewster angle (by analogy with optics) is an incidence angle at which a reflection or transmission coefficient passes through zero or a deep minimum. Unlike critical angles, which depend only on velocity ratios, Brewster angles depend on all six elastic parameters (Vp, Vs, ρ in both media) and arise from destructive interference between displacement potentials at the interface.

Parameters:

rt_coefficientsdict: Output of psv_rt_coefficients — each value is a 1-D array of complex coefficients.
anglesarray_like: Incidence angles (degrees) corresponding to the ray-parameter samples used in rt_coefficients. Must have the same length as the coefficient arrays.
keyslist of str, optional: Which coefficient keys to search (e.g. ['Rps', 'Rss']). By default all eight keys are searched.
thresholdfloat, optional: Only report minima whose absolute value is below this value. Default 0.05.
orderint, optional: Half-window size passed to scipy.signal.argrelmin for local-minimum detection. Default 20.

Returns:

dict[str, list[float]]: Mapping from coefficient key to a list of Brewster angles (degrees). Keys with no detected minima are omitted.

laytracer.transmission_normal(v1, rho1, v2, rho2)[source]#

Normal-incidence displacement transmission coefficient.

\[T = \frac{2\,Z_1}{Z_1 + Z_2}, \qquad Z_i = \rho_i\,v_i\]

Parameters:

v1, rho1float: Velocity (m/s) and density (kg/m³) on the incident side.
v2, rho2float: Velocity and density on the transmitted side.

Returns:

float: Displacement amplitude transmission coefficient.

Visualisation#

laytracer.plot.velocity_profile(vel_df, param='Vp', ax=None, color=None, label=None, xlim=None, ylim=None, unit='m', **kwargs)[source]#

Plot a 1-D model parameter–depth step profile.

Parameters:

vel_dfpandas.DataFrame: Velocity model.
paramstr, optional: 'Vp' (default), 'Vs', 'Rho', 'Qp', or 'Qs'.
axmatplotlib.axes.Axes, optional: Axes to plot on. Created if None.
colorstr, optional: Line colour.
labelstr, optional: Legend label.
xlim, ylimtuple, optional: Axis limits (min, max).
unitstr, optional: 'm' (default) or 'km'. Scales the vertical depth axis.

Returns:

matplotlib.axes.Axes

laytracer.plot.rays_2d(vel_df, rays, vel_type='Vp', sources=None, receivers=None, ax=None, ray_color='k', ray_alpha=0.6, xlim=None, ylim=None, unit='m', plot_model=True, add_colorbar=False, discrete_colorbar=False, model_alpha=1.0, equal_scale=True, colorbar_orientation='vertical', **kwargs)[source]#

Plot ray paths over a 2-D layered velocity cross-section.

Parameters:

vel_dfpandas.DataFrame: Velocity model.
rayslist of numpy.ndarray: Each element is shape (M, 2) or (M, 3). If 3-D, the first two columns are treated as horizontal/depth.
vel_typestr: 'Vp' or 'Vs' — used for layer colouring.
sources, receiversnumpy.ndarray, optional: Coordinate arrays for plotting markers.
axmatplotlib.axes.Axes, optional
ray_colorstr
ray_alphafloat
xlim, ylimtuple, optional
plot_modelbool: If True (default), plot the velocity model background and set axis labels/titles. If False, only plot the rays and markers.
unitstr: 'm' (default) or 'km'. Scales coordinates and labels.
add_colorbarbool: If True (default False), add a colorbar for the velocity model. Only applies if plot_model is True.
discrete_colorbarbool: If True (default False), quantize the colormap to the unique velocity values in the model.
model_alphafloat: Opacity of the velocity model layers (0.0 to 1.0). Default 1.0.
equal_scalebool: If True (default True), force equal scaling for x and y axes using ax.set_aspect('equal').
colorbar_orientationstr: 'vertical' (default) or 'horizontal'.

Returns:

matplotlib.axes.Axes

laytracer.plot.rays_3d(vel_df, rays, vel_type='Vp', sources=None, receivers=None, ray_color='red', opacity=0.3, **kwargs)[source]#

Interactive 3-D ray visualisation using Plotly.

Parameters:

vel_dfpandas.DataFrame: Velocity model.
rayslist of numpy.ndarray: Each element is shape (M, 3).
vel_typestr: 'Vp' or 'Vs'.
sources, receiversnumpy.ndarray, optional: Coordinate arrays for plotting markers.
ray_colorstr: Ray trace colour.
opacityfloat: Layer surface opacity.

Returns:

plotly.graph_objects.Figure