水波仿真

function [s,Tp,fm,B,SK,kx,ky] = sea_surface(x,y,wind_data,type,spreading);

%
% SEA_SURFACE: generates sea surface realizations for a given intensity,
% fetch, and direction of wind velocity.
%
% Usage: [s,Tp,fm,B,Sk,kx,ky] = sea_surface(x,y,wind_data,type,spreading)
%
% where x and y are vectors defining the surface grid, wind_data is a structure containing
% the intensity, direction and fetch of the wind speed, type describes the spectrum and spreading
% is a string defining the angular spreading of the sea surface spectrum
% ('none' for no spreading, 'cos2' for cosine-squared spreading, 'mits'
% for Mitsuyasu spreading and 'hass' for Hasselmann spreading). The
% output is a matrix s (size(s) = [length(y) length(x)]), the peak period
% Tp, the peak frequency fm, the sea state B, the spectrum Sk and the wavenumber
% vector arguments kx and ky, which define the grid of Sk.
%
% Examples:
% nx = 101; xmin = 0; xmax = 100; x = linspace(xmin,xmax,nx);
% ny = 51; ymin = 0; ymax = 50; y = linspace(ymin,ymax,ny);
% wind_data.U = 10; wind_data.thetaU = 0; wind_data.X = 1e6;
% [s,Tp,fm,B,Sk,kx,ky] = sea_surface(x,y,wind_data,'PM','none');
% figure(1)
% subplot(211),mesh(x,y,s),ylabel('y(m)'),xlabel('x(m)')
% subplot(212),mesh(kx,ky,Sk),ylabel('ky (1/m)'),xlabel('kx (1/m)')
% [s,Tp,fm,B,Sk,kx,ky] = sea_surface(x,y,wind_data,'PM','cos2');
% figure(1)
% subplot(211),mesh(x,y,s),ylabel('y(m)'),xlabel('x(m)')
% subplot(212),mesh(kx,ky,Sk),ylabel('ky (1/m)'),xlabel('kx (1/m)')
% [s,Tp,fm,B,Sk,kx,ky] = sea_surface(x,y,wind_data,'PM','mits');
% figure(1)
% subplot(211),mesh(x,y,s),ylabel('y(m)'),xlabel('x(m)')
% subplot(212),mesh(kx,ky,Sk),ylabel('ky (1/m)'),xlabel('kx (1/m)')
% [s,Tp,fm,B,Sk,kx,ky] = sea_surface(x,y,wind_data,'PM','hass');
% figure(1)
% subplot(211),mesh(x,y,s),ylabel('y(m)'),xlabel('x(m)')
% subplot(212),mesh(kx,ky,Sk),ylabel('ky (1/m)'),xlabel('kx (1/m)')
%
% References:
% 1) Directional wave spectra observed during JONSWAP 1973
% D. E. Hasselmann et al. 1980
% 2) Directional wave spectra using cosine-squared and cosine 2s
% spreading functions
% Coastal Engineering Technical Note 1985
% 3) Fourier Synthesis of Ocean Scenes
% Gary A. Mastin et al. 1987
% 4) The generation of a time correlated 2d random process for ocean
% wave motion
% L. M. Linnet et al. 1997
% 5) Acoustic wave scattering from rough sea surface and sea bed
% Chen-Fen Huang, Master Thesis. 1998
% 6) Tutorial 2: Ocean Waves (1)
% http://www.naturewizard.com
% 7) matlabwaves.zip
% http://neumeier.perso.ch/matlab/waves.html

%***************************************************************************************
% First version: 30/07/2008
% First update : 14/10/2009 => match with Beaufort scale
% Second update : 25/01/2010 => returns sea state
% Third update : 30/03/2010 => JONSWAP spectrum
%
% Contact: orodrig@ualg.pt
%
% Any suggestions to improve the performance of this
% code will be greatly appreciated.
%
%***************************************************************************************
s = [];
Sk = s;
Tp = s;
fm = s;
B = s;
kx = s;
ky = s;

U = wind_data.U;
thetaU = wind_data.thetaU;

if U == 0
s = zeros(ny,nx);
return
end

B = fix( ( U/0.836 )^(2/3) );

imunit = sqrt( -1 );

g = 9.80665;
gxg = g*g;
UxU = U*U;
alpha = 8.1e-3;

nx = length( x );
ny = length( y );

%======================================================================
% Surface generation:
% since most expressions for the spectrum are given in the frequency
% domain we need to convert wavenumbers to frequencies, apply the formulas
% and go back to the wavenumber domain:

dx = x(2) - x(1);
kxmax = 1/( 2*dx );
kx = linspace(-kxmax,kxmax,nx);
dy = y(2) - y(1);
kymax = 1/( 2*dy );
ky = linspace(-kymax,kymax,ny);
[Kx,Ky] = meshgrid(kx,ky);
K = sqrt( Kx.^2 + Ky.^2 );
F = sqrt( g*K )/( 2*pi ); % Valid for surface waves over deep oceans
F( F == 0 ) = Inf ;
K( K == 0 ) = Inf ;
dFdK = sqrt( g./K )/( 4*pi );
OMEGA = 2*pi*F;

