function FISH = spherefish(n, R, md)


% H = spherefish(n, R, md) - Random points in a spheric volume
%
% The SPHEREFISH() function works similar to randfish(), and has the
% same output, but it generates points inside a spherical volume of
% radius R. The number of fish generated is n, and md is the minimum
% starting distance (center to center) between the fish (this should 
% be no less than 1 BL, or fish might start out on top of each other).
% Note that the function does not actually generate numbers inside a
% spherical volume, but rather, inside of a cubic volume that circumscribes
% the sphere, and then tosses out any outside the sphere and re-generates,
% until it has n points inside the cube, and also inside the sphere,
% all of which are > md from one another in 3D.
% 

%Created by Jim Work, University of South Carolina, August 13th, 2003


% Input parameters are:
% n = Number of fish
% R = radius of the spherical volume in which the fish start (cm)
% md = minimum distance between fish centroids (cm)
% 
% A 12-column matrix is returned, with one row per fish. The columns
% contain:
% - the frame number (col. 1)
% - the individual number (col. 2)
% - position in x-y-z (col. 3-5)
% - velocity in cm/frame (col. 6-8)
% - Force in dynes (col. 9-11)
% - flag on whether it is an observed or inferred point (0 for models).
%
% Example:
% F = spherefish(4,60,2)
% 
% F =
% 
%   Columns 1 through 7 
% 
%          0    1.0000   19.3743  -49.6504  -26.6437    0.1254    0.4208
%          0    2.0000    1.1361   -6.7531  -48.1014   -0.6193    0.5305
%          0    3.0000  -37.8746   10.4168    6.8734    0.1327    0.4769
%          0    4.0000   33.1907  -47.0825  -10.3573   -0.4138    0.7151
% 
%   Columns 8 through 12 
% 
%    -0.7371    0.1254    0.4208   -0.7371         0
%     0.0786   -0.6193    0.5305    0.0786         0
%     0.6756    0.1327    0.4769    0.6756         0
%    -0.1033   -0.4138    0.7151   -0.1033         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)  
        
        % Start by generating points randomly inside of a cubic
        % volume (this is the easiest way to generate a coordinate
        % triplet).
        
        X = R * (2 * rand - 1);        % generate an X value from 0 to L 
        Y = R * (2 * rand - 1);        % generate a Y value from 0 to  W 
        Z = R * (2 * rand - 1);        % generate a Z value from 0 to H 
        
        % Now we check to make sure that the point is inside the
        % spherical volume (that is, it is <= R radius away from
        % the center point, (0,0,0). If it is not, we throw it out
        % and re-generate the points.
        
        d2c = dist3d(0,0,0,X,Y,Z);
        
        while ( d2c > R )
            % It's inside the cube but outside the sphere
            % circumscribed byt the cube. Thus we must regenerate
            % the points.
            
            X = R * (2 * rand - 1);        % generate an X value from 0 to L 
            Y = R * (2 * rand - 1);        % generate a Y value from 0 to  W 
            Z = R * (2 * rand - 1);        % generate a Z value from 0 to H 
            
            % Check again to make sure it's inside the sphere as well
            % as the cube.
            d2c = dist3d(0,0,0,X,Y,Z);
        end
        
        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;

% Compute force, assuming t = 1. Since a = v / t, if t = 1, then a = v. 
% Since F = m * a, and a = v, F = m * v. If we assume m = 1g for the
% fish, then (at least to start with), F = v. That is, the Force vectors
% basically start the fish off with the same magnitude as the velocity
% vectors and in the same direction.

FORCE = VEL;  

% 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:
FISH = [ TIME PTS VEL FORCE FLAG];