function [P,C,C1] = k2_local_exact(loc,n1,D);

% K2_LOCAL_EXACT compares two point populations, and computes P-values for 
% clustering at EACH POINT of pop1 relative to pop1 + pop2. This is an 
% analytical solution that counts the total points, ni, in each distance ring 
% about point i in pop1 and also counts the number mi of these point from pop1. 
% Then the hypergeomtric probability of getting at least this many pop1 points
% yields the desired p-value for i.
% 
% Written by: TONY E. SMITH, 12/11/07
% 
% INPUTS:
%     (i)   loc  = combined location file [loc(i)=(Xi, Yi),i=1..n]
%                  where n = n1 + n2. It is assumed that points for
%                  pop1 are in the first n1 rows.]
%     (ii)  n1   = size of pop1.
%     (iii) D = vector of distance values at which to measure clustering.

% OUTPUTS: P  = Matrix of P-values (Pij) for each point i at distance j.
%          C  = Total pop counts at each radius
%          C1 = Pop 1 counts at each radius

n = size(loc,1) ;

B = length(D);

X = loc(:,1);

Y = loc(:,2);

C = zeros(n1,B); % Total Counts matrix

C1 = zeros(n1,B); % Pop1 Counts matrix

for i = 1:n1
    
    %Calculate distances    
       
    x = loc(i,1);
    
    y = loc(i,2);
    
    DIST = sqrt(((x - X).^2) + ((y - Y).^2)); 
    
    %Count points
    
    Dmax = max(D); %maximum relevant distance

    Dmin = min(D); %minimum relevant distance
    
    for b = 1:B
        
        d = D(b);
        
        I = find(DIST <= d);
        
        ni = length(I) - 1; % Tot other pop within d of i
        
        I1 = find(I <= n1); %sublist of other pop1 members
        
        mi = length(I1) - 1 ; % Tot other pop1's within d of i
       
        if mi == 0
            
            P(i,b) = 1;
            
        else
            
            P(i,b) = 1 - hygecdf(mi-1,n-1,n1-1,ni);
            
            % chance of getting at least mi of the other n1-1 pop1 members  
            % in a sample of size ni from the other n-1 pop members.
            
        end
               
        C(i,b) = ni;
        
        C1(i,b) = mi;
        
    end
    
end