clear
clf

N=input('2^N+1=number of grid lines in each direction--suggest N=2,3 or 4 ')

b=input('rotation coefficient, suggest +-1.5, +- 1, +-.5, +-1.5 is best ')

Tmax=input('maximum time before computation stops--suggest no more than 100/N')

a=input('bottom friction coefficient--a=0 corresponds to no friction')

xc=input(' initial x coordinate of the center of mass of the fluid--suggest +or-.25,or +or- .5')

yc=input(' initial y coordinate of the center of mass of the fluid--suggest +or-.25, or +or- .5')

mu=input(' eddy viscosity  .1 works ')

a1=input(' ratio of dt/dx; ie dt=a1*dx ')

 m=2^N;

  x=-.5:1/m:.5;
  y=ones(size(1:m+1));

  X=(y'*x)';
  Y=(x'*y)';

  dx=1/m ;
  H=1;     %%%1;
  t=0;


dt=a1*dx;

step=.1/dt;

  d=0;
  s=0;
  w=0;  %1;  %0;  %5;
  T=0;
  u=((d+w)*X+(w+T)*Y)/2;
  v=((T-w)*X+(d-s)*Y)/2;
  X=X+xc;
  Y=Y+yc;

 %%%%%%%%%%%Size of time window


NN1=10/.1; 
NN2=13/.1;

  DxX1(1:m,1:2:m-1)=X(2:m+1,1:2:m-1)-X(1:m,1:2:m-1);
  DxY1(1:m,1:2:m-1)=Y(2:m+1,1:2:m-1)-Y(1:m,1:2:m-1);
  DxX1(1:m,2:2:m)=X(2:m+1,3:2:m+1)-X(1:m,3:2:m+1);
  DxY1(1:m,2:2:m)=Y(2:m+1,3:2:m+1)-Y(1:m,3:2:m+1);

  DxX2(1:m,1:2:m-1)=X(2:m+1,2:2:m)-X(1:m,2:2:m);
  DxY2(1:m,1:2:m-1)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m);
  DxX2(1:m,2:2:m)=X(2:m+1,2:2:m)-X(1:m,2:2:m);
  DxY2(1:m,2:2:m)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m);

  DyX1(2:2:m,1:m)=X(2:2:m,2:m+1)-X(2:2:m,1:m); 
  DyY1(2:2:m,1:m)=Y(2:2:m,2:m+1)-Y(2:2:m,1:m); 
  DyX1(1:2:m-1,1:m)=X(2:2:m,2:m+1)-X(2:2:m,1:m);
  DyY1(1:2:m-1,1:m)=Y(2:2:m,2:m+1)-Y(2:2:m,1:m);

  DyX2(1:2:m-1,1:m)=X(1:2:m-1,2:m+1)-X(1:2:m-1,1:m);
  DyY2(1:2:m-1,1:m)=Y(1:2:m-1,2:m+1)-Y(1:2:m-1,1:m);
  DyX2(2:2:m,1:m)=X(3:2:m+1,2:m+1)-X(3:2:m+1,1:m);
  DyY2(2:2:m,1:m)=Y(3:2:m+1,2:m+1)-Y(3:2:m+1,1:m);


  A1=(DxX1.*DyY1-DyX1.*DxY1)/2;
  A2=(DxX2.*DyY2-DyX2.*DxY2)/2;

  m1=H*A1;
  m2=H*A2;

  h1=m1./A1;
  h2=m2./A2;

%%%%%%%%%%%%%%%%%  Mass Matrix Computation

  mass=sum(sum(sum(A1))+sum(sum(A2)))/m/m;
  M=zeros(m+1);
  M(2:m,2:m)=mass;
  M(1,2:m)=mass/2;
  M(m+1,2:m)=mass/2;
  M(2:m,1)=mass/2;
  M(2:m,m+1)=mass/2;
  M(1,1)=mass/4;
  M(m+1,m+1)=mass/4;
  M(m+1,1)=mass/4;
  M(1,m+1)=mass/4;

  MM=sum(sum(M));


  Xc=xc;
  Yc=yc;
  R=Xc^2+Yc^2;%%% sqrt(xc^2+yc^2);
 

