function [MVMT] = randMovement3d(VECTOR, n, rule, l, w, h, mass, hz);

% movementrrnd(<array>,N, <rule>, l, w, h, mass, hz) -- Compute fish movement based on the rule.
%
% The movementrnd() function is designed to take an <array> of 
% data in the typical 'initpts' type format, including x-y-z coordinates,
% and use one of the valid (i.e. already programmed) movement rules
% to move the individual listed. Movement is restricted within the 
% 3-D bounding box of l x w x h.  We also need to know the
% mass of the fish (in g) and the speed of the film (in hz).
% The <rule> is a rule number, usually asked for by a script. 

% ----------------------------------------------------------------------------
% INITIALIZATION OF VARIABLES
% ----------------------------------------------------------------------------

% Initialize the range of the uniform number generator. This
% is equivalent to setting the max force on the fish (in "virtual dynes").
%
% How to compute this, using the x-components of the Force, Accel, and velocity
% vectors as examples:
% Fx = m * ax
% m = 4 grams
% Thus, ax = F / 4.
% ax also = vx / t.
% Let t = 1 second.  Since both F/4 and vx * t = ax, then we can say:
% F / 4 = vx / 1, or F / 4 = vx.
% Now, we assume max velocity of the fish is (for ease) roughly 60 cm / s. So we
% want to know what is the maximum force, but we know the maximum velocity. We can compute
% it by substitution:
% F / 4 = 60 ==> F = 60 * 4 = 240 dynes
% So, the maximum force that can be exerted on the fish at any given time, in any
% given frame, is 240 dynes.

global Fmax fsd SightMax PARMS drag_const;

if ( isempty(PARMS) == 1 )
	PARMS = [5.51, 7.56, 12.75, 7.85, -102, -44.1, -98.2, 13.1, 8.1];
end

% Next, we initialize all directional vectors. There are a total of
% 6 force vectors that can be used in the model.
% We initialize them to be 0. This way, if they are not included in the model,
% adding them has no effect. Such a system allows us to include whichever
% pieces of the model we wish.

% The forces follow (vaguely) Matuda and Sannomiya (1980) for 2D.


%
% SOCIAL FORCES
% 

% Let s be the social force, or the tendancy to move toward far companions but
% away from those that are too near (aka. attraction/repulsion force). We will call
% sx the social force in x, sy the social force in y, and sz the social force in z.
% To imagine how this works, pretend two fish are side by side, but are too close to
% each other for comfort in the x plane. Then sx would be negative for the fish
% in the low-x side, and positive for the fish on the high-x side, and they would 
% separate.

% Initialization of the social vectors:
sx = 0;
sy = 0;
sz = 0;

%
% ORIENTATION FORCE
%

% Let o be the orientation force, or the tendancy to match swimming speeds and
% orientation angles with those of your neighbor(s). Thus if you are moving at
% 2 cm/s in the x plane and your neighbor is moving 10 cm/s, you will want to
% speed up (within the limits of your ability to accelerate, of course). We will
% call ox the orientation vector in x, oy the orientation vector in y, and oz
% the orientation vector in z.

% Initialize the orientation vectors:

ox = 0;
oy = 0;
oz = 0;

%
% WALL FORCE
%

% Fish tend not to actually collide with the wall, so we introduce w, the 
% force repelling the fish away from the wall, or the intensity with which
% fish prefer not to swim toward the wall. This "force" gets stronger as they
% get closer to the wall. Let wx be the wall repulsion force in x, and so on.

% Initialize the wall vectors:

wx = 0;
wy = 0;
wz = 0;

%
% DIRECTIONAL FORCE
%

% A "directional" force is a constant force that acts in a particular direction,
% such as fish wanting to "move up a gradient" or away from a predator. If we let 
% d represent the directional force, then dx would be the directional force in x,
% dy the directional force in y, etc.

% Initialize the directional force:
dx = 0;
dy = 0;
dz = 0;

%
% RANDOM FORCE
%

