Home > wafo > wsim > seasim.m

# seasim

## PURPOSE

Spectral simulation of a Gaussian sea, 2D (x,t) or 3D (x,y,t)

## SYNOPSIS

[Y,Mv]=seasim(spec,Nx,Ny,Nt,dx,dy,dt,fftdim,plotflag);

## DESCRIPTION

``` SEASIM Spectral simulation of a Gaussian sea, 2D (x,t) or 3D (x,y,t)

CALL:  [Y, Mv]=seasim(spec,Nx,Ny,Nt,dx,dy,dt,fftdim,plotflag)

Y    = output struct with simulated path
Mv   = movie of output (run with movie(Mv))
spec = spectrum struct, must be two-dimensional spectrum
Nx,Ny,Nt = number of points in grid. Put Ny=0 for 2D simulation
(default 2^7)
dx,dy,dt = steps in grid (default dx,dt by the Nyquist freq, dy=dx)
fftdim   = 1 or 2 gives one- or two dimensional FFT. (default 2)
If 1 then FFT over t and looping over x and y.
If 2 then FFT2 over x and y, looping over t.
plotflag = 0,1 or 2. If 0 no plot, if 2 movie (takes a couple of min.)
if 1 single (first) frame of movie (default 0)

The output struct contains:
.Z matrix of size [Ny Nx Nt], but with singleton dimensions removed
.x, .t (and .y if Ny>1)
For the removal of singleton dimensions, see SQUEEZE.

The output movie can also be created separately by calling SEAMOVIE.
Then can also different plot styles be chosen for the movie. The style
used here is a surf-plot.

Limitations: memory demanding! fftdim=2 is in general better.
When fftdim=1 and NX*NY is large (the limit also depends on spectrum, dt
etc.) then a slower version with more looping is used. This version can be
forced by putting fftdim to something greater than 2, then fftdim will equal
the number of extra loops. Try this after an interrupt by 'Out of Memory'.

If the spectrum has a non-empty field .tr, then the transformation is
applied to the simulated data, the result is a simulation of a transformed
Gaussian sea.

Example: % Simulate using demospec, plot the first frame
Y=seasim(demospec('dir'),2^7,2^7,100,10,10,1,2,1);

## CROSS-REFERENCE INFORMATION

This function calls:
 fwaitbar Fast display of wait bar. gravity returns the constant acceleration of gravity k2w Translates from wave number to frequency seamovie Makes a movie of a 2D (x,t) or 3D (x,y,t) simulated sea simpson Numerical integration with the Simpson method spec2mom Calculates spectral moments from spectrum spec2spec Transforms between different types of spectra specinterp Interpolation and zero-padding of spectrum tranproc Transforms process X and up to four derivatives w2k Translates from frequency to wave number axis Control axis scaling and appearance. clear Clear variables and functions from memory. close Close figure. colormap Color look-up table. contour Contour plot. diff Difference and approximate derivative. error Display message and abort function. fft2 Two-dimensional discrete Fourier Transform. fftshift Shift zero-frequency component to center of spectrum. figure Create figure window. gca Get handle to current axis. gcf Get handle to current figure. ifft Inverse discrete Fourier transform. interp Resample data at a higher rate using lowpass interpolation. interp2 2-D interpolation (table lookup). isfield True if field is in structure array. lower Convert string to lowercase. meshgrid X and Y arrays for 3-D plots. nextpow2 Next higher power of 2. plot Linear plot. rmfield Remove fields from a structure array. set Set object properties. shading Color shading mode. squeeze Remove singleton dimensions. strcmp Compare strings. surfl 3-D shaded surface with lighting. view 3-D graph viewpoint specification. xlabel X-axis label. ylabel Y-axis label.
This function is called by:
 Chapter1 % CHAPTER1 demonstrates some applications of WAFO

## SOURCE CODE

```001 function [Y,Mv]=seasim(spec,Nx,Ny,Nt,dx,dy,dt,fftdim,plotflag);
002 %SEASIM Spectral simulation of a Gaussian sea, 2D (x,t) or 3D (x,y,t)
003 %
004 % CALL:  [Y, Mv]=seasim(spec,Nx,Ny,Nt,dx,dy,dt,fftdim,plotflag)
005 %
006 %       Y    = output struct with simulated path
007 %       Mv   = movie of output (run with movie(Mv))
008 %       spec = spectrum struct, must be two-dimensional spectrum
009 %       Nx,Ny,Nt = number of points in grid. Put Ny=0 for 2D simulation
010 %                  (default 2^7)
011 %       dx,dy,dt = steps in grid (default dx,dt by the Nyquist freq, dy=dx)
012 %       fftdim   = 1 or 2 gives one- or two dimensional FFT. (default 2)
013 %                  If 1 then FFT over t and looping over x and y.
014 %                  If 2 then FFT2 over x and y, looping over t.
015 %       plotflag = 0,1 or 2. If 0 no plot, if 2 movie (takes a couple of min.)
016 %                  if 1 single (first) frame of movie (default 0)
017 %
018 % The output struct contains:
019 %          .Z matrix of size [Ny Nx Nt], but with singleton dimensions removed
020 %          .x, .t (and .y if Ny>1)
021 % For the removal of singleton dimensions, see SQUEEZE.
022 %
023 % The output movie can also be created separately by calling SEAMOVIE.
024 % Then can also different plot styles be chosen for the movie. The style
025 % used here is a surf-plot.
026 %
027 % Limitations: memory demanding! fftdim=2 is in general better.
028 % When fftdim=1 and NX*NY is large (the limit also depends on spectrum, dt
029 % etc.) then a slower version with more looping is used. This version can be
030 % forced by putting fftdim to something greater than 2, then fftdim will equal
031 % the number of extra loops. Try this after an interrupt by 'Out of Memory'.
032 %
033 % If the spectrum has a non-empty field .tr, then the transformation is
034 % applied to the simulated data, the result is a simulation of a transformed
035 % Gaussian sea.
036 %
037 % Example: % Simulate using demospec, plot the first frame
038 %          Y=seasim(demospec('dir'),2^7,2^7,100,10,10,1,2,1);
039 %
041
042 % Tested on Matlab 5.3
043 % revised pab June 2005
044 % added fwaitbar and removed disp statements
045 % revised jr 310301, added some status info, displayed when pg is running
046 % revised jr 180101, line 77, g -> spec.g
047 % revised es 200600, removed orient landscape
048 % revised es 130600, corrected two errors,
049 %                    added more plotting alternatives depending on dimension
050 % revised es 060600, added transformation of data
051 % By es 23.05.00
052
053 % Initialization
054 if ~isfield(spec,'g')
055   spec.g=gravity;
056 end
057 if strcmp(lower(spec.type(end-2:end)),'k2d')
058   spec=spec2spec(spec,'dir');
059 elseif ~strcmp(lower(spec.type(end-2:end)),'dir')
060   error('Spectrum must be two-dimensional')
061 end
062 if isfield(spec,'f')
063   spec.w=spec.f*2*pi;
064   spec.S=spec.S/2/pi;
065   spec=rmfield(spec,'f');
066 end
067 if nargin<2|isempty(Nx)
068   Nx=2^nextpow2(length(spec.w));
069 end
070 if nargin<3|isempty(Ny)
071   Ny=2^nextpow2(length(spec.w));
072 end
073 if Ny==0
074   Ny=1;
075 end
076 if nargin<4|isempty(Nt)
077   Nt=2^nextpow2(length(spec.w));
078 end
079 if nargin<5|isempty(dx)
080   dx=pi/(spec.w(end)^2/spec.g)/2;
081 end
082 if nargin<6|isempty(dy)
083   dy=dx;
084 end
085 if nargin<7|isempty(dt)
086   dt=pi/spec.w(end); %directly by the Nyquist frequency
087 end
088 if nargin<8|isempty(fftdim)
089   fftdim=1;
090 end
091 if nargin<9|isempty(plotflag)
092   plotflag=0;
093 end
094 t=(0:Nt-1)'*dt;
095 x=(0:Nx-1)'*dx;
096 y=(0:Ny-1)'*dy;
097
098 if fftdim~=2
099   spec=specinterp(spec,dt);
100   nf= length(spec.w); %number of frequencies
101   np= length(spec.theta); % number of directions
102   Nxy=Nx*Ny;
103   Nxyp=2^10;   % limit for Nxy to do in one step
104   [Xi,Yi]=meshgrid(x,y);
105   %Normalize to get variance back after discretization
106   Sdiscr=(spec.S)*spec.w(end)/nf*diff(spec.theta([1,end]))/np; % size np x nf
107   Sdiscr=(randn(np,nf)+i*randn(np,nf)).*sqrt(Sdiscr);
108   [k1,k2]=w2k(spec.w,spec.theta,spec.h,spec.g);
109   Sxy=zeros(nf,Nxy);
110
111   h = fwaitbar(0,[],' Integration ... ');
112   for ix=1:Nxy,
113     Sxy(:,ix)=simpson(0:np-1,Sdiscr.*exp(i*(Xi(ix)*k1+Yi(ix)*k2))).';
114     fwaitbar(ix/Nxy,h)
115   end
116   close(h)
117   clear Sdiscr k1 k2
118   nfft=2^nextpow2(2*nf-2);
119   Sxy=[Sxy; zeros(nfft-(2*nf)+2,Nxy); conj(Sxy(nf-1:-1:2,:))];
120   % size(Z) is Nxy x nfft
121   Y.Z=zeros(Nxy,nfft);
122   if (Nxy<Nxyp)|fftdim>2
123     Y.Z=real(ifft(Sxy,nfft,1)).'*nfft/(2*nf-2)*nf;
124   else %do it stepwise
125     if fftdim>2 % new Nxyp to get fftdim steps
126       Nxyp=ceil(Nxy/fftdim);
127     end
128     nr=ceil(Nxy/Nxyp);
129
130     h = fwaitbar(0,[],'  ... ');
131     for j=1:nr
132       Y.Z((j-1)*Nyxp+1:min(j*Nyxp,Nxy),:)=real(ifft(Sxy(:,(j-1)*Nxyp+1:min(j*2^10,Nxy)),nfft,1)).'* ...
133         nfft/(2*nf-2)*nf;
134       fwaitbar(j/nr,h)
135     end
136     close(h)
137   end
138   Y.Z=squeeze(reshape(Y.Z(1:Nxy,1:Nt),Ny,Nx,Nt));
139 else  % if fftdim ...
140   nfftx=2^nextpow2(max(Nx,length(spec.w)));
141   nffty=2^nextpow2(max(Ny,length(spec.w)));
142   if nffty<2
143     nffty=2;
144   end
145   dk1=2*pi/nfftx/dx; % necessary wave number lag to satisfy input space lag
146   dk2=2*pi/nffty/dy; % necessary wave number lag to satisfy input space lag
147   S=spec2spec(spec,'k2d');
148   if S.k(1)>=0
149     S.S=[zeros(size(S.S,1),size(S.S,2)-1) S.S];
150     S.k=[-S.k(end:-1:2) S.k];
151   end
152   % add  zeros just above old max-freq, and a zero at new max-freq
153   % to get non-NaN interpolation
154   S.S=[zeros(2,size(S.S,1)+4); zeros(size(S.S,1),2) S.S ...
155        zeros(size(S.S,1),2);zeros(2,size(S.S,1)+4)];
156   dk1old=S.k(2)-S.k(1);
157   dk2old=S.k2(2)-S.k2(1);
158   S.k=[min(dk1*(-nfftx/2), S.k(1)-2*dk1old) S.k(1)-dk1old S.k S.k(end)+dk1old max(dk1*(nfftx/2-1),S.k(end)+2*dk1old)];
159   S.k2=[min(dk2*(-nffty/2),S.k2(1)-2*dk2old); S.k2(1)-dk2old; S.k2; S.k2(end)+dk2old; max(dk2*(nffty/2-1),S.k2(end)+2*dk2old)];
160
161   % Interpolate in spectrum to get right frequency grid to satisfy input
162   disp(' Interpolating in spectrum')
163   S.S=interp2(S.k,S.k2,S.S,dk1*(-nfftx/2:nfftx/2-1),dk2*(-nffty/2:nffty/2-1)');
164   Sdiscr=S.S*dk1*dk2;
165
166   if sum(Sdiscr(:))<0.9*spec2mom(spec,0)
167     disp('WARNING: Too small dx and/or dy, or too small NX and/or NY')
168     disp('         Information in the spectrum is lost')
169     disp('         May be a good idea to interrupt')
170   end
171
172   Sdiscr=fftshift(Sdiscr);
173   % Simulation
174   Z0=sqrt(Sdiscr).*randn(nffty,nfftx)+i*sqrt(Sdiscr).*randn(nffty,nfftx);
175   W=k2w((-nfftx/2:nfftx/2-1)*dk1,(-nffty/2:nffty/2-1)*dk2,S.h,S.g);
176   W(:,nfftx/2+1:nfftx)=-W(:,nfftx/2+1:nfftx);
177   W=fftshift(W);
178   Y.Z=zeros(nffty,nfftx,Nt);
179   %disp(' Loop over t values')
180   h = fwaitbar(0,[],'Loop over t values...');
181   nt = length(t);
182   for j=1:nt
183     Y.Z(:,:,j)=real(fft2(Z0.*exp(i*W*t(j))));
184      fwaitbar(j/nt,h)
185   end
186   close(h)
187   Y.Z=squeeze(Y.Z(1:Ny,1:Nx,:));
188   clear Z0 W Sdiscr
189
190 end % if fftdim ...
191
192 if Nx>1
193   Y.x=x;
194 end
195 if Ny>1
196   Y.y=y;
197 end
198 if Nt>1
199   Y.t=t;
200 end
201
202 % Transformation of data, if given
203 if isfield(spec,'tr') & ~isempty(spec.tr)
204   disp('   Transforming data.')
205   G=fliplr(spec.tr); % the inverse of the transformation
206   Y.Z(:)=tranproc(Y.Z(:),G);
207 end
208
209 % Plotting if requested
210 Mv=[];
211 if plotflag==2
212   Mv=seamovie(Y,1);
213 elseif plotflag==1
214   figure(gcf)
215   if min(Nx,Ny)>2
216     colormap('winter')
217     surfl(Y.x,Y.y,Y.Z(:,:,1),[-30, 45]);
219     view(-37.5,20)
220     axis([Y.x(1) Y.x(end) Y.y(1) Y.y(end) 7*min(Y.Z(:)) 7*max(Y.Z(:))])
221     set(gca,'xtick',[])
222     set(gca,'ytick',[])
223     axis('square')
224     axis('off')
225   elseif Nt>1
226     if Nx>1
227       contour(Y.t,Y.x,Y.Z,[0 0])
228     else
229       contour(Y.t,Y.y,Y.Z,[0 0])
230     end
231     xlabel('[s]')
232     ylabel('[m]')
233   else
234     if Nx>1
235       plot(Y.x,Y.Z)
236     elseif Ny>1
237       plot(Y.y,Y.Z)
238     end
239     xlabel('[m]')
240   end
241 end
242```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003