!-------------------------------------------------------------------------
      !  Program      :                     nbody
      !-------------------------------------------------------------------------
      !
      !  Purpose      : This program simulates two-dimensional, inviscid 
      !                 (zero viscosity) flow using a particle vortex method. 
      !                 The computational intensive part of this algorithm is
      !                 the solution of the N-body problem: the mutual 
      !                 interaction (the induced velocity) of the N particles.
      !
      !  Input        : stdin: number of grid points in x (and y)
      !                 stdin: overlap factor (q).
      !
      !  Input/output : 
      !
      !  Output       : error.<ngrid>.dat: time history of the velocity error 
      !                                    on the particles.
      !
      !                 moments.dat: time history of the vortcity moments.
      !
      !                 compare.dat: the final particle position and velocity
      !                              and the corresponding analytical velocity.
      !
      !                 grid.dat: final velocity field on the grid.
      !
      !  Routines     : velocity
      !                 analytical
      !                 moments
      !                 meshvelocity
      !
      !  Remarks      : The initial velocity error on the particles using the
      !                 O(4) scheme converge approximately as mq, where m is
      !                 the order of the scheme (4) and q is the overlap 
      !                 factor (0.5 < q < 1.0).
      !
      !  References   : M. Perlman 1985, J. Comput. Phys. 59, 200-223
      !                 J.T. Beale & A. Majda 1985, J. Comput. Phys. 58, 188-208
      !
      !  Revisions    :
      !-------------------------------------------------------------------------
      !  $Log$
      !-------------------------------------------------------------------------
      !  Jens Honore Walther
      !  Institute of Computational Science
      !  ETH Zentrum, CH-8092 Zurich, Switzerland
      !-------------------------------------------------------------------------

      PROGRAM nbody
      IMPLICIT NONE
      !-------------------------------------------------------------------------
      !  Includes 
      !-------------------------------------------------------------------------
      INCLUDE 'param.h'
      !-------------------------------------------------------------------------
      !  Local variables 
      !-------------------------------------------------------------------------
      REAL(MK), DIMENSION(:,:), ALLOCATABLE :: xp,vp,xtilde,vtilde,vort,xm,vm
      REAL(MK), DIMENSION(:  ), ALLOCATABLE :: gamma,gtilde
      REAL(MK) :: dx,dy,du,dv,dt,rr2,xt,yt,zt2,zt,dh
      REAL(MK) :: error,m0,m1x,m1y,m2,sigma,q
      INTEGER :: Np,Mp,istep,nstep,iCase,i,j,k,l,ngrid,info,ni,nj
      INTEGER :: schm_convection
      CHARACTER(LEN=256) :: filename

      !-------------------------------------------------------------------------
      !  Time stepping
      !-------------------------------------------------------------------------
      ngrid = 20
      WRITE(*,'(A)',ADVANCE='NO') 'Enter the number of grid points (ngrid,ngrid) '
      READ(*,*) ngrid
      Np    = ngrid**2
      Mp    = Np

      nstep = 25
      dt    = 0.8

      nstep = 100 ! number of time steps
      dt    = 0.2 ! the time step size (arbitrary units)

      schm_convection = 1 ! euler
      schm_convection = 3 ! rh3
      schm_convection = 2 ! heun
     
      !-------------------------------------------------------------------------
      !  Time stepping
      !-------------------------------------------------------------------------
      ALLOCATE(xp(2,Np),vp(2,Np),gamma(Np),xtilde(2,Np),vtilde(2,Np), &
     &         vort(ngrid,ngrid),xm(2,Mp),vm(2,Mp),STAT=info)
      IF (info.NE.0) THEN
         WRITE(*,'(A)') 'ERROR: allocating memory failed!'
         GOTO 9999
      ENDIF 

      !-------------------------------------------------------------------------
      !  Initialize the vorticity field 
      !-------------------------------------------------------------------------
      k   = 0
      l   = 0
      dh  = 2.0_MK/REAL(ngrid,MK)
      DO j=1,ngrid
         DO i=1,ngrid
            !-------------------------------------------------------------------
            !  Place the particle at the centre of the mesh cell
            !-------------------------------------------------------------------
            xt       = REAL(i-1,MK)*dh - 1.0_MK + 0.5_MK*dh
            yt       = REAL(j-1,MK)*dh - 1.0_MK + 0.5_MK*dh
            zt2      = xt**2 + yt**2
            zt       = SQRT(zt2)

            !-------------------------------------------------------------------
            !  Save the mesh position for the diagnostics later
            !-------------------------------------------------------------------
            l        = l + 1
            xm(1,l)  = xt
            xm(2,l)  = yt
            IF (zt.LT.1.0_MK) THEN
               vort(i,j) = (1.0_MK - zt2)**7
               k         = k + 1
               gamma(k)  = vort(i,j)*dh*dh
               xp(1,k)   = xt
               xp(2,k)   = yt
            ELSE
               vort(i,j) = 0.0_MK
            ENDIF 
         ENDDO
      ENDDO
      Np = k

      WRITE(*,'(A,I8)') 'Number of particles: ',Np

      !-------------------------------------------------------------------------
      !  Enter the overlap factor (q)
      !-------------------------------------------------------------------------
      WRITE(*,'(A)',ADVANCE='NO')'Enter 0.5 < q < 1 (delta = dh**q) '
      READ(*,*) q
      sigma = dh**q
      WRITE(*,'(A,E12.4)') 'Overlap                    : ',sigma/dh
      WRITE(*,'(A,E12.4)') 'dh                         : ',dh

      !-------------------------------------------------------------------------
      !  Initial error
      !-------------------------------------------------------------------------
      CALL velocity(xp,gamma,Np,sigma,vp,info)
      CALL analytical(xp,vtilde,np,info)

      !-------------------------------------------------------------------------
      !  Compute the RMS error
      !-------------------------------------------------------------------------
      OPEN(10,FILE='init.dat')
      error = 0.0_MK
