function mit18086_shallowwater
%MIT18086_SHALLOWWATER
%    Solves the shallow water equations by a particle method.

% 03/2007 by Benjamin Seibold
% Feel free to modify for teaching and learning.

n = 200;                            % number of space gridpoints
dt = 2e-3;                                           % time step
tf = 5e-1;                                          % final time

x = linspace(-2,2,n)';
U = [x*0+1,sin(12*x).*exp(-3*x.^2)+2];
ships = linspace(-1.3,.9,4)';
for tn = 1:ceil(tf/dt)
   x = x+dt*U(:,2);
   [mv,mi] = min(diff(x));
   while abs(mv)<.8/n
      x = [x(1:mi-1);(x(mi)+x(mi+1))/2;x(mi+2:end)];
      U = [U(1:mi-1,:);(U(mi,:)+U(mi+1,:))/2;U(mi+2:end,:)];
      [mv,mi] = min(diff(x));
   end
   [mv,mi] = max(diff(x));
   while abs(mv)>8/n
      x = [x(1:mi);(x(mi)+x(mi+1))/2;x(mi+1:end)];
      U = [U(1:mi,:);(U(mi,:)+U(mi+1,:))/2;U(mi+1:end,:)];
      [mv,mi] = max(diff(x));
   end
   DxcU = (U([2:end end],:)-U([1 1:end-1],:))./...
         ((x([2:end end])-x([1 1:end-1]))*[1 1]);
   U = U-dt*[U(:,1).*DxcU(:,2),9.81*DxcU(:,1)];
   DxU = (U(2:end,:)-U(1:end-1,:))./((x(2:end)-x(1:end-1))*[1 1]);
   DxxU = 2*(DxU(2:end,:)-DxU(1:end-1,:))./((x(3:end)-x(1:end-2))*[1 1]);
   U = U+3e-5*[0 0;DxxU;0 0];
   Us = interp1(x,U,ships);
   ships = ships+dt*Us(:,2);
   clf
   patch([-1;x;2],[0;U(:,1);0],'b')
   hold on
%   plot(x,U(:,1),'k.')
   for k = 1:length(ships)
      patch(ships(k)+[-2 -1 1 2 .1 .1 -.2 -.2]*.03,...
             Us(k,1)+[1 0 0 1 1 2 2 1]*.03,'r')
   end
   hold off
   axis([-1 2 0 1.6])
   title(sprintf('shallow water equations  t=%0.2f',tn*dt))
   drawnow
end