function EGCP = unevenGroupCenters2d(ng,R)

% uegrpcntrs(ng) -- Generate random points in the 2D plane, suitable to be 
%  unevenly-sized group centers.
%
% Generates a set of random points as a 2 x N
% matrix, where the first value is the x-coord of that individual,
% and the second value is the y-coord of that individual.  
%  
% Input parameters are:
% ng = Number of points (groups) to be generated.
% R = A matrix containing the radii of the groups to be formed (in parent script).


% Define length of the field; based on a circle of r=3 BL for each fish.  Width 
%  is the same.
L = 30 * sqrt(pi);          
                         

% A matrix DIST is initialized; it will be filled in with the distances among 
%  points.  We start with all zeros. 
DIST = zeros(ng,ng);
% DIST(1,1) = 999;


% Next, we generate an array, 3 columns long, of individual %'s and
% their X- and Y- coordinates. Again, we do this by starting with 
% all zeros.
PTS = zeros(ng, 2);

% Now, a minimum distance matrix is built (MD).  This contains the min dist for
% each set of points, based on the radii which have been generated for them
% (again, by the parent script).

MD = zeros(ng,ng);                      % Initialize MD matrix.
for main = 1:ng                         % Loop over all points...
    for ne = 1:ng                       %  ...by all points.
        mdi = R(1,main) + R(1,ne) + 5;  % Generate md for the ith pair, 
                                        %  calculating from the radii set down 
                                        %  in R.
        MD(main,ne) = mdi;              % Build MD.
    end                                 % close loops.
end


% Now, we generate random points, and check them to make sure the minimum
%  distance between pts is at least the distance given in MD for each pair.

for i = 1:ng                            % Loop over all potential points.            

    % Now, a check for min dist violations is initialized.  The TOOCLOSE matrix
    % is initialized; coordinates are generated. Then,if non-zero values are
    % added in the check, i.e. the point was too close to one or more others,
    % the coordinates are regenerated.
    TOOCLOSE = [1];                         % Initialize matrix.
    while sum(TOOCLOSE) ~= 0                % Start loop, which checks for non-
                                            %  zero vals.
         
        X = (L-R(1,i)) * (2*rand-1);        % Gen. x-coord.  (L-r) is the max
        Y = (L-R(1,i)) * (2*rand-1);        % x and y value, so that the groups 
                                            % will be far enough from the edge  
                                            % of the field.
        
        if i > 1            % This shorts the subsequent loops the first time
                            % through, since the first point cannot be too close
                            % to the others; there are no others.
            
            % New point is checked against previous ones, to look for md 
            % violations.
            for j = 1:i-1                   % Loop over all previous.
                
                DIST(i,j) = dist(X,  Y, PTS(j,1), PTS(j,2)); % Check distances.
                
                if (DIST(i,j) <= MD(i,j))   % If the dist. is too close, a value                                           
                    TOOCLOSE(1,j) = 1;      % of one is returned to TOOCLOSE.
                else
                    TOOCLOSE(1,j) = 0;      % If the dist is ok, return 0.
                end
            end
            
        else
            TOOCLOSE = 0;           % The first time through, return 0, since 
                                    % there cannot be a md violation (see above).
        end
    end         		        
    
    
    % Now we've found a point that's not too close. We can
    % add it to the PTS matrix
 
    PTS(i,1) = X;           % assign X as the x-coord of individual i
    PTS(i,2) = Y;           % assign Y as the y-coord of individual i
        
end                       

EGCP = PTS;

% DIST
% 
% plot(EGCP(:,1),EGCP(:,2),'.')
% axis([-60,60,-60,60])