function [f,h]=kdest3da(x1,x2,x3,A,k,h)

% KDEST3DA    Tri-variate Kernel Density Estimator with generalised smoothing.
%
%             Examples: load kddata
%                       f=kdest3da(0,0,0,lukesh1(:,1:3),'knorm',0) 
% 
%             KDEST3DA(X1,X2,X3,A,K,H) evaluates a tri-variate Kernel Density 
%             Estimate with kernel given by the function in K.
%             'KNORM' is the multivariate normal distribution.
%
%             (X1,X2,X3) contains the point(s) for which the density 
%             estimate is required.  X1, X2 and X3 can be scalars, vectors 
%             or matrices.
%
%             The tri-variate data points are stored row by row in A.
%
%             The parameter H represents the window width:
%             (1) If a single non-zero value of H is given then smoothing is the 
%             same in all directions.
%             (2) If H=[H1, H2, H3] then smoothing is by HI in the xi direction.  
%             If H=0 is specified, automatic 'optimal' values are used. 
%
%             Christian C. Beardah 1995

if nargin==5 | h==0,
  h=[0 0 0];
end;

if length(h)==1,
  h=[h h h];
end;

[n d]=size(A); % Establish dimensions of A, 
               % n=number of data points,
               % d=dimension of the data. 

[r c]=size(x1);
o=ones(r*c,1);

if nargin<5 | nargin>6,
  error('Incorrect number of input arguments');
end;

if min(h)==0,  % If no value of h has been specified by the user then 
               % calculate optimal values.
  if h(1)==0,
    h(1)=hns(A(:,1),k);
  end;
  if h(2)==0,
    h(2)=hns(A(:,2),k);
  end;
  if h(3)==0,
    h(3)=hns(A(:,3),k);
  end;
end;

X=repmat(x1,[1 1 n]);
Y=repmat(x2,[1 1 n]);
Z=repmat(x3,[1 1 n]);

XA=reshape(o*A(:,1)',[r c n]);
YA=reshape(o*A(:,2)',[r c n]);
ZA=reshape(o*A(:,3)',[r c n]);

X=(X-XA)/h(1);
Y=(Y-YA)/h(2);
Z=(Z-ZA)/h(3);

store=feval(k,X,Y,Z);

f=sum(store,3)/(n*prod(h)); % f holds the kernel density estimates for the points in (x,y)