! ----------------------------------------------------------------------- !
!      PROGRAM Estimate V2.1 (RECTANLINEAR BOUNDARY):                     !
!      ESTIMATE CONTINUOUS MACROSCOPIC FIELDS FROM PARTICLE DATA          !
!      USING MAXIMUM LIKELIHOOD INFERENCE, ASSUMING QUASI-EQUILIBRIUM     !
!      PHASE SPACE DISTRIBUTION.                                          !
!                                                                         !
!                                      MAR.13.1997                        !
!                                      DEVELOPED BY JU LI AT MIT          !
! ----------------------------------------------------------------------- !
      
      PROGRAM ESTIMATE
      USE GENERATOR
      USE ESTIMATOR
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: MYXX,TEMP,PV,PV2
      INTEGER MYORDER_T_X, MYORDER_T_Y, MYORDER_T_Z, MYORDER_VX_X, &
              MYORDER_VX_Y, MYORDER_VX_Z, MYORDER_VY_X, MYORDER_VY_Y,&
              MYORDER_VY_Z, MYORDER_VZ_X, MYORDER_VZ_Y, MYORDER_VZ_Z 
      CHARACTER *70 BUF
      DATA LP_TEMP,MESH /25,200/
      
!     =====================================================================
!     READ INSTRUCTIONS FROM INPUT FILE "CON":
      
      READ (*,'(A70)') buf
!     "Input control parameter 1:"
      read *, par1

      READ (*,'(/,A70)') buf
!     "Input control parameter 2:"
      read *, par2
      
      READ (*,'(/,A70)') buf
!     "Input number of particles:"
      read *, np

      READ (*,'(/,A70)') buf
!     "Dimensionality of the system:"
      read *, MYNDIM
      
      READ (*,'(/,A70)') buf
!     "ORDER_T_X    ORDER_T_Y    ORDER_T_Z"
      read *, MYORDER_T_X, MYORDER_T_Y, MYORDER_T_Z
      
      READ (*,'(/,A70)') buf
!     "ORDER_VX_X    ORDER_VX_Y    ORDER_VX_Z"
      read *, MYORDER_VX_X, MYORDER_VX_Y, MYORDER_VX_Z
      
      READ (*,'(/,A70)') buf
!     "ORDER_VY_X    ORDER_VY_Y    ORDER_VY_Z"
      read *, MYORDER_VY_X, MYORDER_VY_Y, MYORDER_VY_Z
      
      READ (*,'(/,A70)') buf
!     "ORDER_VZ_X    ORDER_VZ_Y    ORDER_VZ_Z"
      read *, MYORDER_VZ_X, MYORDER_VZ_Y, MYORDER_VZ_Z
      
      READ (*,'(/,A70)') buf
!     "Give the random number generator a kick:"
      read *, kick
      do i = 1,kick
      a = randsafe()
      ENDDO
!     ======================================================================
      
      ALLOCATE (MYXX(NP*MYNDIM+1))
      ALLOCATE (PV(NP*MYNDIM+1))
      ALLOCATE (PV2(NP))
      ALLOCATE (TEMP(NP))

      XMIN = 0.D0
      XMAX = 1.D0
      YMIN = 0.D0
      YMAX = 1.D0
      ZMIN = 0.D0
      ZMAX = 1.D0
      PI = DACOS(-1.D0)

!     GENERATE PARTICLE COORDINATES:
      DO N = 1,NP*MYNDIM
      MYXX(N) = RANDSAFE()
      ENDDO
!     CALCULATE THE FIELD AND ASSIGN DRIFT VELOCITY TO PARTICLES
      DO N = 1,NP
      CALL ORIG_FIELD (NDIM=MYNDIM, XX=MYXX(N:NP*MYNDIM:NP), T=TEMP(N), & 
                       VV=PV(N:NP*MYNDIM:NP), PARM1=PAR1, PARM2=PAR2) 
      ENDDO

      PRINT *, ' Generator stage 1:'
      PRINT *, ' Avg TEMP =', SUM(TEMP)/NP
      IF (MYNDIM.GE.1) PRINT *, ' Avg VX =', SUM(PV(1:NP))/NP
      IF (MYNDIM.GE.2) PRINT *, ' Avg VY =', SUM(PV(NP+1:2*NP))/NP
      IF (MYNDIM.GE.3) PRINT *, ' Avg VZ =', SUM(PV(2*NP+1:3*NP))/NP
      PRINT *, ' ' 
      
