function tau_diggle(L1,L2,N,poly,s)


%TAU_DIGGLE compares two point populations L1,L2 by generating N random reference 
%points, measuring nearest neighbor distances to each population from these points,
%and correlating these distance pairs using Kendall's tau. The distribution of this 
%statistic is then used to test for significant attraction and/or repulsion between 
%these spatial populations.

%Reference: Diggle, P.E. and T.F. Cox (1981) "On sparse sampling methods and tests of
%           independence for multivariate spatial point patterns", BULLETIN OF THE
%           INTERNATIONAL STATISTICAL INSTITUTE, 49: 213-228.


%Written by: TONY E. SMITH, 12/2/01

% INPUTS:
%     (i)   L1     = locations of pop1
%     (ii)  L2     = locations of pop2
%     (iii) sims   = number of random reference points
%     (iv)  poly   = ([x11,x12],...,[xk-1,1,xk-1,2],[x11,x12])
%                  = (k:2)-matrix of points defining polygon
%                    (Note: last point must be same as first!)
%     (v)      s   = optional specification of random seed. Default seed = 23311.
%                    s = 1 if random seed is to be used.
%                    s = 5_digit positive integer to be used as seed 


% SCREEN OUTPUT:  Value of Kendall's tau and significance of the test.


% Compute Initial Seed

if nargin == 4
   
   s = 23311;  %forest seed (default)   
    
   %s = 54321; %myrtles seed (sims = 999)
   
elseif s == 1
   
   s = round(100000*rand(1))
   
else
   
   s = s;
   
end



%Initial Calculations

n1 = size(L1,1);

n2 = size(L2,1);

% Generate Points

% Define containing rectangle

Mn = min(poly) ; 
Mx = max(poly) ;
   
x0 = Mn(1) ;
y0 = Mn(2) ;
x1 = Mx(1) ;
y1 = Mx(2) ;


% Generate points
   
Pts = zeros(N,2);

for i = 1:N
   
   kk = 0 ;
   
   counter = 1;
   
   while kk == 0 ;
      
      r = rand(1);
      
      x = x0 + r*(x1 - x0) ;
      
      r = rand(1);            
      
      y = y0 + r*(y1 - y0) ;
      
      if pt_in_poly([x,y],poly) == 1
         
         Pts(i,:) = [x,y] ;
         
         kk = 1;
         
      end %end if
      
      % Check to avoid infinite cycles
      
      counter = counter + 1;
      
      if (kk == 0)&(counter == 100)
         
         error('looping error'); 
         
      end	         
      
   end %end point generation
   
end %for
   

%Plot Points

xp = poly(:,1);

yp = poly(:,2);

%plot(xp,yp); hold on;

x = Pts(:,1);

y = Pts(:,2);

%plot(x,y,'.r');

%Compute Distances and Counts

D1 = zeros(N,1);

D2 = zeros(N,1);



for i = 1:N
   
   Diff1 = [L1 - ones(n1,1)*Pts(i,:)]';
   
   DD1 = sqrt(sum(Diff1.^2))';
   
   D1(i) = min(DD1);
   
   Diff2 = [L2 - ones(n2,1)*Pts(i,:)]';
   
   DD2 = sqrt(sum(Diff2.^2))'; 
   
   D2(i) = min(DD2);    
      
end

%Compute observed tau value

[tau,z] = kendall_tau(D1,D2);


%SCREEN OUTPUT

disp(' ');
disp('TEST RESULTS')
disp(' ');
disp(' ');
disp(['Kendall Tau = ',num2str(tau)]);
disp(' ');
disp(' ');

if z > 0
   
   pval = normcdf(-z);
   
   disp(['Attraction P-Value = ',num2str(pval)]);
   
elseif z < 0
   
   pval = normcdf(z);
   
   disp(['Repulsion P-Value = ',num2str(pval)]);
   
else
   
   disp(['No Significance: P-Value = ',num2str(.5)]);
      
end