
function k_function_sim(loc,area,b,extent,sims,poly)

% K_FUNCTION_SIM computes the k-Function for a point pattern
%            plus N random point patterns for a single polygons
%            and plots the normalized L-Function plus Upper and Lower 
%            envelopes (without edge corrections)

% Written by: TONY E. SMITH, 11/26/01

% INPUTS:
%     (i)   loc    = file of locations (xi,yi), i=1..m
%     (ii)  area   = area of region
%     (iii)  b     = number of bins to use in CDF (and plot)
%     (iv) extent  = 1 if max b = half of max pairwise distance (typical case)
%                  = 2 if max h = max pairwise distance to be considered%                  
%     (v)   sims   = number of simulated random patterns
%     (vi)  poly   = polygon boundary file
%
% SCREEN OUTPUTS:  Plot of the L-Function and Simulation Envelope over the 
%                  specified extent.

% Compute Pairwise Distances

DD = dist_vec(loc);

Z = polyform(poly);

[N,junk] = size(Z); %N = number of polygons

% Compute Dmax

Dmax = max(DD);

if extent == 1
   
   Dmax = Dmax/2;
   
end


%Raw Count of Points

C = count(loc,b,Dmax);

C = C'; 


%Normalize Counts

[n,junk] = size(loc);

lam = n/area;

H = Dmax*[1:b]'/b;

K = (1/lam)*(1/n)*C;

L = sqrt(K/pi) - H;

%Now construct envelopes

C_upper = zeros(b,1);

C_lower = ones(b,1)*n^2;

[n,junk] = size(loc);

index = 10;

disp(' ');

for i = 1:sims
   
   %Generate Point Pattern
   
   loc = rand_loc(poly,Z,n);
   
   %Raw Count of Points

	C = count(loc,b,Dmax);

	C = C';

	C_lower = min(C,C_lower); 

	C_upper = max(C,C_upper); 
    

	if 100*(i/sims)>= index
      
      disp(['percent done = ',num2str(index)]);
      
      index = index + 10;
      
   end
   
end %for

%Make envelopes

K_upper = (1/lam)*(1/n)*C_upper;

L_upper = sqrt(K_upper/pi) - H;

K_lower = (1/lam)*(1/n)*C_lower;

L_lower = sqrt(K_lower/pi) - H;

%SCREEN PLOT

plot(H,L,'LineWidth',2);

hold on;

plot(H,zeros(b,1),':k');

hold on;

plot(H,L_upper,'r','LineWidth',2);

hold on;

plot(H,L_lower,'r','LineWidth',2);

xlabel('Radius','Fontsize',12, 'FontWeight', 'bold');

ylabel('L-Value','Fontsize',12, 'FontWeight', 'bold');


