      SUBROUTINE SPHASE(Q,TQ,NQ,HELJD,TITLE,TIME,AFREQ,AMPLTD,PHASES,  
     *                  SIGAMP,SIGPH,IND,NSHIFT,LL,PHCORR,FILPER,FILRES)
C     
C     THIS ROUTINE TAKES THE NQ DATA IN Q AT TIMES TQ AND FITS TO 
C     SUCCESSIVE SECTIONS (OF LENGTH ISHIFT) A SINUSOID OF FREQUENCY    
C     FREQ.  THE RESULTING AMPLITUDES, PHASES AND ASSOCIATED ERRORS ARE 
C     CAPTURED IN THE ARRAYS IN THE PARAMETER LIST    
C     
      DIMENSION Q(*),TIME(*),AFREQ(2,1),AMPLTD(*),PHASES(*) 
      DIMENSION SIGAMP(*),SIGPH(*),IND(*),FILPER(*),FILRES(*)     
C     THESE ARE JUST DUMMY ARRAYS TO KEEP ROUTINE PLOTPH HAPPY    
      DIMENSION EE(1),OO(1)   
      DOUBLE PRECISION HELJD(*),HJDZER,TQ(*)
      CHARACTER TITLE*80
      SAMPLT=TQ(2)-TQ(1)
      READ *,NI,N,ISHIFT,NSHIFT,FREQ,XL   
      IEB=1 
      IF (XL.LT.0.0) THEN     
        IEB=0     
        XL=-XL    
      ENDIF 
      IDUMP=0     
      IF(XL.EQ.0.0)XL=8.0     
      IF(N.EQ.0)N=NQ    
      IPLTPH=0    
      IF(FREQ.LT.0.0)THEN     
        FREQ=-FREQ
        IPLTPH=1  
      END IF
      IF(NI.LE.0)NI=1   
      IF(N.LE.0.OR.N.GT.NQ)THEN     
        PRINT *,'BAD VALUE FOR NUMBER OF POINTS TO BE FITTED'     
        RETURN    
      END IF
      IF(ISHIFT.EQ.0)ISHIFT=N 
      IF(NSHIFT.EQ.0)NSHIFT=(NQ/ISHIFT)-(N/ISHIFT)+1  
      IF(FREQ.LE.0.0)THEN     
        PRINT *,'FITTING FREQUENCY ZERO OR NEGATIVE'  
        STOP
      END IF
      PRINT 10,NSHIFT,N,FREQ,NI,ISHIFT    
   10 FORMAT(12X,'PHASE DETERMINATIONS FROM LEAST SQUARES FITTING',     
     *       ' OF ',I5,' BLOCKS OF ',I5,' POINTS WITH FREQUENCY ',
     *       1PE15.5/12X,'STARTING FROM POINT ',I5,' AND SHIFTING',     
     *       ' ALONG ',I5,' POINTS AT A TIME')  
      RESMAX=1.0  
      CALL FILDET(FILPER,FILRES,FREQ,RESMAX)    
      IF(RESMAX.NE.1.0.OR.PHCORR.NE.0.0)PRINT 11,RESMAX,PHCORR    
   11 FORMAT(/12X,'THE AMPLITUDES HAVE BEEN SCALED UP BY THE FILTER'/   
     *        12X,'RESPONSE OF 1 / ',F6.3,' FOR THIS FREQUENCY'/  
     *        12X,'AND THE PHASES HAVE BEEN CORRECTED BY ',F8.4)  
      PRINT 25    
   25 FORMAT(/12X,'RESULTS : ',3X,'NI',5X,'NF',5X,'TIME OF NI',6X,
     *       'AMPLITUDE',10X,'SIGMA',10X,'PHASE',10X,'SIGMA',5X,  
     *       'FIT FACTOR'/12X,9('='),1X,102('-'))     
      NPHASE=0    
      LL=LL+1     
      DO 50 I=1,NSHIFT  
        NII=NI+(I-1)*ISHIFT   
        NFI=NII+N-1     
   28   IF(NFI.GT.NQ)GO TO 50 
        IF(TQ(NFI)-TQ(NII).GT.5.*N*SAMPLT)THEN  
          NII=NII+1     
          NFI=NFI+1     
          GO TO 28
C         THIS IS DESIGNED TO JUMP THE DAILY GAPS SO THAT NO BLOCK
C         STRADDLES A GAP     
        END IF    
        NFREQ=1   
        MD=0
        CALL LSTSQ(Q,TQ,NII,NFI,NFREQ,FREQ,AMP,PHI,SIGMAA,SIGMAP,A0,MD, 
     *             YMEAN,YSTDEV,RMEAN,RSTDEV)   
        AMP=AMP/RESMAX  
        SIGMAA=SIGMAA/RESMAX  
        PHI=PHI+PHCORR  
        FFACT=(YSTDEV-RSTDEV)*100.0/YSTDEV
        HJDZER=TQ(NII)/86400.0D0+HELJD(1) 
       PRINT 30,NII,NFI,HJDZER,AMP,SIGMAA,PHI,SIGMAP,FFACT  
   30   FORMAT(22X,I5,2X,I5,0PF15.5,1P4E15.6,5X,0PF10.1)    
        NPHASE=NPHASE+1 
        TIME(NPHASE)=(TQ(NII)+TQ(NFI))/2.0
        PHASES(NPHASE)=PHI    
        AMPLTD(NPHASE)=AMP    
        SIGAMP(NPHASE)=SIGMAA 
        SIGPH(NPHASE)=SIGMAP  
        IND(NPHASE)=LL  
        AFREQ(1,NPHASE)=FREQ  
        AFREQ(2,NPHASE)=1.0/FREQ    
   50 CONTINUE    
      NSHIFT=NPHASE     
      PRINT 60    
   60 FORMAT(22X,102('-'))    
      IF(IPLTPH.EQ.1)THEN     
        CALL PLOTPH(XL,TIME,PHASES,SIGPH,EE,OO,NSHIFT,TITLE,IEB) 
        CALL PLOTAM(XL,TIME,AMPLTD,SIGAMP,NSHIFT,TITLE,IEB)
      END IF
      RETURN
      END   
