      program ana_tauint ! Bernd Berg, Oct 3 2008, tau_int calculation.
C Prepares plot for t-dependent integrated autocorrelation time and 
C tries to get an estimate. Calculates the average observable.
C Approach: Copies all nrpt segments into a single time series.
C Also averages and error bars for space- and time-like plaquettes.
      include '../../Libs/Fortran/implicit.sta'
      character cd*1,clx*3,clt*3
      include '../../Libs/Fortran/constants.par'
      include 'lat.par'
      include 'mc.par'
C For autocorrelation calculations:
C nt     Maximum time distance considered.
C nrpt  Number of bins used for jackknife. 
C nskip  Number of data skipped in each step of the analysis. 
C        Reduces the length of the time series by 1/(nskip+1).
C nstep  Stepsize used in the data file for the plot.
C Set observable iobs. Possibly define in code with a new if statement.
      parameter(nt=100,nskip=1,nstep=1,lmean=.true.)
      parameter(iobs=1,nskp1=nskip+1,nhist=102)
      parameter(nall=nrpt*((nmeas+nskip)/nskp1))
      include '../../Libs/Fortran_par/common_lat.f'
      include '../../Libs/Fortran_par/common_tsm.f'
      dimension acorj(nrpt,0:nt),acor(0:nt),work(nrpt),
     &  acintJ(nrpt,0:nt),acint(0:nt),ace(0:nt),hist(nhist), 
     &  tsall(nall),act(nrpt),acpt(nrpt),ws(nrpt),wt(nrpt)
      include 'lat.dat'
      ltest=.true.
      ltest=.false.
