%***********************************************************************
%  This code is written by Chiu-Yen Kao
%  Copyright is reserved  !!

function Tstep=getTimeStep(u,v,phi,tol,eps)
global xmin xmax ymin ymax Mx My
global gx gy nx ny nxs nys nxe nye dx dy x y


%
%  When advance
%
%     S.t = u*S.x + v*S.y 
%
%  the chosing of time steps is governed by the space meshsize;
%  When using central differencing in space and RK in time, we
%  can show using spectrum method the following CFL condition:
% 
%       dt <= Const /( Umax/dx + Vmax/dy )            
%
%  For RK3, Const can be 1.5, for RK4, Const = 2.7 .
%

cfl = 0.5;
umax = 0.0;  
vmax = 0.0;
smax = 0.0;  % maximum velocity

for i = nxs : nxe
   for j = nys : nye
      vel = sqrt( u(i,j)*u(i,j) + v(i,j)*v(i,j) );
      if abs( u(i,j) ) > umax
         umax = abs(u(i,j));
      end
      if abs( v(i,j) ) > vmax
         vmax = abs(v(i,j));
      end
      if abs( vel ) > smax
         smax = vel;
      end
   end
end

%
% make sure each time step moves no more than half grid
%

gLimit = .5 * sqrt( dx*dx + dy*dy ) / ( vel + tol );

%
% stability requarement
%

sLimit = cfl / ( umax/dx + vmax/dy);

if gLimit < sLimit
   Tstep = gLimit;
else
   Tstep = sLimit;
end