function [PVal, D] = k_count(loc,sims,b,M,s,poly)

%K_COUNT computes the raw K-function for a point population 
%in an areal lattice (with boundary file and reference measure given).
%The P-values under the 'Randomness Hypothesis' with respect to the reference 
%measure are both saved and plotted.

%Written by: TONY E. SMITH, 5/11/99

%Programs Called: dist_vec.m, count.m, polyform.m, rand_loc.m
%
% INPUTS:
%     (i)   loc  = population location file [loc(i)=(Xi, Yi),i=1..n]
%     (ii)  sims = number of random patterns generated
%     (iii) b    = number of ticks (bins) to use on the distance axis 
%     (iv)  M    = k-vector of measure values for each of k polygons
%     (v)   s    = seed to be used in random number generator  
%     (vi)   poly = (n:2) matrix describing boundaries of k polygons in the form:
%
%          [1   n1
%           x11 y11 
%           x12 y12
%           .
%           .
%           x11 y11
%           2   n2
%           x21 x22
%           .
%           .
%           .
%           k   nk
%           xk1 yk1
%           .
%           .
%           xk1 yk1]

% DATA OUTPUTS: D = Array of distance values for bins.
%                   [Here the max pairwise distance, DDmax, is computed
%                   and distance, Dmax = DDmax/2, then defines the distance 
%                   range. Hence each bin, k = 1,.., bins, is associated 
%                   with distance value, dist(k) = (k/b)*d.]

%            PVal = P-Value under null hypothesis at each distance in D


% SCREEN OUTPUT: Plots PVal against D


%Compiling assertions

mbintscalar(sims);
mbintscalar(b);

% variable list
% s,n,k,DD,Dmax,Z,cdf0,u,M,P,F,Num,
% out,x,Freq,LL,CDF,d,inc,D0


n = length(loc(:,1)) ; % number of locations

k = length(M) ; % number of polygons


% Compute Pairwise Distances

DD = dist_vec(loc);

% Compute Dmax

DDmax = max(DD);

Dmax = DDmax/2;

stage1 = 'Construct Count for given points in loc'

Z = polyform(poly) ;

C0 = count(loc,b,Dmax);

stage2 = 'Construct Polygon Sampling Distribution'

u = ones(1,k);

m = u*M ;

P = M/m ;

% Form cdf

F = zeros(k,1);

F(1) = P(1) ;

i = 2 ;

while i <= k
   
   F(i) = P(i) + F(i-1) ;
   
   i = i + 1;
   
end

stage = 'Construct Polygon Sample Numbers'

Num = zeros(k,1) ;

i = 1 ;

while i <= n
   
   [x,s] = rand_num(s) ;
        
   j = 1 ;
   
   % Locate Polygon Interval containing x
   
   while j <= k
      
      if x <= F(j)
         
         Num(j) = Num(j) + 1 ;
         
         break
         
      end
      
      j = j + 1 ;
      
   end
   
   i = i + 1 ;
   
end


%Compute Sample Counts plus P-Values

Freq = zeros(b,1);

i = 1 ;

stage = 'Begin Simulations'

index = 10;


while i <= sims
   
   %Form random point sample
   
   [LL,s] = rand_loc(poly,Z,Num,s) ;        
      
   %Construct Count
   
   C = count(LL,b,Dmax) ;
        
   j = 1 ;
   
   while j <= b
      
      if C(j) > C0(j)
         
         Freq(j) = Freq(j) + 1;
         
      elseif C(j) == C0(j)
         
         Freq(j) = Freq(j) + 0.5;
         
      end
      
      j = j + 1 ;
      
   end
   
     
   if 100*(i/sims)>= index
      
      percent_done = index
      
      index = index + 10;
      
   end
        
   i = i + 1;
   
end

PVal = (1 + Freq) ./(sims + 1) ;

D = [1:b].*(Dmax/b);