Neumann_and_PM_and_Jonswap.m

Neumann_and_PM_and_Jonswap.m

clear all;
close all;
clc;
%%
format short

w_spectrum = 0:(pi/2 + pi/10)/2047:(pi/2 + pi/10);
w_size = size(w_spectrum);

g = 9.8;
U10 = 15;
U7_5 = 14.46;
U19_5 = 16.24;

[Sw, sigma2_neumann] = Sw_neumann(w_spectrum);
A2_neumann = 2 * Sw;
[Sw,  m0, m2, m4] = Sw_PM(w_spectrum);
A2_PM = 2 * Sw;
[Sw, sigma2_jonswap] = Sw_jonswap(w_spectrum);
A2_jonswap = 2 * Sw;
gca = subplot(3, 2, 1);
plot(w_spectrum, A2_neumann,'r', 'linewidth', 1.3);
hold on
plot(w_spectrum, A2_PM,'b', 'linewidth', 1.3);
hold on
plot(w_spectrum, A2_jonswap,'g', 'linewidth', 1.3);

xlabel({'ω (s^-^1)', '图1 不同海浪谱的区别和联系'});
ylabel('A^2 (m^2·s)');
set(gca, 'XTick',[pi/10, pi/5, 3*pi/10, 2*pi/5, pi/2]);
set(gca,'TickLabelInterpreter','latex')
set(gca, 'XTicklabel',{'$\frac{\pi}{10}$', '$\frac{\pi}{5}$', '$\frac{3\pi}{10}$', '$\frac{2\pi}{5}$', '$\frac{\pi}{2}$'});
legend('Neumann', 'PM', 'Jonswap')


%%
delta_w = 2*pi/2048
w_plot = (0.0001 + delta_w / 2):delta_w:2*pi;
w_plot_size = size(w_plot)

f_head = w_plot / (2*pi);
% delta_f = 0.05 / (2*pi); % 0.05 %delta_w
delta_f = delta_w / (2*pi); % 0.05 %delta_w
f_head_size = size(f_head);

%%
k_head = (2*pi*f_head).^2 / g; % 4.8.1-2

a_fi_head = sqrt(Sf_PM(f_head).*delta_f);
a_fi_head_size = size(a_fi_head);

x = 1:10;
x_size = size(x);
epsilon = rand(w_plot_size(1), w_plot_size(2)) * 2 * pi - pi;

x_record = [];
i_next = 0;
Fs = 2; % 采样频率为 2 Hz
T = 1 / Fs;
L = 2048; % 信号持续时间为 1024 秒，17 分钟
t = (0:L-1)*T; % 时间向量
t_size = size(t);

for t_index = t
    ksi_array = [];
    for x_index = x
        row = 2 .* (a_fi_head .* cos(k_head*x_index - 2.*pi*f_head*t_index + epsilon)); % P276 4.8.1-1
        ksi = sum(row, 2);
        ksi_array = [ksi_array ksi];
    end
    x_record = [x_record ksi_array(1, x_size(2) / 2)];
    i_next = i_next + 1;
end

x_record_size = size(x_record);
    
subplot(3, 2, 2);
plot(t, x_record);
xlabel('图2 采样序列');
ylabel('波高');

