Intermediate cookbook

This page contains several examples of Capytaine’s features, mostly focused on post-processing features such as reconstruction of the free surface elevation or computation of the Kochin functions. The scripts can be downloaded individually as Python files from the Github folder with examples for version 3.1.dev0.

B1. Plot pressure on hull

This example solves a BEM problem and display the computed pressure on the hull on a 3D view of the mesh.

import numpy as np
import capytaine as cpt
import matplotlib.pyplot as plt

mesh = cpt.mesh_sphere(radius=1.0, faces_max_radius=0.1).immersed_part()
body = cpt.FloatingBody(
    mesh=mesh,
    lid_mesh=mesh.generate_lid(),
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0))
)

pb = cpt.RadiationProblem(wavelength=1.0, body=body, radiating_dof="Surge", water_depth=np.inf)
solver = cpt.BEMSolver()
res = solver.solve(pb, keep_details=True)

body.mesh.show(
        backend="matplotlib",
        color_field=np.real(res.pressure_on_hull),
        cmap=plt.get_cmap("viridis"),  # Colormap
        )

# # Plotting also the (not physically meaningful) pressure on the lid
# body.mesh_including_lid.show(
#         backend="matplotlib",
#         color_field=np.real(res.pressure),
#         cmap=plt.get_cmap("viridis"),  # Colormap
#         )

B2. Haskind’s relation

This example computes the excitation force from the radiation potential using Haskind’s relation. The result is compared with the one obtained by direct integration of the potentials from incident waves and from the diffraction problem.

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_potential, airy_waves_velocity, froude_krylov_force

r = 1.
draft = 0.5
depth = 10.
omega = 1.
rho = 1000

# Define geometry and heave degree of freedom
body = cpt.FloatingBody(
        mesh=cpt.mesh_vertical_cylinder(
            length=2*draft, radius=r, center=(0.,0.,0.),
            resolution=(10, 50, 10)
            ))
body.add_translation_dof(name='Heave')
body = body.immersed_part()

solver = cpt.BEMSolver()

# Define and solve the diffraction and radiation problems
diff_problem = cpt.DiffractionProblem(body=body, water_depth=depth,
                                      omega=omega, wave_direction=0.)
rad_problem = cpt.RadiationProblem(body=body, water_depth=depth,
                                   omega=omega, radiating_dof='Heave')

diff_solution = solver.solve(diff_problem)
rad_solution = solver.solve(rad_problem)


# Read mesh properties
faces_centers = body.mesh.faces_centers
faces_normals = body.mesh.faces_normals
faces_areas = body.mesh.faces_areas

# Computation from the diffraction solution (Capytaine)
FK = froude_krylov_force(diff_problem)['Heave']
diff_force = diff_solution.forces['Heave']

# Get potentials
phi_inc = airy_waves_potential(faces_centers, diff_problem)
v_inc = airy_waves_velocity(faces_centers, diff_problem)
phi_rad = rad_solution.potential

# Indirect computation from the radiation solution, via the Haskind relation
integrand = - (phi_inc * faces_normals[:,2]
              - 1j/omega * phi_rad * np.diag(v_inc@faces_normals.T))
F_hask = 1j * omega * rho * np.sum(integrand*faces_areas)

# Direct integration of the potentials
diff_force_recomputed = - 1j * omega * rho * np.sum(
                        diff_solution.potential*faces_areas*faces_normals[:,2])
FK_recomputed = - 1j * omega * rho * np.sum(
                        phi_inc*faces_areas*faces_normals[:,2])

print('Result from Capytaine = ', FK + diff_force)
print('Result from recomputed direct integration',
        FK_recomputed + diff_force_recomputed)
print('Result from Haskind relation = ', F_hask)

B3. Free surface elevation

This example computes the free surface elevation as a post-processing step of a diffraction problem.

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
import matplotlib.pyplot as plt

mesh = cpt.mesh_sphere(radius=1, resolution=(10, 10)).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs())

