Ray Tracing Canonical Examples¶

# data handling
import numpy as np
import xarray as xr

# plotting
import matplotlib.pyplot as plt
from IPython.display import HTML
import cmocean

# ray tracing
import mantaray

# helper functions
import utils

Introduction¶

This notebook covers how to run Mantaray with the custom depth and current fields that were constructed in idealized_fields.ipynb. It is recommended that users run through this notebook first.

Initial Wave and Model Parameters¶

The initial wave parameters that we need to define for Mantaray are \(k_{x0}\)/\(k_{y_0}\), the initial x/y wavenumber. We will calculate these from the following parameters: - \(T_0\): The initial way period, - \(\theta_0\): The initial wave direction (convention: what direction the waves are going to), - \(n_{rays}\): The number of rays.

All waves will have the same x/y wavenumbers. We will use these initial parameters for the 4 canonical examples.

T0 = 10  # Period [s]
theta0 = 0  # Direction [rad]

# Convert period to wavenumber magnitude
k_0 = utils.period2wavenumber(T0)

# Calculate wavenumber components
k_x0 = k_0 * np.cos(theta0)
k_y0 = k_0 * np.sin(theta0)

# Number of rays
n_rays = 100

# Initialize wavenumber for all rays
K_x0 = k_x0 * np.ones(n_rays)
K_y0 = k_y0 * np.ones(n_rays)

Additionally, we need to define the initial positions for all of our rays. To do this, we will extract the grid from one of our canonical examples. All examples have the same grid, and the initial positions for all rays will be the same for each canonical example.

ds_jet = xr.open_dataset("./forcing/zonal_jet.nc")

# Read x and y from file to get domain size
x = ds_jet.x.values
y = ds_jet.y.values

x0 = 0 * np.ones(n_rays)
y0 = np.linspace(0, y.max(), n_rays)

Additionally, using the grid, we can compute the optimal timestep and duration for each model.

timestep = utils.compute_cfl(x, y, k_0)
duration = utils.compute_duration(x, k_0)

Interactive Animations¶

Any ray tracing example can be animated. For any example, you need to extract: - The x/y grid from the currents/bathymetry dataset - The background field from the currents/bathymetry dataset

X = ds_island.x
Y = ds_island.y
depth = ds_island.depth

Below, we create the frames of the animation.

island_anim = utils.animate_rays(
    X, Y, depth, bundle_island, style="bathymetry", ray_sample=2, time_sample=5
)

Then, we can plot the animation so that it is interactive. WARNING: this cell may take some time to run.

HTML(island_anim.to_jshtml())

Optionally, the animation can be saved.

# anim.save('./output/waves.gif')

Bonus Examples¶

Sloped Beach¶

For this example, we will start the rays from a line on an angle.

theta0 = np.deg2rad(35)

# Calculate wavenumber components
k_x0 = k_0 * np.cos(theta0)
k_y0 = k_0 * np.sin(theta0)

n_rays = 50

# Initialize wavenumber for all rays
Kx0 = k_x0 * np.ones(n_rays)
Ky0 = k_y0 * np.ones(n_rays)

The CFL and duration funtions assume that the rays are moving from the left side of the domain to the right, so they are not a perfect estimate here, but they’ll still work.

step_size = utils.compute_cfl(x, y, k_0)
duration = utils.compute_duration(x, k_0)

We have calculated the x/y wavenumber components for waves traveling at 35 degrees, but we now need to define their initial positions along an angled line.

# Parameters
x_start = 50_000  # Start x coordinate [m]
y_start = 0  # Start y coordinate [m]
line_length = 50_000  # Total length of the line [m]

# Calculate the step size in the x and y directions
dx = -line_length / n_rays * np.cos(np.pi / 2 - theta0)
dy = line_length / n_rays * np.sin(np.pi / 2 - np.pi / 4)

# Generate the coordinates along the line
x0 = x_start + np.linspace(0, dx * (n_rays - 1), n_rays)  # x coordinates
y0 = y_start + np.linspace(0, dy * (n_rays - 1), n_rays)  # y coordinates

beach_path = "./forcing/gentle_slope.nc"

bundle_beach = mantaray.ray_tracing(
    x0, y0, Kx0, Ky0, duration, 10, beach_path, static_cur_path
)

# bundle_beach.to_netcdf("./output/bundle4beach.nc", format="NETCDF3_CLASSIC")

ds_beach = xr.open_dataset(beach_path)

X = ds_beach.x
Y = ds_beach.y
H = ds_beach.depth

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

cf = plt.contourf(X, Y, H, cmap=cmocean.cm.deep, levels=40)
line = plt.plot(x0, y0, label="Line at points (x0, y0)", color="red")

plt.contour(X, Y, H, 0, colors="magenta")


for i in range(bundle_beach.ray.size)[::2]:
    ray = bundle_beach.isel(ray=i)
    plt.plot(ray.x, ray.y, "k", lw=0.78, color="white")

plt.xlabel("x [m]")
plt.ylabel("y [m]")

plt.colorbar(cf, label="depth [m]")
plt.title("Ray Tracing for Gentle Slope")
plt.show()

/var/folders/30/l7lkprrx2c32s7j5y82c67300000gq/T/ipykernel_90471/1316495833.py:15: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "k" (-> color=(0.0, 0.0, 0.0, 1)). The keyword argument will take precedence.
  plt.plot(ray.x, ray.y, "k", lw=0.78, color="white")

Ray Tracing Canonical Examples¶

Introduction¶

Initial Wave and Model Parameters¶

Ray Tracing with Mantaray¶

Zonal Jet¶

Mesoscale Eddy¶

Gaussian Island¶

Interactive Animations¶

Bonus Examples¶

Sloped Beach¶