function XYZ = kernel(XY,CXY)
% This routine computes kernel home range on a single matrix
% of two dimensional locations XY, or of one matrix XY,
% and any number of others stored in a cell array CXY.  
% Result is matrix XYZ with first column grid x, second grid 
% y, and third with estimated probability of animal at that 
% location in third vector. Probabilities all add to 1.0. 
% Single entry matrix results in drawing of contour plot and
% 3D plot. Colors can be selected ahead of time.
% If more than one argument is supplied, the second one CXY MUST
% be a cell array with the ADDITIONAL xy data sets in each cell.
% A cell array is used because lengths of the different arrays
% are likely to differ. Once input, the routine then constructs
% a single grid for all, computes kernel probablities for all on
% that grid, draws contour plot for given % probability for all
% animals, and computes overlaps as % area and % of respective
% probability densities. Output XYZ is then grid axes in first
% 2 columns, and probabilities for each animal in same order as
% input in successive columns. 
hres=20;                  % number of points to check per search for h
np=17;                    % number of pixels to offset successive figures
dyp=50;                   % Relative distance between legend elements
CM1='default';
CM2='copper';
ng=nargin;
if (ng==1)                % Begin single input routine
 	[n,m]=size(XY);
	minx=min(XY(:,1));
	maxx=max(XY(:,1));
	miny=min(XY(:,2));
	maxy=max(XY(:,2));
	dx=maxx-minx;
	dy=maxy-miny;
	disp(' ');
	rsp=input('1. Enter skirt (1), world limits (2), or use default (3): ','s');
	if (rsp=='2')
		tx=num2str(minx);
		minx1=input(['   Enter minimum x value (<=' tx '): ']);    
		tx=num2str(maxx);
		maxx1=input(['   Enter maximum x value (>=' tx '): ']);  
		tx=num2str(miny);
		miny1=input(['   Enter minimum y value (<=' tx '): ']);  
		tx=num2str(maxy);
		maxy1=input(['   Enter maximum y value (>=' tx '): ']);  
	else
		if (rsp=='1')
			skt=input('   Enter skirt as % extension of world axes: ');
		else
			skt=5;         % Use default 5% skirt
		end
		skt=skt/100;
		ddx=skt*dx/2;
		ddy=skt*dy/2;
		minx1=minx-ddx;
		maxx1=maxx+ddx;
		miny1=miny-ddy;
		maxy1=maxy+ddy;
	end
	disp(' ');
	disp('2. Enter grid size as # rows and columns: ');
	NX=input('   N for x axis: ');
	NY=input('   N for y axis: ');
	dxn=(maxx1-minx1)/(NX-1);    % length of each grid cell in x direction
	dyn=(maxy1-miny1)/(NY-1);    % length of each grid cell in y direction
	dA=dxn*dyn;                  % area of each grid cell
	GX=zeros(NX*NY,1);
	GY=zeros(NX*NY,1);
	GZ=zeros(NX*NY,1);
	ii=0;
	disp('3. Getting autocorrelation parameter h: ');
	for i=0:(NX-1)
		x=minx1+(i*dxn);
		for j=0:(NY-1)
			y=miny1+(j*dyn);
			ii=ii+1;
			GX(ii)=x;
			GY(ii)=y;
		end
	end	
	TN=ii;
	h=findh(XY,hres);
	h2=h^2;
	disp('4. Computing kernel matrix... ');
	TT=n*TN;
	ii=0;
	grain=2*h2;
	hh=waitbar(0,'Creating kernel matrix...');
	dummy=0;
	for i=1:TN
		dummy=((XY(:,1)-GX(i)).^2+(XY(:,2)-GY(i)).^2)/grain;
		dummy=sum(exp(-dummy));
		GZ(i)=dummy/(2*pi*n*h2);
		waitbar(i/TN);
	end
	close(hh);
	GZ=GZ/sum(GZ);
	NCC=[];
	txc=[];
	C=[];
	XYZ=[];
	NC=input('5. Enter number of contour layers to plot: ');
	if (NC>0)
		for i=1:NC
			dnc=input('   Contour probability (%): ');
			NCC=[NCC;dnc];
		end
		NCC=sort(NCC);
		for i=1:NC
			tx=num2str(NCC(i));
			if (i<NC)
				tx=[tx ', '];
			else
				tx=['& ' tx ' %'];
			end
			txc=[txc tx];
		end
		OGZ=[GZ GX GY];
		XYZ=[GX GY GZ];
		OGZ=flipud(sortrows(OGZ,1));
		MO=OGZ(:,1);
		MO=100*cumsum(MO);
		NCD=[];
		tx=num2str(dA);
		disp(['6. Results (resolution = ' tx '):']);
		for i=1:NC
			wc=find(MO>NCC(i));  
			cc=wc(1)-1;
			NCD=[NCD;OGZ(cc,1)];
			tx1=['   ' num2str(NCC(i)) '% area: '];
			tx2=num2str(cc*dA);
			disp([tx1 tx2]);
		end; 
		disp(' ');
		disp('7. Ready for 3D Plot Plot: ');
		CM2=changemaps(CM2);
		ti=minx1:dxn:maxx1;
		tj=miny1:dyn:maxy1;
		[xi,yi]=meshgrid(ti,tj);
		zi=griddata(GX,GY,GZ,xi,yi);
		surfl(xi,yi,zi);
		hold on
		shading interp
		colormap(CM2);
		xlabel('x axis');
		ylabel('y axis');
		zlabel('Probability');
		set(gcf,'NumberTitle','off','Name','Kernel Range: 3D Plot');
		hold off
		wp=get(gcf,'Position');
		wp(1)=wp(1)+np;
		wp(2)=wp(2)-np;
		disp('8. Ready for Contour Plot: ');
		CM1=changemaps(CM1);
		figure;
		C=contourf(xi,yi,zi,NCD,'k');
		hold on
		colormap(CM1);
		%clabel(C);
		xlabel('x axis');
		ylabel('y axis');
		tx=['Kernel Contours at ' txc];
		set(gcf,'NumberTitle','off','Position',wp,'Name',tx);
		hold off
	end
	disp(' ');
