# The Clear BSD License
#
# Copyright (c) 2026 Tobias Heibges
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted (subject to the limitations in the disclaimer
# below) provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""
Time Spread Estimation Utilities
Provides geometric estimation of time spread for diverging beams reflecting
off surfaces and reaching a detector. Useful for computing upper/lower bounds
on arrival time distributions without full ray tracing.
Examples
--------
>>> from lsurf.utilities.time_spread import estimate_time_spread
>>> from lsurf.surfaces import GerstnerWaveSurface, GerstnerWaveParams
>>>
>>> # Create a wave surface
>>> wave = GerstnerWaveParams(amplitude=1.0, wavelength=50.0)
>>> surface = GerstnerWaveSurface(wave_params=[wave])
>>>
>>> result = estimate_time_spread(
... source_position=(0, 0, 500),
... beam_direction=(0.98, 0, -0.17),
... divergence_angle=np.radians(1.0),
... detector_position=(32000, 0, 5700),
... surface=surface,
... )
>>> print(f"Time spread: {result.time_spread_ns:.2f} ns")
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import TYPE_CHECKING
import numpy as np
from numpy.typing import NDArray
if TYPE_CHECKING:
from ..surfaces import Surface
# Speed of light in m/s
SPEED_OF_LIGHT = 299792458.0
[docs]
@dataclass
class TimeSpreadResult:
"""
Results from time spread estimation.
Attributes
----------
min_path : float
Shortest path length in meters
max_path : float
Longest path length in meters
path_spread : float
Difference between max and min path in meters
time_spread_s : float
Time spread in seconds
time_spread_ns : float
Time spread in nanoseconds
min_path_point : ndarray
Surface point with shortest path
max_path_point : ndarray
Surface point with longest path
edge_points : ndarray
All beam edge intersection points on surface
path_lengths : ndarray
Path lengths for all edge points
center_point : ndarray
Center ray intersection point
"""
min_path: float
max_path: float
path_spread: float
time_spread_s: float
time_spread_ns: float
min_path_point: NDArray[np.float64]
max_path_point: NDArray[np.float64]
edge_points: NDArray[np.float64]
path_lengths: NDArray[np.float64]
center_point: NDArray[np.float64]
[docs]
def estimate_time_spread(
source_position: tuple[float, float, float],
beam_direction: tuple[float, float, float],
divergence_angle: float,
detector_position: tuple[float, float, float],
surface: Surface | None = None,
n_edge_points: int = 100,
speed_of_light: float = SPEED_OF_LIGHT,
) -> TimeSpreadResult:
"""
Estimate time spread for a diverging beam reflecting to a detector.
Computes geometric path lengths from source through surface edge points
to detector, giving an upper bound estimate on arrival time spread.
Parameters
----------
source_position : tuple
(x, y, z) source position in meters
beam_direction : tuple
Beam center direction (will be normalized)
divergence_angle : float
Beam half-angle divergence in radians
detector_position : tuple
(x, y, z) detector position in meters
surface : Surface, optional
Surface object to intersect with (e.g., GerstnerWaveSurface,
CurvedWaveSurface, PlanarSurface). If None, uses flat surface at z=0.
n_edge_points : int
Number of points around beam edge for sampling
speed_of_light : float
Speed of light in m/s
Returns
-------
TimeSpreadResult
Dataclass containing path lengths, time spread, and geometry info
Examples
--------
>>> # Flat surface estimate
>>> result = estimate_time_spread(
... source_position=(-2800, 0, 500),
... beam_direction=(0.98, 0, -0.17),
... divergence_angle=np.radians(1.0),
... detector_position=(32000, 0, 5700),
... )
>>> print(f"Flat surface: {result.time_spread_ns:.2f} ns")
>>> # Wavy surface estimate with GerstnerWaveSurface
>>> from lsurf.surfaces import GerstnerWaveSurface, GerstnerWaveParams
>>> wave = GerstnerWaveParams(amplitude=1.0, wavelength=50.0, direction=(0, 1))
>>> surface = GerstnerWaveSurface(wave_params=[wave])
>>> result = estimate_time_spread(
... source_position=(-2800, 0, 500),
... beam_direction=(0.98, 0, -0.17),
... divergence_angle=np.radians(1.0),
... detector_position=(32000, 0, 5700),
... surface=surface,
... )
>>> print(f"Wavy surface: {result.time_spread_ns:.2f} ns")
"""
src = np.array(source_position, dtype=np.float64)
det = np.array(detector_position, dtype=np.float64)
# Get beam footprint on surface
footprint = compute_beam_footprint(
source_position,
beam_direction,
divergence_angle,
n_edge_points=n_edge_points,
surface=surface,
)
edge_points = footprint["edge_points"]
center_point = footprint["center_point"]
if len(edge_points) == 0:
# No intersections found
return TimeSpreadResult(
min_path=0.0,
max_path=0.0,
path_spread=0.0,
time_spread_s=0.0,
time_spread_ns=0.0,
min_path_point=center_point,
max_path_point=center_point,
edge_points=edge_points,
path_lengths=np.array([]),
center_point=center_point,
)
# Compute total path lengths: source -> surface -> detector
path_lengths = []
for pt in edge_points:
d_src_to_surface = np.linalg.norm(pt - src)
d_surface_to_det = np.linalg.norm(det - pt)
total_path = d_src_to_surface + d_surface_to_det
path_lengths.append(total_path)
path_lengths = np.array(path_lengths)
min_path = np.min(path_lengths)
max_path = np.max(path_lengths)
path_spread = max_path - min_path
time_spread_s = path_spread / speed_of_light
time_spread_ns = time_spread_s * 1e9
min_idx = np.argmin(path_lengths)
max_idx = np.argmax(path_lengths)
return TimeSpreadResult(
min_path=min_path,
max_path=max_path,
path_spread=path_spread,
time_spread_s=time_spread_s,
time_spread_ns=time_spread_ns,
min_path_point=edge_points[min_idx],
max_path_point=edge_points[max_idx],
edge_points=edge_points,
path_lengths=path_lengths,
center_point=center_point,
)
[docs]
def estimate_time_spread_bounds(
source_position: tuple[float, float, float],
beam_direction: tuple[float, float, float],
divergence_angle: float,
detector_positions: NDArray[np.float64],
surface: Surface | None = None,
n_edge_points: int = 100,
) -> dict[str, NDArray[np.float64]]:
"""
Estimate time spread bounds for multiple detector positions.
Useful for computing time spread maps across a detector surface.
Parameters
----------
source_position : tuple
(x, y, z) source position in meters
beam_direction : tuple
Beam center direction
divergence_angle : float
Beam half-angle divergence in radians
detector_positions : ndarray, shape (N, 3)
Array of detector positions to evaluate
surface : Surface, optional
Surface object to intersect with. If None, uses flat surface at z=0.
n_edge_points : int
Number of edge points for beam footprint
Returns
-------
dict
Dictionary containing:
- time_spread_ns : ndarray, shape (N,) - Time spread at each position
- path_spread : ndarray, shape (N,) - Path spread at each position
- min_path : ndarray, shape (N,) - Minimum path at each position
- max_path : ndarray, shape (N,) - Maximum path at each position
Examples
--------
>>> # Compute time spread for grid of detector positions
>>> lon = np.linspace(-5, 5, 20)
>>> lat = np.linspace(5, 15, 20)
>>> lon_grid, lat_grid = np.meshgrid(lon, lat)
>>> r = 33000 # 33 km
>>> x = r * np.cos(np.radians(lat_grid)) * np.cos(np.radians(lon_grid))
>>> y = r * np.cos(np.radians(lat_grid)) * np.sin(np.radians(lon_grid))
>>> z = r * np.sin(np.radians(lat_grid))
>>> positions = np.stack([x.ravel(), y.ravel(), z.ravel()], axis=1)
>>>
>>> bounds = estimate_time_spread_bounds(
... source_position=(-2800, 0, 500),
... beam_direction=(0.98, 0, -0.17),
... divergence_angle=np.radians(1.0),
... detector_positions=positions,
... )
>>> time_spread_map = bounds['time_spread_ns'].reshape(lat_grid.shape)
"""
n_detectors = len(detector_positions)
time_spreads = np.zeros(n_detectors)
path_spreads = np.zeros(n_detectors)
min_paths = np.zeros(n_detectors)
max_paths = np.zeros(n_detectors)
# Compute footprint once (same for all detectors)
footprint = compute_beam_footprint(
source_position,
beam_direction,
divergence_angle,
n_edge_points=n_edge_points,
surface=surface,
)
src = np.array(source_position, dtype=np.float64)
edge_points = footprint["edge_points"]
if len(edge_points) == 0:
return {
"time_spread_ns": time_spreads,
"path_spread": path_spreads,
"min_path": min_paths,
"max_path": max_paths,
}
# Pre-compute source-to-surface distances (same for all detectors)
src_to_surface = np.linalg.norm(edge_points - src, axis=1)
for i, det_pos in enumerate(detector_positions):
det = np.array(det_pos, dtype=np.float64)
# Surface-to-detector distances
surface_to_det = np.linalg.norm(edge_points - det, axis=1)
# Total paths
total_paths = src_to_surface + surface_to_det
min_paths[i] = np.min(total_paths)
max_paths[i] = np.max(total_paths)
path_spreads[i] = max_paths[i] - min_paths[i]
time_spreads[i] = path_spreads[i] / SPEED_OF_LIGHT * 1e9
return {
"time_spread_ns": time_spreads,
"path_spread": path_spreads,
"min_path": min_paths,
"max_path": max_paths,
}