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

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)