      SUBROUTINE SOMINC(Q,TQ,NQ,HELJD,TITLE,TIME,AFREQ,AMPLTD,PHASES,  
     *                  SIGAMP,SIGPH,NSHIFT,IND,E,O,C,OMINC,LL,PHCORR,
     *                  FILPER,FILRES)    
C     
C     ROUTINE TO CALCULATE OMINC FOR A GIVEN DATA SET AND EPHEMERIS     
C     
      DOUBLE PRECISION HELJD(1),TERM(10),TIMING(5000),CYCLE(5000) 
      DOUBLE PRECISION O(1000),C(1000),EVENT,TBEGIN,TEND
      DOUBLE PRECISION TQ(1),TZERO,TQQ(2000)
      DIMENSION TIME(1),AFREQ(2,1000),AMPLTD(1000),PHASES(1000)
      DIMENSION SIGAMP(1000),SIGPH(1000),IND(1000),E(1000),OMINC(1000) 
      DIMENSION Q(1),QQ(2000),FILPER(1),FILRES(1) 
      CHARACTER TITLE*80     
      NOPLOT=0    
      LL=LL+1     
      READ *,NI,N,ISHIFT,NSHIFT     
      READ *,NTERM,(TERM(I),I=1,NTERM)    
      IF(TERM(2).LT.0.0D0)THEN
        TERM(2)=-TERM(2)
        NOPLOT=1  
      END IF
      IF(NI.LE.0.OR.NI.GT.NQ)NI=1   
      IF(N.LE.0.OR.N.GT.NQ)N=NQ     
      IF(ISHIFT.EQ.0)ISHIFT=N 
      IF(NSHIFT.EQ.0)NSHIFT=(NQ/ISHIFT)-(N/ISHIFT)+1  
      IF(N.GT.2000)THEN 
        PRINT 5,N 
    5   FORMAT(' NUMBER OF FITTING POINTS IS ',I6,' AND EXCEEDS'/ 
     *         ' 2000 THE MAXIMUM ALLOWED.')    
        STOP
      END IF
      IF(NSHIFT.GT.1000)THEN   
        PRINT 6,NSHIFT  
    6   FORMAT(' NUMBER OF SECTIONS TO BE FITTED IS ',I5,/  
     *         ' WHICH EXCEEDS THE MAXIMUM ALLOWED WHICH IS 1000') 
        STOP
      END IF
      FREQ=1.0D0/(TERM(2)*86400.0D0)
      TBEGIN=HELJD(1)+TQ(NI)/86400.0D0    
      TEND=HELJD(1)+TQ(NI+N+(NSHIFT-1)*ISHIFT-1)/86400.0D0  
      PRINT 10,NSHIFT,N,FREQ,NI,ISHIFT,TBEGIN,TEND
   10 FORMAT(12X,'O-C 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'/12X,'SPANNING HJD ',
     *       0PF12.5,' TO ',F12.5)
      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)  
C     
C     CALCULATE THE MOMENTS OF MAXIMUM FROM THE EPHEMERIS.  USE UP
C     TO THE QUADRATIC TERMS TO GET AN INITIAL ESTIMATE OF THE    
C     CYCLE NUMBER.     
C     
      IF(NTERM.GT.2)THEN
        ICYC=0.5-(TERM(2)-DSQRT(TERM(2)**2-4.*TERM(3)*(TERM(1)-TBEGIN)))
     *           /(2.0*TERM(3))     
      ELSE  
        ICYC=0.5+(TBEGIN-TERM(1))/TERM(2) 
      END IF
      CYC=ICYC    
