clear all; close all; % Parameters g = 9.81; % gravitational acceleration dx = 500; % spatial resolution nbin=80; x=linspace(0,dx*nbin/1000,nbin); Hs=1; Tm=6; % Waves fmin = 1/20; % minimum wave frequency fmax = 1/2.5; % maximum wave frequency om1 = 2*pi*fmin; % minimum wave radial frequency om2 = 2*pi*fmax; % maximum wave radial fequency nw = 61; % number of frequency bins dw = (om2-om1)/(nw-1); % integral interval for wave radial frequencies om = om1+(0:nw-1)'*dw; % wave radial frequencies vector T = 2*pi./om; % wave periods vector wlng = g.*T.^2./(2.*pi); % wavelength as a function only of wave period cp = sqrt(g.*wlng./(2.*pi)); % phase speed cg = cp./2; % group speed cgmax = max(cg); % group speed maximum cg(:) = cgmax; % no dispersion = all group speed are the same (maximum) CN=0.5; dt = CN*dx/cgmax; % time interval (temporal resolution) nsteps=ceil(nbin/CN); Ei = jonswap(Tm,Hs,om); % JONSWAP spectrum E=zeros(nsteps,nbin,length(om)); E(1,1,:)=Ei; EE=reshape(E(1,:,:),nbin,length(om)); figure(1) rr = waitforbuttonpress; for n=2:nsteps for w=1:nw % Advection loop over each frequency EE(:,w) = advection(EE(:,w),cg(w,:),dx,dt); % SEE advection routine end % Incident wave spectrum %EE(1,:) = Ei; figure(1) h=mesh(om,x,EE); axis([om1 om2 min(x) max(x) 0 max(Ei)]) xlabel('Frequency [s^{-1}]') zlabel('Energy') ylabel('x [km]') pause(0.1) E(n,:,:)=reshape(EE,1,nbin,length(om)); end for i=1:nsteps E2=reshape(E(i,:,:),nbin,length(om)); m0(i)=sum(sum(E2))/sum(Ei); end figure plot(m0)