water wave particle animation

Question

I have a problem animating the interference of multiple waves. my code receives from the user (by manual input) the number of waves, and for each wave the duration, amplitude, and angle (0 angle means the wave direction is from -x to + x, 180 means from -y to +y).

on the gif generated I'm plotting the particle paths and also the wave crest (eta). I'm generating 2 plots... one is the interference on plane XZ, and one on plane YZ.

and from some reason, the particles do not behave properly. the particles on z=0 do not match the wave crest.

• This is the Code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from sympy import symbols, Eq, nsolve, tanh
import math
import pandas as pd
import time



# Constants
g = 9.81
water_depth = 20
frames_per_second = 10

# Function to get wave data from the user
def get_waves():
    """
    Prompt the user to input data for each wave.
    Returns a list of tuples containing wave parameters.
    """
    global num_waves
    num_waves = int(input("Enter the number of waves: "))
    waves_data = []
    for i in range(num_waves):
        duration = float(input(f"Enter duration for wave {i+1}: "))
        amplitude = float(input(f"Enter amplitude for wave {i+1}: "))
        #phase = float(input(f"Enter start phase for wave {i+1}: "))
        phase = 0
        angle = float(input(f"Enter angle for wave {i+1} (0 for left to right, 180 for right to left): "))
        waves_data.append((duration, amplitude, phase, angle))
    return waves_data

# Function to calculate wavelength
def calculate_wavelength(duration, amplitude, angle):
    """
    Calculate wavelength based on wave parameters.
    """
    T = duration
    wavelength = symbols('wavelength')
    equation = Eq(g * T**2 * tanh(2 * np.pi * water_depth / wavelength) / (2 * np.pi) - wavelength, 0)
    return float(nsolve(equation, wavelength, 1))

# Function to calculate wave parameters for all waves
def calculate_wave_parameters(waves_data):
    """
    Calculate wave parameters including wavelength, wavenumber, and Omega for all waves.
    Returns a list of wave parameters and the maximum wavelength.
    """
    wave_parameters = []
    class MaxVariables:
        def __init__(self):
            self.wavelength = 0
            self.duration = 0
            self.amplitude = 0

    global Max
    Max = MaxVariables()
    for duration, amplitude, phase, angle in waves_data:
        wavelength = calculate_wavelength(duration, amplitude, angle)
        wavenumber = 2 * np.pi / wavelength
        theta = np.radians(angle) #convert the wave angle to radigns and give the name theta
        kx = wavenumber * np.cos(theta)
        ky = wavenumber * np.sin(theta)
        Omega = math.sqrt(g * wavenumber * math.tanh(wavenumber * water_depth))
        wave_parameters.append((duration, amplitude, phase ,angle, wavelength, wavenumber, kx, ky, Omega))
        if wavelength > Max.wavelength:
            Max.wavelength = wavelength
        if duration > Max.duration:
            Max.duration = duration
        if amplitude > Max.amplitude:
            Max.amplitude = amplitude
  
    return wave_parameters, Max

# Function to calculate an array of particle positions for each wave
def calculate_waves_particles(wave_parameters, Max, nh, nv, t_range, plane):
    
    waves_particles = np.zeros((num_waves, len(t_range), nh, nv, 2))    
    
    #creation of a reference points array for all the waves.
    pos2d_init = np.zeros((nh, nv, 2))
    for i in range(nh):
        for j in range(nv):
            pos2d_init[i, j, 0] = j * Max.wavelength * times / (nv - 1)
            pos2d_init[i, j, 1] = -i * water_depth / (nh - 1)
    
    #calculation of particle posotions for all waves in each wave's direction vector and in all time frames            
    for i, (duration, amplitude, phase, angle, wavelength, wavenumber, kx, ky, Omega) in enumerate(wave_parameters):
        #    Pre-calculate particle data for each wave, 
        particle_data = np.zeros((len(t_range),nh, nv, 2))
        theta = angle * np.pi / 180
        particle_data[:,:,:,0] = pos2d_init[:,:,0]
        particle_data[:,:,:,1] = pos2d_init[:,:,1]
        kn = kx   #use default value for xz plane
        if plane == 'yz': # if plane is yz, change default value of k to ky
            kn = ky  
        for n in range(len(t_range)):
            t = t_range[n]
            #Eliptical motion:
            for j in range(nh):
                for k in range(nv):   #"""the problem is here.  something with the removal of the initial positions I think"""
                    e1 = pos2d_init[j, k, 0] #calculation of initial X/y position
                    z = pos2d_init[j, k, 1] #calculation of initial Z position
                    particle_data[n, j, k, 0] = -(g * amplitude * kn / Omega**2) * (math.cosh(kn * (water_depth + z)) / math.cosh(kn * water_depth)) * math.sin(kn*e1 - Omega * t)
                    particle_data[n, j, k, 1] = (g * amplitude * kn / Omega**2) * (math.sinh(kn * (water_depth + z)) / math.cosh(kn * water_depth)) * math.cos(kn*e1 - Omega * t)    
            # Update waves_particles with the calculated particle positions for the current wave 
            #But, it reduces the initial position so that each point stores onl its circular motion. this is important for summing purposes
            waves_particles[i,n,:,:,:] = particle_data[n,:,:,:]
  
    return waves_particles, pos2d_init
                        
