real function fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)

      parameter (NMAX=60*96000)          !Samples per 60 s
      complex cx(npts),cy(npts)
      real a(5)
      complex w,wstep,za,zb,z
      real ss(2600)
      complex csx(0:NMAX/64),csy(0:NMAX/64)
      data twopi/6.283185307/a1,a2,a3/99.,99.,99./
      save

      baud=11025.0/4096.0
      if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then
         a1=a(1)
         a2=a(2)
         a3=a(3)

C  Mix and integrate the complex X and Y signals
         csx(0)=0.
         csy(0)=0.
         w=1.0
         x0=0.5*(npts+1)
         s=2.0/npts
         do i=1,npts
            x=s*(i-x0)
            if(mod(i,100).eq.1) then
               p2=1.5*x*x - 0.5
!               p3=2.5*(x**3) - 1.5*x
!               p4=4.375*(x**4) - 3.75*(x**2) + 0.375
               dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample)
               wstep=cmplx(cos(dphi),sin(dphi))
            endif
            w=w*wstep
            csx(i)=csx(i-1) + w*cx(i)
            csy(i)=csy(i-1) + w*cy(i)
         enddo
      endif

C  Compute 1/2-symbol powers at 1/16-symbol steps.
      fac=1.e-4
      pol=a(4)/57.2957795
      aa=cos(pol)
      bb=sin(pol)
      nsps=nint(fsample/baud)                  !Samples per symbol
      nsph=nsps/2                              !Samples per half-symbol

      ndiv=16                                  !Output ss() steps per symbol
      nout=ndiv*npts/nsps
      dtstep=1.0/(ndiv*baud)                   !Time per output step

      do i=1,nout
         j=i*nsps/ndiv
         k=j-nsph
         ss(i)=0.
         if(k.ge.1) then
            za=csx(j)-csx(k)
            zb=csy(j)-csy(k)
            z=aa*za + bb*zb
            ss(i)=fac*(real(z)**2 + aimag(z)**2)
         endif
      enddo

      ccfmax=0.
      call ccf2(ss,nout,nflip,ccf,lagpk)
      if(ccf.gt.ccfmax) then
         ccfmax=ccf
         dtmax=lagpk*dtstep
      endif
      fchisq=-ccfmax

      return
      end