k=0;
time=0;
Time=time;


%%% Update

while time <=Tmax         
  k=k+1;
  time=time+dt;

%% Pressure and Viscous Stress Computations

%% Substitute u for X and v for Y


  Dxu1(1:m,1:2:m-1)=u(2:m+1,1:2:m-1)-u(1:m,1:2:m-1);
  Dxv1(1:m,1:2:m-1)=v(2:m+1,1:2:m-1)-v(1:m,1:2:m-1);
  Dxu1(1:m,2:2:m)=u(2:m+1,3:2:m+1)-u(1:m,3:2:m+1);
  Dxv1(1:m,2:2:m)=v(2:m+1,3:2:m+1)-v(1:m,3:2:m+1);

  Dxu2(1:m,1:2:m-1)=u(2:m+1,2:2:m)-u(1:m,2:2:m);
  Dxv2(1:m,1:2:m-1)=v(2:m+1,2:2:m)-v(1:m,2:2:m);
  Dxu2(1:m,2:2:m)=u(2:m+1,2:2:m)-u(1:m,2:2:m);
  Dxv2(1:m,2:2:m)=v(2:m+1,2:2:m)-v(1:m,2:2:m);

  Dyu1(2:2:m,1:m)=u(2:2:m,2:m+1)-u(2:2:m,1:m);
  Dyv1(2:2:m,1:m)=v(2:2:m,2:m+1)-v(2:2:m,1:m);
  Dyu1(1:2:m-1,1:m)=u(2:2:m,2:m+1)-u(2:2:m,1:m);
  Dyv1(1:2:m-1,1:m)=v(2:2:m,2:m+1)-v(2:2:m,1:m);

  Dyu2(1:2:m-1,1:m)=u(1:2:m-1,2:m+1)-u(1:2:m-1,1:m);
  Dyv2(1:2:m-1,1:m)=v(1:2:m-1,2:m+1)-v(1:2:m-1,1:m);

  Dyu2(2:2:m,1:m)=u(3:2:m+1,2:m+1)-u(3:2:m+1,1:m);
  Dyv2(2:2:m,1:m)=v(3:2:m+1,2:m+1)-v(3:2:m+1,1:m);

%% Pressure terms in x momentum eqn.

   p1RLx=(h1.^2).*(DyY1)/4;
   p2RLx=(h2.^2).*(DyY2)/4;
   p1UDx=(h1.^2).*(DxY1)/4;
   p2UDx=(h2.^2).*(DxY2)/4;

   P1rxx=0*X;
   P1rxx(1:m,1:m)=p1RLx;

   P1Lxx=0*X;
   P1Lxx(2:m+1,1:m)=p1RLx;

   P2rxx=0*X;
   P2rxx(1:m,1:m)=p2RLx;

   P2Lxx=0*X;
   P2Lxx(2:m+1,1:m)=p2RLx;

%%  x differences in x eqn

%%% for  j=1 and i=1:m+1 take

%%  dt*(PLxx(:,1)-P1rxx(:,1))

%% for  j=3:2:m+1 and i=1:m+1 take

%%      dt*( P1Lxx(:,3:2:m+1) +P1Lxx(:,2:2,m)-P1rxx(:,3:2:m+1)-P1rxx(:,2:2:m))
  
%% for eqs indexed j=2:2:m and i=1:m+1 take

%%   dt*(P2Lxx(:,1:2:m-1)+P2Lxx(:,2:2:m)-P2rxx(:,1:2:m-1)-P2rxx(:,2:2:m))

   P2Uxy=0*Y;
   P2Uxy(1:m,1:m)=p2UDx;

   P2Dxy=0*Y;
   P2Dxy(1:m,2:m+1)=p2UDx;

   P1Uxy=0*Y;
   P1Uxy(1:m,1:m)=p1UDx;

   P1Dxy=0*Y;
   P1Dxy(1:m,2:m+1) = p1UDx;

