% K_FUNCTION_DEMO computes the L-function for a random pattern of points
%            with a buffer to demonstrate the theoretical shape of the
%            L-Function (with no edge corrections)

% Written by: TONY E. SMITH, 11/26/01


% Generate N points in a square of sides = R


N = 1000;  %Number of points

R = 100;  %Length of side for square

b = 20;   %Number of bins in plot

rand('seed',123453);

LOC = rand(N,2);

LOC = R*LOC;


loc = find((R/4 <= LOC(:,1))&(LOC(:,1) <= 3*R/4)&(R/4 <= LOC(:,2))&(LOC(:,2) <= 3*R/4));


Dmax = R/2;

%*************************************

%START Raw Count of Points

ctr = 1 ;

m = length(loc);

C = zeros(b,1);

while ctr <= m
   
   i = loc(ctr); %ctr-th element in loc
   
   for j = 1:N
      
      if j ~= i
         
         Dij = sqrt((LOC(i,1) - LOC(j,1))^2 + ...
            (LOC(i,2) - LOC(j,2))^2);
         
         if Dij <= Dmax
            
            bb = ceil((Dij/Dmax)*b);
            
            bin = max([1,bb]); %allows zero dists
      
            C(bin) = C(bin) + 1 ; %adds (ij) to count
            
         end %end if

      end %end if     
           
   end %end for
   
   ctr = ctr + 1 ;
     
end %end while


%Construct Cumulative Counts

k = 1 ;

while k < b
   
   C(k + 1) = C(k + 1) + C(k) ;
   
   k = k + 1 ;
   
end


%Normalize Counts

area = R^2;

lam = N/area;

n = length(loc);

H = Dmax*[1:b]'/b;

K = (1/lam)*(1/n)*C;

L = sqrt(K/pi) - H;

plot(H,L);

hold on;

plot(H,zeros(b,1),'r');