function tracking
% Program to simulate particle movement in the capture zone of a
% pumping well
% Define velocity grid
nx=15
ny=15
u=1.5  % regional groundwater velocity (constant in x direction)
Q=-20. % well pumping rate
dx=1.  % x grid spacing
dy=1.  % y grid spacing
d=1.   % flow depth
% define grid-related arrays
x=dx*[-(nx-1)/2:(nx-1)/2];
y=dy*[(ny-1)/2:-1:-(ny-1)/2];
xmin=x(1);
xmax=x(nx);
ymax=y(1);
ymin=y(ny);
% bigx and bigy define x and y corrdinates respectively, of every 
% node in grid
bigx=ones(ny,1)*x;
bigy=y'*ones(1,nx);
% compute velocities (from mass balance)
denom=2*pi*d*sqrt(bigx.^2+bigy.^2); % radial velocity denom.
zeroflag=(denom==0); % check for divide by zero
denom(zeroflag)=1.; % set zero elements of denom =1.
vr=Q./denom;  % compute vr 
vr(zeroflag)=0.; % reset vr = 0 where denom was originally zero
vx=u*ones(ny,nx)+vr.*cos(atan2(bigy,bigx));
vy=vr.*sin(atan2(bigy,bigx));
% plot velocity field
close all
figure
quiver(bigx,bigy,vx,vy)
axis([xmin, xmax, ymin, ymax])
% specify number of paricles and initial locations
% y locations distributed randomly
np=100
yminp=ymin+0.10*(ymax-ymin);
ymaxp=ymin+0.90*(ymax-ymin);
xp(1:np)=x(2);
yp(1:np)=yminp*ones(1,np)+(ymaxp-yminp)*rand(1,np);
% specify number and size of tracking time steps
nt=40
dt=0.5
% start time loop
for t=1:nt
   k=1;
   while (k<=np)
      % interpolate vecloities from grid to particle locations
      vxint=interp2(bigx,bigy,vx,xp(k),yp(k));
      vyint=interp2(bigx,bigy,vy,xp(k),yp(k));
      % compute nominal new locations
      ytemp=yp(k)+vyint*dt;
      xtemp=xp(k)+vxint*dt;
      % move particles while they remain within grid
      if(((xtemp~=0.)|(ytemp~=0.))&(xtemp>=xmin)&(xtemp<=xmax)...
         &(ytemp>=ymin)&(ytemp<=ymax))
         xp(k)=xtemp;
         yp(k)=ytemp;
      % remove particles that have left domain
      else
         xp(k)=[];
         yp(k)=[];
         np=np-1
      end
      k=k+1;
   end
   % plot particle locations at selected times
   if((t==1)|(mod(t,2)==0))
      figure
      plot(xp,yp,'o')
      axis([xmin, xmax, ymin, ymax])
   end
end
return