subroutine xcor(ss,ipk,nsteps,nsym,lag1,lag2,ccf,ccf0,lagpk,flip,fdot,nrobust) ! Computes ccf of a row of ss and the pseudo-random array pr. Returns ! peak of the CCF and the lag at which peak occurs. For JT65, the ! CCF peak may be either positive or negative, with negative implying ! the "OOO" message. use jt65_mod parameter (NHMAX=3413) !Max length of power spectra parameter (NSMAX=322) !Max number of half-symbol steps real ss(NSMAX,NHMAX) !2d spectrum, stepped by half-symbols real a(NSMAX) real ccf(-11:540) data lagmin/0/ !Silence g77 warning save df=12000.0/8192. dtstep=0.5/df fac=dtstep/(60.0*df) do j=1,nsteps ii=nint((j-nsteps/2)*fdot*fac)+ipk if( (ii.ge.1) .and. (ii.le.NHMAX) ) then a(j)=ss(j,ii) endif enddo if(nrobust.eq.1) then ! use robust correlation estimator to mitigate AGC attack spikes at beginning ! this reduces the number of spurious candidates overall call pctile(a,nsteps,50,xmed) do j=1,nsteps if( a(j).ge.xmed ) then a(j)=1 else a(j)=-1 endif enddo endif ccfmax=0. ccfmin=0. do lag=lag1,lag2 x=0. do i=1,nsym j=2*i-1+lag if(j.ge.1 .and. j.le.nsteps) x=x+a(j)*pr(i) enddo ccf(lag)=2*x !The 2 is for plotting scale if(ccf(lag).gt.ccfmax) then ccfmax=ccf(lag) lagpk=lag endif if(ccf(lag).lt.ccfmin) then ccfmin=ccf(lag) lagmin=lag endif enddo ccf0=ccfmax flip=1.0 if(-ccfmin.gt.ccfmax) then do lag=lag1,lag2 ccf(lag)=-ccf(lag) enddo lagpk=lagmin ccf0=-ccfmin flip=-1.0 endif return end subroutine xcor