%% y differences in x eqn

%%for i=1 and  j=1:m+1 take 

%%    - dt*(P2Dxy(1,:)-P2Uxy(1,:))

%%for i=3:2:m+1 & j=1:m+1 take 

%%       - dt*(P2Dxy(3:2:m+1,:)+P2Dxy(2:2:m,:)-P2Uxy(3:2:m+1,:)- P2Uxy(2:2:m,:))

%%  for i=2:2:m & j =1:m+1 take 

%%      -dt*(P1Dxy(1:2:m-1,:)+P1Dxy(2:2:m,:)-P1Uxy(1:2:m-1,:)-P1Uxy(2:2:m,:))

%%% Pressure Forces in x eqn 

PPx=0*X;
PPx(:,1)=PPx(:,1)+ dt*(P1Lxx(:,1)-P1rxx(:,1));
PPx(:,3:2:m+1)=PPx(:,3:2:m+1)+dt*( P1Lxx(:,3:2:m+1)-P1rxx(:,3:2:m+1));
PPx(:,3:2:m+1)=PPx(:,3:2:m+1)+dt*(P1Lxx(:,2:2:m)-P1rxx(:,2:2:m));
PPx(:,2:2:m)=PPx(:,2:2:m)+dt*(P2Lxx(:,1:2:m-1)-P2rxx(:,1:2:m-1));
PPx(:,2:2:m)=PPx(:,2:2:m)+dt*(P2Lxx(:,2:2:m)-P2rxx(:,2:2:m));
PPx(1,:)=PPx(1,:)-dt*(P2Dxy(1,:)-P2Uxy(1,:));
PPx(3:2:m+1,:)=PPx(3:2:m+1,:)-dt*(P2Dxy(3:2:m+1,:)-P2Uxy(3:2:m+1,:));
PPx(3:2:m+1,:)=PPx(3:2:m+1,:)-dt*(P2Dxy(2:2:m,:)-P2Uxy(2:2:m,:));
PPx(2:2:m,:)=PPx(2:2:m,:)-dt*(P1Dxy(1:2:m-1,:)-P1Uxy(1:2:m-1,:));
PPx(2:2:m,:)=PPx(2:2:m,:)-dt*(P1Dxy(2:2:m,:)-P1Uxy(2:2:m,:));


%% Pressure terms in y momentum eqn.

   p2udy=(h2.^2).*(DxX2)/4;
   p1udy=(h1.^2).*(DxX1)/4;
   p2RLy=(h2.^2).*(DyX2)/4;
   p1RLy=(h1.^2).*(DyX1)/4;
 
   P2uyy=0*Y;
   P2uyy(1:m,1:m)= p2udy;

   P2dyy=0*Y;
   P2dyy(1:m,2:m+1)=p2udy;

   P1uyy=0*Y;
   P1uyy(1:m,1:m)= p1udy;

   P1dyy=0*Y;
   P1dyy(1:m,2:m+1)= p1udy;

%% y differences in y eqn

%%%   for i=1 and  j=1:m+1   take

%%  dt*(P2dyy(1,:)-P2uyy(1,:))

%% for i=3:2:m+1 and  j=1:m+1   take

%%  dt*(P2dyy(2:2:m,:)+P2dyy(3:2:m+1,:)-P2uyy(2:2:m,:)-P2uyy(3:2:m+1,:))
 
%% for i=2:2:m and  j=1:m+1 take

