%***********************************************************************
%  This code is written by Chiu-Yen Kao
%  Copyright is reserved  !!
%
%  compute the rhs of the following 
%
%      phi_t = -((u*phi)_x + (v*phi)_y)
%
function rhs = getRHS(phi,u,v)                           

global xmin xmax ymin ymax Mx My
global gx gy nx ny nxs nys nxe nye dx dy x y

visx = max(max(abs(u)));
visy = max(max(abs(v)));
vis = max(visx,visy);


phixp = phi(nxs+1:nxe+1,nys:nye);
phixc = phi(nxs:nxe,nys:nye);
phixn = phi(nxs-1:nxe-1,nys:nye);
phiyp = phi(nxs:nxe,nys+1:nye+1);
phiyc = phi(nxs:nxe,nys:nye);
phiyn = phi(nxs:nxe,nys-1:nye-1);
up = u(nxs+1:nxe+1,nys:nye);
uc = u(nxs:nxe,nys:nye);
un = u(nxs-1:nxe-1,nys:nye);
vp = v(nxs:nxe,nys+1:nye+1);
vc = v(nxs:nxe,nys:nye);
vn = v(nxs:nxe,nys-1:nye-1);
fxp = 0.5*(up.*phixp+uc.*phixc-vis*(phixp-phixc));
fxn = 0.5*(uc.*phixc+un.*phixn-vis*(phixc-phixn));
fyp = 0.5*(vp.*phiyp+vc.*phiyc-vis*(phiyp-phiyc));
fyn = 0.5*(vc.*phiyc+vn.*phiyn-vis*(phiyc-phiyn));
rhs(nxs:nxe,nys:nye) = -(fxp-fxn)/dx-(fyp-fyn)/dy ;