Auxiliary Functions

def dist(lat, lon, units="km"):
    """
    Distance between geographic positions and phase angle.
    
    Calculates great circle distance and initial bearing between
    geographic coordinates using spherical Earth approximation.
    
    Parameters:
    - lat: array_like, latitude coordinates [decimal degrees]
    - lon: array_like, longitude coordinates [decimal degrees]  
    - units: str, distance units ("km" or "nm"), default "km"
    
    Returns:
    - tuple: (distance, phase_angle)
        - distance: array_like, great circle distance [km or nm]
        - phase_angle: array_like, initial bearing [degrees]
    """

def f(lat):
    """
    Coriolis parameter (Coriolis frequency).
    
    Calculates the Coriolis parameter which represents the component
    of Earth's rotation affecting horizontal motion at given latitude.
    
    Parameters:
    - lat: array_like, latitude [decimal degrees]
    
    Returns:
    - array_like: Coriolis parameter [s⁻¹]
    """

def satO2(s, t):
    """
    Oxygen solubility in seawater.
    
    Calculates equilibrium dissolved oxygen concentration in seawater
    at atmospheric pressure using Weiss (1970) equations.
    
    Parameters:
    - s: array_like, salinity [psu (PSS-78)]
    - t: array_like, temperature [°C (ITS-90)]
    
    Returns:
    - array_like: oxygen solubility [ml/l]
    """

def satN2(s, t):
    """
    Nitrogen solubility in seawater.
    
    Calculates equilibrium dissolved nitrogen concentration in seawater
    at atmospheric pressure.
    
    Parameters:
    - s: array_like, salinity [psu (PSS-78)]
    - t: array_like, temperature [°C (ITS-90)]
    
    Returns:
    - array_like: nitrogen solubility [ml/l]
    """

def satAr(s, t):
    """
    Argon solubility in seawater.
    
    Calculates equilibrium dissolved argon concentration in seawater
    at atmospheric pressure.
    
    Parameters:
    - s: array_like, salinity [psu (PSS-78)]
    - t: array_like, temperature [°C (ITS-90)]
    
    Returns:
    - array_like: argon solubility [ml/l]
    """

def swvel(length, depth):
    """
    Surface wave velocity using shallow/deep water approximations.
    
    Calculates phase velocity of surface gravity waves based on
    wavelength and water depth using linear wave theory.
    
    Parameters:
    - length: array_like, wavelength [m]
    - depth: array_like, water depth [m]
    
    Returns:
    - array_like: wave phase velocity [m/s]
    """

import seawater as sw
import numpy as np

# Define oceanographic station positions
# Station coordinates (lat, lon in decimal degrees)
stations = {
    'A': (35.0, -75.0),   # Off Cape Hatteras
    'B': (40.0, -70.0),   # Georges Bank
    'C': (45.0, -65.0),   # Nova Scotia shelf
}

# Calculate distances between all station pairs
station_names = list(stations.keys())
coords = list(stations.values())

print("Distance matrix between stations:")
print("     ", end="")
for name in station_names:
    print(f"{name:>8}", end="")
print()

for i, (name1, (lat1, lon1)) in enumerate(stations.items()):
    print(f"{name1:>4} ", end="")
    for j, (name2, (lat2, lon2)) in enumerate(stations.items()):
        if i <= j:
            # Calculate distance
            lats = np.array([lat1, lat2])
            lons = np.array([lon1, lon2])
            distance, bearing = sw.dist(lats, lons, units="km")
            
            if i == j:
                print(f"{'0':>8}", end="")
            else:
                print(f"{distance[1]:.0f}km", end="")
        else:
            print(f"{'':>8}", end="")
    print()

# Calculate bearing from Station A to Station B
lat_ab = np.array([35.0, 40.0])
lon_ab = np.array([-75.0, -70.0])
dist_ab, bearing_ab = sw.dist(lat_ab, lon_ab)
print(f"\\nA→B: {dist_ab[1]:.1f} km at bearing {bearing_ab[1]:.1f}°")

