function m=density(x,y,v,s,f);

% DENSITY     Smooth density plot.
%
%             DENSITY(X,Y,V,S,F) creates a density plot of the data in (X,Y).
%             A square density matrix is formed by moving a square of size S 
%             around the data.
%             V specifies the number of overlaps of each square.
%             If the optional argument F is specified then a movie
%             of F frames is produced. 
%
%             Christian C. Beardah 1993

xmin=min(x);
ymin=min(y);
xmax=max(x);
ymax=max(y);
l=length(x);
xs=(x-xmin)/(xmax-xmin);
ys=(y-ymin)/(ymax-ymin);

inc=s/v;
t=1/inc;

% Size of the final density matrix.

n=(1-s)/inc+1;
m=zeros(n,n);
T=zeros(t,t);

xa=linspace(xmin,xmax,n);
ya=linspace(ymin,ymax,n);

for i=1:l;
  ix=xs(i)*t-rem(xs(i)*t,1)+1;
  if ix==(t+1),
    ix=t;
  end;	
  iy=ys(i)*t-rem(ys(i)*t,1)+1;
  if iy==(t+1),
    iy=t;
  end;	
  T(ix,iy)=T(ix,iy)+1;
end;			

for i=1:n,
  for j=1:n,
    for k1=i:i+v-1,
      for k2=j:j+v-1,
        m(i,j)=m(i,j)+T(k1,k2);
      end;
    end;
  end;
end; 


figure;

subplot(2,2,1);
surf(xa,ya,m');
axis('square');
colormap(brighten(hot,0.75));
shading interp;

subplot(2,2,2);
pcolor(xa,ya,m');
axis('square');
colormap(brighten(hot,0.75));
shading interp;
hold;
contour(xa,ya,m','k');

subplot(2,2,3);
pcolor(xa,ya,m');
axis('square');

subplot(2,2,4);
plot(x,y,'ro');

pause;

if nargin==5,
  figure;
  M=moviein(f+1);
  for j=1:f+1,

%   Rotation 'up and over'.
%   A=viewmtx(90,(j-1)*90/f,10);

%   Roataion 'around and about'.  
    A=viewmtx((j-1)*360/f,45,10);  

    surfl(xa,ya,m',[0 45]);
    shading interp;
    view(A);
    axis('off');
    colormap(brighten(hot,0.5)); 
    M(:,j)=getframe;
  end;
  m=M;
end;

close;