!     WRITE(10,'(A,8A14)') '#','xp(1)','xp(2)','vp(1)','vp(2)','ve(1)','ve(2)','rr2','SQRT(rr2)'
      DO i=1,Np 
         dx    = vp(1,i) - vtilde(1,i)
         dy    = vp(2,i) - vtilde(2,i)
         rr2   = dx*dx + dy*dy
         error = error + rr2
!        WRITE(10,'(8E14.6)') xp(:,i),vp(:,i),vtilde(:,i),rr2,SQRT(rr2)
         WRITE(10,*) xp(:,i),vp(:,i)
      ENDDO
      error = SQRT(error/REAL(Np,MK))
      WRITE(*,'(A,E12.4)') 'The initial error is       : ',error
      CLOSE(10)

      !-------------------------------------------------------------------------
      !  Compute the file name for the time history of the error 
      !-------------------------------------------------------------------------
      WRITE(filename,'(A,I2,A)') 'error.',ngrid,'.dat'

      !-------------------------------------------------------------------------
      !  Open the output files
      !-------------------------------------------------------------------------
      OPEN(10,FILE=filename)
      OPEN(11,FILE='moment.dat')

      !-------------------------------------------------------------------------
      !  Time stepping
      !-------------------------------------------------------------------------
      DO istep=1,nstep
         !----------------------------------------------------------------------
         !  Advance the particles 
         !----------------------------------------------------------------------
         IF     (schm_convection.EQ.1) THEN
            !-------------------------------------------------------------------
            !  Euler time integration
            !-------------------------------------------------------------------
            CALL velocity(xp,gamma,Np,sigma,vp,info)
            DO i=1,Np
               xp(1,i) = xp(1,i) + dt*vp(1,i)
               xp(2,i) = xp(2,i) + dt*vp(2,i)
            ENDDO
         ELSEIF (schm_convection.EQ.2) THEN
            !-------------------------------------------------------------------
            !  Heuns method 
            !-------------------------------------------------------------------
            CALL velocity(xp,gamma,Np,sigma,vp,info)
            DO i=1,Np
               xtilde(1,i) = xp(1,i) + dt*vp(1,i)
               xtilde(2,i) = xp(2,i) + dt*vp(2,i)
            ENDDO
            CALL velocity(xtilde,gamma,Np,sigma,vtilde,info)
            DO i=1,Np
               xp(1,i) = xp(1,i) + 0.5_MK*dt*(vp(1,i) + vtilde(1,i))
               xp(2,i) = xp(2,i) + 0.5_MK*dt*(vp(2,i) + vtilde(2,i))
            ENDDO  
         ELSE
            !-------------------------------------------------------------------
            !  RK3 time integration
            !-------------------------------------------------------------------
            CALL velocity(xp,gamma,Np,sigma,vp,info)
            DO i=1,Np
               vtilde(1,i) = vp(1,i)
               vtilde(2,i) = vp(2,i)
               xp(1,i)  = xp(1,i) + (dt/3.0_MK)*vp(1,i)
               xp(2,i)  = xp(2,i) + (dt/3.0_MK)*vp(2,i)
            ENDDO
            CALL velocity(xp,gamma,Np,sigma,vp,info)
            DO i=1,Np
               vtilde(1,i) = (-5.0_MK/9.0_MK)*vtilde(1,i) + vp(1,i)
               vtilde(2,i) = (-5.0_MK/9.0_MK)*vtilde(2,i) + vp(2,i)
               xp(1,i)  = xp(1,i) + (15.0_MK*dt/16.0_MK)*vtilde(1,i)
               xp(2,i)  = xp(2,i) + (15.0_MK*dt/16.0_MK)*vtilde(2,i)
            ENDDO
            CALL velocity(xp,gamma,Np,sigma,vp,info)
            DO i=1,Np
               vtilde(1,i) = (-153.0_MK/128.0_MK)*vtilde(1,i) + vp(1,i)
               vtilde(2,i) = (-153.0_MK/128.0_MK)*vtilde(2,i) + vp(2,i)
               xp(1,i)  = xp(1,i) + (8.0_MK*dt/15.0_MK)*vtilde(1,i)
               xp(2,i)  = xp(2,i) + (8.0_MK*dt/15.0_MK)*vtilde(2,i)
            ENDDO
         ENDIF 

         !----------------------------------------------------------------------
         !  Compute the vorticity moments
         !----------------------------------------------------------------------
         CALL moments(xp,gamma,Np,m0,m1x,m1y,m2,info)
         WRITE(11,'(4E12.4)') m0,m1x,m1y,m2

         !----------------------------------------------------------------------
         !  Compare on the fly the error
         !----------------------------------------------------------------------
         IF (MOD(istep,5).EQ.0) THEN ! currently every 5 time step
            !-------------------------------------------------------------------
            !  Compute the particle velocity (again, since the particles moved)
            !-------------------------------------------------------------------
            CALL velocity(xp,gamma,Np,sigma,vp,info)

            !-------------------------------------------------------------------
            !  Compute the analytical velocity 
            !-------------------------------------------------------------------
            CALL analytical(xp,vtilde,np,info)

            !-------------------------------------------------------------------
            !  Compute the RMS error
            !-------------------------------------------------------------------
            error = 0.0_MK
            DO i=1,Np 
               error = error + (vp(1,i) - vtilde(1,i))**2 + &
      &                        (vp(2,i) - vtilde(2,i))**2
            ENDDO
            error = SQRT(error/REAL(Np,MK))
            WRITE(10,'(2E12.4)') istep*dt,error

            !-------------------------------------------------------------------
            !  Compute the file name for the particle snapshot
            !-------------------------------------------------------------------
            WRITE(filename,'(A,I4.4)') 'out',istep
            OPEN(12,FILE=filename)

            !-------------------------------------------------------------------
            !  Write the particle position and velocity to the file
            !-------------------------------------------------------------------
            DO i=1,Np
               WRITE(12,'(4E12.4)') xp(:,i),vp(:,i)
            ENDDO
            CLOSE(12)
         ENDIF
      ENDDO 
      CLOSE(10)
      CLOSE(11)

      !-------------------------------------------------------------------------
      !  Compute the error in the velocity on the particles at the end of the
      !  simulation
      !-------------------------------------------------------------------------
      CALL velocity(xp,gamma,Np,sigma,vp,info)
      CALL analytical(xp,vtilde,np,info)

      OPEN(10,FILE='compare.dat')
      error = 0.0_MK
      DO i=1,Np 
         WRITE(10,'(7E12.4)') xp(:,i),vp(:,i),vtilde(:,i),gamma(i)
         error = error + (vp(1,i) - vtilde(1,i))**2 + &
                         (vp(2,i) - vtilde(2,i))**2
      ENDDO
      error = SQRT(error/REAL(Np,MK))
      WRITE(*,'(A,E12.4)') 'Velocity error on particles: ',error
      CLOSE(10)

      !-------------------------------------------------------------------------
      !  Compute the error in the velocity on the mesh at the end of the 
      !  simulation
      !-------------------------------------------------------------------------
      CALL meshvelocity(xp,gamma,Np,sigma,xm,vm,Mp,info)
      CALL analytical(xm,vtilde,Mp,info)

      !-------------------------------------------------------------------------
      !  Write it out
      !-------------------------------------------------------------------------
      k = 0
      OPEN(10,FILE='grid.dat')
      error = 0.0_MK
      DO j=1,ngrid
         DO i=1,ngrid
            k = k + 1
            WRITE(10,'(6E12.4)') xm(:,k),vm(:,k),vtilde(:,k)
            error = error + (vm(1,k) - vtilde(1,k))**2 + &
                            (vm(2,k) - vtilde(2,k))**2
         ENDDO
         WRITE(10,'(A)')
      ENDDO
      error = SQRT(error/REAL(Mp,MK))
      WRITE(*,'(A,E12.4)') 'Velocity error on mesh     : ',error
      CLOSE(10)
 9999 CONTINUE
      END PROGRAM nbody

      !-------------------------------------------------------------------------
      !  Subroutine for solving the n-body problem 
      !-------------------------------------------------------------------------
      SUBROUTINE velocity(xp,gamma,Np,sigma,vp,info)
      IMPLICIT NONE
      INCLUDE 'param.h'
      REAL(MK), DIMENSION(2,*), INTENT(IN ) :: xp
      REAL(MK), DIMENSION(*)  , INTENT(IN ) :: gamma
      REAL(MK), DIMENSION(2,*), INTENT(OUT) :: vp
      INTEGER                 , INTENT(IN ) :: Np
      INTEGER                 , INTENT(OUT) :: info
      REAL(MK) :: dx,dy,r2,rr2,du,dv,scale,sigma,rsig2,expterm
      INTEGER  :: i,j,k
    
      info = 0 
      DO i=1,Np
         vp(1,i) = 0.0_MK
         vp(2,i) = 0.0_MK
      ENDDO 
      scale = 1.0_MK/(2.0_MK*ACOS(-1.0))
      rsig2 = 1.0_MK/sigma**2
      DO j=1,Np
         DO i=1,Np
            IF (i.EQ.j) CYCLE ! skip i = j
            dx  = xp(1,i) - xp(1,j)
            dy  = xp(2,i) - xp(2,j)
            r2  = dx*dx + dy*dy