%%
subplot(3, 2, 3);
Y = fft(x_record);
P2 = abs(Y / L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

size_p1 = size(P1)
size_f = size(f)

plot(f,P1) 
xlabel('图3 频率');
ylabel('fft幅值');

%%
subplot(3, 2, 4);
pp = (P1) .^ 2 * 2 * pi;
pp = pp * max(Sf_PM(f_head)) / max(pp) * 1.1;
plot(f, pp); %smooth()
xlabel('图4 频率');
ylabel('fft幅值，未平滑');
hold on
plot(f, Sf_PM(f));

%%
subplot(3, 2, 5);
plot(f, smooth(pp) * 2); %smooth
xlabel('图5 频率');
ylabel('fft幅值，平滑后');
hold on
plot(f, Sf_PM(f));

%%
H_mean = sqrt(2 * pi * m0) % 4.4.1-2
H_rms = sqrt(8 * m0) % 4.4.1-4

H_1_3 = 1.146 * H_rms % = 4.005 * sqrt(m0)
H_1_10 = 1.80 * H_rms % = 5.091 * sqrt(m0)
H_1_100 = 2.359 * H_rms % = 6.672 * sqrt(m0)

HF_1_100 = H_mean * sqrt(4 / pi * log10(1/(1/100))) % 4.3.4-4
HF_5_100 = H_mean * sqrt(4 / pi * log10(1/(5/100))) % 4.3.4-4
HF_10_100 = H_mean * sqrt(4 / pi * log10(1/(10/100))) % 4.3.4-4
m0, m2, m4

T0_2 = 2 * pi * sqrt(m0 / m2)
T2_4 = 2 * pi * sqrt(m2 / m4)
    
delta = m0 * m4 - m2 ^ 2 % 4.3.2-5
epsilon = sqrt(delta / (m0 * m4)) % 4.3.2-16

%% 计算波高
H = [];

buffer = [];
flag = 1;
for i_next = x_record
    if flag == 1
        if i_next > 0
            buffer = [buffer i_next];
            continue
        elseif i_next < 0
            H = [H, max(buffer)];
            buffer = [];
            buffer = [buffer i_next];
            flag = -1;
            continue
        end
    elseif flag == -1
        if i_next < 0
            buffer = [buffer i_next];
            continue
        elseif i_next > 0
            H = [H, min(buffer)];
            buffer = [];
            buffer = [buffer i_next];
            flag = 1;
            continue
        end
    end
end

% H_delta = abs(H(1:end-1) - H(2:end));
% % H_delta = sort(H_delta);
% subplot(3, 2, 6);
% size_H = size(H_delta);
% plot(1:1:size_H(2), H_delta, 'o-');
% xlabel('时间');
% ylabel('波高 P170');

%%
function [Sw, sigma2] = Sw_neumann(w)
    w_size = size(w)
    w_size2 = w_size(2);

    g = 9.8;
    U10 = 15;
    U7_5 = 14.46;
    U19_5 = 16.24;
    C = 3.05;

    Sw_neumann = C .* pi / 4 ./ (w.^6) .* exp(-2 * g^2 ./ (U7_5^2 .* w.^2));
    sigma2 = nansum(Sw_neumann .* (w(end)/w_size2));
    Sw = Sw_neumann;
end

%%
function [Sw, m0, m2, m4] = Sw_PM(w)
    w_size = size(w);
    w_size2 = w_size(2);
    
    g = 9.8;
    U10 = 15;
    U7_5 = 14.46;
    U19_5 = 16.24;
    a = 8.10*10^(-3);
    b = 0.74;
    
    Sw_PM = a .* g^2 ./ w.^5 .* exp(- b .* (g ./ (U19_5.*w)).^4);
    m0 = nansum(Sw_PM .* w .* (w(end)/w_size2)); % 4.3.2-2
    m2 = nansum(Sw_PM .* (w .^ 2) .* (w(end)/w_size2)); % 4.3.2-2
    m4 = nansum(Sw_PM .* (w .^ 4) .* (w(end)/w_size2)); % 4.3.2-2
    
    Sw = Sw_PM;
end

%%
function [Sf, sigma2] = Sf_PM(f)
    w = 2 * pi .* f;
    w_size = size(w);
    w_size2 = w_size(2);
    
    g = 9.8;
    U10 = 15;
    U7_5 = 14.46;
    U19_5 = 16.24;
    a = 8.10*10^(-3);
    b = 0.74;
    
    Sf_PM = 2 * pi * a .* g^2 ./ w.^5 .* exp(- b .* (g ./ (U19_5.*w)).^4);
    sigma2 = nansum(Sf_PM .* (w(end)/w_size2));
    Sf = Sf_PM;
end

%%
function [Sw, sigma2] = Sw_jonswap(w)
    w_size = size(w);
    w_size2 = w_size(2);

    g = 9.8;
    U10 = 15;
    U7_5 = 14.46;
    U19_5 = 16.24;
    w0 = 0.877 * g / U10;
    x_scope = (0.877 / 22)^(- 1/0.33);
    a_jonswap = 0.076 * (x_scope ^ (-0.23));
    delta_jonswap = repmat(0.07, w_size(1), w_size(2));
    delta_jonswap(w>w0) = 0.09;
    r_jonswap = 3.3;
    r_exp = exp(- (w - w0).^2 ./ (2 * delta_jonswap.^2 .* w0.^2));
    
    Sw_jonswap = a_jonswap .* g^2 ./ w.^5 .* exp(- 5 / 4 .* (w0./w).^4) .* r_jonswap .^ r_exp;
    sigma2 = nansum(Sw_jonswap .* (w(end)/w_size2));
    Sw = Sw_jonswap;
end