%%  dt*(P1duyy(1:2:m-1,:)+P1duyy(2:2:m,:)-P1uyy(1:2:m-1,:)-P1uyy(2:2:m,:))


   P1ryx=0*X;
   P1ryx(1:m,1:m)= p1RLy ;

   P1Lyx=0*X;
   P1Lyx(2:m+1,1:m)=p1RLy; 

   P2ryx=0*X;
   P2ryx(1:m,1:m)=p2RLy;

   P2Lyx=0*X;
   P2Lyx(2:m+1,1:m)= p2RLy;

%% x differences in y eqn

%%% for  j=1 and i=1:m+1 take

%%  - dt*(P1Lryx(:,1)-P1ryx(:,1))

%% for  j=3:2:m+1 and i=1:m+1 take

%%     - dt*( P1Lyx(:,3:2:m+1) +P1Lyx(:,2:2,m)-P1ryx(:,3:2:m+1)-P1ryx(:,2:2:m))

%% for eqs indexed j=2:2:m and i=1:m+1 take

%%   -dt*(P2Lyx(:,1:2:m-1)+P2Lyx(:,2:2:m)-P2ryx(:,1:2:m-1)-P2ryx(:,2:2:m))

%% Pressure Forces in y eqn

PPy=0*Y;
PPy(1,:)=PPy(1,:)+dt*(P2dyy(1,:)-P2uyy(1,:));

PPy(3:2:m+1,:)=PPy(3:2:m+1,:)+dt*(P2dyy(3:2:m+1,:)-P2uyy(3:2:m+1,:));
PPy(3:2:m+1,:)=PPy(3:2:m+1,:)+dt*(P2dyy(2:2:m,:)-P2uyy(2:2:m,:));

PPy(2:2:m,:)=PPy(2:2:m,:)+dt*(P1dyy(1:2:m-1,:)-P1uyy(1:2:m-1,:));
PPy(2:2:m,:)=PPy(2:2:m,:)+dt*(P1dyy(2:2:m,:)-P1uyy(2:2:m,:));

PPy(:,1)=PPy(:,1)-dt*(P1Lyx(:,1)-P1ryx(:,1));

PPy(:,3:2:m+1)=PPy(:,3:2:m+1)-dt*(P1Lyx(:,3:2:m+1)-P1ryx(:,3:2:m+1));
PPy(:,3:2:m+1)=PPy(:,3:2:m+1)-dt*(P1Lyx(:,2:2:m)-P1ryx(:,2:2:m));

PPy(:,2:2:m)=PPy(:,2:2:m)-dt*(P2Lyx(:,2:2:m)-P2ryx(:,2:2:m));
PPy(:,2:2:m)=PPy(:,2:2:m)-dt*(P2Lyx(:,1:2:m-1)-P2ryx(:,1:2:m-1));

%% this ends the pressure computation

%% viscous stresses computation --p=2

   s1=mu*A1.*(h1.^2).*(-(DxX1).*(Dyv1)+(Dxu1).*(DyY1)+(DyX1).*(Dxv1)-(DxY1).*(Dyu1))/dx/dx;
   s2=mu*A2.*(h2.^2).*(-(DxX2).*(Dyv2)+(Dxu2).*(DyY2)+(DyX2).*(Dxv2)-(DxY2).*(Dyu2))/dx/dx;

   T1=mu*A1.*(h1.^2).*((DxX1).*(Dyu1)+(DyY1).*(Dxv1)-(DyX1).*(Dxu1)-(DxY1).*(Dyv1))/dx/dx;
   T2=mu*A2.*(h2.^2).*((DxX2).*(Dyu2)+(DyY2).*(Dxv2)-(DyX2).*(Dxu2)-(DxY2).*(Dyv2))/dx/dx;

   V1rxx=0*X;
   V1rxx(1:m,1:m)=(s1.*DyY1-T1.*DyX1)/2;

   V1Lxx=0*X;
   V1Lxx(2:m+1,1:m)=(s1.*DyY1-T1.*DyX1)/2;

   V2rxx=0*X;
   V2rxx(1:m,1:m)=(s2.*DyY2-T2.*DyX2)/2;

   V2Lxx=0*X;
   V2Lxx(2:m+1,1:m)=(s2.*DyY2-T2.*DyX2)/2 ;

