sea_elevation(x,U);

function [s,fm,Tp,SK,k] = sea_elevation(x,U);

%
% SEA_ELEVATION: generates one-dimensional sea surface realizations for 
% a given value of wind speed v (uses the Pierson-Moskowitz spectrum). 
% 
%
% Usage:  [s,fm,Tp,p,k] = sea_elevation(x,v)
% 
%  where x is a vector defining the surface coordinates, v is the intensity of 
%  the wind speed. The output is an elevation vector s (length(s) = length(x)), 
%  the peak period Tp, the peak frequency fm, the spectrum Sk in the wavenumber 
%  domain and its argument k. The Pierson-Moskowitz spectrum is valid for v >= 10. 
% 
%  Example: 
%  n = 201; xmin = 0; xmax = 200; x = linspace(xmin,xmax,n); 
%  wind_velocity = 10; 
%  [s,fm,Tp,p,k] = sea_elevation(x,wind_velocity);
%  figure(1)
%  subplot(211),plot(x,s),ylabel('Sea surface elevation (m)'),xlabel('x(m)'),box on, grid on
%  subplot(212),plot(k,p),ylabel('Power Spectrum S(k)'),xlabel('k(1/m)'),box on, grid on 
% 
% References:
% 1) Fourier Synthesis of Ocean Scenes 
%    Gary A. Mastin et al. 1987 
% 2) Acoustic wave scattering from rough sea surface and sea bed 
%    Chen-Fen Huang, Master Thesis. 1998 
% 3) Tutorial 2: Ocean Waves (1)
%    http://www.naturewizard.com 
% 4) matlabwaves.zip 
%    http://neumeier.perso.ch/matlab/waves.html

%***************************************************************************************
% First version: 29/07/2008
% 
% Contact: [e-mail]
% 
% Any suggestions to improve the performance of this 
% code will be greatly appreciated. 
% 
%***************************************************************************************

s  = []; 
SK = []; 
Tp = []; 
k  = [];

imunit = sqrt( -1 ); 

g     = 9.80665; 
gg    = g*g;
alpha = 8.1e-3;

UU = U*U; 

n = length( x ); 

fm = 0.13*g/U;
Tp = 1/fm;    

%======================================================================
% Surface generation: 
% most expressions for the spectrum are given for the frequency domain;
% therefore, let us convert wavenumbers to frequencies, apply the 
% formulas and go back to the wavenumber domain:

dx = x(2) - x(1);
kmax = 1/( 2*dx );
k = linspace(-kmax,kmax,n);
K = abs( k );             % A real surface requires a symmetric spectrum...
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 );

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

SF = alpha*gg./( (2*pi)^4*F.^5 ).*exp( -5/4*( fm./F ).^4 );

%======================================================================   
%  Convert spectrum from frequency domain to wavenumber domain:

SK = SF.*dFdK;

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

white_noise     = unifrnd(-127,127,1,n)/127;
WHITE_NOISE     = fft( 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 surface realization through an inverse fft:  
s = ifft( filtered_NOISE );
s = real( s );