75 lines
1.7 KiB
FortranFixed
75 lines
1.7 KiB
FortranFixed
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 |