wave-theory

What this skill does

# Wave Theory SME Skill

Comprehensive ocean wave theory expertise including wave mechanics, spectral analysis, wave statistics, and irregular sea modeling for offshore engineering applications.

## When to Use This Skill

Use wave theory knowledge when:
- **Wave spectra** - JONSWAP, Pierson-Moskowitz, scatter diagrams
- **Wave statistics** - Significant wave height, spectral parameters
- **Irregular seas** - Generate time series from spectra
- **Wave kinematics** - Particle velocities and accelerations
- **Wave transformation** - Shoaling, refraction, diffraction
- **Extreme values** - Design wave estimation

## Core Knowledge Areas

### 1. Regular Wave Theory

**Linear (Airy) Wave Theory:**
```python
import numpy as np

def airy_wave_properties(
    H: float,
    T: float,
    d: float,
    g: float = 9.81
) -> dict:
    """
    Calculate Airy wave properties.

    Valid for: H/L < 0.14, d/L > 0.5 (deep water) or d/L < 0.05 (shallow)

    Args:
        H: Wave height (m)
        T: Wave period (s)
        d: Water depth (m)
        g: Gravity (m/s²)

    Returns:
        Wave properties dictionary
    """
    # Wave frequency
    omega = 2 * np.pi / T

    # Dispersion relation: ω² = gk·tanh(kd)
    # Solve iteratively for wave number k
    from scipy.optimize import fsolve

    def dispersion(k):
        return omega**2 - g * k * np.tanh(k * d)

    k0 = omega**2 / g  # Deep water approximation
    k = fsolve(dispersion, k0)[0]

    # Wave length
    L = 2 * np.pi / k

    # Wave celerity (phase speed)
    C = omega / k

    # Group velocity
    n = 0.5 * (1 + 2*k*d / np.sinh(2*k*d))  # Shoaling coefficient
    Cg = n * C

    # Deep water classification
    if d / L > 0.5:
        regime = "Deep water"
    elif d / L < 0.05:
        regime = "Shallow water"
    else:
        regime = "Intermediate"

    return {
        'height_m': H,
        'period_s': T,
        'depth_m': d,
        'wavelength_m': L,
        'wave_number': k,
        'frequency_rad_s': omega,
        'frequency_hz': omega / (2*np.pi),
        'celerity_m_s': C,
        'group_velocity_m_s': Cg,
        'regime': regime,
        'd_over_L': d / L,
        'H_over_L': H / L,
        'steepness': H / L
    }

# Example
wave = airy_wave_properties(H=8, T=12, d=1500)

print(f"Wave Properties (H={wave['height_m']}m, T={wave['period_s']}s):")
print(f"  Wavelength: {wave['wavelength_m']:.1f} m")
print(f"  Regime: {wave['regime']} (d/L = {wave['d_over_L']:.3f})")
print(f"  Celerity: {wave['celerity_m_s']:.2f} m/s")
print(f"  Steepness: {wave['steepness']:.4f}")
```

**Wave Kinematics:**
```python
def wave_particle_kinematics(
    z: float,
    H: float,
    T: float,
    d: float,
    t: float = 0,
    x: float = 0,
    g: float = 9.81
) -> dict:
    """
    Calculate wave particle velocities and accelerations.

    Args:
        z: Vertical position (0 at SWL, negative below)
        H: Wave height (m)
        T: Wave period (s)
        d: Water depth (m)
        t: Time (s)
        x: Horizontal position (m)
        g: Gravity (m/s²)

    Returns:
        Particle kinematics
    """
    # Wave properties
    wave = airy_wave_properties(H, T, d, g)
    k = wave['wave_number']
    omega = wave['frequency_rad_s']

    # Amplitude
    a = H / 2

    # Hyperbolic functions
    cosh_kz_d = np.cosh(k * (z + d))
    sinh_kz_d = np.sinh(k * (z + d))
    cosh_kd = np.cosh(k * d)
    sinh_kd = np.sinh(k * d)

    # Wave phase
    phase = k * x - omega * t

    # Horizontal velocity
    u = (omega * a * cosh_kz_d / sinh_kd) * np.cos(phase)

    # Vertical velocity
    w = (omega * a * sinh_kz_d / sinh_kd) * np.sin(phase)

    # Horizontal acceleration
    ax = -(omega**2 * a * cosh_kz_d / sinh_kd) * np.sin(phase)

    # Vertical acceleration
    az = (omega**2 * a * sinh_kz_d / sinh_kd) * np.cos(phase)

    # Dynamic pressure
    p_dynamic = g * a * (cosh_kz_d / cosh_kd) * np.cos(phase)

    return {
        'horizontal_velocity': u,
        'vertical_velocity': w,
        'horizontal_acceleration': ax,
        'vertical_acceleration': az,
        'dynamic_pressure': p_dynamic,
        'total_velocity': np.sqrt(u**2 + w**2),
        'total_acceleration': np.sqrt(ax**2 + az**2)
    }

# Example: Surface velocity (z=0)
kinematics = wave_particle_kinematics(z=0, H=8, T=12, d=1500, t=0, x=0)

print(f"Surface Particle Kinematics:")
print(f"  Horizontal velocity: {kinematics['horizontal_velocity']:.2f} m/s")
print(f"  Vertical velocity: {kinematics['vertical_velocity']:.2f} m/s")
print(f"  Total velocity: {kinematics['total_velocity']:.2f} m/s")
```