C     
C     NOW GENERATE 5000 MOMENTS     
C     
      DO 15 J=1,5000    
   12   EVENT=0.0D0     
        DO 14 K=1,NTERM 
          IF(K-1.NE.0)THEN
            EVENT=EVENT+TERM(K)*((CYC+J-1)**(K-1))
          ELSE
            EVENT=EVENT+TERM(K)
          END IF
   14   CONTINUE  
        IF(J.EQ.1.AND.EVENT.GT.TBEGIN) THEN     
          ICYC=ICYC-1   
          CYC=CYC-1.0   
          GO TO 12
        END IF    
        IF(EVENT.LE.TEND)THEN 
          NMOMNT=J
          CYCLE(J)=CYC+J-1    
          TIMING(J)=EVENT     
        END IF    
   15 CONTINUE    
      PRINT 25    
   25 FORMAT(/12X,3X,'NI',5X,'NF',6X,'AMPLITUDE',     
     *       10X,'SIGMA',10X,'PHASE',10X,'SIGMA',5X,'CYCLE',11X,'O',    
     *       11X,'C',9X,'O-C'/12X,118('-'))     
      LAST=1
C     
C     THIS IS THE MAIN LOOP   
C     
      DO 50 I=1,NSHIFT  
        NII=NI+(I-1)*ISHIFT   
        NFI=NII+N-1     
        IF(NFI.GT.NQ)GO TO 50 
        NFREQ=1   
        MD=0
C     
C       EXTRACT THE CURRENT SECTION FOR FITTING.  ADJUST ITS TIME 
C       SO THAT IT IT IS RELATIVE TO TZERO, THE MID POINT OF THE  
C       SECTION   
C     
        TZERO=(TQ(NII)+TQ(NFI))/2.0 
        DO 30 J=1,N     
          QQ(J)=Q(NII+J-1)    
          TQQ(J)=TQ(NII+J-1)-TZERO  
   30   CONTINUE  
        CALL LSTSQ(QQ,TQQ,1,N,NFREQ,FREQ,AMP,PHI,SIGMAA,SIGMAP,A0,MD,   
     *             YMEAN,YSTDEV,RMEAN,RSTDEV)   
C     
C       CALCULATE O(I)  
C     
        O(I)=(-(PHI+PHCORR)/(2.0*3.141592653*FREQ)+TZERO)/86400.0D0     
     *       +HELJD(1)  
C     
C       NOW FIND THE NEAREST C(I) FROM THE EPHEMERIS  
C     
        DO 35 J=LAST,NMOMNT   
          IF(DABS(O(I)-TIMING(J))/TERM(2).LT.0.5)GO TO 40   
   35   CONTINUE  
        PRINT 36,O(I),LAST,(TIMING(J),J=1,NMOMNT)     
   36   FORMAT(' CANT FIND THE NEAREST C: O(I)= ',F15.5,' LAST= ',I5/   
     *         ' WHILE ARRAY TIMING IS '//
     *         (1X,10F13.5))  
        STOP
C     
   40   LAST=J    
        C(I)=TIMING(J)  
        E(I)=CYCLE(J)   
        OMINC(I)=(O(I)-C(I))/TERM(2)
        TIME(I)=(O(I)-HELJD(1))*86400.0D0 
        PHI=PHI+PHCORR  
        PHASES(I)=PHI   
        AMP=AMP/RESMAX  
        AMPLTD(I)=AMP   
        SIGMAA=SIGMAA/RESMAX  
        SIGAMP(I)=SIGMAA
        SIGPH(I)=SIGMAP 
        IND(I)=LL 
        AFREQ(1,I)=FREQ 
        AFREQ(2,I)=1.0/FREQ   
        PRINT 45,NII,NFI,AMP,SIGMAA,PHI,SIGMAP,E(I),O(I),C(I),O(I)-C(I) 
   45   FORMAT(12X,I5,I7,1P4E15.6,0PF10.0,3F12.5)     
   50 CONTINUE    
      PRINT 60    
   60 FORMAT(12X,118('-'))    
      IF(NOPLOT.EQ.1)RETURN   
      XL=8.0
      IEB=-1
      CALL PLOTPH(XL,TIME,PHASES,SIGPH,E,OMINC,NSHIFT,TITLE,IEB) 
      CALL PLOTAM(XL,TIME,AMPLTD,SIGAMP,NSHIFT,TITLE,IEB)  
      RETURN
      END   
