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

%K_ENV computes cdf-version of 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 respec to the reference 
%measure are both saved and plotted.

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

%Programs Called: dist_vec.m, cdf.m, polyform.m
%
% INPUTS:
%     (i)   loc  = population location file [loc(i)=(Xi, Yi),i=1..n]
%     (ii)  sims = number of relabelings
%     (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)   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: Four bin-vectors of the following form:
%     (i)   cdf0 = CDF for population
%     (ii)  U = Upper Envelope of CDF's.
%     (iii) L = Lower Envelope of CDF's
%     (iv)  D = Array of distance values for bins.
%           [Here the max pairwise distance, Dmax, is computed
%           and distance, Diam = Dmax/2, then defines the distance 
%           range. Hence each bin, k = 1,.., bins, is associated 
%           with distance value, dist(k) = (k/(bins - 1))*d.]

% SCREEN OUTPUT: Plots cdf0 - U and cdf0 - L as blue curves 
%                against a red line base. So areas where red
%                is above BOTH blue lines indicates significant
%                clustering at level: alpha = 1/(1+sims)


%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

s = 23311 ; %random seed

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 CDF for given points in loc'

Z = polyform(poly) ;

cdf0 = cdf(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 CDF's plus P-Values

Freq = zeros(b,1);

i = 1 ;

stage = 'Begin Simulations'

while i <= sims
   
   %Form random point sample
   
   [LL,s] = rand_loc(poly,Z,Num,s) ;        
      
   %Construct CDF
   
   CDF = cdf(LL,b,Dmax) ;
        
   j = 1 ;
   
   while j <= b
      
      if CDF(j) > cdf0(j)
         
         Freq(j) = Freq(j) + 1;
         
      end
      
      j = j + 1 ;
      
   end
   
   i = i + 1 
   
end

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

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

plot(D,PVal,'b');