function C = point_counts(L0,L1,D)

%POINT_COUNTS counts the number of points in pattern L1 that are with
%distance d of each reference point in pattern L0, for each of a set 
% of distance values D;

%INPUTS: L0  = n0 x 2 matrix of reference locations (xi,yi)
%        L1  = n1 x 2 matrix of data locations (xi,yi)
%         D  = vector of distances {d1 < d2 < .. <dk}

%OUTPUT:  C  = n0 x k matrix of point counts C(ij) = no. of points
%              in L1 within distance dj of ref point i.

%First merge and sort both sets of points by coordinates

Dmax = max(D);

k = length(D);

n0 = size(L0,1);

n1 = size(L1,1);

C = zeros(n0,k);

for i = 1:n0
    
    dist = sqrt(sum(((L1 - ones(n1,1)*L0(i,:)).^2)')'); %col vect of dist
    
    I = find(dist <= Dmax);
    
    dist = dist(I); %reduce to feasible distances
    
    %Now compute distance vector for i and count points
    
    ni = length(I);
           
    for j = 1:ni
        
       ctr = 1;
       
       flag = 0;
           
       while flag == 0  

           if dist(j) <= D(ctr); %found smallest positive count

               C(i,ctr:k) = C(i,ctr:k) + 1; % add one to all larger dist counts

               flag = 1;  %stop

            else

                ctr = ctr + 1;

            end %End if 

        end %End while
        
    end %End for 1 = 1:ni    
   
    
end %End Main Loop