%http://www.math.lsa.umich.edu/~psmereka/
% Two dimensional level method for normal speed = vn > 0
% Based on Godunov's method
%
clear all
close all
N=200;         % number of grid points in one direction
h=2/(N-1);     % grid spacing
dt=h/5;       % time step
tfin=.6;       % total simulation time
nit=tfin/dt;   % number of time steps
vn=@(x,y)abs(sin(x));           % normal speed this must be positive
x=-1:h:1;
y=x;
[X,Y]=meshgrid(x);
%  Circle of radius R centered at (x0,y0)
%x0=.1;
%y0=.3;
%R=.3;
%phi=((X-x0).*(X-x0)+(Y-y0).*(Y-y0)).^.5-R;
% %
% %  Ellipse with major and minor axes a and b centered at (x0,y0)
% a=.6;
% b=.1;
% x0=.2;
% y0=.2;
% phi=(1./(a*a))*(X-x0).*(X-x0)+(1./(b*b))*(Y-y0).*(Y-y0)-1;
% %  Rotated Ellipse
% %  Ellipse with major and minor axes a and b centered at (x0,y0)
% %  and rotated about (x0,y0) by th
% a=.6;
% b=.1;
% x0=.2;
% y0=.4;
% th=.7;
% XP=(X-x0)*cos(th)+(Y-y0)*sin(th);
% YP=-(X-x0)*sin(th)+(Y-y0)*cos(th);
% phi=(1./(a*a))*XP.*XP+(1./(b*b))*YP.*YP-1;
% %
% %  
% %  Rectangle centered at x0, y0 
% %  with sides of length a and b
% %
 x0=-.2;
 y0=.2;
 a=.5;
 b=1;
 phi1=-((X-x0)+a/2);
 phi2=(X-x0)-a/2; 
 phix=max(phi1,phi2);
 phi3=-((Y-y0)+b/2);
 phi4=(Y-y0)-b/2; 
 phix=max(phi1,phi2);
 phiy=max(phi3,phi4);
 phi=max(phix,phiy);
figure(1)
contour(X,Y,phi,[0,0],'r');
axis([-1 1 -1 1])
axis('square')
%
%      arrays for the periodic boundary conditions
       for i=1:N
         ip(i)=i+1;
         im(i)=i-1;
       end
       im(1)=N;
       ip(N)=1;
%
%      begin simulation loop
       for iter=1:nit
           for i=1:N
             for j=1:N
              dmx=(phi(i,j)-phi(im(i),j))/h;                     % x backward difference
              dpx=(phi(ip(i),j)-phi(i,j))/h;                     % x forward difference
              dmy=(phi(i,j)-phi(i,im(j)))/h;                     % y backward difference
              dpy=(phi(i,ip(j))-phi(i,j))/h;                     % y forward difference
            
              fluxx=max(dmx,0)^2+min(dpx,0)^2;       % Godunov Flux x direction
              fluxy=max(dmy,0)^2+min(dpy,0)^2; 
              phin(i,j)=phi(i,j)-vn(x(i),y(j))*(fluxx+fluxy)^(1/2)*dt;         % advance by dt
           
           
            end
           end 
           phi=phin;                                             % update
%
%         Plotting
          contour(X,Y,phi,[0,0],'r');
          axis([-1 1 -1 1])
          axis('square')
          pause(.001)
       end
