海面几何建模-过滤法-PM谱

clear;clc;
[x,y]=meshgrid(-100:0.5:99.75);
dkx=2.*pi./200;
dky=2.*pi./200;
kx=[(dkx:dkx:(200.*dkx)),-fliplr(dkx:dkx:(200.*dkx))];
ky=kx;
[kx,ky]=meshgrid(kx,ky);
A=dkx.*dky;

B=normrnd(0,1,400,400);%均值，方差,矩阵,正太分布
B2=normrnd(0,1,400,400);
r=(B+1i.*B2)./2.^0.5;
rr=(fliplr(flipud(B))+1i.*fliplr(flipud(B2)))./2.^0.5;
z=conj(rr);
u=5;s=10;phi0=pi./2;
S=(0.5.*0.0081.*((kx.^2+ky.^2).^0.5).^(-3).*exp(-0.74.*9.81.^2.*u.^(-4).*((kx.^2+ky.^2).^0.5).^(-2))).*(2^(2*s-1)*gamma(s+1)^2./(pi*gamma(2*s+1))*cos((atan2(kx,ky)-phi0)./2).^(2*s))./((kx.^2+ky.^2).^0.5);
Q=(0.5.*0.0081.*((kx.^2+ky.^2).^0.5).^(-3).*exp(-0.74.*9.81.^2.*u.^(-4).*((kx.^2+ky.^2).^0.5).^(-2))).*(2^(2*s-1)*gamma(s+1)^2./(pi*gamma(2*s+1))*cos((pi-atan2(kx,ky)-phi0)./2).^(2*s))./((kx.^2+ky.^2).^0.5);

H=r.*(S.*A).^0.5.*exp(1i.*(9.81.*((kx.^2+ky.^2).^0.5.*(1+(kx.^2+ky.^2)./363)).^0.5))+conj(rr).*(Q.*A).^0.5.*exp(-1i.*(9.81.*((kx.^2+ky.^2).^0.5.*(1+(kx.^2+ky.^2)./363)).^0.5));

Z=ifft2(400.*400.*H);
Z=real(Z);

mesh(y,x,Z);

matlab运行结果如封面图所示，如有错误，欢迎指正。