function RPTS = randFish3d(n, L, W, H, md )

% randfish(n, L, W, H, md) -- Build a formatted file of randomly position fish in 3D.
%
% Function for generating a set of random points as a 12 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, and the third
% value is the y-coord of that individual.
% 
% Input parameters are:
% n = Number of points to be generated
% L = length of the field
% W = width of field 
% H = height of the field
% md = minimum distance allowable between individuals ("personal space")
%
%
% Output files need to contain, on a single line:
% - time
% - individual number
% - position (x1,x2,x3)
% - velocity (v1,v2,v3)
% - angular orientation vector (w1,w2,w3)
% - a flag whether it is an observed position or an inferred position (in
%   the expts)
%
% A single line, therefore, looks like this:
% t N x1 x2 x3 v1 v2 v3 w1 w2 w3 flag

% This will be accomplished by generating five matrices, and concatenating them.
%  TIME: will contain the time (0, in this case, since it's the starting position)
%   PTS: will contain N, x1, x2, and x3 for each fish
%   VEL: will contain Vx, Vy, and Vz for each fish
% FORCE: will contain Fx, Fy, Fz for each fish
%  FLAG: will contain the flag (0 for observed, 1 for inferred). 
% 		 In the model, it's always 0

% First, generate the time = 0 column vector:
TIME = zeros(n,1);

% Next, 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, 4 columns long, of individual %'s and
% their X- and Y- coordinates. Again, we do this by starting with 
% all zeros.

% Set up the random number seed so that we have true, rather than pseudo, random
% numbers (note that if you don't do this, the random numbers start off the same
% every time you restart MatLab). We use clock as the variable, since the system 
% clock will be different every time this function is called.

rand('state',sum(100*clock));		


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 you the minimum NND
	% is > md.
	
	while (mindist <= md)  
		X = L * rand;        % generate an X value from 0 to L 
		Y = W * rand;        % generate a Y value from 0 to  W 
		Z = H * rand;        % generate a Z value from 0 to H 
		for j = 1:i-1        % loop over all the previously-generated guys
			DIST(i,j) = dist3d(X,  Y,  Z, PTS(j,2),PTS(j,3), PTS(j,4));
		end
	
		mindist = min(DIST(i,:));
		
	end					    % end k loop    

	% 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 of individual i
	
end	                    % terminate i loop


% Now we build the matrix for velocity. At the start, mean velocity is between 0 and
% 1. Since each time step is equivalent to one 'frame' of video, and we capture
% and film at 30 Hz (frames/sec), a velocity of 1 cm/step = 1 cm/frame = 30 cm/s, so
% that it could cross the tank in 3.3 seconds. Any faster than that would probably
% cause the model fish to swim THROUGH the walls if things were set up wrong. Since
% fish can swim in either direction, this number should be between -1 and 1, actually.

VEL = (2 * rand(n,3)) - 1;

% Penultimately, we need to build the angles in X, Y, and Z. These are represented
% as an n x 3 vector, with the first column being the X direction of orientation,
% the second column being the Y direction for orientation, and the third column being Z.
% Thus, for example, a 45 degree angle in the X-Y plane would be represented as something
% like [1,1], since, if you start at (0,0) and draw a line to (1,1), it's at a 45 degree
% angle to the horizontal. This can be done for all the other planes as well. So, the
% three numbers are ranged from 0 to 1, and represent the 'strength' of the orientation along
% the positive and negative axes of each direction. Thus, these numbers can range from -1 to +1.

ANGLE = VEL;        % assume fish start out facing in the same direction as their velocity

% And finally, the FLAG matrix, which is all 0's because it's known data
% and not interpolated.
FLAG = zeros(n,1);

% concatenate all the pertinent matrices into one big giant one for output:
RPTS = [ TIME PTS VEL ANGLE FLAG];