*  Move particles (fish) in periodic box according to gravity
*  Inputs :
*    nfish = number of fish
*    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
*    tsteps = number of time steps
*  Outputs:
*    plot of fish positions after each time step
*
*  Algorithm: integrate using Euler's method and stepsize dt = .1
*
      program prob2
      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 :: n = 256, m = 256
      integer factor
*     Max. No. of fish
      integer, parameter :: MAXF = 128 	
      complex   fishp(MAXF), fishv(MAXF), force(MAXF), fishpp(MAXF)
CMF$  layout    fishp(:news), fishv(:news), force(:news), fishpp(:news)
      complex   dif(MAXF)
CMF$  layout    dif(:news)
      real      fishm(MAXF), fishmp(MAXF), show(m, n)
CMF$  layout    fishm(:news), fishmp(:news), show(:news,:news)
      integer   x(MAXF), y(MAXF)
CMF$  layout    x(:news), y(:news)
      integer   i, kd, nfish, j, tsteps
      real      dt, range, etime0, etime1, etime2

      call CMF_DESCRIBE_ARRAY(fishp)
      print*, 'Input...'
      print*, 'Give number of fish(max 128) : '
      read (*,*) nfish
      call cmf_randomize(100)
      call cmf_random(fishp,1.0)
      call cmf_random(fishv,1.0)
      fishm = 1
      tsteps = 200
      dt = .1
      kd = 2
      range = 40
*      read (*,*) (fishp(i), i=1,nfish)
*      read (*,*) (fishv(i), i=1,nfish)
*      read (*,*) (fishm(i), i=1,nfish)
*      read (*,*) tsteps
*      read (*,*) dt
*      read (*,*) kd
*      read (*,*) range
      print *, 'nfish =', nfish
      print *, 'fishp =', (fishp(i), i=1,nfish)
      print *, 'fishv =', (fishv(i), i=1,nfish)
      print *, 'fishm =', (fishm(i), i=1,nfish)
      print *, 'tsteps =', tsteps
      print *, 'dt =', dt
      print *, 'kd =', kd
      print *, 'range =', range

      show = 0
* enlarge each pixel for better visibility 
      factor = m / 128
*
* Initialize the timer
      call CM_timer_clear(0)
      call CM_timer_clear(1)
      call CM_timer_clear(2)
      call CM_timer_start(0)
*
* Loop over time steps
      do i=1, tsteps
*
* Compute forces on fish from gravity
	 print *,'i =',i
         call CM_timer_start(1)
         force = 0
         fishpp = fishp
         fishmp = fishm
         do k=1, nfish-1
	     fishpp(1:nfish) = cshift(fishpp(1:nfish), DIM=1, SHIFT=-1)
	     fishmp(1:nfish) = cshift(fishmp(1:nfish), DIM=1, SHIFT=-1)
	     dif = fishpp - fishp
	     force = force + fishmp * fishm * dif / (abs(dif)*abs(dif))
         enddo
*
* Update fish positions and velocities
         fishp = fishp + dt*fishv
         fishv = fishv + dt*force/fishm
         call CM_timer_stop(1)
*
* Display output every kd steps
          if (mod(i-1,kd) .eq. 0) then
* Scale the coordinates according to the range, and shift the origin
             call CM_timer_start(2)
             x = INT((real(fishp)+range)/(range)*(m/2))
             y = INT((aimag(fishp)+range)/(range)*(n/2))
             forall(j=1:nfish) show(x(j),y(j)) = show(x(j),y(j)) + 1
*
* The following enlarges 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 CM_timer_stop(2)
             call display (show,m,n,MINVAL(show),MAXVAL(show)+0.0001)
             call CM_timer_start(2)
             show = 0
             call CM_timer_stop(2)
          endif
      enddo
*
* Read the timer
      call CM_timer_stop(0)
      etime0 = CM_timer_read_elapsed(0)
      etime1 = CM_timer_read_elapsed(1)
      etime2 = CM_timer_read_elapsed(2)
      print *, ' Computation time =', etime1, ' sec'
      print *, ' Histogram time =', etime2, ' sec'
      print *, ' Total Elapsed time =', etime0, ' sec'

      stop
      end