# Function to calculate surface elevation (eta) for all waves
def calculate_waves_eta(wave_parameters, t_range, Max, plane):
    global eta_intervals
    eta_intervals = 100 * times
    #eta_waves = []  # List to store individual wave eta data
    waves_eta = np.zeros((num_waves, len(t_range), eta_intervals))  # Initialize waves_eta with 100 points
        
    #for duration, amplitude, angle, wavelength, wavenumber, Omega in wave_parameters:
    for i, (duration, amplitude, phase, angle,  wavelength, wavenumber, kx, ky, Omega) in enumerate(wave_parameters):    
        kn = kx   #use default value for xz plane
        if plane == 'yz': # if plane is yz, change default value of k to ky
            kn = ky
        theta = angle * np.pi / 180
        global N
        N = np.linspace(0, Max.wavelength * times, eta_intervals ) #create for each wave a linspace in its n direction with size equal to biggest wavelength
        for n in range(len(t_range)):
            for ni in range(eta_intervals):
                t = t_range[n]
                waves_eta[i, n, ni] =  amplitude * np.cos(kn * N[ni] - Omega * t)
               
    return waves_eta


def diffract(waves_particles_xz, waves_particles_yz, waves_eta_xz, waves_eta_yz, t_range, pos2d_init):
    diff_particles_xz = np.zeros((len(t_range),nh, nv, 2)) 
    diff_particles_yz = np.zeros_like(diff_particles_xz)
    diff_particles_xz[:,:,:,:] += pos2d_init
    diff_particles_yz[:,:,:,:] += pos2d_init
    diff_eta_xz = np.zeros((len(t_range), eta_intervals))
    diff_eta_yz = np.zeros((len(t_range), eta_intervals))
    
    Z_yz_y0 = np.zeros((num_waves, len(t_range),nh, nv, 2))
    Z_xz_x0 = np.zeros_like(Z_yz_y0)
   
    # Loop over each wave
    for i in range(num_waves):
        # Loop over each point in the grid    
        for ti in range(len(t_range)):
            diff_eta_xz[ti, :] += waves_eta_xz[i, ti, :] 
            diff_eta_yz[ti, :] += waves_eta_yz[i, ti, :] 
            for ni in range(nh):
                for nj in range(nv):
                    Z_yz_y0[i, ti, ni, nj, 1] = waves_particles_yz[i, ti, ni, 0, 1]
                    Z_xz_x0[i, ti, ni, nj, 1] = waves_particles_xz[i, ti, ni, 0, 1]
                    
            diff_particles_xz[ti, :, :, 0] += waves_particles_xz[i, ti, :, :, 0] 
            diff_particles_xz[ti, :, :, 1] += waves_particles_xz[i, ti, :, :, 1] + Z_yz_y0[i, ti, :, :, 1]
            
            diff_particles_yz[ti, :, :, 0] += waves_particles_yz[i, ti, :, :, 0]
            diff_particles_yz[ti, :, :, 1] += waves_particles_yz[i, ti, :, :, 1] + Z_xz_x0[i, ti, :, :, 1]
    ##### for debug
    # Convert the array to a pa    
    # Convert the array to a pandas DataFrame
    ixz = pd.DataFrame(diff_particles_xz[30,0])
    iyz = pd.DataFrame(diff_particles_yz[30,0])
    etaxz = pd.DataFrame(np.column_stack((N,diff_eta_xz[30])))
    etayz = pd.DataFrame(np.column_stack((N,diff_eta_yz[30])))   
    ixz.to_excel('ixz.xlsx', index=False)
    iyz.to_excel('iyz.xlsx', index=False)
    etaxz.to_excel('etaxz.xlsx', index=False)
    etayz.to_excel('etayz.xlsx', index=False)


    return diff_particles_xz, diff_particles_yz, diff_eta_xz, diff_eta_yz
 
