clear 

H = 4.0;   % altezza onda 
d = 80.0;  % profondita 
L = 150.0; % lunghezza onda
T = 5.0;   % periodo
z = 0.0;   % quota particella 
Nperiodi = 3; % periodi da plottare

k     = 2*pi / L;
sigma = 2*pi / T;

% Coefficienti per zeta (orizzontale)
C1 = -(H/2)*cosh(k*(d+z))/sinh(k*d);
C2 = ((pi*H^2)/(8*L*(sinh(k*d))^2))*(1 - 3*cosh(2*k*(d+z))/(2*(sinh(k*d))^2));
C3 = ((pi*H^2)/(4*L))*cosh(2*k*(d+z))/(sinh(k*d))^2;

% Coefficienti per epsilon (verticale)
D1 = (H/2) * sinh(k*(d+z))/sinh(k*d);
D2 = ((3*pi*H^2)/(16*L))*sinh(2*k*(d+z))/(sinh(k*d))^4;

% Tempo
t = linspace(0, Nperiodi*T, 500);
theta = -sigma*t;  % x0=0

% Posizioni lagrangiane esatte
zeta    = C1*sin(theta) + C2*sin(2*theta) + C3*sigma*t;
epsilon = D1*cos(theta) + D2*cos(2*theta);

% Velocita (---> con derivate numeriche <---)
u = gradient(zeta, t);
w = gradient(epsilon, t);

figure(1)
hold on;
quiver(zeta,epsilon,10*u,10*w)
Npts = round(length(t)/Nperiodi);
for i=1:Nperiodi
    idx = (i-1)*Npts + 1;
    plot(zeta(idx), epsilon(idx), 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');
end
grid on; axis equal;
title('Traiettoria particella (\zeta, \epsilon)'); xlabel('\zeta (m)'); ylabel('\epsilon (m)');


Minimum_Horiz_velocity    = min(u)
Maximum_Horiz_velocity    = max(u)
Minimum_Vertical_velocity = min(w)
Maximum_Vertical_velocity = max(w)

module = sqrt(u.^2 + w.^2);
Minimum_module = min(module)
Maximum_module = max(module)