function h = findh(XY,NS)
% This routine finds the optimal h for the kernel home
% range method using the cross-validation method.
% XY is a two column matrix of x,y locations and NS is
% the number of steps to examine between limits to look
% for h.
[n,m]=size(XY);
vx=var(XY(:,1));
vy=var(XY(:,2));
S=sqrt((vx+vy)/2);
dS=S*exp(-log(n)/6);
d1=0.1*dS;
d2=10*dS;
more=1;
TN=n*(n-1)/2;
DG=zeros(TN,2);
DH=zeros(TN,2);
DD=zeros(TN,1);
ii=0;
MD=[];
MV=[];
hb=waitbar(0,'Finding interpoint distances...');
for i=1:(n-1)
	dm=i*ones(n-i,1);
	MD=[MD;dm];
	dn=[(i+1):n]';
	MV=[MV;dn];
	waitbar(i/(n-1));
end
DG=[XY(MD,1) XY(MV,1)];
DH=[XY(MD,2) XY(MV,2)];	
DD=(diff(DG,1,2).^2)+(diff(DH,1,2).^2);
close(hb);
TT=TN*(NS+1);
lasth=0;
while (more==1)
	dh=(d2-d1)/NS;
	hV=zeros(NS,1);
	CVV=zeros(NS,1);
	ii=0;
	hb=waitbar(0,'Plotting h graph...');
	for k=0:NS
		hh=d1+(dh*k);
		hV(k+1)=hh;
		h2=hh^2;
		coef1=1/(pi*h2*n);
		coef2=1/(4*pi*h2*n*n);
		w=-1/(4*h2);
		LG=DD.*w;
		UG=(DD.*2*w);
		CV=sum(exp(LG)-4*exp(UG));
		CVV(k+1)=coef1+(coef2*CV);
		waitbar(k/NS);
	end;
	close(hb);
	[vh,hi]=min(CVV);
	th=hV(hi);
	if (hi>1) 
		d1=hV(hi-1);
	else
		d1=hi;
	end
	if (hi<(NS+1))
		d2=hV(hi+1);
	else
		d2=NS+1;
	end
	colormap('default');
	plot(hV,CVV,'b');
	xlabel('h value');
	ylabel('test value');
	set(gcf,'NumberTitle','off','Name','H Plot');
	tx=num2str(th);
	dh=abs(th-lasth);
	lasth=th;
	tx1=num2str(dh);
	disp(['   Current h = ' tx '   Change = ' tx1]);
	rsp=input('     Refine further (1), accept (2), use custom h (3), change range (4): ','s');
	switch rsp
	case '2',
		more=0;
	case '3',
		th=input('       Input custom value for h: ');
		more=0;
	case '4',
		d1=input('       Enter lower bound for h: ');
		d2=input('       Enter upper bound for h: ');
	end
	figure(1)
end
delete(gcf);
h=th;