pb = cpt.DiffractionProblem(body=body, wave_direction=np.pi/2, wavelength=3.0)
solver = cpt.BEMSolver()
res = solver.solve(pb)

grid = np.meshgrid(np.linspace(-10.0, 10.0, 100), np.linspace(-10.0, 10.0, 100))
fse = solver.compute_free_surface_elevation(grid, res)
incoming_fse = airy_waves_free_surface_elevation(grid, res)

fig, ax = plt.subplots(1, 1)
ax.pcolormesh(grid[0], grid[1], np.real(fse + incoming_fse))
ax.add_patch(plt.Circle((0, 0), radius=1, color="black"))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect("equal")
ax.set_title("Wave crests and troughs at t=0 seen from above")
plt.show()

B4. Kochin function

This example computes the Kochin function for a surging buoy and plot the results.

import numpy as np
from numpy import pi
import xarray as xr
import matplotlib.pyplot as plt

import capytaine as cpt

cpt.set_logging('INFO')

# SET UP THE MESH
buoy = cpt.FloatingBody(
        mesh=cpt.mesh_vertical_cylinder(
            radius=1.0, length=3.0, center=(0, 0, -1.0),
            resolution=(15, 15, 3)
            ))
buoy.add_translation_dof(name="Surge")

# SOLVE THE BEM PROBLEMS AND COMPUTE THE KOCHIN FUNCTIONS
theta_range = np.linspace(0, 2*pi, 40)
omega_range = np.linspace(1.0, 5.0, 3)

test_matrix = xr.Dataset(coords={
    'omega': omega_range, 'theta': theta_range, 'radiating_dof': ["Surge"],
})
solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, buoy.immersed_part(), hydrostatics=False)

# PLOT THE KOCHIN FUNCTION
plt.figure()
dataset = dataset.swap_dims({'omega': 'wavenumber'})  # Use wavenumber to index data
for k in dataset.coords['wavenumber']:
    (dataset['kochin_radiation']
     .sel(radiating_dof="Surge", wavenumber=k)
     .real
     .plot(x='theta', label=f"k={float(k):.2f} m¯¹", ax=plt.gca()))
plt.legend()
plt.title("Real part of the Kochin function as a function of theta")
plt.tight_layout()
plt.show()

B5. Plot velocity in domain

This example reconstruct the velocity field in the domain and display it using matplotlib.

import numpy as np
import capytaine as cpt

# Solve the potential flow problems
radius = 1.0
mesh = cpt.mesh_horizontal_cylinder(radius=radius, length=1.0, center=(0, 0, 0), resolution=(5, 20, 5)).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs())
solver = cpt.BEMSolver()
results = solver.solve_all([cpt.RadiationProblem(body=body, omega=1.0, radiating_dof=dof) for dof in body.dofs])

# Define a set of points in the fluid on which the fluid velocities will be computed
Y, Z = np.meshgrid(np.linspace(-3.0, 3.0, 40), np.linspace(-3.0, -0.05, 20))
points = np.vstack([np.zeros(Y.size), Y.ravel(), Z.ravel()]).T
# Filter out points that are inside the floating body
points = np.array([[x, y, z] for (x, y, z) in points if y**2 + z**2 > radius**2])

# Compute the velocity field for each of the result object, that is each of the radiating dof
velocities = {}
for res in results:
    velocities[res.radiating_dof] = solver.compute_velocity(points, res)

import matplotlib.pyplot as plt
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
for ax, dof in zip(axs.ravel(), ["Sway", "Heave"]):
    # Display the arrows of the velocity
    y_velocities = -velocities[dof][:, 1].imag
    z_velocities = -velocities[dof][:, 2].imag
    ax.quiver(points[:, 1], points[:, 2], y_velocities, z_velocities, angles="xy")
    ax.axis('equal')

    # Display the boundary of the floating body
    theta = np.linspace(-np.pi, 0.0, 100)
    ax.plot(radius*np.cos(theta), radius*np.sin(theta), color="red")
    ax.set_title(dof)
    ax.set_xlabel("y")
    ax.set_ylabel("z")

