Wave Surface Modeling

L-SURF implements sophisticated ocean wave surface models for realistic light-water interactions.

Gerstner Waves

Overview

Gerstner waves (also called trochoidal waves) are a more physically accurate model than simple sinusoidal waves. They produce sharper crests and flatter troughs, matching real ocean waves.

Mathematical Formulation

A single Gerstner wave component displaces a point \((x_0, y_0)\) to:

\[\begin{split}x &= x_0 - \frac{A k_x}{k} \sin(k_x x_0 + k_y y_0 - \omega t + \phi) \\ y &= y_0 - \frac{A k_y}{k} \sin(k_x x_0 + k_y y_0 - \omega t + \phi) \\ z &= A \cos(k_x x_0 + k_y y_0 - \omega t + \phi)\end{split}\]

where:

  • \(A\) = wave amplitude

  • \(\vec{k} = (k_x, k_y)\) = wave vector

  • \(k = |\vec{k}| = 2\pi/\lambda\) = wavenumber

  • \(\omega\) = angular frequency

  • \(\phi\) = phase offset

  • \(t\) = time

Physical Parameters

Wavelength: \(\lambda\) (meters)

Typical ocean wavelengths:

  • Capillary waves: < 1 cm

  • Wind waves: 1-100 m

  • Swell: 100-1000 m

Wave steepness: \(s = A k\)

Physical constraint: \(s < 1\) (typically \(s < 0.5\) for stable waves)

Dispersion relation: \(\omega^2 = gk \tanh(kh)\)

For deep water (\(h \to \infty\)):

\[\omega = \sqrt{gk} = \sqrt{\frac{2\pi g}{\lambda}}\]

where \(g = 9.81\) m/s² is gravitational acceleration.

Multi-Component Waves

Realistic ocean surfaces combine multiple wave components:

\[\vec{r}(x_0, y_0, t) = (x_0, y_0, 0) + \sum_{i=1}^{N} \vec{\Delta}_i(x_0, y_0, t)\]

where each \(\vec{\Delta}_i\) is a Gerstner wave component.

Implementation:

class GerstnerWaveSurface:
    def __init__(self, wave_components):
        """
        Parameters
        ----------
        wave_components : List[GerstnerWaveParams]
            List of wave parameters (A, k, omega, phi, direction)
        """
        self.waves = wave_components

    def position(self, x0, y0, t=0):
        """Compute displaced position."""
        x, y, z = x0, y0, 0.0
        for wave in self.waves:
            A, k, omega, phi = wave.amplitude, wave.k, wave.omega, wave.phase
            kx, ky = wave.k_vector

            phase = kx * x0 + ky * y0 - omega * t + phi
            sin_p = sin(phase)
            cos_p = cos(phase)

            x -= (A * kx / k) * sin_p
            y -= (A * ky / k) * sin_p
            z += A * cos_p

        return x, y, z

Surface Normal

The normal vector is computed from the Jacobian:

\[\vec{n} = \frac{\partial\vec{r}}{\partial x_0} \times \frac{\partial\vec{r}}{\partial y_0}\]

For Gerstner waves:

\[\begin{split}\frac{\partial x}{\partial x_0} &= 1 - \sum_i A_i k_{x,i} \cos(\ldots) \\ \frac{\partial y}{\partial x_0} &= -\sum_i A_i k_{y,i} \cos(\ldots) \\ \frac{\partial z}{\partial x_0} &= -\sum_i A_i k_{x,i} \sin(\ldots)\end{split}\]

Implementation:

