function [RIPK_OUT] = ripleysK2d(M);
% RipleysK(M) ---  Will calculate Ripley's K stat, of a given array M, for a given 
% distance h.  Note that no edge correction is used here, since M should be an 
% array within the field of interest (omega).  Also, M should have x-coords in 
% column 2, y-coords in column 3.  The L transform is used on the final K.

% First, the two parameters needed are set.  N is the number of points, and R is the 
% area of the field.
N = size(M,1);
R = 3600*pi;

% 'r_field' is the radius of the field being used (taken as a circle).
r_field = sqrt(R/pi);

% Initialize matrices.
K_VALS = [];
FUNC = [];
                         

for i = 1:size(M,1)               % From each individual...
    
    DIST = [];
    for j = 1:size(M,1)           % ...to each other...
        
        % ...the distance is calculated.  
        dist = dist2d(M(i,2),M(i,3),M(j,2),M(j,3));
        
        % These distances are stored in the 'DIST' matrix.
        DIST = [DIST dist];
    end
    
    % The DIST matrix for each individual is stored in the i_Dist structure.
    i_Dist(i).all = DIST;
end

% The 'h' values are distances; at each distance h from each individual, 
% the K(h) values is found.  The h values go from 1.0 to 1/4 the radius of 
% the field (BL).  
for h = linspace(1,(r_field/4),100)
    
    I_VALS = [];                % Initialize I_VALS matrix.
    
    % Now the distances to the other individuals will be examined for each
    % individual i.  
    for i = 1:size(M,1)
        DISTS = i_Dist(i).all;
        
        % Each distance is checked and given an 'I' value.  This is simply
        % 1 if the distance to 'other' is less than h, and 0 otherwise.  
        % Distance to oneself is not counted, so 0 values are thrown out.
        for d2nbr_ind = 1:size(M,1)
            
            dist = DISTS(1,d2nbr_ind);
            
            % Check and assign 'I'.
            if dist == 0
                I = [];
            elseif dist <= h
                I = 1;
            else
                I = 0;
            end
            
            % The I values are stored in a matrix.
            I_VALS = [I_VALS I];
        end
    end
       
    
    % The K stat rests on the sum of the I func values; it is now calculated.
    K = (R/(N^2)) * sum(I_VALS);
    
    % Lastly, the L transform is done; this L value should be less than h for Poisson 
    % random arrays, and greater than h when there is clustering (theoretically).
    L = sqrt(K / pi);
    
    % K and L values are placed in matrices.
    FUNC = [FUNC L];
    K_VALS = [K_VALS; K];
    
      
end


% Output for the function contains h, K, and L values.
RIPK_OUT = [[linspace(1,(r_field/4),100)]' K_VALS [FUNC]'];


% Plot commands for testing purpose.  Not necessary to the function.
%
% plot([linspace(1,(r_field/4),100)],FUNC,'+')
% x = linspace(0,(r_field/4),1000);
% y = x;
% hold on
% plot(x,y)
% hold off