      PROGRAM TAU_INT1 ! Bernd Berg Oct 26 2007.
C Calculates tau_int(t) from time series on formatted input file data.d.
C Output on tau_int.d. Plot with tau_int.plt.
      INCLUDE '../../Libs/Fortran/implicit.sta'
      INCLUDE '../../Libs/Fortran/constants.par'
C NDAT  Number of data
C NT    Maximum time distance considered.
C NSTEP Stepsize in the data file for the plot.
C LMEAN Calculation of mean value from data. TRUE works generically. 
C       For the case at hand FALSE would actually be better, as we 
C       know the mean to be zero.
C     PARAMETER(NDAT=2**14,NT=2**05,NBINS=2**5,NSTEP=2**0,LMEAN=.FALSE.)
      PARAMETER(NDAT=2**14,NT=2**05,NBINS=2**5,NSTEP=2**0,LMEAN=.TRUE.)
C NSW   Allows to set multiple stepsize in multiples of some unit.
      PARAMETER(NSW=1, IUO=6,IUD1=11,IUD2=12)
      DIMENSION ACORJ(NBINS,0:NT),ACOR(0:NT),WORK(NBINS)
      DIMENSION ACINTJ(NBINS,0:NT),ACINT(0:NT),ACE(0:NT)
      DIMENSION TSA(NDAT)
      LTEST=.TRUE.
      LTEST=.FALSE.
C
      WRITE(IUO,'(/," tau_int1: ndat =",I11,/,/,2X,
     & "Analyzes one time series. Gives jackknife error bars.")') NDAT
      OPEN(IUD1,FILE="data.d",FORM='formatted',STATUS='old')
      IDAT=0
1     IDAT=IDAT+1
        READ(IUD1,*,END=2) JDAT,TSA_IN
        IF(IDAT.GT.NDAT) STOP "TAU_INT1: Error IDAT.GT.NDAT."
        TSA(IDAT)=TSA_IN
      GOTO 1
2     CONTINUE
      CLOSE(IUD1)
      IDAT=IDAT-1
      IF(LTEST) PRINT*,"IDAT,JDAT,TSA(NDAT) =",IDAT,JDAT,TSA(NDAT)
      IF(NDAT.NE.IDAT) STOP "TAU_INGT: Adjust NDAT to data.d size."
      IF(LTEST) STOP "TAU_INT1: LTEST."
C
      DO IT=0,NT
        ACOR(IT)=AUTCORF(IT,ndat,tsa,LMEAN)
        CALL AUTCORJ(IT,ndat,NBINS,TSA,WORK,ACORJ(1,IT),LMEAN)
        CALL STEBJ0(NBINS,ACORJ(1,IT),ACM,ACV,ACE(IT))
        CALL WRITE_PROGRESS(IUD2,"IT,NT,ACM:",IT,NT,ACM)
      END DO
C
      OPEN(IUD1,FILE="tau_int.d",FORM="formatted",STATUS="unknown")
      CALL AC_INT(NT,ACOR,ACINT)
      CALL AC_INTJ(NT,NBINS,ACORJ,ACINTJ)
      TAU_INT=HALF*ACINT(0)
      LMONOTON=.TRUE.
      DO IT=0,NT,NSTEP
        CALL STEBJ0(NBINS,ACINTJ(1,IT),ACM,ACV,ACE(IT))
        BIAS=(NBINS-1)*(ACINT(IT)-ACM) ! Bias.
        WRITE(IUD1,'(I6,3G15.6)') (IT*NSW),(ACINT(IT)*NSW),
     &                            (ACE(IT)*NSW),(BIAS*NSW)
        IF(LMONOTON) THEN
          IF(ACINT(IT).GE.TAU_INT) THEN
            TAU_INT=ACINT(IT)
            TAU_INTE=ACE(IT)
            EBIAS=BIAS ! Estimate of bias at ITEST.    
            ITEST=IT
          ELSE
            LMONOTON=.FALSE. 
          ENDIF
        ENDIF
        CALL WRITE_PROGRESS(IUD2,"IT,NT,ACINT:",IT,NT,ACINT(IT))
      END DO
      CLOSE(IUD1)
      WRITE(IUO,'(/,2X,"Autocorrelation times are on tau_int.d.")')
      WRITE(IUO,'(2X,
     &  "Bias is in 3rd column. Estimates are NOT bias-corrected.")')
      WRITE(IUO,'(2X,"Plot with tau_int.plt and extract estimate.")')
      WRITE(IUO,'(/,2X,"A somewhat unreliable automatic pick:",/)')
      WRITE(IUO,'(5X,"ITEST",4X,"TAU_INT",7X,"TAU_INTE",10X,"BIAS")')
      WRITE(IUO,'(/,I10,3G15.6,/)') ITEST,TAU_INT,TAU_INTE,EBIAS
      WRITE(IUO,'(10X,1G15.6," bias-corrected.",/)') (TAU_INT+EBIAS)
C
      STOP
      END

      INCLUDE '../../Libs/Fortran/ac_int.f'
      INCLUDE '../../Libs/Fortran/ac_intj.f'
      INCLUDE '../../Libs/Fortran/autcorf.f'
      INCLUDE '../../Libs/Fortran/autcorj.f'
      INCLUDE '../../Libs/Fortran/datjack.f'
      INCLUDE '../../Libs/Fortran/stebj0.f'
      INCLUDE '../../Libs/Fortran/stmean.f'
      INCLUDE '../../Libs/Fortran/write_progress.f'
