function m=dens3d(x,y,z,n,r)

% DENS3D      Density plots of 3D data.
% 
%             DENS3D(X,Y,Z,N,R) outputs various graphical representations 
%             of the data (X,Y,Z).  
%             N is the dimension of the square density matrices formed at each of 
%             R increments through the three dimensional data.  
%
%             Christian C. Beardah 1993

xmin=min(x);
xmax=max(x);
xa=linspace(xmin,xmax,n);
ymin=min(y);
ymax=max(y);
ya=linspace(ymin,ymax,n);
zmin=min(z);
zmax=max(z);
za=linspace(zmin,zmax,n);

l=length(x);

xs=(x-xmin)/(xmax-xmin);
ys=(y-ymin)/(ymax-ymin);
zs=(z-zmin)/(zmax-zmin);

figure;

for j=1:r;

m=zeros(n,n);

% Find the x and y index of each data point

  for i=1:l;
    ix=n*xs(i)-rem(n*xs(i),1)+1;
    if ix==(n+1),
       ix=n;
    end;	
    iy=n*ys(i)-rem(n*ys(i),1)+1;
    if iy==(n+1),
       iy=n;
    end;	
    iz=r*zs(i)-rem(r*zs(i),1)+1;
    if iz==(r+1),
       iz=r;
    end; 
    if iz==j,
      m(ix,iy)=m(ix,iy)+1;
    end;
  end;

% 3D plot with boxed outline of region under consideration

  subplot(2,2,1);

  zl=zmin+(j-1)*(zmax-zmin)/r;
  zh=zmin+j*(zmax-zmin)/r;

  plot3(x,y,z,'.');
% xlabel('Component 1');
% ylabel('Component 2');
% zlabel('Frequency');  
  hold;
  plot3([xmin xmax],[ymin ymin],[zl zl],':c');
  plot3([xmax xmax],[ymin ymax],[zl zl],':c');
  plot3([xmin xmax],[ymax ymax],[zl zl],':c');
  plot3([xmin xmin],[ymin ymax],[zl zl],':c');
  plot3([xmin xmax],[ymin ymin],[zh zh],':c');
  plot3([xmax xmax],[ymin ymax],[zh zh],':c');
  plot3([xmin xmax],[ymax ymax],[zh zh],':c');
  plot3([xmin xmin],[ymin ymax],[zh zh],':c');
  hold;

% 3D density plot of region under consideration

  subplot(2,2,2);

  surf(xa,ya,m');
  shading interp;
  colormap(hsv(64));
% xlabel('Component 1');
% ylabel('Component 2');
% zlabel('Frequency');
  view(-37.5,70);

% 2D interpolated shaded density plot

  subplot(2,2,3);

  pcolor(xa,ya,m');
  shading interp;
% xlabel('Component 1');
% ylabel('Component 2');
  colormap(hsv(64));

% Side on view

  subplot(2,2,4);

  zl=zmin+(j-1)*(zmax-zmin)/r;
  zh=zmin+j*(zmax-zmin)/r;

  plot(x,z,'.');
% xlabel('Component 1');
% ylabel('Component 2');
% zlabel('Frequency');  
  hold;
  plot([xmin xmax],[zl zl],':c');
  plot([xmin xmax],[zh zh],':c');
  hold;

  pause;

end;

close;