### 2. Wave Spectra

**JONSWAP Spectrum:**
```python
def jonswap_spectrum(
    frequencies: np.ndarray,
    Hs: float,
    Tp: float,
    gamma: float = 3.3,
    alpha: float = None
) -> np.ndarray:
    """
    Calculate JONSWAP wave spectrum.

    S(f) = α g² (2π)^-4 f^-5 exp[-5/4(f/fp)^-4] γ^exp[-(f-fp)²/(2σ²fp²)]

    Args:
        frequencies: Frequency array (Hz)
        Hs: Significant wave height (m)
        Tp: Peak period (s)
        gamma: Peak enhancement factor (3.3 for North Sea)
        alpha: Phillips constant (calculated if None)

    Returns:
        Spectral density S(f) (m²/Hz)
    """
    g = 9.81
    fp = 1 / Tp  # Peak frequency (Hz)

    # Calculate alpha if not provided
    if alpha is None:
        # Relationship: Hs = 4*sqrt(m0)
        # For JONSWAP: alpha ≈ 5.061 * Hs² / Tp⁴ * (1 - 0.287*ln(γ))
        alpha = 5.061 * Hs**2 / Tp**4 * (1 - 0.287 * np.log(gamma))

    # Sigma parameter
    sigma = np.where(frequencies <= fp, 0.07, 0.09)

    # Pierson-Moskowitz spectrum
    S_PM = alpha * g**2 * (2*np.pi)**(-4) * frequencies**(-5) * \
           np.exp(-1.25 * (frequencies / fp)**(-4))

    # Peak enhancement
    r = np.exp(-(frequencies - fp)**2 / (2 * sigma**2 * fp**2))
    gamma_factor = gamma ** r

    # JONSWAP spectrum
    S = S_PM * gamma_factor

    return S

# Example: Generate JONSWAP spectrum
freq = np.linspace(0.01, 0.5, 500)
S = jonswap_spectrum(freq, Hs=8.5, Tp=12.0, gamma=3.3)

# Verify Hs
m0 = np.trapz(S, freq)
Hs_calc = 4 * np.sqrt(m0)

print(f"JONSWAP Spectrum:")
print(f"  Input Hs: 8.5 m")
print(f"  Calculated Hs: {Hs_calc:.2f} m")
print(f"  Peak frequency: {1/12:.4f} Hz")
```

**Pierson-Moskowitz Spectrum:**
```python
def pierson_moskowitz_spectrum(
    frequencies: np.ndarray,
    Hs: float,
    Tp: float = None,
    U19_5: float = None
) -> np.ndarray:
    """
    Calculate Pierson-Moskowitz spectrum (fully developed sea).

    Args:
        frequencies: Frequency array (Hz)
        Hs: Significant wave height (m)
        Tp: Peak period (s) - optional
        U19_5: Wind speed at 19.5m height (m/s) - optional

    Returns:
        Spectral density S(f) (m²/Hz)
    """
    g = 9.81

    if Tp is not None:
        # Use peak period
        fp = 1 / Tp
    elif U19_5 is not None:
        # Calculate from wind speed
        fp = 0.877 * g / (2 * np.pi * U19_5)
    else:
        raise ValueError("Must provide either Tp or U19_5")

    # Phillips constant
    alpha = 0.0081  # For fully developed seas

    # P-M spectrum
    S = alpha * g**2 * (2*np.pi)**(-4) * frequencies**(-5) * \
        np.exp(-1.25 * (frequencies / fp)**(-4))

    return S

# Example
S_PM = pierson_moskowitz_spectrum(freq, Hs=8.5, Tp=12.0)

m0_PM = np.trapz(S_PM, freq)
Hs_PM = 4 * np.sqrt(m0_PM)

print(f"P-M Spectrum Hs: {Hs_PM:.2f} m")
```

### 3. Wave Statistics

**Spectral Parameters:**
```python
def calculate_spectral_parameters(
    S: np.ndarray,
    frequencies: np.ndarray
) -> dict:
    """
    Calculate spectral wave parameters.

    Args:
        S: Wave spectrum (m²/Hz)
        frequencies: Frequency array (Hz)

    Returns:
        Spectral parameters
    """
    # Spectral moments
    m0 = np.trapz(S, frequencies)
    m1 = np.trapz(S * frequencies, frequencies)
    m2 = np.trapz(S * frequencies**2, frequencies)
    m4 = np.trapz(S * frequencies**4, frequencies)

    # Significant wave height
    Hs = 4 * np.sqrt(m0)

    #