随机波浪Jonswap谱

随机海浪往往具有统计特征，组成频率会呈现出某一频率集中的特征。由此而衍生出的海浪谱多种多样。其中较为著名的一种海浪谱Jonswap被广泛应用在海洋科学、海洋工程领域。
以合田改进的Jonswap谱（1999）为例：
$S(f)={\beta }_j{H}_{1/3}^2{T}_P^{-4}{f}^{-5}\exp [-\frac{5}{4}(T_{P}f{)}^{-4}]{\gamma }^{\exp [-(\frac{f}{f_P}-1{)}^2/2{\sigma }^2]}$
其中，
${\beta }_j=\frac{0.06238}{0.230+0.0336\gamma -0.185(1.9+\gamma {)}^{-1}}[1.094-0.01915\ln \gamma ]$,
$T_P=\frac{T_{H_{1/3}}}{1-0.132(\gamma +0.2{)}^{-0.559}}$,

对于平均Jonswap谱来说：
$\gamma =3.3,{\sigma }_a=0.07,{\sigma }_b=0.09,$
$\alpha =0.076{\bar{X}}^{-0.22},\bar{X}=10^{-1}to10^5,$
${\omega }_m=22(g/\bar{U}){\bar{X}}^{-0.33}$
或$f_m=3.5(g/\bar{U}){\bar{X}}^{-0.33}，$
$\bar{X}=gX/U^2.$
峰型参数$\sigma ={\sigma }_a$(当$\omega <={\omega }_m$时)，$\sigma ={\sigma }_b$(当$\omega >{\omega }_m$)时。

% Improved Jonswap Spectral
% Designed by: JN-Cui 
% Modified on 12/09/2019
%% DEFINITIONS
% alpha - energy scale factor; gama - spectral peak elevation factor;
% omega_m - spectral peak circular frequency; f_m - spectral peak frequency;
% U - wind speed at 10 m above sea surface; H_s - significant wave height;
% g - gravity acceleration;
%% FOR AVERAGE JONSWAP SPECTRAL
% gama=3.3; k=83.7; sigma_a=0.07; sigma_b=0.09;
% alpha=0.076*(X_bar)^(-0.22);
% X_bar=10^(-1)~20^(5); omega_m=22(g/U)*(X_bar)^(-0.33);
% f_m=3.5(g/U)(X_bar)^(-0.33);
%% IMPUT PARAMETERS
% H_s - significant wave height; T_s -  wave period at 1/3 wave height
% dm - calculation interval of omega
%% FUNCTION

function [W_s,t,x,Omega,S]=random_wave_v2_0(H_s,T_p,dm,dt,T,dx,X)
d=2000;
g=9.8;
gama=3.3; 
sigma_a=0.07;
sigma_b=0.09;
T_s=T_p*(1-0.132*(gama+0.2)^(-0.559));
f_p=1/T_p;
omega_p=f_p*2*pi;
i=1;
S_o=zeros(1,length(0:dm:1/T_s*2*pi*4));
Omega1=zeros(1,length(0:dm:1/T_s*2*pi*4));
for omega=0:dm:1/T_s*2*pi*4
    if omega<omega_p
        sigma=sigma_a;
    else
        sigma=sigma_b;
    end
    S_o(i)=320*H_s^2/(1/f_p)^4*omega^(-5)*exp(-1950/(1/f_p)^4*omega^(-4))*gama^exp(-0.5*((omega-omega_p)/(sigma*omega_p))^2.0);
    Omega1(i)=omega;
    i=i+1;
end
Omega=Omega1(2:end);
S=S_o(2:end);
N=length(S);
M=T/dt;
e=rand(1,N)*2*pi;
w=Omega;
f=zeros(1,length(w));
T_w=zeros(1,length(w));
K=zeros(1,length(w));
k=zeros(1,length(w));
a=zeros(1,length(w));
L=zeros(1,length(w));
for i=1:N
    f(i)=w(i)/2/pi;
    T_w(i)=1/f(i);
    K(i)=g*T_w(i)^2/(2*pi);
    fun=@(L1) L1/tanh(2*pi/L1*d)-K(i);
    L(i)=fzero(fun,K(i));
    k(i)=(2*pi/L(i));
    a(i)=sqrt(2*S(i)*(w(2)-w(1)));
end
t=zeros(1,M);
for i=1:M
    t(i)=(i-1)*dt;
end
NN=floor(X/dx)+1;
x=0:dx:X;
W_s=zeros(M,NN-1);
for j=1:M
    for ii=1:NN
        W_s(j,ii)=sum(a.*cos(k.*x(ii)-w.*t(j)+e));
    end
end
end