%% will take following difference in x in x mom eqn. :dt*(V1rxx-V1Lxx+V2rxx-V2Lxx)

%%  x differences in x eqn

%%% for  j=1 and i=1:m+1 take

%%  (V1rxx(:,1)-V1Lxx(:,1))

%% for  j=3:2:m+1 and i=1:m+1 take

%% ( V1rxx(:,3:2:m+1) +V1rxx(:,2:2,m)-V1Lxx(:,3:2:m+1)-V1Lxx(:,2:2:m))

%% for eqs indexed j=2:2:m and i=1:m+1 take

%%  (V2rxx(:,1:2:m-1)+V2rxx(:,2:2:m)-V2Lxx(:,1:2:m-1)-V2Lxx(:,2:2:m))


   V1uxy=0*X;
   V1uxy(1:m,1:m)=(T1.*DxX1-s1.*DxY1)/2;

   V1dxy=0*X;
   V1dxy(1:m,2:m+1)=(T1.*DxX1-s1.*DxY1)/2;

   V2uxy=0*X;
   V2uxy(1:m,1:m)=(T2.*DxX2-s2.*DxY2)/2;

   V2dxy=0*X;
   V2dxy(1:m,2:m+1)=(T2.*DxX2-s2.*DxY2)/2;

%% y differences in x eqn. 

%%   for i=1 and  j=1:m+1   take

%%  (V2uxy(1,:)-V2dxy(1,:))

%% for i=3:2:m+1 and  j=1:m+1   take

%%  (V2uxy(2:2:m,:)+V2uxy(3:2:m+1,:)-V2dxy(2:2:m,:)-V2dxy(3:2:m+1,:))

%% for i=2:2:m and  j=1:m+1 take

%%  (V1uxy(1:2:m-1,:)+V1uxy(2:2:m,:)-V1dxy(1:2:m-1,:)-V1dxy(2:2:m,:))

%% Viscous Forces in x eqn

VVx=0*X;

VVx(:,1)=VVx(:,1)+(V1rxx(:,1)-V1Lxx(:,1));

VVx(:,3:2:m+1)=VVx(:,3:2:m+1)+(V1rxx(:,3:2:m+1)-V1Lxx(:,3:2:m+1));
VVx(:,3:2:m+1)=VVx(:,3:2:m+1)+(V1rxx(:,2:2:m)-V1Lxx(:,2:2:m));

VVx(:,2:2:m)=VVx(:,2:2:m)+ (V2rxx(:,2:2:m)-V2Lxx(:,2:2:m));
VVx(:,2:2:m)=VVx(:,2:2:m)+ (V2rxx(:,1:2:m-1)-V2Lxx(:,1:2:m-1));

VVx(1,:)=VVx(1,:)+(V2uxy(1,:)-V2dxy(1,:));

VVx(2:2:m,:)=VVx(2:2:m,:)+(V1uxy(2:2:m,:)-V1dxy(2:2:m,:));
VVx(2:2:m,:)=VVx(2:2:m,:)+(V1uxy(1:2:m-1,:)-V1dxy(1:2:m-1,:));

VVx(3:2:m+1,:)=VVx(3:2:m+1,:)+(V2uxy(3:2:m+1,:)-V2dxy(3:2:m+1,:));
VVx(3:2:m+1,:)=VVx(3:2:m+1,:)+(V2uxy(2:2:m,:)-V2dxy(2:2:m,:));


   V1ryx=0*X;
   V1ryx(1:m,1:m)=(T1.*DyY1+s1.*DyX1)/2;

   V1Lyx=0*X;
   V1Lyx(2:m+1,1:m)=(T1.*DyY1+s1.*DyX1)/2;

   V2ryx=0*X;
   V2ryx(1:m,1:m)=(T2.*DyY2+s2.*DyX2)/2;

   V2Lyx=0*X;
   V2Lyx(2:m+1,1:m)=(T2.*DyY2+s2.*DyX2)/2; 

