function [xa,ya,z,w]=pccont2(p,x,y,A,k,sym,h)

% PCCONT2     p% contour plot for bi-variate data.
%
%             PCCONT2(P,X,Y,A,K,SYM,H) splits the bi-variate data (X,Y) 
%             into groups of equal size and displays a contour (using
%             SYM as a plot style) which encloses P% of the data.  
%             Groups are formed by ranking the data by kernel density estimate.
%
%             The rows of A contain the whole data-set.
% 
%             The kernel is given by the function in K.
%             'KNORM' is the multivariate normal distribution.
%             'KEPAN' is the Epanechnikov kernal.
%
%             The optional parameter H represents the window width.
%             If omitted, an 'optimal' value is used.
%
%             Christian C. Beardah 1993

if nargin==6,
  h=0;
end;

w=zeros(length(p),2);

p=100-p;

[n d]=size(A);

xrange=max(x)-min(x);
yrange=max(y)-min(y);

xa=linspace(min(x)-xrange/4,max(x)+xrange/4,40);
ya=linspace(min(y)-yrange/4,max(y)+yrange/4,40);

[X,Y]=meshgrid(xa,ya);

z=kdest2d(X,Y,A,k,h); 
r=kdest2d(x,y,A,k,h);

[order index]=sort(r);

m=ceil(n*p/100);

for i=1:length(p),
  if m(i)==0,
    v=0.9*order(1); 
  end;

  if m(i)>0 & m(i)<n,
    v=(order(m(i))+order(m(i)+1))/2;
  end;

  if m(i)==n,
    v=(order(n-1)+order(n))/2;
  end;

  w(i,:)=[v Inf];

  contour(xa,ya,z,w,sym);  
  
  if i==1 & length(p)>1,
    hold on;
  end;
  
end;

r=order;