% Move particles (fish) in periodic box according to gravity
% Inputs (for a sample set of inputs, run fish2init):
% fishp = column vector of initial fish positions (stored as complex numbers)
% fishv = column vector of initial fish velocities (stored as complex numbers)
% fishm = column vector 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, clg
% Set initial time step
dt = .03;
% Set up permutation array for computing gravity
nfish = length(fishp);
perm = [nfish,(1:nfish-1)];
% Initialize array of mean square velocities and times
mnsqvel=[];
time=[];
%
% loop over time steps
t = 0;
i = 0;
while t < tfinal;
% compute forces on fish from gravity
force = 0;
fishpp=fishp;
fishmp=fishm;
for k=1:nfish-1,
fishpp = fishpp(perm);
fishmp = fishmp(perm);
dif = fishpp - fishp;
force = force + fishmp .* fishm .* dif ./ max(abs(dif).^2,1e-6) ;
end
% update time, fish positions and fish velocities
t = t + dt;
i = i + 1;
fishp = fishp + dt * fishv;
accel = force ./ 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,'+w');
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 versus time
plotoutput