!     RANDOMIZE PARTICLE VELOCITIES WITH SPECIFIED TEMPERATURES:
      IP = 1 
      V2AVG = 0.D0
 78   S1=DSQRT(-2.D0*LOG(RANDSAFE()))
      S2=2*PI*RAN1(KICK)
      DV1 = SQRT(TEMP(MOD(IP-1,NP)+1))*S1*COS(S2)
      PV(IP) = PV(IP)+DV1
      DV2 = SQRT(TEMP(MOD(IP,NP)+1))*S1*SIN(S2)
      PV(IP+1) = PV(IP+1)+DV2
      V2AVG = V2AVG+DV1**2+DV2**2
      IP = IP+2
      IF (IP.LE.MYNDIM*NP) GOTO 78

      PRINT *, ' Generator stage 2:'
      PRINT *, ' Avg TEMP =', V2AVG/MYNDIM/NP
      IF (MYNDIM.GE.1) PRINT *, ' Avg VX =', SUM(PV(1:NP))/NP
      IF (MYNDIM.GE.2) PRINT *, ' Avg VY =', SUM(PV(NP+1:2*NP))/NP
      IF (MYNDIM.GE.3) PRINT *, ' Avg VZ =', SUM(PV(2*NP+1:3*NP))/NP
      PRINT *, ' ' 
            
      CALL ESTIMATOR_INITIALIZE (NP, DIMENSION=MYNDIM,	 & 
       ORDER_T_X = MYORDER_T_X,  ORDER_VX_Y = MYORDER_VX_Y)
      CALL ESTIMATE_FIELD (PXX=MYXX, PVV=PV)
      
!     GET THE TEMPERATURE AT TWO ENDS AND COMPARE FITTING WITH ORIGINAL
      CALL GIVE_ESTIMATE (X=XMIN, T=T_XMIN)
      CALL GIVE_ESTIMATE (X=XMAX, T=T_XMAX)
      CALL ORIG_FIELD (X=XMIN, T=ORIG_T_XMIN)
      CALL ORIG_FIELD (X=XMAX, T=ORIG_T_XMAX)
      PRINT *,'T_XMIN: ORIGINAL = ', ORIG_T_XMIN, 'FITTED = ', T_XMIN
      PRINT *,'T_XMAX: ORIGINAL = ', ORIG_T_XMAX, 'FITTED = ', T_XMAX
      PRINT *,'T_XMID: ORIGINAL = ', (ORIG_T_XMAX+ORIG_T_XMIN)/2,   & 
               'FITTED = ', (T_XMAX+T_XMIN)/2
      PRINT *, ' ' 
      
!     GET THE VELOCITY AT TWO ENDS AND COMPARE FITTING WITH ORIGINAL
      CALL GIVE_ESTIMATE (Y=YMIN, VX=VX_YMIN)
      CALL GIVE_ESTIMATE (Y=YMAX, VX=VX_YMAX)
      CALL ORIG_FIELD (Y=YMIN, VX=ORIG_VX_YMIN)
      CALL ORIG_FIELD (Y=YMAX, VX=ORIG_VX_YMAX)
      PRINT *,'VX_YMIN: ORIGINAL = ', ORIG_VX_YMIN, 'FITTED = ', VX_YMIN
      PRINT *,'VX_YMAX: ORIGINAL = ', ORIG_VX_YMAX, 'FITTED = ', VX_YMAX
      PRINT *,'VX_YMID: ORIGINAL = ', (ORIG_VX_YMAX+ORIG_VX_YMIN)/2,  & 
              'FITTED = ', (VX_YMAX+VX_YMIN)/2
      PRINT *, ' ' 
 
      OPEN (UNIT=LP_TEMP, STATUS='UNKNOWN', FORM='FORMATTED',FILE='temp.out')
!     GIVE THE AVERAGE DRIFT AND FLUCTUATION OF THE TEMPERATURE FIELD
      XDEL = (XMAX-XMIN)/(MESH-1)
      XC = XMIN
      T_SUM = 0.D0
      T_VAR = 0.D0
      DO I = 1,MESH
      CALL ORIG_FIELD (X=XC,T=T_ORIG)
      CALL GIVE_ESTIMATE (X=XC,T=T_FIT)
      WRITE (LP_TEMP,'(3(1X,F10.5))') XC,T_ORIG,T_FIT
      DIFF = T_FIT - T_ORIG
      T_SUM = T_SUM + DIFF
      T_VAR = T_VAR + DIFF*DIFF
      XC = XC+XDEL
      ENDDO
      CLOSE (LP_TEMP)
      
      PRINT *,'estimator: THE AVERAGE DRIFT IN T IS',T_SUM/MESH
      PRINT *,'estimator: THE AVERAGE VARIATION IN T IS',DSQRT(T_VAR/MESH)   
      
      STOP
      END PROGRAM ESTIMATE