!     --------------------------------------------------------------      
!     Draw the flowlines of a pipe flow scenario with obstacle in 
!     the mid-stream. Compile as:  
!     f90 -O3 -n32 -r10000 flowline.f90 CG_MINIMIZER.o ESTIMATOR.o
!     --------------------------------------------------------------      
      PROGRAM FLOWLINE          
      USE ESTIMATOR
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (IDROP=4,MAXSTEP=20000)
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: X,Y,VX,VY
      INTEGER ORDER_T_X, ORDER_T_Y, ORDER_VX_X, ORDER_VX_Y, &
              ORDER_VY_X, ORDER_VY_Y
      DATA LP_DAT,LP_OUT,NLINE /25,26,2/      
      DOUBLE PRECISION DROPX(IDROP),DROPY(IDROP)
!     DATA DROPX /78,78, 80,80,     80,80 /
!     DATA DROPY /34,16, 31.7,18.3, 21,29 /
      DATA DROPX /78,78, 80,80 /
      DATA DROPY /34,16, 21,29 /
      
      XMIN = 20.D0
      XMAX = 140.D0
      YMIN = 0.D0
!     YMAX = 25.D0
      YMAX = 50.D0

      NIN  = 0
!     FIRST COUNT HOW MANY ATOMS ARE INSIDE THE REGION:
      OPEN (UNIT=LP_DAT,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="par.dat")
      READ (LP_DAT,*) NP
      DO I=1,NP
      READ (LP_DAT,*) XI,YI,VXI,VYI
!     SELECT A REGION
      IF (XI.GT.XMIN.AND.XI.LT.XMAX) THEN   
      NIN = NIN+1
      ENDIF
      ENDDO
      CLOSE (LP_DAT)

      ALLOCATE (X(NIN),Y(NIN),VX(NIN),VY(NIN))

      IP = 0
      VYSUM = 0.
      OPEN (UNIT=LP_DAT,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="par.dat")
      READ (LP_DAT,*) NP
      DO I=1,NP
      READ (LP_DAT,*) XI,YI,VXI,VYI
      IF (XI.GT.XMIN.AND.XI.LT.XMAX) THEN 
      IP = IP+1
      X(IP) = XI 
      Y(IP) = YI 
      VX(IP) = VXI 
      VY(IP) = VYI
!     IF (Y(IP).GT.25.) THEN 
!     Y(IP) = 50.-Y(IP)
!     VY(IP) = -VY(IP)
!     ENDIF
      VYSUM = VYSUM+VY(IP)
      ENDIF
      ENDDO
      CLOSE (LP_DAT)

      PRINT *, 'Average VY = ',VYSUM/NIN
      DO IP=1,NIN
      VY(IP) = VY(IP)-VYSUM/NIN
      ENDDO
      
      ORDER_T_X = 1
      ORDER_T_Y = 1
      ORDER_VX_X = 11
      ORDER_VX_Y = 11
      ORDER_VY_X = 11
      ORDER_VY_Y = 11
      
      CALL ESTIMATOR_INITIALIZE (NIN, DIMENSION=2,     & 
      XMIN=XMIN, XMAX=XMAX, YMIN=YMIN, YMAX=YMAX,      &
      ORDER_T_X = ORDER_T_X,  ORDER_T_Y = ORDER_T_Y,   &
      ORDER_VX_X = ORDER_VX_X, ORDER_VX_Y = ORDER_VX_Y, & 
      ORDER_VY_X = ORDER_VY_X, ORDER_VY_Y = ORDER_VY_Y)
      
      CALL ESTIMATE_FIELD &
      (PX=X(1:NIN),PY=Y(1:NIN),PVX=VX(1:NIN),PVY=VY(1:NIN), &
       TOL=1D-2,MAXI=1)      

      OPEN (UNIT=LP_OUT,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="traj.out")
      DELT = 0.1
!     INTEGRATE THE FLOW LINE
      DO 20 IY = 1,NLINE+IDROP
      IF (IY.LE.NLINE) THEN 
      XNOW = XMIN
      YNOW = (YMAX-YMIN)/NLINE*(IY-0.5)
      ELSE
      XNOW = DROPX(IY-NLINE)
      YNOW = DROPY(IY-NLINE)
      ENDIF
      CALL GIVE_ESTIMATE (X=XNOW,Y=YNOW,T=TNOW,VX=VXNOW1,VY=VYNOW1)
      CALL GIVE_ESTIMATE (X=XNOW,Y=50-YNOW,T=TNOW,VX=VXNOW2,VY=VYNOW2)
      VXNOW = (VXNOW1+VXNOW2)/2.
      VYNOW = (VYNOW1-VYNOW2)/2.
!     CALL GIVE_ESTIMATE (X=XNOW,Y=YNOW,T=TNOW,VX=VXNOW,VY=VYNOW)
      XOLD = XNOW-DELT*VXNOW
      YOLD = YNOW-DELT*VYNOW
      ISTEP = 0
  10  CONTINUE
      WRITE (LP_OUT,*) XNOW,YNOW
      XNEW = XOLD+2*DELT*VXNOW
      YNEW = YOLD+2*DELT*VYNOW
      XOLD = XNOW
      YOLD = YNOW
      XNOW = XNEW
      YNOW = YNEW
      CALL GIVE_ESTIMATE (X=XNOW,Y=YNOW,T=TNOW,VX=VXNOW1,VY=VYNOW1)
      CALL GIVE_ESTIMATE (X=XNOW,Y=50-YNOW,T=TNOW,VX=VXNOW2,VY=VYNOW2)
      VXNOW = (VXNOW1+VXNOW2)/2.
      VYNOW = (VYNOW1-VYNOW2)/2.
!     CALL GIVE_ESTIMATE (X=XNOW,Y=YNOW,T=TNOW,VX=VXNOW,VY=VYNOW)
      IF (YNOW.LT.0) YNOW = 0
!     IF (YNOW.GT.25) YNOW = 25
      IF (YNOW.GT.50) YNOW = 50
      ISTEP = ISTEP+1
      IF (XNOW.GT.XMAX.OR.ISTEP.GE.MAXSTEP) GOTO 20
      GOTO 10
  20  CONTINUE
      CLOSE (LP_OUT)

      DEALLOCATE (X,Y,VX,VY)
      STOP
      END PROGRAM FLOWLINE