C
      write(iuo,'(/," ana_tauint: nt,nrpt,nskip,nstep,lmean =",
     &              I9,3I4,L2)') nt,nrpt,nskip,nstep,lmean
      if(iobs.eq.1) write(iuo,'(/," Plaquette Action iobs =",I2)') iobs
      if(iobs.eq.2) write(iuo,'(/," Spacelike Plq. iobs =",I2)') iobs
      if(iobs.eq.3) write(iuo,'(/," Timelike Plq. iobs =",I3)') iobs
      write(cd,'(I1.1)') nd ! Prepare to read the data.
      write(clx,'(I3.3)') nla(1)
      write(clt,'(I3.3)') nla(nd)
      write(cfile,'(A,"d",A,"x",A,"t")') cd,clx,clt
      open(iud1,file=cfln//cfile//".D",form="unformatted",status="old")
      read(iud1) beta_in,nd_in,nla,nlink,nequi_in,nrpt_in,nmeas_in,
     &           nsw_in
      write(iuo,'(/," beta,nlink,nd,nla:",F15.10,I10,I3,2X,4I5)')
     &                beta,nlink,nd_in,nla
      write(iuo,'(" nequi,nrpt,nmeas: ",I15,2I10)')
     &              nequi_in,nrpt_in,nmeas_in
      if(nd.ne.nd_in) stop "nd.ne.nd_in."
      if(nequi.ne.nequi_in) stop "nequi.ne.nequi_in."
      if(nrpt.ne.nrpt_in) stop "nrpt.ne.nrpt_in."
      if(nmeas.ne.nmeas_in) stop "nmeas.ne.nmeas_in."
      if(nsw.ne.nsw_in) stop "nsw.ne.nsw_in."
      if(ltest) stop "ana_tauint: ltest."
C
      ndat=0
      do irpt=1,nrpt
        read(iud1) tsa,tsws,tswt,tspr,tspi,acpt(irpt),a_min,a_max,i_in
        act(irpt)=stmean(nmeas,tsa)
        ws(irpt)=stmean(nmeas,tsws)
        wt(irpt)=stmean(nmeas,tswt)
        if(irpt.ne.i_in) stop "irpt.ne.i_in incomplete read."
        if(iobs.ne.1) then ! Redefine tsa to desirec observable.
          do imeas=1,nmeas
            if(iobs.eq.2) tsa(imeas)=tsws(imeas)
            if(iobs.eq.3) tsa(imeas)=tsws(imeas)
          enddo
        endif
        do imeas=1,nmeas,nskp1 ! Stepsize nskip+1.
          ndat=ndat+1
          if(ndat.gt.nall) stop "ana_tauint: ndat.gt.nall."
          tsall(ndat)=tsa(imeas)
        enddo
      enddo
      close(iud1)
      write(iuo,'(" ndat =",I16)') ndat
C
      call steb0(nrpt,act,actm,av,actme)
      write(iuo,'(" actm = ",F15.10," +/-",F15.10)') actm,actme
      call steb0(nrpt,acpt,acptm,av,acptme)
      write(iuo,'(" acptm =",F15.10," +/-",F15.10)') acptm,acptme
      call steb0(nrpt,ws,wsm,wsv,wsme)
      write(iuo,'(" wsm =  ",F15.10," +/-",F15.10)') wsm,wsme
      call steb0(nrpt,wt,wtm,wtv,wtme)
      write(iuo,'(" wtm =  ",F15.10," +/-",F15.10,/)') wtm,wtme
C
      do it=0,nt
        acor(it)=autcorf(it,ndat,tsall,lmean)
        call autcorj(it,ndat,nrpt,tsall,work,acorj(1,it),lmean)
        call stebj0(nrpt,acorj(1,it),acm,acv,ace(it))
        call write_PROGRESS(IUD2,"it,nt,acm:",it,nt,acm)
      enddo
C
      open(iud1,file="1a"//cd//"d"//clx//"x"//clt//"t.d",
     &          form="formatted",status="unknown")
      call ac_int(nt,acor,acint)
      call ac_intj(nt,nrpt,acorj,acintJ)
      lmonoton=.true.
      tau_int=zero
      do it=0,nt,nstep
        call stebj0(nrpt,acintJ(1,it),acm,acv,ace(it))
        bias=(nrpt-1)*ABS(acint(it)-acm) ! Bias: Compare with
        itt=nskp1*nsw*it                  ! signal and error bar.
        acint(it)=nskp1*nsw*acint(it)
        ace(it)=nskp1*nsw*ace(it)
        bias=nskp1*nsw*bias
        write(iud1,'(I6,3G15.6)') itt,acint(it),ace(it),bias
        if(lmonoton) then
          if(acint(it).GE.tau_int) then
            tau_int=acint(it)
            tau_inte=ace(it)
            ebias=bias ! Estimate of bias at ittest.    
            ittest=nskp1*nsw*it
          else
            lmonoton=.false. 
          endif
        endif
      enddo
      close(iud1)
      write(iuo,*)"Autocorrelation times are on the 1a*.d file. Bias is"
      write(iuo,'(1X,"in 3rd column. Estimates in ",
     &               "file are NOT bias-corrected.")')
      write(iuo,'(1X,"Plot with tau_int.plt and extract estimate.")')
      write(iuo,'(/,2X,"A somewhat unreliable automatic pick: ittmax =",
     &           I10,/)') (nskp1*nsw*nt)
      write(iuo,'(5X,"ittest",4X,"tau_int",7X,"tau_inte",10X,"bias")')
      write(iuo,'(/,I10,3G15.6)') ittest,tau_int,tau_inte,Ebias
      write(iuo,'(10X,1G15.6," bias-corrected.")') (tau_int+ebias)
      xmin=two
      xmax=one
      miud=-iud1
      call hist_gnu(miud,nhist,ndat,hist,tsall,xmin,xmax)
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/bias.f'
      include '../../Libs/Fortran/bining.f'
      include '../../Libs/Fortran/datjack.f'
      include '../../Libs/Fortran/hist_gnu.f'
      include '../../Libs/Fortran/steb0.f'
      include '../../Libs/Fortran/stebj0.f'
      include '../../Libs/Fortran/stmean.f'
      include '../../Libs/Fortran/write_progress.f'
