function waterwave
% WATER WAVE
% 2D Shallow Water Model
%
% Lax-Wendroff finite difference method.
% Reflective boundary conditions.
% Random water drops initiate gravity waves.
% Surface plot displays height colored by momentum.
% Plot title shows t = simulated time and tv = a measure of total variation.
% An exact solution to the conservation law would have constant tv.
% Lax-Wendroff produces nonphysical oscillations and increasing tv.
%
% See:
%    http://en.wikipedia.org/wiki/Shallow_water_equations
%    http://www.amath.washington.edu/~rjl/research/tsunamis
%    http://www.amath.washington.edu/~dgeorge/tsunamimodeling.html
%    http://www.amath.washington.edu/~claw/applications/shallow/www


% Parameters
clf
clear

n = (100);               % grid size
g = 9.8;                 % gravitational constant
dt = 0.09;               % hardwired timestep
dx = 1.0;
dy = 1;
nplotstep = 8;           % plot interval
ndrops = 1;              % maximum number of drops
dropstep = 500;          % drop interval


G = zeros(n+2,n+2);
for i=1:n+2 
  for j=1:n+2
      G(i,j) = 0.5; % Enter you ocean floor elevation here
  end
end

% Initialize graphics

[surfplot,top] = initgraphics(n);

H = ones(n+2,n+2);   U = zeros(n+2,n+2);  V = zeros(n+2,n+2);
Hx = zeros(n+1,n+1); Ux = zeros(n+1,n+1); Vx = zeros(n+1,n+1);
Hy = zeros(n+1,n+1); Uy = zeros(n+1,n+1); Vy = zeros(n+1,n+1);
ndrop = ceil(rand*ndrops);
nstep = 0;

% Set inital conditions
for i=2:n+1 
   for j=2:n+1
      H(i,j) = H(i,j)+7*exp(-0.2*i^2);
   end
end

% loop over time steps.
while (nstep<12000)   

   nstep = nstep + 1
      
   for i=1:1000
   if nstep==i*100
   
 for i=2:n+1 
    for j=2:n+1
       H(i,j) = H(i,j)+7*exp(-0.2*i^2);
    end
 end

   end
   end
      
   
  for i=1:1000
   if nstep==i*100-3
    for i=1:(.24*n);
    for j=2:n+1;
       U(i,j) = 0;
       V(i,j) = 0;
       H(i,j) = 1;
end
end
end
end 


  for i=1:1000
   if nstep==i*100-40
    for i=.25*n:(.4*n);
    for j=.10*n:.35*n
       U(i,j) = 0;
       V(i,j) = 0;
       H(i,j) = 1;
end
end
end
end
   
   
  for i=1:1000
   if nstep==i*100-40
    for i=.25*n:(.4*n);
    for j=.65*n:.9*n
       U(i,j) = 0;
       V(i,j) = 0;
       H(i,j) = 1;
end
end
end
end
   
   
if mod(i+174,175) == 0  
        f=f+1;
        surf(X,Y,T)
        axis([0 1 0 1 -1 3])
        saveas(gcf, ['wave_propagation',sprintf('%03d',f),'.png'])
        pause(0.01)
