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);  
 
  See also  seamovie, spec2sdat

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %
040 % See also  seamovie, spec2sdat
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]);
218     shading interp
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