function DAT = k12_count(L1,L2,n1,n2,len,hgt,D,sx,sy)

% This program computes k12-counts for k12_shift in (len x hgt) box

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

% INPUTS:
%     (i)   L1,L2  = locations,(xi,yi), of pops 1 and 2
%     (ii)  n1,n2  = sizes of pops 1 and 2
%     (iii) len  = length of box
%     (iv)  hgt  = height of box
%     (v)     D  = vector of distance values (ordered low to high)
%     (vi)   sx  = fixed x-shift (optional)
%     (vii)  sy  = fixed y-shift (optional)


% OUTPUTS: DAT = (D x 1) vector of k12-counts


%compiling initialization

mbrealvector(D);

n1 = length(L1(:,1)) ;

n2 = length(L2(:,1)) ;


%random shift of pop 2

if nargin == 7 %Use random shifts

	sx = (2*len)*rand(1) - len; %x-shift

	sy = (2*hgt)*rand(1) - hgt; %y-shift

end

	L2(:,1) = mod(L2(:,1) + sx,len); %x-wrap

	L2(:,2) = mod(L2(:,2) + sy,hgt); %y-wrap  


%Compute pairwise distances

D12 = zeros(n1,n2) ; %Matrix of loc-dist pairs

for i = 1:n1
   
   for j = 1:n2      
     
      D12(i,j) = sqrt((L1(i,1) - L2(j,1))^2 + ...
         (L1(i,2) - L2(j,2))^2);            
            
   end   
      
end

%Count point frequencies

D12_min = min(min(D12));

B = length(D);

C = zeros(n1,B) ;

Dmax = D(B); %maximum relevant distance

mbrealscalar(Dmax);

Dmin = D(1); %minimum relevant distance

mbrealscalar(Dmin);


for i = 1:n1
   
   for j = 1:n2
      
      dist = D12(i,j);
      
      mbrealscalar(dist);
      
      if dist <= Dmax
         
         if dist <= Dmin %j in all rings
            
            C(i,:) = C(i,:) + 1;
            
         else
            
            b = B;
            
            while D(b) >= dist
               
               C(i,b) = C(i,b) + 1;
               
               b = b -1; %b = 1 already done
               
            end
            
         end %if <= Dmin
         
      end %if <= Dmax 
      
   end % for j
   
end %for i


%Compute Output

DAT = sum(C) ; %sum columns of C to get total counts