%  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, clf
%  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,'+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 versus time
      plotoutput