plt.show()

B6. Animated free surface elevation

This example solves a diffraction problem, it computes the free surface elevation and shows it as a 3D animation.

"""This example shows how to animate the free surface elevation and pressure on a floating body using vedo."""

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import (
    airy_waves_free_surface_elevation,
    airy_waves_pressure,
)
from capytaine.ui.vedo_animations import Animation


# Parameters
omega = 3.0
T = 2 * np.pi / omega


## RESOLUTION WITH CAPYTAINE
positions = [
    (np.cos(angle) * 5.0, np.sin(angle) * 5.0, 0.0)
    for angle in [0.0, 2 * np.pi / 3, -2 * np.pi / 3]
]
meshes = [
    cpt.mesh_vertical_cylinder(length=4.0, radius=2.0, center=p, faces_max_radius=0.3)
    for p in positions
]
mesh = cpt.Mesh.join_meshes(*meshes)
sphere = cpt.FloatingBody(mesh=mesh.immersed_part(), lid_mesh=mesh.generate_lid())
solver = cpt.BEMSolver()
diffraction_problem = cpt.DiffractionProblem(body=sphere, wave_direction=0.0, omega=omega)
print("Solving the diffraction problem...")
diffraction_result = solver.solve(diffraction_problem)
print("Diffraction problem solved.")


## POST-PROCESSING THE PRESSURE AND FREE SURFACE ELEVATION
free_surface_mesh = cpt.mesh_rectangle(
    size=(100.0, 100.0), center=(0.0, 0.0, 0.0), resolution=(300, 300)
)
amplitude = 0.2
print("Computing the free surface elevation and pressure...")
diffraction_elevation_at_faces = amplitude * (
    solver.compute_free_surface_elevation(free_surface_mesh.vertices, diffraction_result)
    + airy_waves_free_surface_elevation(free_surface_mesh.vertices, diffraction_problem)
)
pressure_on_sphere = amplitude * (
    solver.compute_pressure(sphere.mesh.vertices, diffraction_result)
    + airy_waves_pressure(sphere.mesh.vertices, diffraction_problem)
)
print("Free surface elevation and pressure computed.")


## ANIMATION WITH VEDO
anim = Animation(loop_duration=T, fps=24)
anim.add_body(sphere, vertices_pressure=pressure_on_sphere)
anim.add_free_surface(free_surface_mesh, vertices_elevation=diffraction_elevation_at_faces)
anim.save("diffraction.mp4", size=(1024, 512))

B7. Animation of the RAO

This script generates the animation of the RAO motion for a wave incoming in front of a ship, such as the one used on the main page of this documentation. This script requires the mesh of the ship boat_200.mar. It can be downloaded from: https://raw.githubusercontent.com/capytaine/capytaine/master/docs/examples/src/boat_200.mar

from numpy import pi

import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
from capytaine.ui.vedo_animations import Animation

from vedo import Light


bem_solver = cpt.BEMSolver()


def generate_boat():
    boat_mesh = cpt.load_mesh("docs/examples/src/boat_200.mar", file_format="mar")
    boat = cpt.FloatingBody(
            mesh=boat_mesh,
            dofs=cpt.rigid_body_dofs(rotation_center=boat_mesh.center_of_buoyancy),
            center_of_mass = boat_mesh.center_of_buoyancy,
            )
    boat.inertia_matrix = boat.compute_rigid_body_inertia() / 40 # Artificially lower to have a more appealing animation
    return boat


