      program ana_scalar1 ! Bernd Berg, Feb 19 2012.
C Analysis of scalar field data, 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. 
      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 iobs=1 action, iobs=2 field.
C ntau   Maximum time distance considered for tau_int calculation.
C        Normally nbins=nrpt. Can be diferent.
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.
c     parameter(iobs=1,ntau=2**5,nskip=2**0-1,nstep=1,lmean=.true.)
      parameter(iobs=2,ntau=2**6,nskip=2**0-1,nstep=1,lmean=.false.)
C nbins  Number of bins used for jackknife. 
      parameter(nbins=nrpt,nskp1=nskip+1,nhist=102)
      parameter(nall=nrpt*((nmeas+nskip)/nskp1))
      include '../../Libs/Fortran_par/common_lat.f'
      include '../../Libs/Scalar_par/common_scalar.f'
C act(nrpt)    average action.
C fld(nrpt)    average field.
      dimension acorj(nbins,0:ntau),acor(0:ntau),work(nbins),
     &  acintJ(nbins,0:ntau),acint(0:ntau),ace(0:ntau),hist(nhist), 
     &  tsall(nall),ann(nrpt),ac2(nrpt),act(nrpt),fld(nrpt),acpt(nrpt)
      include 'lat.dat'
      ltest=.true.
      ltest=.false.
C     if(ltest) stop "ana_scalar: ltest."
C
      write(iuo,'(/," ana_scalar: ntau,nbins,nskip,nstep,lmean =",
     &              I9,3I4,L2)') ntau,nbins,nskip,nstep,lmean
      if(iobs.eq.1) write(iuo,'(/," Scalar Action iobs =",I2)') iobs
      if(iobs.eq.2) write(iuo,'(/," Field tau_int iobs =",I2)') 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) xm2_in,xla_in,nd_in,nla,nlink,nequi_in,nrpt_in,
     &           nmeas_in,nsw_in
      write(iuo,'(/," xm2,xla,nlink,nd:",2F15.10,I10,I3,2X)')
     &                xm2,xla,nlink,nd_in
      write(iuo,'(/," nla:",4I5)') nla
      write(iuo,'(" nequi,nrpt,nmeas,nsw: ",I15,2I10,I6)')
     &              nequi_in,nrpt_in,nmeas_in,nsw_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."
C
      ndat=0
      do irpt=1,nrpt
        read(iud1) tsann,tsa2,tsphi,acpt(irpt),phi0sum,
     &             cor0,cor1r,cor1i,i_in
        if(irpt.ne.i_in) stop "irpt.ne.i_in incomplete read."
        do imeas=1,nmeas,nskp1 ! Stepsize nskip+1.
          ndat=ndat+1
          if(ndat.gt.nall) stop "ana_scalar: ndat.gt.nall."
          if(iobs.eq.1) tsall(ndat)=tsann(imeas) ! -tsa2(imeas)
          if(iobs.eq.2) tsall(ndat)=tsphi(imeas)
        enddo
        ann(irpt)=stmean(nmeas,tsann)
        ac2(irpt)=stmean(nmeas,tsa2)
        act(irpt)=half+ann(irpt)-ac2(irpt)
        fld(irpt)=stmean(nmeas,tsphi)
      enddo
      close(iud1)
      write(iuo,'(" ndat =",I16)') ndat
C
      call steb0(nrpt,ann,annm,av,annme)
      write(iuo,'(" annm = ",F15.10," +/-",F15.10)') annm,annme
      call steb0(nrpt,ac2,ac2m,av,ac2me)
      write(iuo,'(" ac2m = ",F15.10," +/-",F15.10)') ac2m,ac2me
      call steb0(nrpt,act,actm,av,actme)
      write(iuo,'(" actm = ",F15.10," +/-",F15.10)') actm,actme
      call steb0(nrpt,fld,fldm,fv,fldme)
      write(iuo,'(" fldm = ",F15.10," +/-",F15.10)') fldm,fldme
      call steb0(nrpt,acpt,acptm,av,acptme)
      write(iuo,'(" acptm =",F15.10," +/-",F15.10)') acptm,acptme
C
      do it=0,ntau
        acor(it)=autcorf(it,ndat,tsall,lmean)
        call autcorj(it,ndat,nbins,tsall,work,acorj(1,it),lmean)
        call stebj0(nbins,acorj(1,it),acm,acv,ace(it))
        call write_PROGRESS(IUD2,"it,ntau,acm:",it,ntau,acm)
      enddo
C
      open(iud1,file="1a"//cd//"d"//clx//"x"//clt//"t.d",
     &          form='formatted',status='unknown')
      call ac_int(ntau,acor,acint)
      call ac_intj(ntau,nbins,acorj,acintJ)
      lmonoton=.true.
      tau_int=zero
      do it=0,ntau,nstep
        call stebj0(nbins,acintJ(1,it),acm,acv,ace(it))
        bias=(nbins-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,"With ittmax =",I10," a somewhat unreliable ",
     &                 "automatic pick:",/)') (nskp1*nsw*ntau)
      write(iuo,'(/,2X,"A somewhat unreliable automatic pick:",/)')
      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/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'