%%  x diff in y  eqn 

%% %%% for  j=1 and i=1:m+1 take

%%  (V1ryx(:,1)-V1Lyx(:,1))

%% for  j=3:2:m+1 and i=1:m+1 take

%% ( V1ryx(:,3:2:m+1) +V1ryx(:,2:2,m)-V1Lyx(:,3:2:m+1)-V1Lyx(:,2:2:m))

%% for eqs indexed j=2:2:m and i=1:m+1 take

%% (V2ryx(:,1:2:m-1)+V2ryx(:,2:2:m)-V2Lyx(:,1:2:m-1)-V2Lyx(:,2:2:m))
 
   V1uyy=0*X;
   V1uyy(1:m,1:m)=(T1.*DxY1+s1.*DxX1)/2;

   V1dyy=0*X;
   V1dyy(1:m,2:m+1)=(T1.*DxY1+s1.*DxX1)/2; 

   V2uyy=0*Y;
   V2uyy(1:m,1:m)=(T2.*DxY2+s2.*DxX2)/2;

   V2dyy=0*X;
   V2dyy(1:m,2:m+1)=(T2.*DxY2+s2.*DxX2)/2;

%%  y diff in y eqn: dt*(V1dyy-V1uyy+V2dyy-V2uyy);

%%   for i=1 and  j=1:m+1   take

%% (V2uyy(1,:)-V2dyy(1,:))

%% for i=3:2:m+1 and  j=1:m+1   take

%% (V2uyy(2:2:m,:)+V2uyy(3:2:m+1,:)-V2dyy(2:2:m,:)-V2dyy(3:2:m+1,:))

%% for i=2:2:m and  j=1:m+1 take

%% (V1uyy(1:2:m-1,:)+V1uyy(2:2:m,:)-V1dyy(1:2:m-1,:)-V1dyy(2:2:m,:))

%% Viscous Forces in y eqn

VVy=0*Y;

VVy(:,1)=VVy(:,1)+(V1ryx(:,1)-V1Lyx(:,1));

VVy(:,3:2:m+1)=VVy(:,3:2:m+1)+(V1ryx(:,3:2:m+1) -V1Lyx(:,3:2:m+1));
VVy(:,3:2:m+1)=VVy(:,3:2:m+1)+(V1ryx(:,2:2:m)-V1Lyx(:,2:2:m));

VVy(:,2:2:m)=VVy(:,2:2:m)+(V2ryx(:,2:2:m)-V2Lyx(:,2:2:m));
VVy(:,2:2:m)=VVy(:,2:2:m)+(V2ryx(:,1:2:m-1)-V2Lyx(:,1:2:m-1));

VVy(1,:)=VVy(1,:)-(V2uyy(1,:)-V2dyy(1,:));

VVy(3:2:m+1,:)=VVy(3:2:m+1,:)-(V2uyy(3:2:m+1,:)-V2dyy(3:2:m+1,:));
VVy(3:2:m+1,:)=VVy(3:2:m+1,:)-(V2uyy(2:2:m,:)-V2dyy(2:2:m,:));

VVy(2:2:m,:)=VVy(2:2:m,:)-(V1uyy(2:2:m,:)-V1dyy(2:2:m,:));
VVy(2:2:m,:)=VVy(2:2:m,:)-(V1uyy(1:2:m-1,:)-V1dyy(1:2:m-1,:));


%%%  Computation of x and y Forces


Fx=PPx+VVx;
Fy=PPy+VVy;

%% this concludes the force computation

%%%%  Main update


     u1=(1-a*dt)*u+(Fx./M-dt*(X-b*v));
     v1=(1-a*dt)*v+(Fy./M-dt*(Y+b*u));


