function [MVMT] = randMovement3d(VECTOR, n, rule, l, w, h, mass, hz); % movementrrnd(,N, , l, w, h, mass, hz) -- Compute fish movement based on the rule. % % The movementrnd() function is designed to take an 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 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) ];