%  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