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
Fast display of wait bar. | |
returns the constant acceleration of gravity | |
Translates from wave number to frequency | |
Makes a movie of a 2D (x,t) or 3D (x,y,t) simulated sea | |
Numerical integration with the Simpson method | |
Calculates spectral moments from spectrum | |
Transforms between different types of spectra | |
Interpolation and zero-padding of spectrum | |
Transforms process X and up to four derivatives | |
Translates from frequency to wave number | |
Control axis scaling and appearance. | |
Clear variables and functions from memory. | |
Close figure. | |
Color look-up table. | |
Contour plot. | |
Difference and approximate derivative. | |
Display message and abort function. | |
Two-dimensional discrete Fourier Transform. | |
Shift zero-frequency component to center of spectrum. | |
Create figure window. | |
Get handle to current axis. | |
Get handle to current figure. | |
Inverse discrete Fourier transform. | |
Resample data at a higher rate using lowpass interpolation. | |
2-D interpolation (table lookup). | |
True if field is in structure array. | |
Convert string to lowercase. | |
X and Y arrays for 3-D plots. | |
Next higher power of 2. | |
Linear plot. | |
Remove fields from a structure array. | |
Set object properties. | |
Color shading mode. | |
Remove singleton dimensions. | |
Compare strings. | |
3-D shaded surface with lighting. | |
3-D graph viewpoint specification. | |
X-axis label. | |
Y-axis label. |
% CHAPTER1 demonstrates some applications of WAFO |
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
Comments or corrections to the WAFO group