% PLOT of Dispersion Relations and eta(x,t)
% knowing d and T ---> L and C can be computed
clear  
g = 9.8;        % gravit. acc.

%Chose d and T
d = 2.3;            % water deph  (d at the deep water limit kd \sim 1 - d=L/(2*pi)) 
T = 10;            % wave period
H = 2;             % wave height

%% FIRST SECTION - plot tanh

k = 0:0.01:100; % wavenumber

% visualize tanh
f     = tanh(k*d);
xaxis = d.*k./(2*pi);
figure(1)
%
subplot(2,1,1)
loglog(k,f)
hold on
title("visualize tanh(kd)")
xlabel("k")
ylabel("tanh(kd)")
grid on
%
subplot(2,1,2)
loglog(k,xaxis)
title("convert k to meters")
xlabel("k")
ylabel("d/L")
grid on

%% SECOND SECTION - compute C,L 

% Initial (guess) value for L
L     = 0.1;
niter = 100; % Number of max iterations 
for i=1:niter
    C = (g*T/(2*pi))*tanh(2*pi*d/L);
    % Compute L based on C and T
    L = C*T;
end

L
C
% Check L
%check = (g*T*T/(2*pi))*tanh(2*pi*d/L)
relative_depth = d/L
steepness = H/L
particle_vel = pi*H/T


%% THIRD SECTION - particle orbit

% chose z < d 
z = -2;

t = 0:0.01:T;            % time window
nmt = length(t);

x   = 0:0.1:100;
sigma = (2*pi)/T;  % Compute the angular frequency

for j=1:10:nmt
    time = t(j);
    % position
    zeta  = -0.5*H*( cosh((2*pi/L)*(z+d))/sinh((2*pi/L)*d) )*sin((2*pi/L)*x(2) - sigma*time);
    eps   =  0.5*H*( sinh((2*pi/L)*(z+d))/sinh((2*pi/L)*d) )*cos((2*pi/L)*x(2) - sigma*time);
    % velocity
    u   =  (pi*H/T)*( cosh((2*pi/L)*(z+d))/sinh((2*pi/L)*d) )*cos((2*pi/L)*x(2) - sigma*time);
    w   =  (pi*H/T)*( sinh((2*pi/L)*(z+d))/sinh((2*pi/L)*d) )*sin((2*pi/L)*x(2) - sigma*time);
    figure(4)
    scatter(zeta,eps)
    hold on
    quiver(zeta,eps,u,w)
    grid on
    maxrad = 0.5*H*( sinh((2*pi/L)*(z+d))/sinh((2*pi/L)*d) ) + 0.5;
    xlim([-maxrad maxrad])
    ylim([-maxrad maxrad])
    drawnow
    title("visualize particle position varying in time")
    xlabel("zeta (meters)")
    ylabel("eps  (meters)")
end

%% THIRD PART - visualize eta

to = 0:0.01:T*0.2;        % time window
nmt = length(to);

figure(2)
for j=1:nmt
    time = to(j);
    eta  = 0.5*H*cos((2*pi/L).*x - sigma*time);
    plot(x,eta)
    grid on
    plotex = 2.0*H;
    ylim([-d d+plotex])
    drawnow
    %pause(0.0001)
    title("visualize eta(x,t) varying in time")
    xlabel("length (meters)")
    ylabel("eta(x,t)")
end

%close all
%% FOURTH PART - visualize eta - Stokes

% compute selected k
k = 2*pi/L;

ampStokes = cosh( k*d*( 2 + cosh(2*k*d) ) )/sinh(k*d)^3;

figure(3)
for j=1:nmt
    time = to(j);
    etaL  = 0.5*H*cos((2*pi/L).*x - sigma*time);
    etaS  = (pi*H*H/(8*L))*ampStokes*cos(2*((2*pi/L).*x - sigma*time));
    plot(x,etaL,x,etaS)
    plotex = 2.0*H;
    ylim([-d d+plotex])
    grid on
    drawnow
    title("visualize eta(x,t) varying in time")
    xlabel("length (meters)")
    ylabel("eta(x,t)")
end


%% Group velocity - wave sovrapposition
dt   = 0.001;
endt = 10*pi;
t  = 0:dt:endt;
f1 = 3;
f2 = 7;

fun1 = sin(2*pi*f1*t);
fun2 = sin(2*pi*f2*t);

sumf = fun1 + fun2;


% fourier
n     = length(sumf);
nhalf = n*0.5;
trasf = abs(fft(sumf))/nhalf;
freq  = (1/endt):((1/dt) - (1/endt))/(n -1):(1/dt);


figure(5)
subplot(2,1,1)
plot(t,fun1)
hold on
plot(t,fun2)
plot(t,sumf)
xlim([0 pi])
subplot(2,1,2)
loglog(freq(1:end/2),trasf(1:end/2))
grid on