% quiver demo *** must use Matlab 6.0 ***
 % examples of gradient, curl and divergence
 del=.1;
 x=-1:del:1; y=-1:del:1; 
 % cartesian coordinate meshgrid
 [xx,yy] = meshgrid(x,y);
 % cylindrical coordinates at meshgrid locations
 phi=atan2(yy,xx);
 rho=sqrt(+xx.^2+yy.^2);
 % ---------------------------------------------------------------------
 % EXAMPLE 1: gradient of a function zz
  figure(1)
  subplot(221)
  zz = xx.*exp(-xx.^2-yy.^2);
  [Ax,Ay] = gradient(zz,del,del);
  quiver(x,y,Ax,Ay,4);
  xlabel('x'),ylabel('y')
  axis square
  axis([min(x)-2*del,max(x)+2*del,min(y)-2*del,max(y)+2*del])
  title('quiver plot of A(x,y)')
  subplot(222)
  div=divergence(xx,yy,Ax,Ay);
  meshc(div)
  axis([0,length(x)+del,0,length(y),-5,5])
  view(20,30)
  xlabel('x'),ylabel('y')
  title('plot of div[A(x,y)]')
  Az=zeros(size(Ax));
  Cz=curl(xx,yy,Ax,Ay);
  subplot(223)
  meshc(Cz);
  view(20,30)
  xlabel('x'),ylabel('y')
  axis([0,length(x)+del,0,length(y),-1,1])
  title('quiver plot of curl[A(x,y)]_z')
   
 % ---------------------------------------------------------------------
 % EXAMPLE 2: phi-hat vector (has no divergence)
  figure(2)
  subplot(221)
  Ax=-sin(phi); Ay=cos(phi);                 % phi-hat vector
  quiver(x,y,Ax,Ay,1);
  xlabel('x'),ylabel('y')
  axis square
  axis([min(x)-2*del,max(x)+2*del,min(y)-2*del,max(y)+2*del])
  title('quiver plot of A(x,y)=phi-hat')
  subplot(222)
  div=divergence(xx,yy,Ax,Ay);
  meshc(div)
  axis([0,length(x)+del,0,length(y),-5,5])
  view(20,30)
  xlabel('x'),ylabel('y')
  title('plot of div[A(x,y)]')
  Cz=curl(xx,yy,Ax,Ay);
  subplot(223)
  meshc(Cz);
  view(20,30)
  xlabel('x'),ylabel('y')
  axis([0,length(x)+del,0,length(y),-1,1])
  title('quiver plot of curl[A(x,y)]_z')
 
 % ---------------------------------------------------------------------
 % EXAMPLE 3: phi-hat vector (has no divergence)
  figure(3)
  subplot(221)
  Ax=rho.*cos(phi); Ay=rho.*sin(phi);        % rho-hat vector
  quiver(x,y,Ax,Ay,3);
  xlabel('x'),ylabel('y')
  axis square
  axis([min(x)-2*del,max(x)+2*del,min(y)-2*del,max(y)+2*del])
  title('quiver plot of A(x,y)=rho*rho-hat')
  subplot(222)
  div=divergence(xx,yy,Ax,Ay);
  meshc(div)
  axis([0,length(x)+del,0,length(y),-5,5])
  view(20,30)
  xlabel('x'),ylabel('y')
  title('plot of div[A(x,y)]')
  Cz=curl(xx,yy,Ax,Ay);
  subplot(223)
  meshc(Cz);
  view(20,30)
  xlabel('x'),ylabel('y')
  MIN=min([min(min(Cz)),-1]);
  MAX=max([max(max(Cz)),1]);
  axis([0,length(x)+del,0,length(y),MIN,MAX])
  title('quiver plot of curl[A(x,y)]_z')
 
 
 
 
 
 
 
 % ndgrid example
 % figure(2)
 % [x1,x2,x3] = ndgrid(-2:.2:2, -2:.25:2, -2:.16:2);
 % z = x2 .* exp(-x1.^2 - x2.^2 - x3.^2);
 % slice(x2,x1,x3,z,[-1.2 .8 2],2,[-2 -.2])