      SUBROUTINE RESIDL(Q,TQ,NQ,KD,XNU,FF,FREMAX,NTERMS,
     *                  FREQ,AMP,PHI,SIGMAA,SIGMAP)
      DIMENSION Q(NQ),XNU(*),FREQ(*),AMP(*),PHI(*),SIGMAA(*)   
      DIMENSION FF(*),SIGMAP(*)   
      DOUBLE PRECISION TQ(NQ),DF,DFDF
      DOUBLE PRECISION XNUMAX,F,PERIOD
      X=NQ
      X=X*X
      AMAX=0.0    
      AMIN=1.0E30
      MODE=0
      IF(NTERMS.LT.0)THEN
        NTERMS=-NTERMS
        MODE=1
      END IF
C     
C     THIS ROUTINE CALCULATES THE RESIDUAL SUM OF SQUARES AGAINST
C     FREQUENCY FOR AN M-ORDER FOURIER SERIES FIT TO THE LIGHT CURVE
C     NTERMS IS THE NUMBER OF FOURIER COMPONENTS (>=1). THE QUANTITY
C     RETURNED IS 1/RESIDUAL VARIANCE
C     
      DF=XNU(2)-XNU(1)  
      DO K=1,KD     
        DO L=1,NTERMS
          FREQ(L)=L*XNU(K)
        END DO
        MD=0
        CALL LSTSQ(Q,TQ,1,NQ,NTERMS,FREQ,AMP,PHI,SIGMAA,
     *             SIGMAP,AFIT,MD,YMEAN,YSTDEV,RMEAN,RSTDEV)
        FF(K)=1.0/(RSTDEV*RSTDEV)
        IF(MODE.EQ.1)WRITE(11,4)K,XNU(K),FF(K)
    4   FORMAT(I10,1P3E20.8)
 
        IF(FF(K).GT.AMAX)THEN 
          IMAX=K  
          AMAX=FF(K)    
        END IF    
      END DO      

      IF(FREMAX.NE.-1.0)THEN
        PRINT 8   
    8   FORMAT(/14X,'DETERMINATION OF PRECISE MAXIMUM FREQUENCY', 
     *         ' FROM A FINER GRID AROUND THE ROUGH MAXIMUM'//    
     *         16X,'I',7X,'FREQUENCY',8X,'PERIOD',8X,'AMPLIT',5X, 
     *         '    PHASE'/12X,63('='))   
        AMAX=0.0  
        DFDF=DF/10.0
        XNUMAX=XNU(IMAX)
        DO 15 K=1,21    
          F=XNUMAX-(11.0D0-K)*DFDF
          DO L=1,NTERMS
            FREQ(L)=L*F
          END DO
          MD=0
          CALL LSTSQ(Q,TQ,1,NQ,NTERMS,FREQ,AMP,PHI,SIGMAA,
     *               SIGMAP,AFIT,MD,YMEAN,YSTDEV,RMEAN,RSTDEV)
          FFFF=1.0/(RSTDEV*RSTDEV)
          IF(FFFF.GT.AMAX)FREMAX=F  
          IF(FFFF.GT.AMAX)AMAX=FFFF 
          PERIOD=1000000.0D0
          IF(F.NE.0.0)PERIOD=1.0D0/F
          PRINT 12,K,F,PERIOD,FFFF,PHI(1)     
   12     FORMAT(13X,I4,1PE16.8,0PF14.6,1P2E14.5)     
   15   CONTINUE  
        PRINT 16  
   16   FORMAT(12X,63('='))   
        PRINT 18,FREMAX 
   18   FORMAT(/12X,'FREQUENCY OF MAXIMUM IS ',1PE14.6/)    
      END IF
      RETURN
      END   