import seawater as sw
import numpy as np
import matplotlib.pyplot as plt

# Calculate Coriolis parameter across latitudes
latitudes = np.linspace(-90, 90, 181)
coriolis = sw.f(latitudes)

# Key oceanographic latitudes
key_lats = [-60, -30, 0, 30, 60]  # Southern Ocean, Subtropics, Equator, etc.
key_f = sw.f(key_lats)

print("Coriolis parameter at key latitudes:")
lat_names = ["60°S", "30°S", "Equator", "30°N", "60°N"]
for name, lat, f_val in zip(lat_names, key_lats, key_f):
    print(f"{name:>8}: f = {f_val:.2e} s⁻¹")

# Rossby radius calculation example (first baroclinic mode)
# Typical values for mid-latitude ocean
N = 5e-3  # Brünt-Väisälä frequency [rad/s]
H = 1000  # depth scale [m]
f_30N = sw.f(30)  # Coriolis at 30°N

rossby_radius = N * H / abs(f_30N) / 1000  # convert to km
print(f"\\nTypical Rossby radius at 30°N: {rossby_radius:.0f} km")

# Plot Coriolis parameter
plt.figure(figsize=(10, 6))
plt.plot(latitudes, coriolis * 1e4, 'b-', linewidth=2)
plt.axhline(0, color='k', linestyle='--', alpha=0.5)
plt.xlabel('Latitude (°)')
plt.ylabel('Coriolis Parameter (×10⁻⁴ s⁻¹)')
plt.title('Coriolis Parameter vs Latitude')
plt.grid(True, alpha=0.3)
plt.xlim(-90, 90)

# Mark key latitudes
for lat, f_val, name in zip(key_lats, key_f, lat_names):
    plt.plot(lat, f_val * 1e4, 'ro', markersize=8)
    plt.annotate(name, (lat, f_val * 1e4), xytext=(5, 5), 
                textcoords='offset points', fontsize=9)

plt.tight_layout()
plt.show()

import seawater as sw
import numpy as np
import matplotlib.pyplot as plt

# Temperature and salinity ranges for gas solubility analysis
temperatures = np.linspace(0, 30, 31)  # 0-30°C
salinities = [0, 20, 35]  # Fresh, brackish, seawater

plt.figure(figsize=(15, 5))

# Plot oxygen solubility
plt.subplot(1, 3, 1)
for s in salinities:
    o2_sol = sw.satO2(s, temperatures)
    plt.plot(temperatures, o2_sol, label=f'S = {s} psu', linewidth=2)
plt.xlabel('Temperature (°C)')
plt.ylabel('O₂ Solubility (ml/l)')
plt.title('Oxygen Solubility')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot nitrogen solubility
plt.subplot(1, 3, 2)
for s in salinities:
    n2_sol = sw.satN2(s, temperatures)
    plt.plot(temperatures, n2_sol, label=f'S = {s} psu', linewidth=2)
plt.xlabel('Temperature (°C)')
plt.ylabel('N₂ Solubility (ml/l)')
plt.title('Nitrogen Solubility')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot argon solubility
plt.subplot(1, 3, 3)
for s in salinities:
    ar_sol = sw.satAr(s, temperatures)
    plt.plot(temperatures, ar_sol, label=f'S = {s} psu', linewidth=2)
plt.xlabel('Temperature (°C)')
plt.ylabel('Ar Solubility (ml/l)')
plt.title('Argon Solubility')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Calculate gas ratios at standard conditions
std_temp, std_sal = 20.0, 35.0  # Standard seawater conditions
o2_std = sw.satO2(std_sal, std_temp)
n2_std = sw.satN2(std_sal, std_temp)
ar_std = sw.satAr(std_sal, std_temp)

