2018美赛A题海浪仿真

function s=bopu(fengji,duanshu)
u=[3,5,7,9,11,13,15,17];
wavewmax = [16.444115 9.866469 7.047480 5.481373 4.484760 3.794799 3.288826 2.90190];
if fengji>8
    fengji=8;
end
if fengji<1
    fengji=1;
end
fi=fengji;
u=u(fi);
wmin=0;
wmax=wavewmax(fi);
m=duanshu;
wavemn=(wmax-wmin)/m;
w=[wmin:wavemn:wmax];
s=0.81*exp(-7400./(w*u+eps).^4)./(w.^5+eps);

plot(w,s);

function [z]=erweihailangboxing(fengji,pinpushu,jiaodushu)


% 3 2.438306 16.444115  4.053570
% 5 1.462983 9.866469   2.432142
% 7 1.044989 7.047480   1.737244
% 9 0.812770 5.481373   1.351190
% 11 0.664988 4.484760  1.105519
% 13 0.562683 3.794799  0.935439
% 15 0.487659 3.288826  0.810714
% 17 0.430288 2.901905  0.715336


wavewmin = [2.438306 1.462983 1.044989 0.812770 0.664988 0.562683 0.487659 0.430288];
wavewmax = [16.444115 9.866469 7.047480 5.481373 4.484760 3.794799 3.288826 2.90190];
wavewp=[4.053570 2.432142 1.737244 1.351190 1.105519 0.935439 0.810714 0.715336];


%-----------------------------------------------------
u=[3,5,7,9,11,13,15,17];
%---------------------------------------------------


if fengji>8
    fengji=8;
end
if fengji<1
    fengji=1;
end
fi=fengji;
wmin=wavewmin(fi);
wmax=wavewmax(fi);
wp=wavewp(fi);
ui=u(fi);
M=pinpushu;
N=jiaodushu;
wavewn=(wmax-wmin)/M;
thetawn=pi/N;
dx=1;
dy=1;
x=[0:dx:500];
y=[0:dy:300];
[x,y]=meshgrid(x,y);
z=zeros(size(x));
for wi=1:M
    for ki=1:N
        theta=-pi/2+(ki-1)*thetawn;
         epsin=rand*2*pi;
         w=wmin+(wi-1)*wavewn+wavewn/2;
         swi=0.81*exp(-7400/(w*ui+eps).^4)*2*(cos(theta)).^2/(pi*(w.^5+eps));
        an=sqrt(2*swi*wavewn*theta);
       z1=w*w*x*cos(theta)/9.8+w*w*y*sin(theta)/9.8+epsin;
       z=an*cos(z1)+z;
    end
end
surfl(x,y,z);
shading interp;
lightangle(-45,30);
set(findobj(gca,'type','surface'),'FaceLighting','phong','AmbientStrength',.3,'DiffuseStrength',.8,...
    'SpecularStrength',.9,'SpecularExponent',200)



function [y1]=hailangboxing(fengji,duanshu)


% 3 2.438306 16.444115  4.053570
% 5 1.462983 9.866469   2.432142
% 7 1.044989 7.047480   1.737244
% 9 0.812770 5.481373   1.351190
% 11 0.664988 4.484760  1.105519
% 13 0.562683 3.794799  0.935439
% 15 0.487659 3.288826  0.810714
% 17 0.430288 2.901905  0.715336


wavewmin = [2.438306 1.462983 1.044989 0.812770 0.664988 0.562683 0.487659 0.430288];
wavewmax = [16.444115 9.866469 7.04748 5.481373 4.484760 3.794799 3.288826 2.901905];
wavewp=[4.053570 2.432142 1.737244 1.351190 1.105519 0.935439 0.810714 0.715336];


%-----------------------------------------------------
u=[3 5 7 9 11 13 15 17];
%---------------------------------------------------


if fengji>8
    fengji=8;
end
if fengji<1
    fengji=1;
