function OUT = Radius_4000(Incin,Lung,Larynx);

%RADIUS_4000 identifies the set of lung and larynx cases
%within a distance of 4000 from the incinerator

%This script assumes that Larynx.mat is the given workspace.
%Here: 

%	Incin  = location coords of incinerator
%	 Lung  = location coords of Lung cases
%	Larynx = location coords of Larynx cases

%OUTPUTS: OUT is an output structure with: 

%     OUT.L1 = coords of Larnyx cases within distance 4000
%     OUT.D1 = Distances of L1 to the Incinerator
%     OUT.L2 = coords of Lung cases within distance 4000
%     OUT.D2 = Distances of L2 to the Incinerator  
%    OUT.LOC = [L1;L2] = stacked file of coords
%   OUT.DIST = matrix DIST for JMP with 
%              DIST(:,1) = [ones(N1,1);2*ones(N2,1)];
%              DIST(:,2) = [D1;D2];

[n1,junk] = size(Larynx);
[n2,junk] = size(Lung);


% Construct L1
  
DD1 = zeros(n1,1);

for i = 1:n1 
   
   DD1(i) = sqrt((Incin(1) - Larynx(i,1))^2 + ...
         (Incin(2) - Larynx(i,2))^2);            
            
end   

I = find(DD1 <= 4000);

L1 = Larynx(I,:); %Larynx cases <= 4000

D1 = DD1(I); %Distances to the Incinerator


% Construct L2

DD2 = zeros(n2,1);

for i = 1:n2 
   
   DD2(i) = sqrt((Incin(1) - Lung(i,1))^2 + ...
         (Incin(2) - Lung(i,2))^2);            
            
end   

I = find(DD2 <= 4000);

L2 = Lung(I,:); %Lung cases <= 4000

D2 = DD2(I);  %Distances to the Incinerator

LOC = [L1;L2];

%Output for JMPIN

N1 = length(D1);
N2 = length(D2);

DIST = zeros(N1+N2,2);

DIST(:,1) = [ones(N1,1);2*ones(N2,1)];
DIST(:,2) = [D1;D2];

OUT.L1 = L1;
OUT.D1 = D1;
OUT.L2 = L2;
OUT.D2 = D2;  
OUT.LOC = [L1;L2];
OUT.DIST = DIST;