def setup_animation(body, fs, omega, wave_amplitude, wave_direction):
    # SOLVE BEM PROBLEMS
    radiation_problems = [cpt.RadiationProblem(omega=omega, body=body.immersed_part(), radiating_dof=dof) for dof in body.dofs]
    radiation_results = bem_solver.solve_all(radiation_problems)
    diffraction_problem = cpt.DiffractionProblem(omega=omega, body=body.immersed_part(), wave_direction=wave_direction)
    diffraction_result = bem_solver.solve(diffraction_problem)

    dataset = cpt.assemble_dataset(radiation_results + [diffraction_result])
    rao = cpt.post_pro.rao(dataset, wave_direction=wave_direction)

    # COMPUTE FREE SURFACE ELEVATION
    # Compute the diffracted wave pattern
    incoming_waves_elevation = airy_waves_free_surface_elevation(fs.vertices, diffraction_result)
    diffraction_elevation = bem_solver.compute_free_surface_elevation(fs.vertices, diffraction_result)

    # Compute the wave pattern radiated by the RAO
    radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs.vertices, res) for res in radiation_results}
    radiation_elevation = sum(rao.sel(omega=omega, radiating_dof=dof).data * radiation_elevations_per_dof[dof] for dof in body.dofs)

    # SET UP ANIMATION
    # Compute the motion of each face of the mesh for the animation
    dofs_vertices_motions = {dof: body.dofs[dof].evaluate_motion_at_points(body.mesh.vertices) for dof in body.dofs}
    rao_vertices_motion = sum(rao.sel(omega=omega, radiating_dof=dof).data * dofs_vertices_motions[dof] for dof in body.dofs)

    # Set up scene
    animation = Animation(loop_duration=2*pi/omega)
    animation.add_body(body, vertices_motion=wave_amplitude*rao_vertices_motion)
    animation.add_free_surface(fs, wave_amplitude * (incoming_waves_elevation + diffraction_elevation + radiation_elevation))
    return animation


if __name__ == '__main__':
    body = generate_boat()
    fs = cpt.mesh_rectangle(size=(200.0, 200.0), resolution=(100, 100))

    anim = setup_animation(body, fs, omega=1.5, wave_amplitude=1.0, wave_direction=pi)

    params = dict(
        camera={"pos": (70, 70, 100)},
        lights=[Light([0, 0, 100], intensity=0.8)],
        resolution=(800, 600)
    )
    # anim.run(**params)
    anim.save("animated_boat.mp4", **params)

B8. Pressure field for zero or infinite frequency

The computations for zero and infinite frequencies are tricky, because several magnitudes such as the pressure are zero or infinity, but they can be expressed as \(\omega\) or \(\omega^2\) times some finite value. The example below shows an example of extracting this value.

import numpy as np
import capytaine as cpt

# Initialize an example mesh and test case.
mesh = cpt.mesh_sphere().immersed_part()
body = cpt.FloatingBody(mesh, dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)))

pb = cpt.RadiationProblem(body=body, omega=np.inf, radiating_dof="Heave")
# Zero and infinity are accepted values for "omega", as well as for "period", "wavelength" and "wavenumber"

solver = cpt.BEMSolver()
res = solver.solve(pb)

print(f"Heave-heave added mass: {res.added_mass['Heave']:.2e}")

p = np.real(res.pressure)
# p is a custom Python object storing the pressure on each panel

i_panel = 0  # Some panel for testing

# All pressures are actually infinite
print(f"Pressure on panel {i_panel}:", float(p[i_panel]))

# However p/omega^2 is finite (hence the finite added mass)
p_over_omega2 = p/pb.omega**2
print(f"Pressure/omega^2 on panel {i_panel}: {p_over_omega2[i_panel]:.2e}")

B9. Mean drift force

This example computes the mean drift force on a floating box with both the near field and the far field formulations, and compares them. The near field mean drift force is still an experimental feature, whose interface may change in the future.

import numpy as np
import xarray as xr

import capytaine as cpt
from capytaine.post_pro.mean_drift_force import far_field_mean_drift_force, near_field_mean_drift_force

mesh = cpt.mesh_parallelepiped(resolution=(30, 30, 30)).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs(), center_of_mass=(0,0,0))

wave_direction = [0, np.pi/4]
test_matrix = xr.Dataset(coords={
    'wavenumber': [1.0, 2.0],
    'wave_direction': wave_direction,
    'theta': np.linspace(-np.pi/8, 2*np.pi, 100),  # For far-field mean-drift-force.
    'radiating_dof': list(body.dofs.keys())
})

solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, body, keep_details=True)
rao = cpt.post_pro.rao(dataset)
nf_mdf = near_field_mean_drift_force(rao, dataset, body, solver=solver)
ff_mdf = far_field_mean_drift_force(rao, dataset)

with np.printoptions(precision=2, suppress=True):
    print('Surge')
    print("Far field: ", ff_mdf['drift_force_surge'].values)
    print("Near field: ", nf_mdf['near_field_mean_drift_force'].sel(influenced_dof='Surge').values)

    print()
    print('Yaw')
    print("Far field: ", ff_mdf['drift_force_yaw'].values)
    print("Near field: ", nf_mdf['near_field_mean_drift_force'].sel(influenced_dof='Yaw').values)

B10. Mean drift pressure field

This example extracts the pressure fields from the near field mean drift force computation, for a sea made of two waves with the same frequency: the first-order pressure on the hull, the mean drift pressure on the hull, and the waterline term spread along a semi-elliptic profile around the mean waterline. They are plotted on the mesh of the whole box, and the first-order pressure is animated over a wave period. To save the animation as a video, uncomment the call to anim.save (requires ffmpeg).

"""Example of extracting mean drift force pressure field from mean drift computation in Capytaine.
The whole file is heavily using DataArray, such that the meaning of each array dimension are explicit, for teaching purpose.
"""
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.spatial import cKDTree

import capytaine as cpt
from capytaine.post_pro.mean_drift_force import near_field_mean_drift_force

full_mesh = cpt.mesh_parallelepiped(resolution=(30, 30, 30))
mesh = full_mesh.immersed_part()
body = cpt.FloatingBody(
    mesh=mesh,
    lid_mesh=mesh.generate_lid(),
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)),
    center_of_mass=(0,0,0)
)
test_matrix = xr.Dataset(coords={
    'wavelength': [2.0],
    'wave_direction': [0.0, np.pi/3],
    'radiating_dof': list(body.dofs.keys())
})

solver = cpt.BEMSolver()
ds = solver.fill_dataset(test_matrix, body, keep_details=True)
rao = cpt.post_pro.rao(ds)
ds.update(near_field_mean_drift_force(rao, ds, body, solver=solver, output_pressure=True))

# The sea is made of the two waves, with the same frequency.
# Their complex amplitudes set their heights and their phases (at the origin, at t=0).
amplitude = xr.DataArray(
    np.array([[0.1, 0.05*np.exp(1j*np.pi/3)]]),
    dims=["wavelength", "wave_direction"],
    coords={"wavelength": ds.coords["wavelength"], "wave_direction": ds.coords["wave_direction"]}
)

# Pressures on the immersed hull, transferred to the faces of the full mesh
# without interpolation (nearest neighbors), setting zero on the emerged part.
distance, immersed_face = cKDTree(mesh.faces_centers).query(full_mesh.faces_centers)
immersed_face = xr.DataArray(immersed_face, dims=["full_mesh_face"])
is_immersed = xr.DataArray(distance < 1e-8, dims=["full_mesh_face"])
def remap_from_mesh_to_full_mesh(data):
    return data.isel(hull_face=immersed_face).where(is_immersed, 0.0)


# ============================ FIRST ORDER PRESSURE ============================
first_order_pressure = (
    amplitude * (
        ds["pressure_diffraction"] + ds["pressure_Froude_Krylov"]
        + (ds["pressure_radiation"] * rao).sum("radiating_dof")
        )).sum(dim=["wave_direction"])
# dims = (wavelength, hull_face)
first_order_pressure_full_mesh = remap_from_mesh_to_full_mesh(first_order_pressure)
# dims = (wavelength, full_mesh_face)


# ========================= SECOND ORDER HULL PRESSURE =========================
# Sum over all pairs of waves of the difference-frequency QTF
amplitude_kl = (amplitude.rename({"wave_direction": "wave_direction_k"})
                    * np.conjugate(amplitude.rename({"wave_direction": "wave_direction_l"})))
