subroutine afc65b(cx,npts,fsample,nflip,mode65,a,ccfbest,dtbest)
! Find delta f, f1, f2 ==> a(1:3)
complex cx(npts)
real a(5),deltaa(5)
a(1)=0.
a(2)=0.
a(3)=0.
a(4)=0.
deltaa(1)=2.0*mode65
deltaa(2)=2.0*mode65
deltaa(3)=1.0
nterms=2 !Maybe 2 is enough?
! Start the iteration
chisqr=0.
chisqr0=1.e6
do iter=1,100 !How many iters is enough?
do j=1,nterms
chisq1=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
fn=0.
delta=deltaa(j)
10 a(j)=a(j)+delta
chisq2=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
if(chisq2.eq.chisq1) go to 10
if(chisq2.gt.chisq1) then
delta=-delta !Reverse direction
a(j)=a(j)+delta
tmp=chisq1
chisq1=chisq2
chisq2=tmp
endif
20 fn=fn+1.0
chisq3=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
if(chisq3.lt.chisq2) then
chisq2=chisq3
go to 20
! Find minimum of parabola defined by last three points
delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
a(j)=a(j)-delta
deltaa(j)=deltaa(j)*fn/3.
! write(*,4000) iter,j,a(1:2),-chisq2
!4000 format(2i2,4f9.4)
enddo
chisqr=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
fdiff=chisqr/chisqr0-1.0
! write(*,4000) 0,0,a(1:2),-chisqr,fdiff
if(abs(fdiff).lt.0.0001) exit
chisqr0=chisqr
ccfbest=ccfmax * (1378.125/fsample)**2
dtbest=dtmax
return
end subroutine afc65b