def normal(self, x0, y0, t=0):
    """Compute surface normal at (x0, y0)."""
    # Initialize with flat surface Jacobian
    dx_dx0 = 1.0
    dy_dx0 = 0.0
    dz_dx0 = 0.0
    dx_dy0 = 0.0
    dy_dy0 = 1.0
    dz_dy0 = 0.0

    for wave in self.waves:
        A, kx, ky, k, omega, phi = ...
        phase = kx * x0 + ky * y0 - omega * t + phi
        cos_p = cos(phase)
        sin_p = sin(phase)

        # Accumulate derivatives
        dx_dx0 -= A * kx * cos_p
        dy_dx0 -= A * ky * cos_p
        dz_dx0 -= A * kx * sin_p

        dx_dy0 -= A * kx * cos_p
        dy_dy0 -= A * ky * cos_p
        dz_dy0 -= A * ky * sin_p

    # Cross product for normal
    tangent_x = (dx_dx0, dy_dx0, dz_dx0)
    tangent_y = (dx_dy0, dy_dy0, dz_dy0)
    normal = cross(tangent_x, tangent_y)
    return normalize(normal)

Curved Earth Geometry

For large-scale simulations (> 10 km), Earth’s curvature becomes significant.

Spherical Coordinate System

L-SURF uses a local spherical coordinate system centered on Earth:

\[\begin{split}x &= (R_E + h) \cos\phi \cos\lambda \\ y &= (R_E + h) \cos\phi \sin\lambda \\ z &= (R_E + h) \sin\phi\end{split}\]

where:

  • \(R_E = 6371\) km = Earth radius

  • \(\phi\) = latitude

  • \(\lambda\) = longitude

  • \(h\) = height above reference ellipsoid

Local Tangent Plane

For computational efficiency, waves are computed in a local tangent plane and then mapped to the curved surface:

  1. Define local patch: Small region (~100 km × 100 km)

  2. Compute waves in 2D: Use Gerstner waves in local (x, y) coordinates

  3. Map to sphere: Apply spherical transformation

  4. Adjust normals: Account for curvature

Curvature correction:

The true surface normal includes Earth curvature:

\[\vec{n}_{\text{true}} = \vec{n}_{\text{local}} + \frac{\vec{r}}{R_E}\]

where \(\vec{r}\) is the position vector from Earth center.

CurvedWaveSurface Class

class CurvedWaveSurface:
    """Wave surface with Earth curvature."""

    def __init__(self, wave_params, earth_radius=6371e3):
        self.waves = wave_params
        self.R_E = earth_radius

    def position(self, x, y, t=0):
        """Position accounting for curvature."""
        # Compute local wave displacement
        dx, dy, dz = self._local_wave_displacement(x, y, t)

        # Map to sphere
        r = sqrt(x**2 + y**2)
        if r > 0:
            # Angle from center
            theta = r / self.R_E

            # Radial direction
            cos_theta = cos(theta)
            sin_theta = sin(theta)

            # Curved position
            x_curved = x * (self.R_E + dz) / r * sin_theta + dx * cos_theta
            y_curved = y * (self.R_E + dz) / r * sin_theta + dy * cos_theta
            z_curved = (self.R_E + dz) * cos_theta - r * sin_theta
        else:
            x_curved, y_curved, z_curved = dx, dy, self.R_E + dz

        return x_curved, y_curved, z_curved

Wave Spectra

Realistic Ocean Spectra

Real ocean waves follow statistical distributions. Common models:

Pierson-Moskowitz Spectrum:

\[S(\omega) = \frac{\alpha g^2}{\omega^5} \exp\left[-\beta\left(\frac{\omega_0}{\omega}\right)^4\right]\]

where:

  • \(\alpha = 0.0081\) (Phillips constant)

  • \(\beta = 0.74\)

  • \(\omega_0\) = peak frequency

  • \(g\) = 9.81 m/s²

JONSWAP Spectrum (for fetch-limited seas):

\[S(\omega) = S_{PM}(\omega) \cdot \gamma^{\exp[-(\omega-\omega_0)^2/(2\sigma^2\omega_0^2)]}\]

where \(\gamma \approx 3.3\) is peak enhancement factor.

Directional Spread

Waves propagate in various directions with angular distribution:

\[D(\theta) = \frac{2}{\pi} \cos^{2s}(\theta - \theta_0)\]