mean_drift_hull_pressure_per_pair_of_waves = amplitude_kl * ds["pressure_hull_mean_drift"]
mean_drift_hull_pressure = mean_drift_hull_pressure_per_pair_of_waves.sum(dim=["wave_direction_k", "wave_direction_l"])
# dims = (wavelength, hull_face)
assert np.max(np.abs(mean_drift_hull_pressure.values.imag)) < 1e-10
mean_drift_hull_pressure = mean_drift_hull_pressure.real  # Time independent pressure
mean_drift_hull_pressure_full_mesh = remap_from_mesh_to_full_mesh(mean_drift_hull_pressure)
# dims = (wavelength, full_mesh_face)


# ====================== SECOND ORDER WATERLINE PRESSURE =======================
def mean_of_product(a, b):
    # Difference-frequency QTF of the product of the first order fields a and b,
    # 1/4 (a_k conj(b_l) + conj(a_l) b_k), summed over all pairs of waves (k, l).
    a_k = a.rename({"wave_direction": "wave_direction_k"})
    a_l = a.rename({"wave_direction": "wave_direction_l"})
    b_k = b.rename({"wave_direction": "wave_direction_k"})
    b_l = b.rename({"wave_direction": "wave_direction_l"})
    qtf = (a_k * np.conjugate(b_l) + np.conjugate(a_l) * b_k) / 4
    return qtf.sum(["wave_direction_k", "wave_direction_l"]).real  # The imaginary part is zero.


def remap_waterline_pressure(amplitude, pressure_waterline_total, relative_elevation_waterline, mesh_waterline, box_mesh):
    """Mean waterline pressure of the sea made of all the waves, spread on the faces of a box around the mean waterline.

    The waterline term of the mean drift force is a lineic force <p w> per unit length of waterline, where p is the
    pressure of the waterline strip and w = h / sqrt(1 - n_z²) the width of the strip along the hull, with h the relative
    elevation of the free surface and n the normal of the hull. The walls of the box are vertical (n_z = 0), so w = h
    here: the lineic force is <p h>. It is spread in the vertical
    coordinate z around the mean waterline along a semi-elliptic profile of half-width |h| (the amplitude of h):

        <p>(z) = 2 <p h> / (pi |h|) * sqrt(1 - (z/|h|)^2)    for |z| < |h|,

    whose integral over z is <p h>. Each face of the box gets the value of this profile at its center, using the
    waterline segment with the same horizontal position.

    The waves have the same frequency: the mean value in the sea made of all of them is the sum over all pairs of
    waves (k, l) of the difference-frequency QTF, which includes the cross terms between different waves.

    Parameters
    ----------
    amplitude: DataArray with dims (..., wave_direction), complex
        Complex amplitude of each wave, possibly with other dimensions such as `wavelength`.
    pressure_waterline_total, relative_elevation_waterline: DataArrays with dims (..., wave_direction, waterline_segment), complex
        Outputs of `near_field_mean_drift_force(..., output_pressure=True)`, for a unit amplitude.
    mesh_waterline: DataArray with dims (waterline_segment, 2, 3)
        Output of `near_field_mean_drift_force(..., output_pressure=True)`.
    box_mesh: Mesh
        Mesh of the whole box, including its part above the waterline, with the same horizontal resolution as the
        immersed mesh.

    Returns
    -------
    DataArray with dims (..., full_mesh_face), real
        The other dimensions (such as `wavelength`) of the inputs are kept.
    """
    p = amplitude * pressure_waterline_total
    h = amplitude * relative_elevation_waterline
    ph = mean_of_product(p, h)  # dims = (..., waterline_segment)
    abs_h = np.sqrt(2*mean_of_product(h, h))  # dims = (..., waterline_segment), since <h h> = |h|^2/2

    _, segment = cKDTree(mesh_waterline.mean(axis=1)[:, :2]).query(box_mesh.faces_centers[:, :2])
    segment = xr.DataArray(segment, dims=["full_mesh_face"])
    u = xr.DataArray(box_mesh.faces_centers[:, 2], dims=["full_mesh_face"]) / abs_h.isel(waterline_segment=segment)

    # Value at the center of the faces. When |h| is not much larger than the faces, their average
    # (with the antiderivative u sqrt(1 - u^2) + arcsin(u) of 2 sqrt(1 - u^2)) is more accurate.
    return (2 * ph.isel(waterline_segment=segment) / (np.pi * abs_h.isel(waterline_segment=segment))
              * np.sqrt(np.clip(1 - u**2, 0.0, None)))  # zero for |u| > 1


