clear all; close all;
n=40; h=2/n; x=-1:h:1; y=x;     % Domain: [-1, 1] x [-1, 1]

%%%%% % Define a level set function for the geometry
%%
%%              x^2/0.5^2 + y^2/0.4^2 = 1,   an ellipse

for i=1:n+1,
  for j=1:n+1,
%     phi(i,j) = sqrt( x(i)*x(i)/0.5^2 + y(j)*y(j)/0.5^2)-1; 
      ph1 = sqrt( (x(i)+0.5)^2 + y(j)^2) - 0.3;
      phi2 = sqrt( (x(i)-0.5)^2 + y(j)^2) - 0.25;
      if y(j) >=0
         phi3 = y(j)-0.1;
      else
          phi3 = -( y(j) + 0.1);
      end
      phvt = [ph1,phi2,phi3];
      phi(i,j) = min(phvt);
  end                                    
end

z=zeros(n+1,n+1); contour(x,y,phi,[0 0]);   axis('square')  % Plot the ellipse
phn = z; curv=z;

for k=1:55;

for i=1:n+1,
  for j=1:n+1,
     curv(i,j) = curvature(n,n,i,j,h,x,y,phi);
%    phi(i,j) = phi(i,j) + dt *curv*phn;
  end
end

%dt = 0.25*h/max(max(abs(curv)));
 dt = 0.25*h;

for i=1:n+1,
  for j=1:n+1,
    phi(i,j) = phi(i,j) - dt *curv(i,j)*phn(i,j);
%     phi(i,j) = phi(i,j) + dt*phn(i,j);
  end
end

contour(x,y,phi,[0 0]); axis('square')
pause(1)

end