function RPTS = randpoints(N)

% randpoints(N) -- Generate random points in the 3D plane.
%
% Function for generating a set of random points as a 3 x N
% matrix, where the first value is the point # (fish #1, bird #3, etc),
% the second value is the x-coord of that individual, the third
% value is the y-coord of that individual, and the fourth value is the z-coord.
% 
% Input parameter is:
% N = Number of points to be generated

L = 30 * sqrt(pi);             % Length (and width) of the field; based 
% on a circle of r=3 BL for each fish.    
% Min. distance set at 1 body length (BL).
% if exist('md') == 0
md = 1;
    % end

% First, we generate a matrix n x n wide, and fill it with numbers
% astronomically large. These will then be replaced with actual
% distances.

DIST = 9999 * ones(N,N);


% 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(N, 4);


% Now, we generate random points, and check them to make sure
% the minimum distance between pts is at least 'md' units

for i = 1:N            % loop over all members of the group
    
    mindist = 0;
    
    % keep generating randing numbers until the minimum NND
    % is > md
    
    while (mindist <= md)  
        X = L * (2*rand-1);        % generate an X value from -L to L
        Y = L * (2*rand-1);        % generate a Y value from -L to L
	Z = L * (2*rand-1);
        
        if i > 1                % Shorts the loop the first time through, since
                                % the first point cannot violate md.

            for j = 1:i-1              % loop over all the previously-generated guys
                %DIST(i,j) = dist(X,  Y, PTS(j,2), PTS(j,3));
		DIST(i,j) = dist3d(X,  Y, Z, PTS(j,2), PTS(j,3), PTS(j,4));
            end
            
            mindist = min(DIST(i,:));
        else
            mindist = 999;
        end
    end         		        
    
    % Now we've found a point that's not too close. We can
    % add it to the PTS matrix
    
    PTS(i,1) = round(i);    % assign the index number to individual i
    PTS(i,2) = X;           % assign X as the x-coord of individual i
    PTS(i,3) = Y;           % assign Y as the y-coord of individual i
    PTS(i,4) = Z;           % assign Z as the z-coord
end                       

 RPTS = PTS;
 
%  plot(RPTS(:,2),RPTS(:,3),'.')


