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 :math:`(x_0, y_0)` to:

.. math::

   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)

where:

- :math:`A` = wave amplitude
- :math:`\vec{k} = (k_x, k_y)` = wave vector
- :math:`k = |\vec{k}| = 2\pi/\lambda` = wavenumber
- :math:`\omega` = angular frequency
- :math:`\phi` = phase offset
- :math:`t` = time

Physical Parameters
~~~~~~~~~~~~~~~~~~~

**Wavelength:** :math:`\lambda` (meters)

Typical ocean wavelengths:

- Capillary waves: < 1 cm
- Wind waves: 1-100 m
- Swell: 100-1000 m

**Wave steepness:** :math:`s = A k`

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

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

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

.. math::

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

where :math:`g = 9.81` m/s² is gravitational acceleration.

Multi-Component Waves
~~~~~~~~~~~~~~~~~~~~~

Realistic ocean surfaces combine multiple wave components:

.. math::

   \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 :math:`\vec{\Delta}_i` is a Gerstner wave component.

**Implementation:**

.. code-block:: python

   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:

.. math::

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

For Gerstner waves:

.. math::

   \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)

**Implementation:**

.. code-block:: python

   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:

.. math::

   x &= (R_E + h) \cos\phi \cos\lambda \\
   y &= (R_E + h) \cos\phi \sin\lambda \\
   z &= (R_E + h) \sin\phi

where:

- :math:`R_E = 6371` km = Earth radius
- :math:`\phi` = latitude
- :math:`\lambda` = longitude
- :math:`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:

.. math::

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

where :math:`\vec{r}` is the position vector from Earth center.

CurvedWaveSurface Class
~~~~~~~~~~~~~~~~~~~~~~~

.. code-block:: python

   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:**

.. math::

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

where:

- :math:`\alpha = 0.0081` (Phillips constant)
- :math:`\beta = 0.74`
- :math:`\omega_0` = peak frequency
- :math:`g` = 9.81 m/s²

**JONSWAP Spectrum** (for fetch-limited seas):

.. math::

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

where :math:`\gamma \approx 3.3` is peak enhancement factor.

Directional Spread
~~~~~~~~~~~~~~~~~~

Waves propagate in various directions with angular distribution:

.. math::

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

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

Sampling Wave Components
~~~~~~~~~~~~~~~~~~~~~~~~

To generate realistic waves:

1. **Sample frequencies** from spectrum :math:`S(\omega)`
2. **Sample directions** from directional distribution :math:`D(\theta)`
3. **Compute amplitudes:**

   .. math::

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

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

**Example:**

.. code-block:: python

   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:

.. math::

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

Phase velocity:

.. math::

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

Group velocity:

.. math::

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

Computational Considerations
----------------------------

Memory Requirements
~~~~~~~~~~~~~~~~~~~

For :math:`N` wave components:

- **Per-wave storage:** 32 bytes (A, kx, ky, omega, phi as float32)
- **Total:** :math:`32N` bytes
- **Typical:** 20 components = 640 bytes (negligible)

Computation Cost
~~~~~~~~~~~~~~~~

**Surface evaluation:**

- Per point: :math:`O(N)` trigonometric evaluations
- With :math:`M` rays: :math:`O(MN)` total
- **Optimization:** Vectorize with NumPy/CUDA

**GPU acceleration:**

.. code-block:: cuda

   __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