test_adv.m 1.84 KB
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.7;
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));



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
      
    E(n,:,:)=reshape(EE,1,nbin,length(om));
    
    
  
end



figure(1)
   rr = waitforbuttonpress;
for n=1:nsteps
     figure(1)
    
   
     E3=reshape(E(n,:,:),nbin,length(om));
    h=mesh(om,x,E3);
    axis([om1 om2 min(x) max(x) 0 max(Ei)])
         xlabel('Frequency [s^{-1}]')
        zlabel('Energy')
   
    ylabel('x [km]')
  pause(0.1)
   
    
end


for i=1:nsteps
    
E2=reshape(E(i,:,:),nbin,length(om));
m0(i)=sum(sum(E2))/sum(Ei);
end
figure
plot(m0)