%%% Leap Frog

     X=X+dt*u1;
     Y=Y+dt*v1;


     u=u1;
     v=v1;


  DxX1(1:m,1:2:m-1)=X(2:m+1,1:2:m-1)-X(1:m,1:2:m-1);
  DxY1(1:m,1:2:m-1)=Y(2:m+1,1:2:m-1)-Y(1:m,1:2:m-1);
  DxX1(1:m,2:2:m)=X(2:m+1,3:2:m+1)-X(1:m,3:2:m+1);
  DxY1(1:m,2:2:m)=Y(2:m+1,3:2:m+1)-Y(1:m,3:2:m+1);

  DxX2(1:m,1:2:m-1)=X(2:m+1,2:2:m)-X(1:m,2:2:m);
  DxY2(1:m,1:2:m-1)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m);
  DxX2(1:m,2:2:m)=X(2:m+1,2:2:m)-X(1:m,2:2:m);
  DxY2(1:m,2:2:m)=Y(2:m+1,2:2:m)-Y(1:m,2:2:m);

  DyX1(2:2:m,1:m)=X(2:2:m,2:m+1)-X(2:2:m,1:m);
  DyY1(2:2:m,1:m)=Y(2:2:m,2:m+1)-Y(2:2:m,1:m);
  DyX1(1:2:m-1,1:m)=X(2:2:m,2:m+1)-X(2:2:m,1:m);
  DyY1(1:2:m-1,1:m)=Y(2:2:m,2:m+1)-Y(2:2:m,1:m);

  DyX2(1:2:m-1,1:m)=X(1:2:m-1,2:m+1)-X(1:2:m-1,1:m);
  DyY2(1:2:m-1,1:m)=Y(1:2:m-1,2:m+1)-Y(1:2:m-1,1:m);

  DyX2(2:2:m,1:m)=X(3:2:m+1,2:m+1)-X(3:2:m+1,1:m);
  DyY2(2:2:m,1:m)=Y(3:2:m+1,2:m+1)-Y(3:2:m+1,1:m);




  A1=(DxX1.*DyY1-DyX1.*DxY1)/2;
  A2=(DxX2.*DyY2-DyX2.*DxY2)/2;

  h1=m1./A1;
  h2=m2./A2;

 
%% sampling and graphing

 jj=k-step*floor(k/step);

 if jj==0
  SS=s1.*s1+T1.*T1+s2.^2 +T2.^2;
  SS1=[(10^13)*max(max(SS)),max(max([A1(1,1),A2(1,1)]))/dx/dx, max(max([A1 A2]))/dx/dx ,min(min([A1 A2]))/dx/dx,  N,time ]


  HH=0*X;
  AA=0*X;

AA(2:2:m,2:2:m)=(A2(2:2:m,2:2:m)+A2(2:2:m,1:2:m-1)+A2(1:2:m-1,2:2:m)+A2(1:2:m-1,1:2:m-1))/4;
AA(2:2:m,2:2:m)=AA(2:2:m,2:2:m)+(A1(2:2:m,2:2:m)+A1(2:2:m,1:2:m-1)+A1(1:2:m-1,2:2:m)+A1(1:2:m-1,1:2:m-1))/4;
AA(2:2:m,3:2:m-1)=(A1(2:2:m,3:2:m-1)+A1(1:2:m-1,3:2:m-1)+A1(1:2:m-1,2:2:m-2)+A1(2:2:m,2:2:m-2))/2;

AA(3:2:m-1,2:2:m)=(A2(3:2:m-1,2:2:m)+A2(3:2:m-1,1:2:m-1)+A2(2:2:m-2,2:2:m)+A2(2:2:m-2,1:2:m-1))/2;
AA(3:2:m-1,3:2:m-1)=(A2(3:2:m-1,3:2:m-1)+A2(2:2:m-2,3:2:m-1)+A2(3:2:m-1,2:2:m-2)+A2(2:2:m-2,2:2:m-2))/4;
AA(3:2:m-1,3:2:m-1)=AA(3:2:m-1,3:2:m-1)+(A1(3:2:m-1,3:2:m-1)+A1(1:2:m-2,3:2:m-1)+A1(3:2:m-1,2:2:m-2)+A1(2:2:m-2,2:2:m-2))/4;