print(f"Gas solubilities at {std_temp}°C, {std_sal} psu:")
print(f"O₂: {o2_std:.2f} ml/l")
print(f"N₂: {n2_std:.2f} ml/l")
print(f"Ar: {ar_std:.3f} ml/l")
print(f"N₂/Ar ratio: {n2_std/ar_std:.1f}")

import seawater as sw
import numpy as np
import matplotlib.pyplot as plt

# Wave analysis for different water depths and wavelengths
wavelengths = np.logspace(0, 3, 100)  # 1 m to 1000 m
depths = [5, 20, 100, 1000]  # Shallow to deep water

plt.figure(figsize=(12, 8))

# Plot wave velocity vs wavelength for different depths
plt.subplot(2, 2, 1)
for depth in depths:
    velocities = sw.swvel(wavelengths, depth)
    plt.loglog(wavelengths, velocities, label=f'h = {depth} m', linewidth=2)

# Add deep water limit (c = sqrt(g*L/(2*pi)))
g = 9.81
c_deep = np.sqrt(g * wavelengths / (2 * np.pi))
plt.loglog(wavelengths, c_deep, 'k--', label='Deep water limit', alpha=0.7)

plt.xlabel('Wavelength (m)')
plt.ylabel('Phase Velocity (m/s)')
plt.title('Wave Velocity vs Wavelength')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot wave velocity vs depth for fixed wavelengths
plt.subplot(2, 2, 2)
depths_fine = np.logspace(0, 3, 100)  # 1 m to 1000 m depth
wave_lengths = [10, 50, 200]  # Different wavelengths

for wl in wave_lengths:
    velocities = sw.swvel(wl, depths_fine)
    plt.semilogx(depths_fine, velocities, label=f'λ = {wl} m', linewidth=2)

plt.xlabel('Water Depth (m)')
plt.ylabel('Phase Velocity (m/s)')
plt.title('Wave Velocity vs Water Depth')
plt.legend()
plt.grid(True, alpha=0.3)

# Calculate wave parameters for typical wind waves
plt.subplot(2, 2, 3)
periods = np.linspace(5, 20, 16)  # Wave periods 5-20 seconds
g = 9.81

# Deep water wavelength: L = g*T²/(2π)
deep_wavelengths = g * periods**2 / (2 * np.pi)
deep_velocities = g * periods / (2 * np.pi)

plt.plot(periods, deep_wavelengths, 'b-', linewidth=2, label='Wavelength')
plt.plot(periods, deep_velocities, 'r-', linewidth=2, label='Phase velocity')
plt.xlabel('Wave Period (s)')
plt.ylabel('Wavelength (m) / Velocity (m/s)')
plt.title('Deep Water Wave Characteristics')
plt.legend()
plt.grid(True, alpha=0.3)

# Calculate relative depth parameter (kh = 2πh/L)
plt.subplot(2, 2, 4)
h_values = [10, 50, 200]  # Water depths
for h in h_values:
    kh = 2 * np.pi * h / deep_wavelengths
    plt.plot(periods, kh, label=f'h = {h} m', linewidth=2)

plt.axhline(np.pi, color='k', linestyle='--', alpha=0.7, label='Deep water (kh > π)')
plt.axhline(np.pi/10, color='k', linestyle=':', alpha=0.7, label='Shallow water (kh < π/10)')
plt.xlabel('Wave Period (s)')
plt.ylabel('Relative Depth (kh)')
plt.title('Wave Dispersion Parameter')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Example: Tsunami wave speed calculation
print("\\nTsunami wave analysis:")
tsunami_wavelength = 200000  # 200 km typical tsunami wavelength
ocean_depths = [1000, 3000, 5000]  # Typical ocean depths

for depth in ocean_depths:
    tsunami_speed = sw.swvel(tsunami_wavelength, depth)
    print(f"Depth {depth:4d} m: Tsunami speed = {tsunami_speed:.0f} m/s ({tsunami_speed*3.6:.0f} km/h)")