function [cdf0,U,L,D] = k_env_s(loc,DD,sims,M,poly)

%K_ENV_S is the same as K_ENV except that it assumes DD has been
%        sorted in ascending order (to save time)

%Written by: TONY E. SMITH, 3/30/99

%Programs Called: mySort.m, cdf.m, polyform.m
%
% INPUTS:
%     (i)    loc  = population location file [loc(i)=(Xi, Yi),i=1..n]
%     (ii)   DD = vector of pairwise distances
%     (iii)  sims = number of relabelings
%     (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);

% variable list
% s,n,k,DD,kk,Dmax,Diam,bins,b,Z,cdf0,u,M,P,F,Num,
% out,x,U,L,D,LL,CDF,d,inc,D0

s = 23311 ; %random seed

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

k = length(M) ; % number of polygons


% Compute Dmax

Dmax = max(DD);

Diam = Dmax/2;

% Construct Bin Numbers

bins = DD;

b = length(bins);

kk = 1;

while kk <= b  %reduce bins to < Diam
   
   if bins(kk) >= Diam
      
      break
      
   end
   
   kk = kk + 1 ;
   
end

b = kk;

binsR = bins(1:b); % reduced bins

stage = 'Construct CDF for given points in loc'

Z = polyform(poly) ;

cdf0 = cdf(loc,binsR,Dmax);

stage = '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 Upper and Lower Envelopes

U = zeros(1,b) ; %Upper envelope

L = ones(1,b) ;  %Lower envelopes

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,binsR,Dmax) ;
        
   j = 1 ;
   
   while j <= b
      
      U(j) = max([U(j),CDF(j)]);
   
      L(j) = min([L(j),CDF(j)]);
      
      j = j + 1 ;
      
   end
   
   i = i + 1 
   
end

D = binsR';



D0 = zeros(1,length(D));

plot(D,cdf0 - U,'b',D,cdf0 - L,'b',D,D0,'r')