% This is the tendancy of the fish to move randomly. There might be no tendancy
% to do so (as in Selfish Herd theory), or some slight tendancy. If the Wall
% and Random forces are the only things acting, this should pretty much be
% the equivalent of bounded random movement. We will call r the random force, 
% so that rx is the random movement in the x direction, and so on.

% Initialize the random force:

rx = 0;
ry = 0;
rz = 0;

% DRAG FORCE
% The force of drag (opposite to the current velocity vector).
% ----------------------------------------------------------------------------
% IMPLEMENTATION OF THE MOVEMENT RULE
% ----------------------------------------------------------------------------


% switch rule									    % check the rule # requested
% 
% case 1,											% RANDOM MOVEMENT is rule #1
% 	[rx, ry, rz] = randmove(fsd); 		  		% compute random x-y-x components
% 	% Zero out the velocity from last time (there's no memory).
% 	VECTOR(n,6:8) = 0;
% 	
% case 2,												% RANDOM MOVEMENT with MEMORY is rule #2
% 	[rx, ry, rz] = randmove(fsd); 		  			% call the randmove function
% 	
% case 3,												% the TENDENCY DISTANCE model
	global frontdeg backdeg drag_const AR_N;
	
	[rx, ry, rz] = randmove(fsd); 		  							% call the randmove function
	
% otherwise,										% There is no such rule # coded
% 	error('Invalid movement rule number!');
% end

% Wall repulsion ALWAYS occurs. Note that we must check the velocity without
% the wall force present, and make sure the wall force away from the wall (if 
% you're too close) is greater than the vector toward the wall. (This is a
% potential problem if you have a bunch of forces that all happen to line up 
% -- likely to be very rare but still possible.). Ultimately what we will want to
% do is scale the wall force to the non-wall force. So, if you have a force of 5
% and you're close to the wall, wall-repulsion would be O(5) (where O(5) is 
% read as "on the order of 5", or "approximately 5").



% All we need to do is find the largest of the three components and use that.
% Forces without wall force: (note, Kinert was already applied in the intertia function).

pre_vx = sx + ox + dx + rx;		
pre_vy = sy + oy + dy + ry;		
pre_vz = sz + oz + dz + rz;        

PRE_V = [pre_vx, pre_vy, pre_vz, Fmax];	% Fmax is the smallest "maximum force" you can have

pmv = max(PRE_V);						% Fmax is the "order of magnitude" we're
% going to have to deal with.

[wx, wy, wz] = wallforce(VECTOR(n,3:5), l, w, h, pmv); 

% Now compute the x, y, and z components of the Force vector:
Fx = pre_vx + wx;
Fy = pre_vy + wy;
Fz = pre_vz + wz;

% Now, we cap Force at Fmax dynes, which corresponds to the max accel. of 2 cm/s/frame

% compute the absolute value of the velocity using the 3d distance formula:

F3D = dist3d(0,0,0,Fx, Fy, Fz);

% Check to see if it's > 60 cm/s, the arbitrary imposed limit. If it is,
% we will have to reduce each component by a given amount.

if ( F3D > Fmax )
	% The percentage by which V3D exceeds 60:
	force_excess = F3D / Fmax; 
	
	% The reciprocal of that, which is the amount by which we will have to 
	% multiply each component vector to reduce it to the proper size:
	reduce_factor = 1 / force_excess;
	
	% Now multiply the compnents of velocity by reduce_vector 
	Fx = reduce_factor * Fx;
	Fy = reduce_factor * Fy;
	Fz = reduce_factor * Fz;
end


% This represents the "force" on the fish in "dynes". Realistically these are 
% "vitrual" dynes, since there's no real physical force pushing on the fish.
% However, it is a FORCE, and what we need is a velocity, so we must convert.

% ----------------------------------------------------------------------------
% CONVERT FORCE TO VELOCITY
% ----------------------------------------------------------------------------

%
% STEP 1: Convert Force to Acceleration
%

