function C = k2ct_loc(loc,bins,index)

% This program computes k-counts for k2_local

% Written by: TONY E. SMITH, 6/18/99

% INPUTS:
%     (i)   loc    = file of locations (xi,yi), i=1..m
%     (ii)  bins   = vector of distance values (ordered low to high)
%     (iii) index  = 0 for initial population k-counts
%                    1 for simulated counts (using only first point)

% OUTPUTS: DAT = (bins x 1) vector of counts

%compiling initializations

%mbrealvector(bins);
%mbintscalar(index);

m = length(loc(:,1)) ;


if index == 0 %k-counts for all pop1 points

	%Compute pairwise distances

	D = zeros(m*(m-1)/2, 3) ; %Matrix of loc-dist pairs

	i = 1 ;

	k = 1 ;  %counter for D.

	while i < m
   
   	j = i + 1 ;
   
   	while j <= m
      
      	D(k,1) = sqrt((loc(i,1) - loc(j,1))^2 + ...
         	(loc(i,2) - loc(j,2))^2);

	      D(k,2) = i ;
      
   	   D(k,3) = j ;  
                
      	j = j + 1 ;
      
      	k = k + 1 ;
      
   	end
   
   	i = i + 1 ;
   
	end

	%Count point frequencies

	B = length(bins);

	C = zeros(m,B) ;

	[M,junk] = size(D) ;

	Dmax = max(bins); %maximum relevant distance

	mbrealscalar(Dmax);

	Dmin = min(bins); %minimum relevant distance

	mbrealscalar(Dmin);

	
	k = 1 ;

	while k <= M 
   
  
  	 if D(k,1) <= Dmax
                        
   	   i = D(k,2) ;
      
      	j = D(k,3) ; 
      
      	if D(k,1) <= Dmin
         
         	C(i,:) = C(i,:) + 1;
         
         	C(j,:) = C(j,:) + 1;
         
      	else
         
         	dist = max([D(k,1),Dmin]);
         
         	mbrealscalar(dist);
         
         	b = B ;
         
         	while bins(b) >= dist
         
         		C(i,b) = C(i,b) + 1 ;
         
         		C(j,b) = C(j,b) + 1 ; %add to both i and j 
         
         		b = b - 1 ;
         
         	end
         
      	end %end if >=
      
   	end %end if <=
   
   	k = k + 1 ;
   
	end

	   
   
else  %k-count for only first point in each simulation
   
   %Compute pairwise distances

	D = zeros(m-1) ; %Matrix of loc-dist pairs

	   
	for j = 2:m
   
      	D(j) = sqrt((loc(1,1) - loc(j,1))^2 + ...
         	(loc(1,2) - loc(j,2))^2);
	            
   end
   
   
	%Count point frequencies

	B = length(bins);

	C = zeros(1,B) ;

	Dmax = max(bins); %maximum relevant distance

	mbrealscalar(Dmax);

	Dmin = min(bins); %minimum relevant distance

	mbrealscalar(Dmin);

	%Now compute k-counts
	
	for j = 2:m    
  
  	 if D(j) <= Dmax
                        
   	        
      	if D(j) <= Dmin
         
         	C = C + 1;         
                  
      	else
         
         	dist = max([D(j),Dmin]);
         
         	mbrealscalar(dist);
         
         	b = B ;
         
         	while bins(b) >= dist
         
         		C(b) = C(b) + 1 ;         
         	        
         		b = b - 1 ;
         
         	end
         
      	end %end if >=
      
   	end %end if <=
      	   
	end %end main loop

end %end outer 'if' loop