%%%  Graphics version of h

  HH(2:m,2:m)=M(2:m,2:m)./AA(2:m,2:m);

Xc=sum(sum(M.*X))/MM;
Yc=sum(sum(M.*Y))/MM;
uc=sum(sum(M.*u))/MM;
vc=sum(sum(M.*v))/MM;

if max(size(xc))<NN2
   xc=[xc,Xc];
   yc=[yc,Yc];
else
   xc=[xc(2:NN2),Xc];
   yc=[yc(2:NN2),Yc];
end

if max(size(R))<NN1;
   R=[R,Xc^2+Yc^2];    %%sqrt(Xc^2+Yc^2)];
   Time=[Time,time];
else
   R=[R(2:NN1),Xc^2+Yc^2];  %sqrt(Xc^2+Yc^2)];
   Time=[Time(2:NN1),time];
end


ff=['velocity fields u-uc and v-vc as a function of x-xc and y-yc '];

ff1=['surface plots of water column height h as a function of x-xc and y-yc--i.e. centered at (0,0)'];

ff2=['surface plots of water elevation h+(x^2+y^2)/2 and the parabolic boundary (x^2+y^2)/2'];

ff3=['xc^2+yc^2 vs time and time  ' num2str(time)];



clf

%subplot(2,2,1)
%ZZ1=(((X-Xc).^2+(Y-Yc).^2)/2) ;
%ZZ2=HH+((X-Xc).^2+(Y-Yc).^2)/2;
%surf(X-Xc,Y-Yc,ZZ1)
%shading interp
%hold on
%surf(X-Xc,Y-Yc,ZZ2)
%shading interp
%hold on
%contour3(X-Xc,Y-Yc,ZZ2,20,'k')
%hold on
%contour(X-Xc,Y-Yc,ZZ2,20,'k')
%axis([-2 2 -2 2 0 2.5])
%axis square
%title(ff1)




subplot(2,2,1 )
surf(X-Xc,Y-Yc,HH)
shading interp
hold on
contour3(X-Xc,Y-Yc,HH,20,'k')
hold on
axis([-2 2 -2 2 0 2])
axis square
title(ff1)

subplot(322)
surf(X-Xc,Y-Yc,u-uc)
shading interp
hold on
surf(X-Xc,Y-Yc,v-vc)
shading interp
contour3(X-Xc,Y-Yc,u-uc,20,'k')
title(ff)
axis([-1.5 1.5 -1.5 1.5 -1.5 1.5])% -.5 .5])
axis square

subplot(3,2,4)
plot(xc,yc,'r',Xc,Yc,'+k',0,0,'+k',[-2 2],[0 0],'k',[0 0],[-2 2],'k')
axis([-1 1 -1 1])
title('trajectory of the center of mass-initial position (1/sqrt(2) , 0) ')
axis square

subplot(223)
ZZ1=((X.^2+Y.^2)/2) ;
ZZ2=HH+(X.^2+Y.^2)/2;
surf(X,Y,ZZ1)
shading interp
hold on
surf(X,Y,ZZ2)
shading interp
hold on
contour3(X,Y,ZZ2,20,'k')
hold on
contour(X,Y,ZZ2,20,'k')
axis([-2 2 -2 2 0 2.5])
axis square
title(ff2)

subplot(3,2,6)
plot(Time,R)
axis([min(Time),max(Time), min(R),max(R)])
%title('xc^2+yc^2 vs time')

title(ff3)

drawnow
%%%%%%%%%%%%%%%%%

 else
  end
end