function M=PlotAll(w,r)
% PLOTALL plots a movie of the velocity and vorticity evolution
%   M=PlotAll(w,r) takes the three 3-dimensional array w of snapshots 
%   of a vorticity field and plots a color coded sequence using
%   global scaling if r=1 and local else, including the velocity vector 
%   field, if the resolution is low enough. It returns a Matlab movie M.

if nargin==1,
  r=0;                                        % default local scaling
end;
vecMax=16;
N=size(w,3);
[n,m]=size(w(:,:,1));
if n<=vecMax,
  Dx=repmat(i*[0:n/2-1 -n/2:-1]',1,m);        % derivative operators
  Dy=repmat(i*[0:m/2-1 -m/2:-1],n,1); 
  DS=Dx.^2+Dy.^2;                             % Laplacian in Fourier
  DS(1,1)=1; DD=1./DS; DD(1,1)=0;             % Laplace solver, constant=0
end;
x=0:2*pi/n:2*pi*(n-1)/n;

ma1=max(max(abs(w(:,:,1))));
if r==1,
  m1=max(max(max(w(:,:,:))));
  m2=min(min(min(w(:,:,:))));
  nocol=64;
end;
for j=1:N,
  if r==1,
    image(x,x,nocol*(w(:,:,j)-m2)/(m1-m2));         % real values
  else
    image(x,x,w(:,:,j),'CDataMapping','scaled');    % scaled values
  end;
  if n<=vecMax,
    wh=fft2(w(:,:,j));
    u=real(ifft2(-Dy.*wh.*DD));   
    v=real(ifft2(Dx.*wh.*DD));
    hold on;
    quiver(x,x,v,u);
    hold off;
  end;
  xlabel('y');ylabel('x');
  ma=max(max(abs(w(:,:,j))));
  title(['Image no ' num2str(j) '    vorticity decay = ' num2str(ma/ma1)])
  M(j)=getframe;
end;