!/***************************************************************************
! *   Copyright (C) 2004 by Claudio Attaccalite                             *
! *   claudio@freescience.info                                              *
! *                                                                         *
! *   This program is free software; you can redistribute it and/or modify  *
! *   it under the terms of the GNU General Public License as published by  *
! *   the Free Software Foundation; either version 2 of the License, or     *
! *   (at your option) any later version.                                   *
! *                                                                         *
! *   This program is distributed in the hope that it will be useful,       *
! *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
! *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
! *   GNU General Public License for more details.                          *
! *                                                                         *
! *   You should have received a copy of the GNU General Public License     *
! *   along with this program; if not, write to the                         *
! *   Free Software Foundation, Inc.,                                       *
! *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
! ***************************************************************************/
!
!***********************************************************************
!     
!     Program for calculating autocorrelation functions
!
!***********************************************************************

program auto
      
implicit none
      
integer i,j,ncol,ndat,ndiv,skip
character(30) :: ifile,ofile
real(8) dtim,time,sum,tmp
real(8), DIMENSION(:), ALLOCATABLE :: c,u,t

 write(*,*) ' Autocorrelation Program'
 write(*,'("Input file name :")',ADVANCE='NO')
 read(5,*) ifile
 write(*,'("Output file name :")',ADVANCE='NO')
 read(5,*) ofile
 write(*,'("Initial data to skip (skip>=0) :")',ADVANCE='NO')
 read(5,*) skip

 write(*,'("Select the coloumn (col>1 first=time):")',ADVANCE='NO')
 read(*,*)ncol

 if(ncol.le.1) then
   write(*,*) ' Error! ncol <= 1 ' 
    stop
  endif

! input file
      open(7,file=ifile,form='formatted',status='old')

        read(7,*,end=100) 
      do i=1,100000
        read(7,*,end=100) 
      enddo

100  write(*,*) ' Total number of points = ',i
      ndiv=i-skip-2
      ALLOCATE(c(ndiv),u(ndiv),t(ndiv))

      rewind(7)

      read(7,*,end=100) 
      do i=1,skip
        read(7,*)
      enddo

      do i=1,ndiv
        read(7,*) t(i),(u(i),j=1,ncol-1)
      enddo
      
      close(7)

!    time interval
      
      dtim=t(2)-t(1)
      
!     calculate average
      
      sum=0.d0
      do i=1,ndiv
        sum=sum+u(i)
      enddo
      sum=sum/dble(ndiv)
      write(*,*)'Average :',sum

!     subtract average
      
      do i=1,ndiv
        u(i)=u(i)-sum
      enddo

!     calculate correlation function
      
      do i=1,ndiv
        c(i)=0.d0
        do j=1,ndiv+1-i
          c(i)=c(i)+u(j)*u(i+j-1)
        enddo
      enddo

!     normalise the correlation function

      u(1)=c(1)/dble(ndiv)

      write(*,*) 'Normalisation :',u(1)
      do i=2,ndiv
        u(i)=(c(i)/dble(ndiv+1-i))/u(1)
      enddo
      u(1)=1.d0

!	 correlation file
      
      open(8,file=ofile,form='formatted',status='unknown')

!     write out correlation function
      
! Error estimated using formula of
! Phys Review Vol 22 N 1 Page 280 Eq. 24
! for a discussion see Phys Rev A Vol 36 Number 2 Page 958

      do i=1,ndiv
        time=dble(i-1)*dtim
        write(8,'(3e)') time,u(i),dsqrt(2.d0*i*(1-u(i))**2/dble(ndiv))
      enddo
      close(8)
      write(*,*) ' done '

      end