Utilities Module
Core Utilities
Utilities Module
Core utility functions and data structures for the ray tracing framework. This module contains foundational components used throughout the simulation.
Submodules
- ray_data
Ray batch data structures for GPU processing (RayBatch, create_ray_batch).
- fresnel
Fresnel equations for reflection/refraction calculations.
- propagation
Ray propagation engine (straight-line and curved-path).
- recording_sphere
(Deprecated) Spherical detector for capturing rays at altitude. Use lsurf.detectors instead for RecordingSphereDetector.
Examples
>>> from lsurf.utilities import RayBatch, create_ray_batch
>>> rays = create_ray_batch(1000)
>>>
>>> from lsurf.utilities import fresnel_coefficients
>>> R, T = fresnel_coefficients(1.0, 1.5, cos_theta)
- class lsurf.utilities.RayBatch(positions, directions, wavelengths, intensities, optical_path_lengths, geometric_path_lengths, accumulated_time, generations, domain_ids, active, polarization_s=None, polarization_p=None, polarization_vector=None, phase=None, optical_depth=None)[source]
Bases:
objectStructure-of-Arrays (SoA) layout for efficient GPU processing of ray batches.
All arrays have shape (N,) where N is the number of rays, except positions and directions which have shape (N, 3).
This layout maximizes memory coalescing on GPU and enables efficient SIMD operations on CPU.
- positions
Ray origin positions, shape (N, 3), units: meters [x0, y0, z0; x1, y1, z1; …]
- Type:
Float32Array
- directions
Ray direction unit vectors, shape (N, 3), dimensionless [dx0, dy0, dz0; dx1, dy1, dz1; …]
- Type:
Float32Array
- wavelengths
Wavelengths, shape (N,), units: meters
- Type:
Float32Array
- intensities
Current ray intensities/powers, shape (N,), units: W or arbitrary
- Type:
Float32Array
- optical_path_lengths
Accumulated optical path length ∫n(s)ds, shape (N,), units: meters
- Type:
Float32Array
- geometric_path_lengths
Accumulated geometric path length ∫ds, shape (N,), units: meters
- Type:
Float32Array
- accumulated_time
Accumulated propagation time, shape (N,), units: seconds
- Type:
Float32Array
- generations
Scattering generation (0=primary, 1=first scatter, etc.), shape (N,)
- Type:
Int32Array
- domain_ids
Current domain ID for each ray, shape (N,)
- Type:
Int32Array
- active
Whether ray is still active (not terminated), shape (N,)
- Type:
BoolArray
- polarization_s
S-polarization component (optional), shape (N,)
- Type:
Optional[Float32Array]
- polarization_p
P-polarization component (optional), shape (N,)
- Type:
Optional[Float32Array]
- polarization_vector
3D polarization vector (electric field direction), shape (N, 3) Unit vector perpendicular to ray direction representing E-field orientation. Used for tracking polarization state through reflections/refractions.
- Type:
Optional[Float32Array]
- phase
Phase in radians (for coherent simulations), shape (N,)
- Type:
Optional[Float32Array]
- optical_depth
Accumulated optical depth τ = ∫α·ds, shape (N,), dimensionless Used for Beer-Lambert absorption tracking. τ = -ln(I/I₀)
- Type:
Optional[Float32Array]
Notes
The SoA layout enables efficient GPU memory access patterns. For example, all x-coordinates are contiguous in memory, allowing coalesced reads.
References
- compact()[source]
Remove inactive rays to reduce memory and computation.
- Returns:
New batch containing only active rays
- Return type:
Notes
This operation is called “stream compaction” in GPU programming. Useful for maintaining efficiency as rays terminate.
- __init__(positions, directions, wavelengths, intensities, optical_path_lengths, geometric_path_lengths, accumulated_time, generations, domain_ids, active, polarization_s=None, polarization_p=None, polarization_vector=None, phase=None, optical_depth=None)
- class lsurf.utilities.RayStatistics(total_rays, active_rays, mean_intensity, total_power, mean_optical_path, mean_generation, max_generation)[source]
Bases:
NamedTupleStatistical summary of ray batch state.
- lsurf.utilities.create_ray_batch(num_rays, enable_polarization=False, enable_polarization_vector=False, enable_phase=False, enable_optical_depth=False)[source]
Create an empty ray batch with zero-initialized arrays.
- Parameters:
num_rays (int) – Number of rays to allocate
enable_polarization (bool, optional) – Whether to allocate scalar polarization arrays (s and p components)
enable_polarization_vector (bool, optional) – Whether to allocate 3D polarization vector array
enable_phase (bool, optional) – Whether to allocate phase array
enable_optical_depth (bool, optional) – Whether to allocate optical depth array for absorption tracking
- Returns:
Initialized ray batch with all fields set to zero/false
- Return type:
- lsurf.utilities.compute_statistics(batch)[source]
Compute statistical summary of a ray batch.
- Parameters:
batch (RayBatch) – Ray batch to analyze
- Returns:
Statistical summary
- Return type:
- lsurf.utilities.merge_ray_batches(batches)[source]
Merge multiple ray batches into a single batch.
- Parameters:
- Returns:
Combined ray batch containing all rays from input batches
- Return type:
Notes
This is useful for combining rays from different sources or for collecting rays after ray splitting at interfaces.
- lsurf.utilities.fresnel_coefficients(n1, n2, cos_theta_i, polarization='unpolarized')[source]
Compute Fresnel reflection and transmission coefficients.
- Parameters:
n1 (float or ndarray) – Refractive index of incident medium
n2 (float or ndarray) – Refractive index of transmitted medium
cos_theta_i (ndarray, shape (N,)) – Cosine of incident angle (dot product of direction and normal)
polarization (str, optional) – Polarization state: ‘s’, ‘p’, or ‘unpolarized’ (default)
- Returns:
R (ndarray, shape (N,)) – Reflection coefficient (fraction of intensity reflected)
T (ndarray, shape (N,)) – Transmission coefficient (fraction of intensity transmitted)
- Return type:
tuple[ndarray[tuple[Any, …], dtype[float32]], ndarray[tuple[Any, …], dtype[float32]]]
Notes
For unpolarized light, we average s and p polarizations. Total internal reflection occurs when n1 > n2 and angle exceeds critical.
The Fresnel equations for intensity (not amplitude) are: - s-polarization: R_s = |r_s|², T_s = (n2*cos_theta_t)/(n1*cos_theta_i) * |t_s|² - p-polarization: R_p = |r_p|², T_p = (n2*cos_theta_t)/(n1*cos_theta_i) * |t_p|²
Examples
>>> # Air to glass at 45 degrees >>> cos_theta_i = np.cos(np.radians(45)) >>> R, T = fresnel_coefficients(1.0, 1.5, cos_theta_i) >>> print(f"Reflection: {R:.3f}, Transmission: {T:.3f}")
- lsurf.utilities.compute_reflection_direction(incident, normal)[source]
Compute reflected ray direction using law of reflection.
- Parameters:
incident (ndarray, shape (N, 3)) – Incident ray directions (should be normalized)
normal (ndarray, shape (N, 3)) – Surface normals at intersection points (should be normalized)
- Returns:
reflected – Reflected ray directions (normalized)
- Return type:
ndarray, shape (N, 3)
Notes
Reflection formula: r = d - 2(d·n)n where d is incident direction and n is surface normal.
The normal should point toward the incident side.
Examples
>>> # Reflect ray at 45° off horizontal surface >>> incident = np.array([[1/np.sqrt(2), 0, -1/np.sqrt(2)]]) >>> normal = np.array([[0, 0, 1]]) >>> reflected = compute_reflection_direction(incident, normal) >>> print(reflected) # Should be [1/√2, 0, 1/√2]
- lsurf.utilities.compute_refraction_direction(incident, normal, n1, n2)[source]
Compute refracted ray direction using Snell’s law.
- Parameters:
- Returns:
refracted (ndarray, shape (N, 3)) – Refracted ray directions (normalized) For TIR cases, returns zero vector
tir_mask (ndarray, shape (N,)) – Boolean mask indicating total internal reflection
- Return type:
tuple[ndarray[tuple[Any, …], dtype[float32]], ndarray[tuple[Any, …], dtype[bool]]]
Notes
Snell’s law: n1*sin(θ1) = n2*sin(θ2)
Vector form of Snell’s law: t = (n1/n2)[d - (d·n)n] - n*sqrt(1 - (n1/n2)²*[1-(d·n)²])
Total internal reflection occurs when (n1/n2)²*[1-(d·n)²] > 1
Examples
>>> # Air to glass at normal incidence >>> incident = np.array([[0, 0, -1]]) >>> normal = np.array([[0, 0, 1]]) >>> refracted, tir = compute_refraction_direction(incident, normal, 1.0, 1.5) >>> print(refracted) # Should be [0, 0, -1] (straight through) >>> print(tir) # Should be False
- class lsurf.utilities.PropagationMode(*values)[source]
Bases:
EnumRay propagation mode selector.
- STRAIGHT_LINE = 0
- CURVED_PATH = 1
- ADAPTIVE = 2
- class lsurf.utilities.PropagationConfig(mode=PropagationMode.ADAPTIVE, max_step_size=0.001, min_step_size=1e-06, gradient_threshold=1e-06, adaptive_stepping=True, rk4_integration=True)[source]
Bases:
objectConfiguration for ray propagation.
- mode
Propagation mode
- Type:
Notes
For gradient media, step size should satisfy: Δs ≪ R_c where R_c = n/|∇n| is the radius of curvature
- __init__(mode=PropagationMode.ADAPTIVE, max_step_size=0.001, min_step_size=1e-06, gradient_threshold=1e-06, adaptive_stepping=True, rk4_integration=True)[source]
Initialize propagation configuration.
- Parameters:
mode (PropagationMode, optional) – Propagation mode, default ADAPTIVE
max_step_size (float, optional) – Maximum step size in meters, default 1 mm
min_step_size (float, optional) – Minimum step size in meters, default 1 μm
gradient_threshold (float, optional) – Gradient threshold in m^-1, default 1e-6
adaptive_stepping (bool, optional) – Enable adaptive stepping, default True
rk4_integration (bool, optional) – Use RK4 integration, default True
- class lsurf.utilities.StraightLinePropagator(threads_per_block=256, prefer_gpu=True)[source]
Bases:
objectPropagator for homogeneous media using straight-line propagation.
This is the simplest and fastest propagation mode, appropriate when the refractive index gradient is negligible.
Automatically falls back to CPU if CUDA is not available.
- propagate_step(rays, step_size, material)[source]
Propagate rays in straight lines.
- Parameters:
rays (RayBatch) – Ray batch to propagate (modified in-place)
step_size (float) – Step size in meters
material (MaterialField) – Material field (not used for straight-line propagation)
- class lsurf.utilities.GradientPropagator(use_rk4=None, method=None, adaptive_stepping=True, min_step_size=1e-06, max_step_size=0.001, threads_per_block=256)[source]
Bases:
objectPropagator for inhomogeneous media using ray equation integration.
This propagator handles materials with spatially-varying refractive indices by numerically integrating the ray equation (eikonal equation). Supports both Euler and RK4 integration methods.
References
- __init__(use_rk4=None, method=None, adaptive_stepping=True, min_step_size=1e-06, max_step_size=0.001, threads_per_block=256)[source]
Initialize gradient propagator.
- Parameters:
use_rk4 (bool, optional) – Use RK4 instead of Euler (deprecated, use method instead)
method (str, optional) – Integration method: “euler” or “rk4”, default “rk4”
adaptive_stepping (bool, optional) – Enable adaptive step sizing, default True
min_step_size (float, optional) – Minimum step size in meters, default 1 μm
max_step_size (float, optional) – Maximum step size in meters, default 1 mm
threads_per_block (int, optional) – CUDA threads per block (legacy), default 256
- compute_gradient_threshold(wavelength)[source]
Compute gradient threshold for switching to curved-path mode.
- Parameters:
wavelength (float) – Wavelength in meters
- Returns:
Gradient magnitude threshold in m^-1
- Return type:
Notes
Threshold is chosen such that ray curvature becomes significant over a distance of ~wavelength: |∇n|/n > 1/λ
- propagate_step_cpu(positions, directions, active, material, step_size, wavelength, geometric_path_lengths=None, optical_path_lengths=None, accumulated_time=None)[source]
Propagate rays through gradient medium using CPU (vectorized NumPy).
This method modifies arrays in-place using Euler or RK4 integration.
- Parameters:
positions (ndarray) – Ray positions, shape (N, 3), modified in-place.
directions (ndarray) – Ray directions, shape (N, 3), modified in-place.
active (ndarray) – Active ray mask, shape (N,).
material (MaterialFieldProtocol) – Material with get_refractive_index and get_refractive_index_gradient.
step_size (float) – Step size in meters.
wavelength (float) – Wavelength in meters.
geometric_path_lengths (ndarray, optional) – Geometric path lengths, shape (N,), updated in-place.
optical_path_lengths (ndarray, optional) – Optical path lengths, shape (N,), updated in-place.
accumulated_time (ndarray, optional) – Accumulated time, shape (N,), updated in-place.
- propagate_step(rays, step_size, wavelength=5.32e-07, material=None)[source]
Propagate rays by a single integration step.
This is the unified interface method compatible with GPU propagators. Delegates to propagate_step_cpu internally.
- Parameters:
rays (RayBatch) – Ray batch containing positions, directions, and other ray properties. Modified in-place.
step_size (float) – Integration step size in meters.
wavelength (float, optional) – Wavelength in meters, default 532 nm.
material (MaterialFieldProtocol, optional) – Material to propagate through. If not provided, must be set via a separate mechanism (e.g., stored during initialization).
- propagate(rays, total_distance, step_size, wavelength=5.32e-07, material=None)[source]
Propagate rays through a total distance.
This is the unified interface method compatible with GPU propagators.
- Parameters:
rays (RayBatch) – Ray batch to propagate (modified in-place)
total_distance (float) – Total distance to propagate in meters
step_size (float) – Integration step size in meters
wavelength (float, optional) – Wavelength in meters, default 532 nm
material (MaterialFieldProtocol, optional) – Material to propagate through
- class lsurf.utilities.IPropagator(*args, **kwargs)[source]
Bases:
ProtocolProtocol (interface) for ray propagation implementations.
All propagators must implement this interface to ensure compatibility with the simulation framework.
- propagate_step(rays, step_size, material)[source]
Propagate rays forward by one step.
- Parameters:
rays (RayBatch) – Ray batch to propagate (modified in-place)
step_size (float) – Step size in meters
material (MaterialField) – Material field defining refractive index
- compute_gradient_threshold(wavelength)[source]
Compute threshold for switching between propagation modes.
- __init__(*args, **kwargs)
- class lsurf.utilities.RecordingSphere(altitude=33000.0, earth_center=(0, 0, -6371000.0), earth_radius=6371000.0)[source]
Bases:
objectSpherical detection surface at a specified altitude above Earth.
Records all rays that intersect the sphere, capturing full ray state for later analysis.
- Parameters:
Notes
The recording sphere has radius = earth_radius + altitude, centered at earth_center.
- detect(rays, compute_travel_time=True, speed_of_light=299792458.0, max_propagation_distance=None)[source]
Detect rays intersecting the recording sphere.
- Parameters:
rays (RayBatch) – Rays to detect
compute_travel_time (bool) – If True, add travel time to intersection to ray’s accumulated time
speed_of_light (float) – Speed of light for time computation
max_propagation_distance (float, optional) – Maximum distance rays can propagate before detection (meters). If None, no limit is applied. Use this to prevent detecting rays that would hit the sphere from the opposite hemisphere.
- Returns:
Recorded ray data for all intersecting rays
- Return type:
- class lsurf.utilities.LocalRecordingSphere(radius=33000.0, center=(0, 0, 0))[source]
Bases:
objectSimple spherical detection surface centered at origin.
Records all rays that intersect the sphere from inside, useful for local-scale simulations without Earth curvature.
- Parameters:
- class lsurf.utilities.RecordedRays(positions, directions, times, intensities, wavelengths, generations, polarization_vectors=None, ray_indices=None)[source]
Bases:
objectContainer for recorded ray data at the detection sphere.
All arrays have shape (N,) or (N, 3) where N is the number of recorded rays.
- positions
Intersection positions on the recording sphere (meters)
- Type:
ndarray, shape (N, 3)
- directions
Ray directions at intersection (unit vectors)
- Type:
ndarray, shape (N, 3)
- times
Time of arrival at recording sphere (seconds)
- Type:
ndarray, shape (N,)
- intensities
Ray intensity at recording sphere
- Type:
ndarray, shape (N,)
- wavelengths
Ray wavelength (meters)
- Type:
ndarray, shape (N,)
- generations
Ray generation (number of surface interactions)
- Type:
ndarray, shape (N,)
- polarization_vectors
3D polarization vectors (electric field direction) at intersection. Unit vectors perpendicular to ray direction representing E-field orientation.
- Type:
ndarray, shape (N, 3), optional
- Derived quantities (computed on demand or at save time)
- elevation_angles
Elevation angle above local horizon (radians)
- Type:
ndarray, shape (N,)
- azimuth_angles
Azimuth angle in local tangent plane (radians)
- Type:
ndarray, shape (N,)
- zenith_angles
Angle from local zenith (radians)
- Type:
ndarray, shape (N,)
- local_positions
Position in local tangent coordinates (meters)
- Type:
ndarray, shape (N, 2)
- compute_angular_coordinates(earth_center=None)[source]
Compute angular coordinates for all recorded ray intersection points.
Computes spherical coordinates (latitude/longitude) of intersection points on the detection sphere relative to Earth’s center.
- Parameters:
earth_center (ndarray, optional) – Earth center position, default (0, 0, -EARTH_RADIUS)
- Returns:
Dictionary with: - ‘elevation’: Latitude angle above equator (radians, -π/2 to π/2) - ‘azimuth’: Longitude angle (radians, -π to π) - ‘zenith’: Zenith angle from north pole (radians, 0 to π) - ‘incidence’: Angle between ray direction and outward radial (radians)
- Return type:
- compute_viewing_angle_from_origin(origin=None)[source]
Compute viewing angle from horizontal at specified origin.
Calculates the angle above the horizontal plane (XY plane) when viewing each intersection point from the origin position.
- Parameters:
origin (ndarray, optional) – Observer position, default (0, 0, 0) - Earth surface at reference point
- Returns:
Viewing angle from horizontal in radians (-π/2 to π/2) Positive angles are above horizontal, negative below
- Return type:
ndarray
- compute_ray_direction_angles()[source]
Compute elevation and azimuth angles of ray directions.
Calculates the angles of the ray direction vectors themselves, not the position coordinates. Useful for understanding the angular distribution of ray propagation.
- Returns:
Dictionary with: - ‘elevation’: Angle above horizontal plane in radians (-π/2 to π/2)
0 = horizontal, π/2 = straight up, -π/2 = straight down
- ’azimuth’: Azimuth angle in horizontal plane in radians (-π to π)
0 = +X direction, π/2 = +Y direction
- Return type:
- __init__(positions, directions, times, intensities, wavelengths, generations, polarization_vectors=None, ray_indices=None)
- lsurf.utilities.save_recorded_rays_hdf5(recorded_rays, filepath, metadata=None, compression='gzip')[source]
Save recorded rays to HDF5 file.
- Parameters:
recorded_rays (RecordedRays) – Ray data to save
filepath (str) – Output file path
metadata (dict, optional) – Additional metadata to store (simulation parameters, etc.)
compression (str) – Compression algorithm (‘gzip’, ‘lzf’, or None)
- lsurf.utilities.save_recorded_rays_numpy(recorded_rays, filepath, metadata=None)[source]
Save recorded rays to numpy .npz file.
- Parameters:
recorded_rays (RecordedRays) – Ray data to save
filepath (str) – Output file path
metadata (dict, optional) – Additional metadata to store
- lsurf.utilities.load_recorded_rays_numpy(filepath)[source]
Load recorded rays from numpy .npz file.
- class lsurf.utilities.TimeSpreadResult(min_path, max_path, path_spread, time_spread_s, time_spread_ns, min_path_point, max_path_point, edge_points, path_lengths, center_point)[source]
Bases:
objectResults from time spread estimation.
- min_path_point
Surface point with shortest path
- Type:
ndarray
- max_path_point
Surface point with longest path
- Type:
ndarray
- edge_points
All beam edge intersection points on surface
- Type:
ndarray
- path_lengths
Path lengths for all edge points
- Type:
ndarray
- center_point
Center ray intersection point
- Type:
ndarray
- __init__(min_path, max_path, path_spread, time_spread_s, time_spread_ns, min_path_point, max_path_point, edge_points, path_lengths, center_point)
- lsurf.utilities.estimate_time_spread(source_position, beam_direction, divergence_angle, detector_position, surface=None, n_edge_points=100, speed_of_light=299792458.0)[source]
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:
Dataclass containing path lengths, time spread, and geometry info
- Return type:
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")
- lsurf.utilities.estimate_time_spread_bounds(source_position, beam_direction, divergence_angle, detector_positions, surface=None, n_edge_points=100)[source]
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:
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
- Return type:
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)
- lsurf.utilities.compute_beam_footprint(source_position, beam_direction, divergence_angle, n_edge_points=100, surface=None)[source]
Compute where the edges of a diverging beam hit a surface.
- Parameters:
source_position (tuple) – (x, y, z) position of the source in meters
beam_direction (tuple) – Unit vector of beam center direction
divergence_angle (float) – Half-angle divergence of the beam in radians
n_edge_points (int) – Number of points around the beam edge cone
surface (Surface, optional) – Surface object to intersect with (e.g., GerstnerWaveSurface, CurvedWaveSurface, PlanarSurface). If None, uses flat surface at z=0.
- Returns:
Dictionary containing: - edge_points : ndarray, shape (N, 3) - Surface intersection points - center_point : ndarray, shape (3,) - Center ray intersection - edge_distances : ndarray - Distances from source to each edge point - center_distance : float - Distance from source to center point
- Return type:
- lsurf.utilities.find_peak_energy_density(recorded_rays, detector_radius, detector_center=None, n_bins=50)[source]
Scan the detector sphere to find the location with highest energy density.
- Parameters:
recorded_rays (RecordedRays) – Recorded rays on detection sphere
detector_radius (float) – Radius of the detector sphere (m)
detector_center (array-like, optional) – Center of the detector sphere (default origin)
n_bins (int) – Number of bins in each angular dimension
- Returns:
Dictionary containing: - peak_lon, peak_lat: Peak location in radians - peak_lon_deg, peak_lat_deg: Peak location in degrees - peak_position: Peak location in Cartesian (x, y, z) - peak_irradiance: Irradiance at peak (W/m²) - total_power: Total detected power - histogram: 2D irradiance histogram - lon_edges, lat_edges: Bin edges - lon_centers, lat_centers: Bin centers
- Return type:
Notes
This function uses lat/lon binning which has anisotropic resolution (bins are smaller near poles). For more uniform resolution, use find_peak_irradiance_local() which uses local tangent plane coordinates.
- lsurf.utilities.find_peak_irradiance_local(recorded_rays, detector_radius, detector_center=None, n_bins=50, bin_size_meters=None)[source]
Find peak irradiance using local tangent plane coordinates.
Uses a local East-North coordinate system centered on the peak region, providing constant resolution in physical units (meters). This avoids the polar singularities and anisotropic bin sizes of lat/lon binning.
- Parameters:
recorded_rays (RecordedRays) – Recorded rays on detection sphere
detector_radius (float) – Radius of the detector sphere (m)
detector_center (array-like, optional) – Center of the detector sphere (default origin)
n_bins (int) – Number of bins in each dimension (default 50)
bin_size_meters (float, optional) – Physical size of each bin in meters. If None, auto-computed to cover the data extent with n_bins.
- Returns:
Dictionary containing: - peak_east, peak_north: Peak location in local coordinates (m) - peak_lon, peak_lat: Peak location in radians - peak_lon_deg, peak_lat_deg: Peak location in degrees - peak_position: Peak location in Cartesian (x, y, z) - peak_irradiance: Irradiance at peak (W/m²) - total_power: Total detected power (W) - histogram: 2D irradiance histogram - east_edges, north_edges: Bin edges in local coords (m) - bin_size: Physical bin size (m)
- Return type:
- lsurf.utilities.compute_pareto_front(recorded_rays, detector_radius, detector_center=None, n_bins=30, min_rays_per_bin=5)[source]
Compute Pareto front for irradiance vs time spread trade-off.
- Parameters:
recorded_rays (RecordedRays) – Recorded rays on detection sphere
detector_radius (float) – Radius of the detector sphere (m)
detector_center (array-like, optional) – Center of the detector sphere (default origin)
n_bins (int) – Number of bins in each angular dimension
min_rays_per_bin (int) – Minimum rays required per bin for reliable statistics
- Returns:
Dictionary containing: - bin_data: List of dicts with (lon, lat, irradiance, time_spread, n_rays) - pareto_indices: Indices of Pareto-optimal bins - pareto_front: List of Pareto-optimal bin data - all_irradiances: Array of all bin irradiances (W/m²) - all_time_spreads: Array of all bin time spreads
- Return type:
- lsurf.utilities.analyze_healpix_detector(recorded_rays, detector_radius, nside=128, min_rays_per_pixel=5)[source]
Analyze detected rays using HEALPix equal-area pixelization.
HEALPix provides equal-area pixels on the sphere, avoiding the polar singularities and anisotropic bin sizes of lat/lon binning.
- Parameters:
recorded_rays (RecordedRays) – Recorded ray data from detection sphere
detector_radius (float) – Radius of detection sphere in meters
nside (int) – HEALPix resolution parameter (default 128 → ~13,000 pixels) nside=64 → ~3,400 pixels nside=128 → ~13,000 pixels nside=256 → ~50,000 pixels
min_rays_per_pixel (int) – Minimum rays required per pixel for valid statistics
- Returns:
Dictionary containing: - peak_pixel_id: HEALPix pixel with highest irradiance - peak_lon_deg, peak_lat_deg: Peak location in degrees - peak_irradiance: Irradiance at peak (W/m²) - total_power: Total detected power (W) - pixel_area: Area of each HEALPix pixel (m²) - pixel_data: List of dicts with per-pixel data:
pixel_id, lon_deg, lat_deg
irradiance (W/m²)
time_spread_ns (90th-10th percentile)
ray_count
healpix_data: Full HEALPixData object
- Return type:
- Raises:
ImportError – If astropy-healpix is not installed
- lsurf.utilities.compute_healpix_pareto_front(recorded_rays, detector_radius, nside=64, min_rays_per_pixel=10)[source]
Compute Pareto front for irradiance vs time spread using HEALPix.
A pixel is Pareto-optimal if no other pixel has BOTH higher irradiance AND lower time spread.
- Parameters:
recorded_rays (RecordedRays) – Recorded ray data from detection sphere
detector_radius (float) – Radius of detection sphere in meters
nside (int) – HEALPix resolution parameter (default 64 → ~3,400 pixels)
min_rays_per_pixel (int) – Minimum rays required per pixel for valid statistics
- Returns:
Dictionary containing: - pixel_data: List of all valid pixel dicts - pareto_indices: Indices of Pareto-optimal pixels - pareto_front: List of Pareto-optimal pixel data - all_irradiances: Array of all pixel irradiances - all_time_spreads: Array of all pixel time spreads
- Return type:
- Raises:
ImportError – If astropy-healpix is not installed
- lsurf.utilities.process_surface_interaction(rays, surface, wavelength=5e-07, generate_reflected=True, generate_refracted=True, polarization='unpolarized', track_polarization_vector=False)[source]
Process rays intersecting a surface.
Computes intersections, applies Fresnel equations, and generates reflected and/or refracted ray bundles.
- Parameters:
rays (RayBatch) – Input rays to test for intersection
surface (Surface) – Surface to intersect with
wavelength (float or ndarray, optional) – Wavelength for computing refractive indices (default: 500 nm) If rays have multiple wavelengths, use rays.wavelengths
generate_reflected (bool, optional) – Whether to generate reflected rays (default: True)
generate_refracted (bool, optional) – Whether to generate refracted rays (default: True)
polarization (str, optional) – Polarization state: ‘s’, ‘p’, or ‘unpolarized’ (default) Used for Fresnel coefficient calculation.
track_polarization_vector (bool, optional) – Whether to track 3D polarization vectors through the interaction. If True, polarization vectors are initialized (if not present) and transformed through reflection/refraction. (default: False)
- Returns:
reflected_rays (RayBatch or None) – Reflected rays (None if generate_reflected=False or no hits)
refracted_rays (RayBatch or None) – Refracted rays (None if generate_refracted=False or no hits)
- Return type:
Notes
Only active rays are tested for intersection. Input rays are not modified.
When track_polarization_vector=True, the function: 1. Initializes polarization vectors if rays.polarization_vector is None 2. Transforms polarization vectors through reflection/refraction 3. Stores the transformed vectors in the output ray batches
Examples
>>> # Create surface and rays >>> surface = PlanarSurface( ... point=(0, 0, 1), ... normal=(0, 0, -1), ... material_front=BK7_GLASS, ... material_back=AIR_STP ... ) >>> rays = ... # Create rays >>> >>> # Process interaction with polarization tracking >>> reflected, refracted = process_surface_interaction( ... rays, surface, generate_reflected=True, generate_refracted=True, ... track_polarization_vector=True ... )
- lsurf.utilities.reflect_rays(rays, surface, wavelength=5e-07, polarization='unpolarized', in_place=False)[source]
Apply reflection to rays at surface.
Updates ray directions and intensities based on Fresnel reflection. Rays that don’t hit the surface are deactivated.
- Parameters:
- Returns:
Reflected rays (same object if in_place=True)
- Return type:
- lsurf.utilities.refract_rays(rays, surface, wavelength=5e-07, polarization='unpolarized', in_place=False)[source]
Apply refraction to rays at surface.
Updates ray directions and intensities based on Fresnel transmission. Rays undergoing total internal reflection are deactivated.
- Parameters:
- Returns:
Refracted rays (same object if in_place=True)
- Return type:
- lsurf.utilities.trace_rays_multi_bounce(rays, surface, max_bounces=10, bounding_radius=10000.0, wavelength=5.32e-07, min_intensity=1e-10, track_refracted=True)[source]
Trace rays through multiple surface interactions until termination.
Rays are traced until they: - Exit the bounding sphere - Reach maximum number of bounces - Have intensity below threshold
- Parameters:
rays (RayBatch) – Initial ray batch to trace
surface (Surface) – Surface to interact with
max_bounces (int, optional) – Maximum number of surface interactions (default: 10)
bounding_radius (float, optional) – Radius of bounding sphere in meters (default: 10000)
wavelength (float, optional) – Wavelength for Fresnel calculations (default: 532nm)
min_intensity (float, optional) – Minimum intensity threshold (default: 1e-10)
track_refracted (bool, optional) – Whether to track refracted rays (default: True)
- Returns:
final_reflected (RayBatch) – Final state of all reflected rays that exited bounding sphere
final_refracted (RayBatch) – Final state of all refracted rays (combined from all bounces)
ray_paths (dict) – Dictionary containing: - ‘reflected_paths’: list of paths, each path is Nx3 array of positions - ‘refracted_paths’: list of paths for refracted rays (start from refraction point) - ‘reflected_final_dirs’: list of final direction vectors for each reflected path - ‘refracted_final_dirs’: list of final direction vectors for each refracted path
- Return type:
Notes
The function tracks reflected rays through multiple bounces on the wavy surface. Each ray’s complete path (all positions) is stored for visualization. Refracted rays are collected but not further traced (they go into the water and don’t re-interact with the air-water interface from below in typical scenarios).
Examples
>>> # Trace rays with up to 5 bounces >>> final_refl, final_refr, paths = trace_rays_multi_bounce( ... rays, surface, max_bounces=5, bounding_radius=5000.0 ... ) >>> # Plot a ray's path >>> path = paths['reflected_paths'][0] # First ray's path >>> plt.plot(path[:, 0], path[:, 2]) # x-z projection
- lsurf.utilities.trace_rays_with_splitting(rays, surfaces, max_bounces=10, bounding_radius=10000.0, bounding_center=None, wavelength=5.32e-07, min_intensity=1e-10, polarization='unpolarized')[source]
Trace rays through multiple surfaces with proper Fresnel ray splitting.
At each surface interaction, every ray is split into two child rays: - A reflected ray with intensity scaled by the Fresnel reflection coefficient R - A refracted ray with intensity scaled by the Fresnel transmission coefficient T
This creates a ray tree where both reflected and refracted paths are traced until termination conditions are met.
- Parameters:
rays (RayBatch) – Initial ray batch to trace
surfaces (list of Surface) – List of surfaces to interact with (tested in order for closest hit)
max_bounces (int, optional) – Maximum number of surface interactions per ray tree branch (default: 10)
bounding_radius (float, optional) – Radius of bounding sphere in meters (default: 10000)
bounding_center (tuple of float, optional) – Center of bounding sphere (x, y, z) in meters. Default is (0, 0, 0). For curved Earth simulations, use (0, 0, -EARTH_RADIUS).
wavelength (float, optional) – Wavelength for Fresnel calculations (default: 532nm)
min_intensity (float, optional) – Minimum intensity threshold for ray termination (default: 1e-10)
polarization (str, optional) – Polarization state: ‘s’, ‘p’, or ‘unpolarized’ (default)
- Returns:
final_rays (RayBatch) – All terminal rays (rays that exited the scene or hit max bounces) Each ray has its intensity weighted by the product of all Fresnel coefficients along its path.
trace_info (dict) – Dictionary containing tracing statistics: - ‘total_rays_created’: Total number of rays created during tracing - ‘max_depth_reached’: Maximum tree depth reached - ‘terminated_by_intensity’: Number of rays terminated due to low intensity - ‘terminated_by_bounds’: Number of rays that exited bounding sphere - ‘terminated_by_max_bounces’: Number of rays that hit max bounce limit
- Return type:
Notes
This function implements proper physical ray splitting where each ray at an interface creates two child rays:
Reflected ray: direction from law of reflection, intensity = I_parent * R
Refracted ray: direction from Snell’s law, intensity = I_parent * T
where R and T are the Fresnel reflection and transmission coefficients computed from the refractive indices and incident angle.
For total internal reflection (TIR), T=0 and R=1, so only a reflected ray is created.
The number of rays grows exponentially with depth (up to 2^depth), so the min_intensity threshold is critical for pruning weak ray branches.
Examples
>>> from surface_roughness.surfaces import PlanarSurface >>> from surface_roughness.materials import Glass, Air >>> >>> # Create a glass slab >>> surface1 = PlanarSurface( ... point=(0, 0, 0.01), ... normal=(0, 0, -1), ... material_front=Air(), ... material_back=Glass() ... ) >>> surface2 = PlanarSurface( ... point=(0, 0, 0.02), ... normal=(0, 0, -1), ... material_front=Glass(), ... material_back=Air() ... ) >>> >>> # Trace rays with splitting >>> final_rays, info = trace_rays_with_splitting( ... rays, [surface1, surface2], max_bounces=5 ... ) >>> print(f"Created {info['total_rays_created']} rays total")