where \(s\) controls spreading (typical: \(s = 2-10\)).

Sampling Wave Components

To generate realistic waves:

  1. Sample frequencies from spectrum \(S(\omega)\)

  2. Sample directions from directional distribution \(D(\theta)\)

  3. Compute amplitudes:

    \[A_i = \sqrt{2 S(\omega_i) \Delta\omega}\]

  4. Random phases: \(\phi_i \sim U(0, 2\pi)\)

Example:

def create_ocean_wave_surface(wind_speed, fetch, num_components=20):
    """Create realistic ocean wave surface."""
    # Compute spectrum parameters
    omega_peak = 0.855 * 9.81 / wind_speed

    # Sample frequencies (logarithmically spaced)
    omega_min = 0.1 * omega_peak
    omega_max = 10 * omega_peak
    omegas = np.logspace(log10(omega_min), log10(omega_max), num_components)

    # Compute amplitudes from JONSWAP spectrum
    spectrum = jonswap_spectrum(omegas, wind_speed, fetch)
    delta_omega = omegas[1] - omegas[0]
    amplitudes = np.sqrt(2 * spectrum * delta_omega)

    # Sample directions
    wind_dir = 0.0  # degrees
    spreading = 5
    directions = sample_directional_spread(num_components, wind_dir, spreading)

    # Random phases
    phases = np.random.uniform(0, 2*np.pi, num_components)

    # Create wave components
    waves = []
    for A, omega, theta, phi in zip(amplitudes, omegas, directions, phases):
        k = omega**2 / 9.81  # Deep water dispersion
        kx = k * cos(theta)
        ky = k * sin(theta)
        waves.append(GerstnerWaveParams(A, kx, ky, omega, phi))

    return GerstnerWaveSurface(waves)

Time Evolution

Static vs Dynamic Waves

L-SURF supports both modes:

Static (t=0):

  • Faster computation

  • Suitable for snapshot simulations

  • Good for parameter studies

Dynamic (t variable):

  • Realistic wave motion

  • Required for time-dependent effects

  • Animation capabilities

Wave Propagation

Waves evolve according to dispersion relation:

\[\vec{r}(x_0, y_0, t) = \vec{r}_0 + \sum_i \vec{\Delta}_i(x_0, y_0, t)\]

Phase velocity:

\[v_p = \frac{\omega}{k} = \sqrt{\frac{g}{k}} = \sqrt{\frac{g\lambda}{2\pi}}\]

Group velocity:

\[v_g = \frac{d\omega}{dk} = \frac{1}{2}v_p\]

Computational Considerations

Memory Requirements

For \(N\) wave components:

  • Per-wave storage: 32 bytes (A, kx, ky, omega, phi as float32)

  • Total: \(32N\) bytes

  • Typical: 20 components = 640 bytes (negligible)

Computation Cost

Surface evaluation:

  • Per point: \(O(N)\) trigonometric evaluations

  • With \(M\) rays: \(O(MN)\) total

  • Optimization: Vectorize with NumPy/CUDA

GPU acceleration:

__global__ void evaluate_wave_height(float* x, float* y, float* z,
                                     WaveParams* waves, int N_waves,
                                     int N_points) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N_points) {
        float x0 = x[idx];
        float y0 = y[idx];
        float height = 0.0f;

        for (int i = 0; i < N_waves; i++) {
            float phase = waves[i].kx * x0 + waves[i].ky * y0 + waves[i].phi;
            height += waves[i].A * cosf(phase);
        }

        z[idx] = height;
    }
}

Achieves ~1000x speedup on modern GPUs for large ray counts.

Validation

The wave implementation is validated against:

  1. Analytical solutions: Single Gerstner wave profiles

  2. Oceanographic data: Statistical wave height distributions

  3. Energy conservation: Total wave energy vs. spectrum integral

  4. Geometric constraints: Steepness limits, no surface self-intersection