*  Move particles (fish) according to external force field (current)
*
*  Inputs :
*    nfish  = number of fish
*    fishp  = array of initial fish positions (stored as complex numbers)
*    fishv  = array of initial fish velocities (stored as complex numbers)
*    fishm  = array masses of fish
*    tfinal = final time for simulation
*    zoom   = range of display field
*
*  Outputs:
*    time intervals and root mean square velocity at each step
*    plot of fish positions after each time step
*
*  Algorithm: integrate using Euler's method
*
      program fish1
      include '/usr/cm5/cmf/include/cm/timer-fort.h'
      include '/usr/cm5/cmf/include/cm/CMF_defs.h'

*     Define the size of the display window
      integer, parameter :: m = 256
      integer, parameter :: disp_interval = 10
      include 'color_display.inc'

*     Max. No. of fish
      integer, parameter :: nfish = 10000

      complex 	fishp(nfish), fishv(nfish), force(nfish), accel(nfish)
      real      fishm(nfish), fishspeed(nfish)
CMF$  layout    fishp(:news),fishv(:news),fishm(:news)
CMF$  layout    force(:news),accel(:news)

      integer 	step,i
      real 	dt, t, tfinal, mnsqvel, zoom
      complex, parameter :: ii = (0,1) 	
C
C Initialize fishes, time, and display parameters
C
      force = (0.,0.)
      accel = (0.,0.)
      forall (i=1:nfish) fishp(i) = (float(i)*2.0/nfish)-1.0
      fishv = ii*fishp
      forall (i=1:nfish) fishm(i) = 1.0+float(i)/nfish

      tfinal = 30.
      zoom   = 40.
*
* Loop over time steps
      step = 0
      dt   = .01
      t    = 0.0
      do while (t < tfinal) 
         t = t + dt
*
* Update fish positions and velocities
         fishp   = fishp + dt*fishv
         call compute_current(force,fishp)
         accel   = force/fishm
         fishv   = fishv + dt*accel
* Update time step
         fishspeed = abs(fishv)
         mnsqvel = sqrt(sum(fishspeed*fishspeed)/nfish)
         print 100,t,mnsqvel
         dt      = .1*maxval(fishspeed) / maxval(abs(accel))
*
* Show fish positions
         if (mod(step,disp_interval) .eq. 0) then
            call show_fish(fishp,nfish,zoom,m)
         endif
         step = step + 1
      enddo
 100  format(2f20.10)
      stop
      end

      subroutine compute_current(force,fishp)
      complex force(:),fishp(:)
      complex, parameter :: ii = (0,1) 	
CMF$  layout force(:news),fishp(:news)

      force = (3,0)*(fishp*ii)/(max(abs(fishp),0.01)) - fishp
      end


      subroutine show_fish(fishp,nfish,zoom,m)
      complex fishp(nfish)
      real zoom
      integer 	x(nfish), y(nfish), show(m, m)
      integer 	j

* Scale the coordinates according to the zoom, and shift the origin 
      x = INT(( real(fishp)+zoom)/(zoom)*(m/2)) 
      y = INT((aimag(fishp)+zoom)/(zoom)*(m/2))
      show = 0 
      forall(j=1:nfish) show(x(j),y(j)) = 1
*
* enlarge each pixel for better visibility
      show(:,1:m-1) = show(:,1:m-1) + show(:,2:m)
      show(1:m-1,:) = show(1:m-1,:) + show(2:m,:)
      call display(show,m,m)
      end