%% JONSWAP Spectrum - arbitrary input W and F
clear 

% Input parameters
W = 8;        % wind speed [m/s], constant
F = 10000;     % Fetch [m]



gamma = 3.3;    % Peak factor - typical value
g = 9.81;       

% Compute JONSWAP parameters from wind & fetch
alpha = 0.076*(g*F/W^2)^(-0.22);
fp    = (3.5*g/W)*(g*F/W^2)^(-0.33);

f  = 0.01:0.001:1;   % Frequency [Hz]

sigma = 0.07*ones(size(f));
sigma(f > fp) = 0.09;

% JONSWAP spectrum formula
S = alpha * g^2 ./ (2*pi).^4 ./ f.^5 .* exp(-5/4*(fp./f).^4) ...
    .* gamma.^exp(-(f-fp).^2 ./ (2*sigma.^2*fp^2));


figure(1);
plot(f, S, 'LineWidth', 2);
xlabel('Frequency [Hz]');
ylabel('Spectral Density S(f) [m^2/Hz]');
title('JONSWAP Spectrum');
grid on;
hold on;

% Integrate to get Hs = 4*sqrt(m0)
moment0  = trapz(f, S);       % integral S(f) df
Hs_check = 4*sqrt(moment0);  
Tp       = 1/fp;

fprintf('Hs from integrated spectrum = %.2f m\n', Hs_check);
fprintf('Tp peak period = %.2f s\n', Tp);

%% PDF based on spec
N    = 1000;
Hmin = sqrt(2*min(S))
Hmax = sqrt(2*max(S));
Hran = (Hmax-Hmin)/N;
H    = Hmin:Hran:Hmax;
Hrms = Hs_check/1.416;
pRay = (2*H/(Hrms*Hrms)).*exp(-(H./Hrms).^2);

figure(2)
plot(H,pRay, 'LineWidth', 2); grid on;
xlabel('H [m]');
ylabel('Probability Density Function');
title('Rayleigh distribution');


%% visualize waves
nf   = length(S);
df   = f(2) - f(1);
dt   = 1/(10*f(end));
Per  = 1/f(1);
t    = dt:dt:Per;
nmt  = length(t);

phi    = 2*pi*rand(1,nf);  % fasi casuali
etaSum = zeros(nmt,1);
eta    = zeros(nmt,nf);

for i = 1:10:nf
    etaSum(:,1) = etaSum(:,1) + sqrt(2*S(i))*cos(2*pi*f(i)*t(:) + phi(i));
    eta(:,i)    = sqrt(2*S(i))*cos(2*pi*f(i)*t + phi(i));    
end

figure(3)
subplot(2,1,1)
plot(t, eta); grid on;
hold on
xlim([0 20])
%ylim([-0.05 0.05])
title("Resulting waves composition - from JONSWAP spectrum")
xlabel("Time [s]")
ylabel("Height [m]")

subplot(2,1,2)
plot(t, etaSum, "LineWidth", 2); grid on;
xlim([0 20])
title("Resulting surface η(t) - from JONSWAP spectrum")
xlabel("Time [s]")
ylabel("Height [m]")