generate1.m

通过PM函数产生而为PM随机海面

  %    GENERATE 2D Pierson-Moskowitz Distribution

  %  Parameters
  clear;
  NX=512;    % number in x direction
  NY=512;    % number in y direction

  avg_x=0.1;  % sampling size in x
  avg_y=0.1;  % sampling size in y

%  NX=600;    % number in x direction
%  NY=600;    % number in y direction
%  avg_x=0.005;  % sampling size in x
%  avg_y=0.005;  % sampling size in y

  
  TOTAL_NUM=NX*NY;

  hrms=input('Input the rms height of the surface (suggest: 0.01 to 0.2) =   ');
  U=sqrt( sqrt(hrms*hrms/8.1e-3*4*0.74*9.81*9.81)  );
  fprintf('Equivalent Wind speed (corresponding to the rms) = %f(m/s)',U);
   
  %------generating white noise
    %randn('state',0);
    randn('state',sum(10*clock))
    noisetmp=randn(1,TOTAL_NUM);  
 % plot(noise)
    
    for j=1:NY
        for i=1:NX
            noise(j,i)=(noisetmp((j-1)*NX+i));
        end
    end    
 
    %----- statistic of the nosie
    rms=0.0;
    mean=0.0;

    for j=1:NY
        for i=1:NX
            mean=mean+noise(j,i);
        end
    end    
    mean=mean/TOTAL_NUM;
      
    for j=1:NY
        for i=1:NX
            rms=rms+(abs(noise(j,i)-mean) )^2;
        end
    end    
    rms=sqrt(rms/TOTAL_NUM);
    fprintf('\n')
  %  fprintf('Mean of the input noise=%f,  RMS=%f\n',mean,rms)
    %---spectrum of the input noise
    
    noise=noise/sqrt(TOTAL_NUM);
    spec_noise=fft2(noise);  % transform into spectrum domain
    
    delkx=2.0*pi/(avg_x*NX);
    delky=2.0*pi/(avg_y*NY);
    
    for j=0:NY-1
        for i=0:NX-1
            if(i<=NX/2) ii=i;
            else ii=i-NX;
            end
            if(j<=NY/2) jj=j;
            else jj=j-NY;
            end
           
            kx=ii*delkx;
            ky=jj*delky;
            K=sqrt(kx*kx+ky*ky);
            
            filterout(j+1,i+1)=spec_noise(j+1,i+1)*sqrt( Two_PM(U,K)*delkx/2.0/pi*delky/2.0/pi);
        end
    end    
  
        
    rough=ifft2(filterout);
    outsurface=real(rough);   %output surface
    outsurface=outsurface*TOTAL_NUM;
    surf(outsurface)
    shading interp
   %----- statistic of the output surface
    rms=0.0;
    mean=0.0;

    for j=1:NY
        for i=1:NX
            mean=mean+outsurface(j,i);
        end
    end    
    mean=mean/TOTAL_NUM;
    meanofroughsurface=mean;
    
    for j=1:NY
        for i=1:NX
            rms=rms+(abs(outsurface(j,i)-mean) )^2;
%            pause
        end
    end    
    
    rmsofroughsurface=sqrt(rms/TOTAL_NUM);
    fprintf('Output surface mean=%f,  RMS=%f\n',meanofroughsurface,rmsofroughsurface);
    fprintf('The Presumed RMS is %f\n',hrms);
%    fprintf('END OF GENERATING PERSON-MOSKWITZ DISTRIBUTION')