end
fi=fengji;
wmin=wavewmin(fi);
wmax=wavewmax(fi);
%wp=wavewp(fi);
ui=u(fi);
M=duanshu;
%epsin=rand*2*pi;
wavewn=(wmax-wmin)/M;
dx=0.5;
dz=0.5;
x=[0:dx:125];
z=[0:dz:125];
[x,z]=meshgrid(x,z);
y=zeros(size(x));
for wi=1:M
    epsin=rand*2*pi;
    w=wmin+(wi-1)*wavewn;
    swi=0.81*exp(-7400/(w*ui+eps).^4)/(w.^5+eps);
    an=sqrt(2*swi*wavewn);
    y1=w*w*x/9.8+epsin;
    y=an*cos(y1)+y;
end
y1=y(25,:);

plot(y1);

function [data1]=SDwave(N,initWave) 
%输入N为：方阵的维数2^N+1,在方阵的四个角的顶点上放置的初值，maxWave为波动的程度 
%这里定为四个初值一样，这是可拼接的必要条件 
% Example: 
%       tic; 
%       d=SquareDiamond2(5,0,10); 
%       d=d-mod(d,1); 
%       colormap(pink); 
% 
%       surf(d); 
%       shading faceted 
%       axis equal 
%       toc 


%n=2^N; 
% data=zeros(n+1); 
% data(1,1)=initvalue; 
% data(1,n+1)=initvalue; 
% data(n+1,1)=initvalue; 
% data(n+1,n+1)=initvalue; 
data=initWave;
r0=0.03;
for i1=1:N
   w=(N-i1)/(N-1);
   r=r0*w*w;
   data=mytry(data,r); 
   
   nd=size(data,1);
for idn=2:nd-1
    for idm=2:nd-1
        data1(idn,idm)=(1/6)*(data(idn,idm-1)+data(idn,idm+1)+data(idn+1,idm)+data(idn-1,idm)+data(idn,idm)+data(idn,idm));
    end
end


end
colormap(winter);
surf(data1);






    


function [x]=rnd(absvalue) 
%扩展的随机函数生成器，产生绝对值小于absvalue的随机实数 
x=(rand(1)-0.5)*2*absvalue; 


function [data]=mytry(initdata,r) 
%square========================================= 
%x---x 
%----- 
%--0--  由四个x定中间的0 
%----- 
%x---x 
m=2;
n0=size(initdata,1);
for i1=1:n0
    for i2=1:n0
    data(i1*2-1,i2*2-1)=initdata(i1,i2);
    end
end
  n=n0*2-2;
for i=1:m:n 
  for j=1:m:n 
     data((i+i+m)/2,(j+j+m)/2)=(data(i,j)+data(i,j+m)+data(i+m,j)+data(i+m,j+m))/4+rnd(r); 
  end 
end 
%diamond======================================== 
%---x-- 
%----- 
%x-0-x  由四个x定中间的0 
%----- 
%---x-- 


%钻石步骤的横向部分 
%line No.1 and last 
for j=1+m/2:m:n 
  data(1,j)=(data(1,j+m/2)+data(1+m/2,j)+data(1,j-m/2)+data(n+1-m/2,j))/4+rnd(r); 
  data(n+1,j)=data(1,j); 
end 
%middle 
for i=1+m:m:n 
  for j=1+m/2:m:n 
     data(i,j)=(data(i,j+m/2)+data(i+m/2,j)+data(i,j-m/2)+data(i-m/2,j))/4+rnd(r); 
  end 
end 


%钻石步骤的纵向部分 
%line No.1 and last 
for i=1+m/2:m:n 
  data(i,1)=(data(i,1+m/2)+data(i+m/2,1)+data(i,n+1-m/2)+data(i-m/2,1))/4+rnd(r); 
  data(i,n+1)=data(i,1); 
end 
%middle 
for i=1+m/2:m:n 
  for j=1+m:m:n 
     data(i,j)=(data(i,j+m/2)+data(i+m/2,j)+data(i,j-m/2)+data(i-m/2,j))/4+rnd(r); 
  end 
end 


% if (m>2) 
%   data=mytry(data,m/2,r/3,n); 
% end 
%