%======================================================================
% Calculate the spectrum in the frequency domain:

if type == 'PM'
fm = 0.13*g/U;
Tp = 1/fm;
SF = alpha*gxg/( (2*pi)^4 )*( F.^(-5) ).*exp( -5/4*( fm./F ).^4 );
elseif type == 'JS'
X = wind_data.X;
OMEGAp = 7*pi*( g/U )*( g*X/UxU )^(-0.33);
fm = OMEGAp/( 2*pi );
Tp = 2*pi/OMEGAp;
cgamma = 3.3;
csigma0 = 0.07*ones( size( OMEGA ) );
indexes = find( OMEGA > OMEGAp );
csigma0( indexes ) = 0.09;
calpha = 0.076*( g*X/U )^(-0.22);
delta = exp( -( ( OMEGA - OMEGAp ).^2 )./( 2*csigma0.^2*OMEGAp^2 ) );
SF = calpha*gxg*OMEGA.^(-5).*( exp( -(5/4)*( OMEGA/OMEGAp ).^(-4) ) ).*( cgamma.^( delta ) );
else
disp('Unknown spectrum type...')
end
%======================================================================
% Convert spectrum from frequency domain to wavenumber domain:

SK = SF.*dFdK;

%======================================================================
% A real sea surface requires a symmetric spectrum in the wavenumber
% domain; thus, wherever required, additional calculations will ensure
% that the spreading matrix is indeed symmetrical:

THETA = angle( Kx+imunit*Ky ) - thetaU;

switch spreading

case 'none' % no spreading

D = ones(size(K));

case 'cos2' % cosinus-squared spreading

D = cos( THETA ).^2;

case 'mits' % Mitsuyasu spreading

ss = 9.77*( F/fm ).^(-2.5);
indexes1 = ( F < fm );
ss( indexes1 ) = 6.97*( F(indexes1)/fm ).^5;
Nss = gamma( ss + 1 )./( 2*sqrt(pi)*gamma( ss + 0.5 ) );
D = ( ( cos( THETA/2 ).^2 ).^(ss) ).*Nss;
D = D + fliplr( flipud( D ) );

case 'hass' % Hasselmann spreading

Mu = 4.06*ones( size(K) );
indexes1 = ( F > fm );
Mu( indexes1 ) = -2.34;
pp = 9.77*( F/fm ).^Mu;
Npp = pi*2.^( 1 - 2*pp ).*gamma( 2*pp + 1 )./( gamma( pp + 1 ) ).^2;
D = cos( THETA/2 ).^(2*pp)./Npp;
D = ( ( cos( THETA/2 ).^2 ).^pp )./Npp;
D = D + fliplr( flipud( D ) );

otherwise

disp('Unknown sea surface spreading.')

end

%======================================================================
% Get the power spectrum:

D = D/max( D(:) )*2/pi; % spreading normalization
SK = SK.*D; % power spectrum(k,theta) = spectrum(k)*spreading(theta)

%======================================================================
% Get the surface realization from the spectrum:

white_noise = unifrnd(-127,127,ny,nx)/127;
WHITE_NOISE = fft2( white_noise );
NOISE_amplitude = abs( WHITE_NOISE );

NOISE_energy = sum( WHITE_NOISE(:).^2 );
WHITE_NOISE = WHITE_NOISE/NOISE_energy;

centered_WHITE_NOISE = fftshift( WHITE_NOISE );
NOISE_phase = angle( centered_WHITE_NOISE );

% Modulate noise amplitude with the power spectrum:
NOISE_amplitude = NOISE_amplitude .* SK;

% Randomize modulated noise in the wavenumber space combining
% modulated amplitudes with original phases:
NOISE_ipart = NOISE_amplitude .* sin( NOISE_phase );
NOISE_rpart = NOISE_amplitude .* cos( NOISE_phase );

filtered_NOISE = NOISE_rpart + imunit*NOISE_ipart;
filtered_NOISE = fftshift( filtered_NOISE );

% Get the 2D surface through an inverse fft:
s = ifft2( filtered_NOISE );
s = real( s );
average_height = ( max(s(:)) - min(s(:)) )/2;

%Beaufort scale (according to Wikipedia):
velocities = [0 0.3 1.5 3.3 5.5 8.0 11.0 14.0 17.0 20.0 24 28.0 32];
wave_heights = [0 0 0.2 0.5 1.0 2.0 3.0 4.0 5.5 7.5 10 12.5 16];

if U > max(velocities)
wave_height = 16; % I guess this really should be a very large wave...
elseif U <= 0.3;
wave_height = 0;
else
wave_height = interp1(velocities,wave_heights,U);
end

if average_height > 0
s = wave_height*s/average_height;
else
s = 0*s;
end