% Move particles (fish) in periodic box according to external force field (current)
% Inputs (for a sample set of inputs, run fish1init):
% fishp = array of initial fish positions (stored as complex numbers)
% fishv = array of initial fish velocities (stored as complex numbers)
% fishm = array of masses of fish
% tfinal = final time for simulation (0 = initial time)
% zoom = size of plotting window
% Outputs:
% plot of fish positions after each time step
% plot of mean-square velocity of particles at each time step
% plot of time step of simulation at each step
%
% Algorithm: integrate using Euler's method with varying step size
%
% Set up axes for plotting
hold off, clf
% Set initial time step
dt = .01;
i = 0;
% Initialize array of mean square velocities and times
mnsqvel=[];
time=[];
%
% loop over time steps
t = 0;
while t < tfinal,
% update time, fish positions and fish velocities
t = t + dt;
i = i + 1;
fishp = fishp + dt*fishv;
accel = current(fishp)./fishm;
fishv = fishv + dt*accel;
% update time step
dt = min(.1*max(abs(fishv))/max(abs(accel)),1);
% update list of mean square velocities and times
mnsqvel=[mnsqvel,sqrt(sum(abs(fishv).^2)/length(fishp))];
time = [time,t];
% plot fish positions every 2 steps
if (rem(i,2) == 0),
plot(fishp,'+k');
axis([-zoom,zoom,-zoom,zoom]);axis('square');
title(['Time = ',num2str(t)]); pause(.75)
end
end
disp('Hit return to continue'), pause
%
% plot mean square velocities and timesteps
plotoutput