end   
   

   % Reflective boundary conditions
  
   H(:,n+2) = 1;  U(:,n+2) = 0;   V(:,n+2) = 0;
   H(1,:) = 1;      U(1,:) = 0;      V(1,:) = 0;
   H(n+2,:) = 1;  U(n+2,:) = 0;  V(n+2,:) = 0;
   
     
       H(.15*n,1:.45*n) = 1;
       
       U(.15*n,1:.45*n) = 0;  
       
       V(.15*n,1:.45*n) = 0;
      
       
       %%%%%%%%%%%%
       
       H(.15*n,.55*n:n) = 1;
       
       U(.15*n,.55*n:n) = 0;  
       
       V(.15*n,.55*n:n) = 0;
       
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
       H(.45*n,1:.25*n) = 1;
       
       U(.45*n,1:.25*n) = 0;  
       
       V(.45*n,1:.25*n) = 0;
       
     
       %%%%%%%%%%%%
       H(.45*n,.3*n:.7*n) = 1;
       
       U(.45*n,.3*n:.7*n) = 0;  
       
       V(.45*n,.3*n:.7*n) = 0;
       
       %%%%%%%%%%%%    
       
       H(.45*n,.75*n:n) = 1;
       
       U(.45*n,.75*n:n) = 0;  
       
       V(.45*n,.75*n:n) = 0;

       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
       H(.7*n,1:.2*n) = 1;
       
       U(.7*n,1:.2*n) = 0;  
       
       V(.7*n,1:.2*n) = 0;
       
     
       %%%%%%%%%%%%
       H(.7*n,.25*n:.45*n) = 1;
       
       U(.7*n,.25*n:.45*n) = 0;  
       
       V(.7*n,.25*n:.45*n) = 0;
       
       %%%%%%%%%%%%    
       
       H(.7*n,.5*n:.7*n) = 1;
       
       U(.7*n,.5*n:.7*n) = 0;  
       
       V(.7*n,.5*n:.7*n) = 0;
       
       %%%%%%%%%%%%
       
       H(.7*n,.75*n:n) = 1;
       
       U(.7*n,.75*n:n) = 0;  
       
       V(.7*n,.75*n:n) = 0;
       
       
   % First half step

   % x direction
   i = 1:n+1;
   j = 1:n;

   Horg = H;
   H = H - G;
   % height
   Hx(i,j) = (H(i+1,j+1)+H(i,j+1))/2 - dt/(2*dx)*(U(i+1,j+1)-U(i,j+1));
   Gx(i,j) = (G(i+1,j+1)+G(i,j+1))/2;

   % x momentum
   Ux(i,j) = (U(i+1,j+1)+U(i,j+1))/2 -  ...
             dt/(2*dx)*((U(i+1,j+1).^2./H(i+1,j+1) + g.*(H(i,j+1)+H(i+1,j+1))./2.*Horg(i+1,j+1)) - ...
                        (U(i  ,j+1).^2./H(i  ,j+1) + g.*(H(i,j+1)+H(i+1,j+1))./2.*Horg(i  ,j+1)));

   % y momentum
   Vx(i,j) = (V(i+1,j+1)+V(i,j+1))/2 - ...
             dt/(2*dx)*((U(i+1,j+1).*V(i+1,j+1)./H(i+1,j+1)) - ...
                        (U(i,j+1).*V(i,j+1)./H(i,j+1)));

   % y direction
   i = 1:n;
   j = 1:n+1;

   % height
   Hy(i,j) = (H(i+1,j+1)+H(i+1,j))/2 - dt/(2*dy)*(V(i+1,j+1)-V(i+1,j));
   Gy(i,j) = (G(i+1,j+1)+G(i+1,j))/2;

   % x momentum
   Uy(i,j) = (U(i+1,j+1)+U(i+1,j))/2 - ...
             dt/(2*dy)*((V(i+1,j+1).*U(i+1,j+1)./H(i+1,j+1)) - ...
                        (V(i+1,j).*U(i+1,j)./H(i+1,j)));
   % y momentum
   Vy(i,j) = (V(i+1,j+1)+V(i+1,j))/2 - ...
             dt/(2*dy)*((V(i+1,j+1).^2./H(i+1,j+1) + g.*(H(i+1,j+1)+H(i+1,j))./2.*Horg(i+1,j+1)) - ...
                        (V(i+1,j  ).^2./H(i+1,j  ) + g.*(H(i+1,j+1)+H(i+1,j))./2.*Horg(i+1,j  )));

   % Second half step
   i = 2:n+1;
   j = 2:n+1;

   % height
   H(i,j) = H(i,j) - (dt/dx)*(Ux(i,j-1)-Ux(i-1,j-1)) - ...
                     (dt/dy)*(Vy(i-1,j)-Vy(i-1,j-1));
   % x momentum
   U(i,j) = U(i,j) - (dt/dx)*((Ux(i  ,j-1).^2./Hx(i  ,j-1) + g.*(Hx(i-1,j-1)+Hx(i,j-1))./2.*(Hx(i  ,j-1)+Gx(i  ,j-1)) ) - ...
                              (Ux(i-1,j-1).^2./Hx(i-1,j-1) + g.*(Hx(i-1,j-1)+Hx(i,j-1))./2.*(Hx(i-1,j-1)+Gx(i-1,j-1)))) ...
                   - (dt/dy)*((Vy(i-1,j  ).*Uy(i-1,j  )./Hy(i-1,j  )) - ...
                              (Vy(i-1,j-1).*Uy(i-1,j-1)./Hy(i-1,j-1)));
   % y momentum
   V(i,j) = V(i,j) - (dt/dx)*((Ux(i  ,j-1).*Vx(i  ,j-1)./Hx(i  ,j-1)) - ...
                              (Ux(i-1,j-1).*Vx(i-1,j-1)./Hx(i-1,j-1))) ...
                   - (dt/dy)*((Vy(i-1,j  ).^2./Hy(i-1,j  ) + g.*(Hy(i-1,j)+Hy(i-1,j-1))./2.*(Hy(i-1,j  )+Gy(i-1,j  )) ) - ...
                              (Vy(i-1,j-1).^2./Hy(i-1,j-1) + g.*(Hy(i-1,j)+Hy(i-1,j-1))./2.*(Hy(i-1,j-1)+Gy(i-1,j-1))));

   H = H + G;



   % Update plot
   if mod(nstep,nplotstep) == 0
      C = abs(U(i,j)) + abs(V(i,j));  % Color shows momemtum
      C(.25*n:n,:) = 2.5 .* C(.25*n:n,:);
      t = nstep*dt;
      tv = norm(C,'fro');
      set(surfplot,'zdata',H(i,j),'cdata',C);
      set(top,'string',sprintf('t = %6.2f,  tv = %6.2f',t,tv))
      drawnow
   end
   
   if all(all(isnan(H))), break, end  % Unstable, quit
end

%close(gcf)

function [surfplot,top] = initgraphics(n);
% INITGRAPHICS  Initialize graphics for waterwave.
% [surfplot,top,start,stop] = initgraphics(n)
% returns handles to a surface plot, its title, and two uicontrol toggles.
    
   clf
   shg
   set(gcf,'numbertitle','off','name','Wave_Diffraction_in_Water')
   x = (0:n-1)/(n-1);
   surfplot = surf(x,x,ones(n,n),zeros(n,n));
   grid off
   axis([0 1 0 1 -1 3])
   caxis([-1 1])
   shading faceted
   c = (1:64)'/64;
   cyan = [0*c c c];
   colormap(cyan)
   top = title('Click restart');
   