c autocorrelation function
c
      subroutine autocorr(nsamp,xvar,yvar,step,rcorr,ncorr)
      implicit none
      integer nsamp,ncorr
      double precision xvar(nsamp),yvar(nsamp),step(nsamp),
     &                 rcorr(nsamp)
      
      double precision ymean,sum,hmin,sig0,crit
      integer i,j,k,n
      
      sum=0.0d0
      do i=1,nsamp
        sum=sum+yvar(i)
      enddo
      ymean=sum/dble(nsamp)
      
      step(1)=0.0d0
      rcorr(1)=1.0d0
      sum=0.0d0
      do i=1,nsamp
        sum=sum+(yvar(i)-ymean)*(yvar(i)-ymean)
      enddo
      sig0=sum/dble(nsamp)
      if(sig0.eq.0.0d0)then
        do i=1,nsamp
          step(i)=-9999.0d0
          rcorr(i)=-9999.0d0
        enddo
        return
      endif
            
      hmin=xvar(2)-xvar(1)
      do i=2,nsamp-1
        if((xvar(i+1)-xvar(i)).lt.hmin)then
          hmin=xvar(i+1)-xvar(i)
        endif
      enddo
      crit=0.01d0*hmin
 
      
      i=2
      step(i)=hmin
           
            
10    j=1
      n=0
      sum=0.0d0
100   do k=j+1,nsamp
        if((xvar(k)-xvar(j)).lt.(step(i)+crit).and.
     &     (xvar(k)-xvar(j)).gt.(step(i)-crit))then
          n=n+1
          sum=sum+(yvar(k)-ymean)*(yvar(j)-ymean)
        endif
      enddo
      if(j.lt.nsamp-1)then
        j=j+1
        goto 100
      endif
      if(n.le.1)then
        step(i)=step(i)+hmin
      else
! this form of autocorrelation has less bias
        rcorr(i)=sum/(dble(n)*sig0)
     
! but this form is more common in statistic literature because it has certain
! desirable properties
!        rcorr(i)=sum/(dble(nsamp)*sig0)
        i=i+1    
        step(i)=step(i-1)+hmin     
      endif
      if((xvar(1)+step(i)).lt.xvar(nsamp))goto 10
      ncorr=i-1
      return
      end