else
	MXY=XY;                 % Begin multi-input overlap routines
	[n,m]=size(XY);
	NM=max(size(CXY));      % Get number of additional animals provided
	NN=NM+1;
	NNM=[n];
	for i=1:NM
		MXY=[MXY;CXY{i}];   % Make master list of points to get world
		[n,m]=size(CXY{i});
		NNM=[NNM;n];        % Create vector with matrix lengths
	end
	minx=min(MXY(:,1));
	maxx=max(MXY(:,1));
	miny=min(MXY(:,2));
	maxy=max(MXY(:,2));
	dx=maxx-minx;
	dy=maxy-miny;
	disp(' ');
	rsp=input('1. Enter skirt (1), world limits (2), or use default (3): ','s');
	if (rsp=='2')
		tx=num2str(minx);
		minx1=input(['   Enter minimum x value (<=' tx '): ']);    
		tx=num2str(maxx);
		maxx1=input(['   Enter maximum x value (>=' tx '): ']);  
		tx=num2str(miny);
		miny1=input(['   Enter minimum y value (<=' tx '): ']);  
		tx=num2str(maxy);
		maxy1=input(['   Enter maximum y value (>=' tx '): ']);  
	else
		if (rsp=='1')
			skt=input('   Enter skirt as % extension of world axes: ');
		else
			skt=5;
		end
		skt=skt/100;
		ddx=skt*dx/2;
		ddy=skt*dy/2;
		minx1=minx-ddx;
		maxx1=maxx+ddx;
		miny1=miny-ddy;
		maxy1=maxy+ddy;
	end
	disp(' ');
	disp('2. Enter grid size as # rows and columns: ');
	NX=input('   N for x axis: ');
	NY=input('   N for y axis: ');
	disp(' ');
	dxn=(maxx1-minx1)/(NX-1);    % length of each grid cell in x direction
	dyn=(maxy1-miny1)/(NY-1);    % length of each grid cell in y direction
	dA=dxn*dyn;                  % area of each grid cell
	GX=zeros(NX*NY,1);
	GY=zeros(NX*NY,1);
	GZ=zeros(NX*NY,(NM+1));
	ii=0;
	hm=[];
	for i=0:(NX-1)              % Create grid point matrix
		x=minx1+(i*dxn);
		for j=0:(NY-1)
			y=miny1+(j*dyn);
			ii=ii+1;
			GX(ii)=x;
			GY(ii)=y;
		end
	end	
	TN=ii;
	disp('3. Create kernels for each file: ');
	for k=1:NN
		tx=num2str(k);
		disp([' File ' tx ':']);
		if (k==1)
			XYY=XY;
		else
			XYY=CXY{(k-1)};
		end
		h=findh(XYY,hres);
		h2=h^2;
		TT=n*TN;
		ii=0;
		grain=2*h2;
		hh=waitbar(0,'Creating kernel matrix...');
		dummy=0;
		for i=1:TN
			dummy=((XYY(:,1)-GX(i)).^2+(XYY(:,2)-GY(i)).^2)/grain;
			dummy=sum(exp(-dummy));
			GZ(i,k)=dummy/(2*pi*n*h2);
			waitbar(i/TN);
		end
		close(hh);
		GZ(:,k)=GZ(:,k)/sum(GZ(:,k));
	end
	NCC=[];
	txc=[];
	C=[];
	XYZ=[];
	disp(' ');
	NC=input('4. Enter number of contour layers to plot: ');
	if (NC>0)
		NCC=zeros(NC,NN);
		for i=1:NC
			dnc=input('   Contour probability (%): ');
			NCC(i,1)=dnc;
		end
		NCC=sort(NCC);      % Get sorted list of contour values
		NCA=NCC;
		for i=1:NC
			tx=num2str(NCC(i,1));
			if (i<NC)
				tx=[tx ', '];
			else
				tx=['& ' tx ' %'];
			end
			txc=[txc tx];
		end
		XYZ=[GX GY GZ];
		tx=num2str(dA);
		disp(' ');
		disp(['5. Resulting Areas (resolution = ' tx '):']);
		disp(' ');
		for i=1:NC
			for j=1:NN
				MO=XYZ(:,(j+2));
				MO=flipud(sort(MO));
				VO=100*cumsum(MO);
				wc=find(VO>NCC(i,1));  
				cc=wc(1)-1;
				NCC(i,(j+1))=MO(cc);  % Save cutoff prob for this contour for each file
				NCA(i,(j+1))=cc*dA;   % Save area for this contour for each file
			end
		end;
		tbs=[];
		for i=1:(NC+1)
			tbs=[tbs;char(9)];
		end
		bl1=char('Area %',num2str(NCC(:,1)));
		tx=[tbs bl1];
		for i=1:NN
			tx1=num2str(i);
			tx1=strcat(' File',tx1);
			if (i==1)
				txlb=[tx1];
			else
				txlb=char(txlb,tx1);
			end
			bl=char(tx1,num2str(NCA(:,(i+1))));
			tx=[tx tbs bl];
		end
		disp(tx);
		disp(' ');
		disp('6. Now creating overlap contour plots: ');
		ti=minx1:dxn:maxx1;
		tj=miny1:dyn:maxy1;
		[xi,yi]=meshgrid(ti,tj);
		CM1=changemaps(CM1);
		colormap(CM1);
		map=colormap;
		[g,h]=size(map);
		ss=round(g/NN);
		for i=1:NC
			if (i>1)
				figure;
			end
			ccc=[];
			for j=1:NN
				v=NCC(i,(j+1));
				v=[v v];
				mp5=map((1+(j-1)*ss),:);
				ccc=[ccc;mp5];                      % Store line colors
				zi=griddata(GX,GY,GZ(:,j),xi,yi);
				[c,h]=contour(xi,yi,zi,v);
				set(h,'EdgeColor',mp5);
				if (j==1)
					hold on
				end
			end
			wp=get(gcf,'Position');
			wp(1)=wp(1)+((i-1)*np);
			wp(2)=wp(2)-((i-1)*np);
			txx=[num2str(NCC(i,1)) '%'];
			tx=['Kernel Contours at ' txx];
			set(gcf,'NumberTitle','off','Position',wp,'Name',tx);
			xlabel('x axis');
			ylabel('y axis');
			input('   Click mouse to place UL of legend and hit <return>..');
			p=get(gca,'CurrentPoint');
			p=p(1,1:2);
			yl=get(gca,'YLim');
			legdy=(yl(2)-yl(1))/dyp;
			for j=1:NN
				tf=text(p(1),p(2),txlb(j,:));
				set(tf,'Color',ccc(j,:));
				R=get(tf,'Extent');
				p(2)=R(2)-legdy;
			end
			hold off
		end
		colormap('hot');
		disp(' ');
		disp('7. Overlap computations:');
		more=1;
		while (more==1)
			vv=[];
			ff=[];
			tfx='   Overlap for Files ';
			v1=zeros(NN,1);
			disp('   Select files to include (1) or exclude (2): ');
			nf=0;
			for i=1:NN
				tq=['     ' txlb(i,:) ': '];         
				rsp=input(tq,'s');
				if (rsp=='1')
					vx=v1;
					nf=nf+1;
					vx(i)=1;
					vv=[vv vx];
					ff=[ff;txlb(i,:) '%']; 
					tx=num2str(i);
					if (nf>1)
						tfx=[tfx ', 'tx];
					else
						tfx=[tfx tx];
					end
				end
			end
			tfx=[tfx ':'];
			[OA,OV]=overlap(vv,GZ,NCC,NCA,dA);
			disp(' ');
			disp(tfx);
			disp('      By Overlap Area:');
			[u,v]=size(OA);
			bl1=char('Cutoff(%)',num2str(OA(:,1)));
			bl2=char('Overlap A',num2str(OA(:,4)));
			tx=[tbs tbs bl1 tbs bl2];
			for i=3:v
				bl=char(ff((i-2),:),num2str(OA(:,i),4));
				tx=[tx tbs tbs bl];
			end
			disp(tx);
			disp(' ');
			disp(tfx);
			disp('      By Overlap Volume:');
			[u,v]=size(OV);
			bl1=char('Cutoff(%)',num2str(OV(:,1),4));
			tx=[tbs tbs bl1];
			for i=2:v
				bl=char(ff(i-1,:),num2str(OV(:,i),4));
				tx=[tx tbs tbs bl];
			end
			disp(tx);
			disp(' ');
			rsp=input('  Do other combinations (1) or quit (2): ','s');
			if (rsp=='2')
				more=0;
			end
		end
	end
	disp(' ');
end