# 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.
"""
Gerstner Wave Surface (CPU-Only)
Flat-earth Gerstner wave ocean surface. Uses ray marching for intersection
since the surface shape is too complex for closed-form GPU signed distance.
For curved Earth surfaces with waves, see CurvedWaveSurface.
"""
from dataclasses import dataclass, field
from typing import Any
import numpy as np
from numpy.typing import NDArray
from ..protocol import Surface, SurfaceRole
from ..registry import register_surface_type
from .wave_params import GerstnerWaveParams
[docs]
@dataclass
class GerstnerWaveSurface(Surface):
"""
Flat ocean surface with Gerstner wave physics (CPU-only).
Implements the Gerstner (trochoidal) wave model where water particles
move in circular orbits, producing realistic ocean wave shapes with
sharp crests and flat troughs.
Parameters
----------
wave_params : list of GerstnerWaveParams
List of wave components to superimpose.
role : SurfaceRole
What happens when a ray hits (typically OPTICAL).
name : str
Human-readable name.
time : float, optional
Time for wave animation in seconds. Default is 0.0.
reference_z : float, optional
Mean sea level z-coordinate in meters. Default is 0.0.
material_front : MaterialField, optional
Material above surface (air).
material_back : MaterialField, optional
Material below surface (water).
max_distance : float, optional
Maximum ray marching distance in meters. Default is 10000.0.
Examples
--------
>>> from lsurf.surfaces import GerstnerWaveSurface, GerstnerWaveParams, SurfaceRole
>>> from lsurf.materials import AIR_STP, WATER
>>>
>>> wave = GerstnerWaveParams(amplitude=1.0, wavelength=50.0)
>>> surface = GerstnerWaveSurface(
... wave_params=[wave],
... role=SurfaceRole.OPTICAL,
... name="ocean",
... material_front=AIR_STP,
... material_back=WATER,
... )
"""
wave_params: list[GerstnerWaveParams]
role: SurfaceRole
name: str = "gerstner_wave"
time: float = 0.0
reference_z: float = 0.0
material_front: Any = None
material_back: Any = None
max_distance: float = 10000.0
# CPU-only surface
_gpu_capable: bool = field(default=False, init=False, repr=False)
_geometry_id: int = field(
default=0, init=False, repr=False
) # CPU-only, no GPU geometry
# Precomputed wave arrays (set in __post_init__)
_amplitudes: NDArray = field(default=None, init=False, repr=False)
_wave_numbers: NDArray = field(default=None, init=False, repr=False)
_frequencies: NDArray = field(default=None, init=False, repr=False)
_directions: NDArray = field(default=None, init=False, repr=False)
_phases: NDArray = field(default=None, init=False, repr=False)
_steepness: NDArray = field(default=None, init=False, repr=False)
def __post_init__(self) -> None:
self._precompute_wave_params()
def _precompute_wave_params(self) -> None:
"""Precompute wave parameters as arrays for fast evaluation."""
n_waves = len(self.wave_params)
if n_waves == 0:
self._amplitudes = np.array([], dtype=np.float64)
self._wave_numbers = np.array([], dtype=np.float64)
self._frequencies = np.array([], dtype=np.float64)
self._directions = np.zeros((0, 2), dtype=np.float64)
self._phases = np.array([], dtype=np.float64)
self._steepness = np.array([], dtype=np.float64)
return
self._amplitudes = np.array(
[w.amplitude for w in self.wave_params], dtype=np.float64
)
self._wave_numbers = np.array(
[w.wave_number for w in self.wave_params], dtype=np.float64
)
self._frequencies = np.array(
[w.angular_frequency for w in self.wave_params], dtype=np.float64
)
self._directions = np.array(
[w.direction_normalized for w in self.wave_params], dtype=np.float64
)
self._phases = np.array([w.phase for w in self.wave_params], dtype=np.float64)
self._steepness = np.array(
[w.steepness for w in self.wave_params], dtype=np.float64
)
@property
def gpu_capable(self) -> bool:
"""This surface does NOT support GPU acceleration."""
return False
@property
def geometry_id(self) -> int:
"""GPU geometry type ID (0 = CPU-only)."""
return 0
[docs]
def get_materials(self) -> tuple | None:
"""Return materials for Fresnel calculation."""
if self.role == SurfaceRole.OPTICAL:
return (self.material_front, self.material_back)
return None
[docs]
def get_max_wave_height(self) -> float:
"""Get maximum possible wave height above reference_z."""
if len(self._amplitudes) == 0:
return 0.0
return float(np.sum(self._amplitudes))
def _compute_displacement(
self,
x: NDArray[np.float64],
y: NDArray[np.float64],
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
"""Compute Gerstner wave displacement at positions (x, y)."""
x = np.atleast_1d(x)
y = np.atleast_1d(y)
n_points = len(x)
dx = np.zeros(n_points, dtype=np.float64)
dy = np.zeros(n_points, dtype=np.float64)
dz = np.zeros(n_points, dtype=np.float64)
for i in range(len(self.wave_params)):
A = self._amplitudes[i]
k = self._wave_numbers[i]
omega = self._frequencies[i]
dir_x, dir_y = self._directions[i]
phi = self._phases[i]
Q = self._steepness[i]
phase = k * (dir_x * x + dir_y * y) - omega * self.time + phi
cos_phase = np.cos(phase)
sin_phase = np.sin(phase)
dx -= Q * A * dir_x * sin_phase
dy -= Q * A * dir_y * sin_phase
dz += A * cos_phase
return dx, dy, dz
def _compute_normal(
self,
x: NDArray[np.float64],
y: NDArray[np.float64],
) -> NDArray[np.float64]:
"""Compute surface normal at positions (x, y)."""
n_points = len(x)
nx = np.zeros(n_points, dtype=np.float64)
ny = np.zeros(n_points, dtype=np.float64)
nz = np.ones(n_points, dtype=np.float64)
for i in range(len(self.wave_params)):
A = self._amplitudes[i]
k = self._wave_numbers[i]
omega = self._frequencies[i]
dir_x, dir_y = self._directions[i]
phi = self._phases[i]
Q = self._steepness[i]
phase = k * (dir_x * x + dir_y * y) - omega * self.time + phi
cos_phase = np.cos(phase)
sin_phase = np.sin(phase)
WA = k * A
nx += dir_x * WA * sin_phase
ny += dir_y * WA * sin_phase
nz -= Q * WA * cos_phase
norms = np.sqrt(nx**2 + ny**2 + nz**2)
norms = np.maximum(norms, 1e-10)
normals = np.stack([nx / norms, ny / norms, nz / norms], axis=-1)
return normals.astype(np.float32)
def _surface_z(
self,
x: NDArray[np.float64],
y: NDArray[np.float64],
) -> NDArray[np.float64]:
"""Compute surface height z at positions (x, y)."""
is_scalar = np.isscalar(x) and np.isscalar(y)
_, _, dz = self._compute_displacement(x, y)
result = self.reference_z + dz
if is_scalar:
return float(result[0])
return result
[docs]
def intersect(
self,
origins: NDArray[np.float32],
directions: NDArray[np.float32],
min_distance: float = 1e-6,
) -> tuple[NDArray[np.float32], NDArray[np.bool_]]:
"""
Find ray-surface intersections using ray marching.
Parameters
----------
origins : ndarray, shape (N, 3)
Ray origin positions.
directions : ndarray, shape (N, 3)
Ray direction unit vectors.
min_distance : float
Minimum valid intersection distance.
Returns
-------
distances : ndarray, shape (N,)
Distance to intersection (inf if no hit).
hit_mask : ndarray, shape (N,), dtype=bool
True for rays that hit the surface.
"""
origins = origins.astype(np.float64)
directions = directions.astype(np.float64)
n_rays = len(origins)
distances = np.full(n_rays, np.inf, dtype=np.float64)
hit_mask = np.zeros(n_rays, dtype=bool)
max_wave_height = self.get_max_wave_height()
z_max = self.reference_z + max_wave_height
z_min = self.reference_z - max_wave_height
tolerance = 1e-4
max_iterations = 100
for i in range(n_rays):
ox, oy, oz = origins[i]
dx, dy, dz = directions[i]
# Skip rays parallel to surface
if abs(dz) < 1e-10:
if oz < z_min or oz > z_max:
continue
# Initial guess: intersection with mean plane
if abs(dz) > 1e-10:
t = (self.reference_z - oz) / dz
else:
t = 0.0
t = max(t - max_wave_height / max(abs(dz), 0.1), min_distance)
for _ in range(max_iterations):
px = ox + t * dx
py = oy + t * dy
pz = oz + t * dz
z_surf = self._surface_z(
np.array([px], dtype=np.float64), np.array([py], dtype=np.float64)
)[0]
signed_dist = pz - z_surf
if abs(signed_dist) < tolerance:
if t > min_distance:
distances[i] = t
hit_mask[i] = True
break
if abs(dz) > 0.01:
step = signed_dist / abs(dz) * 0.8
else:
step = signed_dist * 0.5
step = np.clip(step, -self.max_distance * 0.1, self.max_distance * 0.1)
t += step
if t < 0 or t > self.max_distance:
break
return distances.astype(np.float32), hit_mask
[docs]
def normal_at(
self,
positions: NDArray[np.float32],
incoming_directions: NDArray[np.float32] | None = None,
) -> NDArray[np.float32]:
"""
Compute surface normal at given positions.
Parameters
----------
positions : ndarray, shape (N, 3)
Points on the surface.
incoming_directions : ndarray, shape (N, 3), optional
Incoming ray directions.
Returns
-------
normals : ndarray, shape (N, 3)
Unit normal vectors.
"""
x = positions[:, 0].astype(np.float64)
y = positions[:, 1].astype(np.float64)
normals = self._compute_normal(x, y)
if incoming_directions is not None:
dot_products = np.sum(incoming_directions * normals, axis=1)
flip_mask = dot_products > 0
normals[flip_mask] *= -1
return normals
[docs]
def set_time(self, time: float) -> None:
"""Update the wave animation time."""
self.time = time
# Register class with registry
register_surface_type("gerstner_wave", "cpu", cls=GerstnerWaveSurface)