% Recall from Physics 101, that, if F = Force (in dynes, which are measured
% in the fun units of g * cm * s^(-2), m is mass (in g) and a is  acceleration
% (in cm * s^(-2), then F = m * a, which of course also means that a = F / m.
% Fish mass, therefore, is needed in grams. We will assume the fish is approximately
% the same density as water, which is 1 g / cm^3. Assume futher that the fish is
% 1 cm high, 1 cm across, and 4 cm long. This makes him 1x1x4 = 4 cm^3, which means
% each fish weighs (roughly) 4 g.  Thus, a = F / 4. We can change this later on if
% we decide that's not a good mass. Since F is a vector [Fx, Fy, Fz], we must
% perform this operation on each component.

m = mass;						% mass of fish, in grams.

ax = Fx / m;
ay = Fy / m;
az = Fz / m;

% Thus, for example, maximum acceleration, if Fx (max) is 240 dynes for a 1 s 
% duration, would be 240 / 4 or 60 cm/s^2.

%
% STEP 2: Convert acceleration to change in velocity.
%

% Recall from physics 101 that, if we let v represent velocity (in cm / 2) and t distancedistance
% represent time (in s), then since acceleration is in cm / s^2, a = v / t 
% (and the units cancel).  Since the acceleration is in seconds and cm, and the distancedistance
% distances are all in cm, we must convert time to s. Thus if we are "running the
% film" at 30 Hz, t = 1/30 s. Let hz stand for the number of frames in a second.
% Since a = v / t, and thus a * t  = v; and since t = 1 / hz, then a / hz = v, or
% in computer jargon, v = a / hz. This is the change in velocity PER FRAME.

vx = ax / hz;
vy = ay / hz;
vz = az / hz;

% Thus, if we are at maximum acceleration (60 cm/s^2, see STEP 1 above), 
% maximum velocity increment would be + 60 / 30 (assuming 30 Hz) or,
% + 2 cm / frame. This makes sense. If the fish is moving at, say, 1 cm/s
% and is frightened into a 60 cm/s^2 accel for a 1-second burst, 
% he would increase his velocity by +2 cm/s every frame for 30 frames,
% resulting in a final velocity of 61 cm/s, which corresponds (over
% that 1 second) to a 60 cm/s^2 acceleration.

% This allows us to update the coordinates:
DELTA_VEL = [vx vy vz];
FORCE = [Fx Fy Fz];		% Angle stores the force vector for now.

% ----------------------------------------------------------------------------
% UPDATE THE FISH'S MOVEMENT
% ----------------------------------------------------------------------------

% First, update the velocity by adding vectors:
OLD_VEL = VECTOR(n,6:8);

% Remeber, the velocity in VECTOR is in cm/s. Therefore, we add the
% vectors, but then we will divide by "hz" to find out the change in
% position.

VEL = OLD_VEL + DELTA_VEL;

% Now, we cap VEL at 60 cm/s, which is the maximum observed velocity of the
% fish (at least until we do motion analysis).

% compute the absolute value of the velocity using the 3d distance formula:

V3D = dist3d(0,0,0,VEL(1),VEL(2),VEL(3));

% Check to see if it's > 60 cm/s, the arbitrary imposed limit. If it is,
% we will have to reduce each component by a given amount.

if ( V3D > 60 )
	% The percentage by which V3D exceeds 60:
	speed_excess = V3D / 60; 
	
	% The reciprocal of that, which is the amount by which we will have to 
	% multiply each component vector to reduce it to the proper size:
	reduce_factor = 1 / speed_excess;
	
	% Now multiply the compnents of velocity by reduce_vector 
	VEL = reduce_factor * VEL;
end


OLD_COORDS = VECTOR(n,3:5);

COORDS = OLD_COORDS + ( VEL / hz );

% Now generate the output vector as a horizontal array 12 columns long
% 	column 1 = time step (already generated above)
%	column 2 = individual #
%	column 3-5 = coordinates in x-y-z
% 	column 6-8 = velocity x-y-z
%	column 9-11 = angle x-y-z
%	column 12 = 0 for 'not interpolated position' 
%			    (it's a model, so all positions are real)

% fsd
% Fmax

MVMT = [ VECTOR(n,1:2) COORDS VEL FORCE VECTOR(n,12) ];