mean_drift_waterline_pressure_full_mesh = remap_waterline_pressure(
        amplitude,
        ds["pressure_waterline_total"],
        ds["relative_elevation_waterline"],
        ds["mesh_waterline"],
        full_mesh,
        )
# dims = (wavelength, full_mesh_face)


# ================================== DISPLAY ==================================
# Display only a single frequency
freq_display = ds.coords["wavelength"].values[0]

first_order_pressure_displayed = first_order_pressure_full_mesh.sel(wavelength=freq_display)
mean_drift_hull_pressure_displayed = mean_drift_hull_pressure_full_mesh.sel(wavelength=freq_display)
mean_drift_waterline_pressure_displayed = mean_drift_waterline_pressure_full_mesh.sel(wavelength=freq_display)

omega = float(ds["omega"].sel(wavelength=freq_display))
period = 2*np.pi/omega

cmap = plt.get_cmap('coolwarm')
first_order_vabs = np.abs(first_order_pressure_displayed).max()
first_order_norm = plt.Normalize(vmin=-first_order_vabs, vmax=first_order_vabs)
mean_vabs = max(np.abs(mean_drift_hull_pressure_displayed).max(), np.abs(mean_drift_waterline_pressure_displayed).max())
mean_norm = plt.Normalize(vmin=-mean_vabs, vmax=mean_vabs)

fields = [
        ("First order pressure", first_order_pressure_displayed.real, first_order_norm),
        ("Mean pressure on hull", mean_drift_hull_pressure_displayed, mean_norm),
        ("Remapped mean waterline pressure", mean_drift_waterline_pressure_displayed, mean_norm),
        ]

fig = plt.figure(figsize=(15, 5))
axs = [fig.add_subplot(1, len(fields), i+1, projection='3d') for i in range(len(fields))]
for ax, (title, field, norm) in zip(axs, fields):
    full_mesh.show(
            backend="matplotlib",
            ax=ax,
            facecolors=cmap(norm(field)),
            linewidth=0.1,
            )
    ax.set_title(title)
    ax.view_init(elev=30, azim=-135)  # Show the faces x = -0.5 and y = -0.5, on which the incoming waves hit
first_order_collection = axs[0].collections[-1]
fig.suptitle(", ".join(f"wave from {np.degrees(beta):.0f}°: amplitude {np.abs(a):.2f}m, phase {np.degrees(np.angle(a)):.0f}°"
                       for beta, a in zip(ds.coords["wave_direction"].values, amplitude.sel(wavelength=freq_display).values)))
fig.colorbar(plt.cm.ScalarMappable(norm=first_order_norm, cmap=cmap), ax=axs[0], shrink=0.6)
fig.colorbar(plt.cm.ScalarMappable(norm=mean_norm, cmap=cmap), ax=axs[1:], shrink=0.6)

n_frames = 30
animation_duration = 3.0  # seconds, wall-clock time for one wave period


def update(frame):
    t = frame/n_frames * period
    phase = np.exp(-1j*omega*t)
    first_order_collection.set_facecolor(cmap(first_order_norm((first_order_pressure_displayed * phase).real)))
    return first_order_collection,


anim = animation.FuncAnimation(fig, update, frames=n_frames, interval=1000*animation_duration/n_frames, blit=False)

# anim.save("second_order_pressure_field.mp4", fps=n_frames/animation_duration)
plt.show()