def animate1(diff_particles_xz, diff_eta_xz, diff_particles_yz, diff_eta_yz, Max):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) # Create a single subplot with two axes
    def update1(frame):
        ax1.cla()  # Clear previous plots in the first axis
        ax2.cla()  # Clear previous plots in the second axis

        # Define data structure to hold coordinates for both 'xz' and 'yz' planes
        planes_data = [
            (diff_particles_xz[:, :, :, 0], diff_particles_xz[:, :, :, 1], diff_eta_xz[frame, :], ax1),
            (diff_particles_yz[:, :, :, 0], diff_particles_yz[:, :, :, 1], diff_eta_yz[frame, :], ax2)
        ]

        # Plot eta and scatter plot for current positions for both planes
        for n_coords, z_coords, eta, ax in planes_data:
            ax.plot(np.linspace(0, Max.wavelength * times, eta_intervals), eta, c='b', label='Wave crest')
            ax.scatter(n_coords[frame], z_coords[frame], s=10, alpha=0.8, label="Current positions")

        # Plot eta, scatter plot for current positions, and particle paths for both planes
        for n_coords, z_coords, eta, ax in planes_data:
            # Plot lines between particles in the first row
            for j in range(nv-1):
                ax.plot(
                    [n_coords[frame,0, j], n_coords[frame,0, j + 1]],  # X coordinates of the line
                    [z_coords[frame,0, j], z_coords[frame,0, j + 1]],  # Z coordinates of the line
                    color="red",  # Example color
                    alpha=0.5,    # Example alpha value
                    linewidth=2,  # Example line width
                )
                # plots particle paths using all frames
                for i in range(nh):
                    ax.plot(
                        n_coords[:,i, j],
                        z_coords[:,i, j],
                        color="C{}".format(j),
                        alpha=0.3,
                        label="Particle {}".format(i * nv + j + 1),
                    )
        
        ax1.set_title("Interference projection of water waves on XZ plane")
        ax1.grid(True)
        ax2.set_title("Interference projection of water waves on YZ plane")
        ax2.grid(True)
        ax1.set_xlim(-0.1 * Max.wavelength,times * Max.wavelength)
        ax1.set_ylim(-1.1 * water_depth, 0.2 * water_depth)
        ax2.set_xlim(-0.1 * Max.wavelength,times * Max.wavelength)
        ax2.set_ylim(-1.1 * water_depth, 0.2 * water_depth)
        x_range = times * Max.wavelength - (-0.1 * Max.wavelength)
        y_range = 0.2 * water_depth - (-1.1 * water_depth)
        aspect_ratio = x_range / y_range
        # Set the aspect ratio manually
        #ax1.set_aspect(aspect_ratio)
        #ax2.set_aspect(aspect_ratio)
        ax1.autoscale(False)
        ax2.autoscale(False) 
 
    ani = FuncAnimation(fig, update1, frames=len(t_range), interval=500 / frames_per_second)
    # Generate unique run stamp
    run_stamp = time.strftime("%Y%m%d_%H%M%S")
    gif_filename = f'Interference_XZ_YZ_{run_stamp}.gif'
   
    ani.save(gif_filename, writer='pillow', fps=frames_per_second)
    plt.close(fig)

    
# Main program
waves_data = get_waves()
global nh
global nv
global times
times = 2
nh = 15
nv = 15

wave_parameters, Max = calculate_wave_parameters(waves_data)
t_range = np.linspace(0, 2 * Max.duration, int(2 * Max.duration * frames_per_second))  # Create time range

waves_particles_xz, pos2d_init = calculate_waves_particles(wave_parameters, Max, nh, nv, t_range, 'xz')
waves_particles_yz, pos2d_init = calculate_waves_particles(wave_parameters, Max, nh, nv, t_range, 'yz')

waves_eta_xz = calculate_waves_eta(wave_parameters, t_range, Max, 'xz')
waves_eta_yz = calculate_waves_eta(wave_parameters, t_range, Max, 'yz')

diff_particles_xz, diff_particles_yz, diff_eta_xz, diff_eta_yz =  diffract(waves_particles_xz, waves_particles_yz, waves_eta_xz, waves_eta_yz, t_range, pos2d_init)

animate1(diff_particles_xz, diff_eta_xz, diff_particles_yz, diff_eta_yz, Max)

I tried to plot a line between the upper most particles (z=0) (in red color) and see if it resembles the eta line (blue). It matches only in node points.