#ifdef __BLOB
            !-------------------------------------------------------------------
            !  Vortex blob approximation
            !-------------------------------------------------------------------
#ifdef __SECOND_ORDER
            !-------------------------------------------------------------------
            !  Second order (BM p. 192)
            !-------------------------------------------------------------------
            expterm = EXP(-r2*rsig2)
            rr2 = scale*(1.0_MK - expterm)/r2
#else
            !-------------------------------------------------------------------
            !  Fourth order (BM p. 192)
            !-------------------------------------------------------------------
            expterm = EXP(-0.5_MK*r2*rsig2)
            rr2 = scale*(1.0_MK - expterm)*(1.0_MK + 2.0_MK*expterm)/r2 
#endif
#else
            !-------------------------------------------------------------------
            !  Point vortex
            !-------------------------------------------------------------------
            rr2 = scale/r2 ! point vortex
#endif
            du  = -dy*rr2
            dv  =  dx*rr2

            vp(1,j) = vp(1,j) - du*gamma(i)
            vp(2,j) = vp(2,j) - dv*gamma(i)
         ENDDO
      ENDDO
      RETURN
      END SUBROUTINE velocity

      SUBROUTINE analytical(xp,vp,np,info)
      IMPLICIT NONE
      INCLUDE 'param.h'
      REAL(MK), DIMENSION(2,*), INTENT(IN ) :: xp
      REAL(MK), DIMENSION(2,*), INTENT(OUT) :: vp
      INTEGER                 , INTENT(IN ) :: Np
      INTEGER                 , INTENT(OUT) :: info
      REAL(MK) :: dx,dy,z2,fz
      INTEGER :: i

      info = 0
      DO i=1,Np
         dx = xp(1,i)
         dy = xp(2,i)
         z2 = dx*dx + dy*dy
         IF (z2.LT.1.0_MK) THEN
            IF (z2.EQ.0.0_MK) THEN
               fz =  0.0_MK
            ELSE
               fz = -1.0_MK/(16.0_MK*z2)*(1.0_MK - (1.0_MK - z2)**8)
            ENDIF 
         ELSE
            fz = -1.0_MK/(16.0_MK*z2)
         ENDIF 
         vp(1,i) =  fz*xp(2,i)
         vp(2,i) = -fz*xp(1,i)
      ENDDO 
      RETURN
      END SUBROUTINE analytical

      SUBROUTINE moments(xp,gamma,Np,m0,m1x,m1y,m2,info)
      IMPLICIT NONE
      INCLUDE 'param.h'
      REAL(MK), DIMENSION(2,*), INTENT(IN ) :: xp
      REAL(MK), DIMENSION(*)  , INTENT(IN ) :: gamma
      REAL(MK)                , INTENT(OUT) :: m0,m1x,m1y,m2
      INTEGER                 , INTENT(IN ) :: Np
      INTEGER                 , INTENT(OUT) :: info
      INTEGER :: k

      info = 0
      m0  = 0.0_MK
      m1x = 0.0_MK
      m1y = 0.0_MK
      m2  = 0.0_MK
      DO k=1,Np 
         m0       = m0  + gamma(k)
         m1x      = m1x + gamma(k)*xp(1,k)
         m1y      = m1y + gamma(k)*xp(2,k)
         m2       = m2  + gamma(k)*(xp(1,k)**2 + xp(2,k)**2)
      ENDDO
      RETURN
      END SUBROUTINE moments

      SUBROUTINE meshvelocity(xp,gamma,Np,sigma,xm,vm,Mp,info)
      IMPLICIT NONE
      INCLUDE 'param.h'
      REAL(MK), DIMENSION(2,*), INTENT(IN ) :: xp
      REAL(MK), DIMENSION(*)  , INTENT(IN ) :: gamma
      REAL(MK), DIMENSION(2,*), INTENT(IN ) :: xm
      REAL(MK), DIMENSION(2,*), INTENT(OUT) :: vm
      INTEGER                 , INTENT(IN ) :: Np,Mp
      INTEGER                 , INTENT(OUT) :: info
      REAL(MK) :: dx,dy,r2,rr2,du,dv,scale,sigma,rsig2,expterm
      INTEGER  :: i,j
    
      info = 0 
      DO i=1,Mp
         vm(1,i) = 0.0_MK
         vm(2,i) = 0.0_MK
      ENDDO 
      scale = 1.0_MK/(2.0_MK*ACOS(-1.0))
      rsig2 = 1.0_MK/sigma**2
      DO j=1,Mp
         DO i=1,Np
            dx  = xp(1,i) - xm(1,j)
            dy  = xp(2,i) - xm(2,j)
            r2  = dx*dx + dy*dy
#ifdef __BLOB
            !-------------------------------------------------------------------
            !  Vortex blob approximation
            !-------------------------------------------------------------------
#ifdef __SECOND_ORDER
            !-------------------------------------------------------------------
            !  Second order (BM p. 192)
            !-------------------------------------------------------------------
            expterm = EXP(-r2*rsig2)
            rr2 = scale*(1.0_MK - expterm)/r2
#else
            !-------------------------------------------------------------------
            !  Fourth order (BM p. 192)
            !-------------------------------------------------------------------
            expterm = EXP(-0.5_MK*r2*rsig2)
            rr2 = scale*(1.0_MK - expterm)*(1.0_MK + 2.0_MK*expterm)/r2 
#endif
#else
            !-------------------------------------------------------------------
            !  Point vortex
            !-------------------------------------------------------------------
            rr2 = scale/r2 ! point vortex
#endif
            du  = -dy*rr2
            dv  =  dx*rr2
            vm(1,j) = vm(1,j) - du*gamma(i)
            vm(2,j) = vm(2,j) - dv*gamma(i)
         ENDDO
      ENDDO
      RETURN
      END SUBROUTINE meshvelocity