      SUBROUTINE EPHEM(EPOCH,EPER,ESEC,ETHIRD,CYCLE,OMINC,N,TITLE)     
C     
C     ROUTINE TO CALCULATE AN EPHEMERIS OF THE FORM   
C     EPOCH + E*P + E*E*ESEC + E*E*E*ETHIRD     
C
C     NUM IS THE NUMBER OF TERMS. IF NUM<0, A SINUSOIDAL COMPONENT
C     IS INCLUDED IN THE FIT
C     
      DOUBLE PRECISION CM(6,6),SIG,CALC(1000)
      dimension work2(6)                     
      DOUBLE PRECISION EPOCH,EPER,ESEC,ETHIRD,TIMING(1000),DCYC(1000)
      DOUBLE PRECISION COV(6,6),COVINV(6,6),UNIT(6,6)  
      DOUBLE PRECISION DERIV(6),Y(6),SOLN(6),RESID,E,C
      DOUBLE PRECISION EPLONG,PLONG,TCOS,PHICOS,TWOPI,ARG   
      DIMENSION JUNK(6),OMINC(1),CYCLE(1),CYCCYC(100),CALCAL(100) 
      CHARACTER FNAME*45,TITLE*80,FORMAT*80
      FORMAT=' '
      IDEBUG=0
      TWOPI=2.0D0*3.141592653D0     
      EPOCH=0.0D0 
      EPER=0.0D0  
      ESEC=0.0D0  
      ETHIRD=0.0D0
      READ 5,FNAME
    5 FORMAT(A)
      IF(FNAME.EQ.'FORMAT'.OR.FNAME.EQ.'format')THEN
C
C       READ IN THE DATA FILE INPUT FORMAT
C
        READ 5,FORMAT
C
C       NOW READ THE DATA FILE NAME (AGAIN!)
C
        READ 5,FNAME
      END IF
      READ *,NUM  
      IF(NUM.EQ.99)GO TO 400
      IPOLY=IABS(NUM)   
      IF(NUM.LT.0)THEN  
        NUM=-NUM+2
        READ *,PLONG,EPLONG
      END IF
      IF(NUM.LE.0.OR.NUM.GE.6)NUM=4 
      PRINT 10,FNAME,IPOLY-1    
   10 FORMAT(12X,'DATA TO BE READ FROM FILE ',A/
     *       12X,'POLYNOMIAL OF ORDER ',I2,' FITTED') 
      IF(NUM.GT.IPOLY)PRINT 11,PLONG,EPLONG
   11 FORMAT(12X,'A COSINE OF PERIOD ',F15.7,' D IS INCLUDED IN FIT'/
     *       12X,'WITH PHASING RELATIVE TO EPOCH ',F15.6)   
      PRINT 13    
   13 FORMAT(//)  
      PRINT 14    
   14 FORMAT(12X,'     CYCLE                TIME           PREDICTED',  
     *       5X,'       O-C'/ 
     *       12X,'----------        ------------        ------------',  
     *       5X,'----------') 
      CLOSE(11)   
      OPEN(11,FILE=FNAME,STATUS='OLD')     
      READ(11,*)N 
      IF(FORMAT.EQ.' ')THEN
        READ(11,*)(TIMING(J),DCYC(J),J=1,N)
      ELSE
        READ(11,FORMAT)(TIMING(J),DCYC(J),J=1,N)
      END IF
      DO J=1,N
        DCYC(J)=NINT(DCYC(J))
      END DO

C  20 FORMAT(5X,F13.7,85X,F10.0)   OLD FORMAT FOR READING DQ Her TIMINGS
      REWIND 11   
      CLOSE(11)   

C     
C     CALCULATE THE NORMAL EQUATION MATRIX
C     
      DO 30 J=1,NUM     
        Y(J)=0.0  
        DO 30 K=1,NUM   
          COV(J,K)=0.0  
          CM(J,K)=0.0   
   30 CONTINUE    
      DO 50 I=1,N 
        E=DCYC(I)
        DERIV(1)=1.0    
        DO 32 J=2,NUM   
          IF(J.LE.IPOLY)THEN  
            DERIV(J)=DERIV(J-1)*E   
          ELSE    
            ARG=DMOD(TWOPI*(TIMING(I)-EPLONG)/PLONG,TWOPI)    
            IF(J.EQ.IPOLY+1)DERIV(J)=DCOS(ARG)  
            IF(J.EQ.IPOLY+2)DERIV(J)=DSIN(ARG)  
          END IF  
   32   CONTINUE  
        DO 40 J=1,NUM   
          DO 35 K=1,NUM 
            COV(J,K)=COV(J,K)+DERIV(J)*DERIV(K) 
   35     CONTINUE
          Y(J)=Y(J)+DERIV(J)*TIMING(I)    
   40   CONTINUE  
   50 CONTINUE    
C     
C     INVERT THE NORMAL MATRIX
C     
      DO 60 J=1,NUM     
        DO 60 K=1,NUM   
          COVINV(J,K)=COV(J,K)
   60 CONTINUE    
C     
C     INVERT THE NORMAL MATRIX USING DINV.FOR
C     
      call dinv(num,covinv,6,work2,ifail)
      if (ifail.ne.0) then
        print *,'MATRIX INVERSION FAILURE'
        print *,'ifail = ',ifail
        ierror=1
        idebug=1
      endif
C     
C     FROM HERE DOWN TO ABOUT LABEL 135 ARE DEBUG CALCULATIONS & PRINTS 
C     
   65 IF(IDEBUG.EQ.1)THEN     
        PRINT 70  
   70   FORMAT(//' COVARIANCE MATRIX IS : '/)   
        DO 90 I=1,NUM   
          PRINT 80,(COV(I,J),J=1,NUM)     
   80     FORMAT(1X,1P10E13.4)
          PRINT 80
   90   CONTINUE  
        PRINT 95  
   95   FORMAT(///' INVERSE IS : '/)
        DO 100 I=1,NUM  
          PRINT 80,(COVINV(I,J),J=1,NUM)  
          PRINT 80
  100   CONTINUE  
        PRINT 105 
  105   FORMAT(///' THE Y VECTOR IS : '/) 
        PRINT 80,(Y(J),J=1,NUM)     
        PRINT 80  
C     
C       TEST FOR GOOD MATRIX INVERSION BY PRODUCING A UNIT MATRIX 
C       AND PRINTING IT AND THEN CALCULATE THE CONTRIBUTION TO THE
C       SOLUTIONS BY CALCULATING A MATRIX GIVEN BY COVINV(I,J)*Y(J)     
C     
        PRINT 120 
        DO 125 I=1,NUM  
          DO 115 K=1,NUM
            UNIT(I,K)=0.0     
            DO 110 J=1,NUM    
              UNIT(I,K)=UNIT(I,K)+COVINV(I,J)*COV(J,K)
  110       CONTINUE    
  115     CONTINUE
  120     FORMAT(///' PRODUCT OF COVARIANCE AND ITS INVERSE IS : '/)    
          PRINT 80,(UNIT(I,K),K=1,NUM)    
          PRINT 80
  125   CONTINUE  
        PRINT 128 
  128   FORMAT(///' PRODUCT OF COVARIANCE INVERSE AND Y VECTOR :'/)     
        DO 135 I=1,NUM  
          DO 130 J=1,NUM
            UNIT(I,J)=COVINV(I,J)*Y(J)    
  130     CONTINUE
          PRINT 80,(UNIT(I,J),J=1,NUM)    
          PRINT 80
  135   CONTINUE  
        IF(IERROR.EQ.1)STOP   
      END IF
C     
C     CALCULATE THE SOLUTION  
C     
      DO 140 I=1,NUM    
        SOLN(I)=0.0     
        DO 140 J=1,NUM  
          SOLN(I)=SOLN(I)+COVINV(I,J)*Y(J)
  140 CONTINUE    
      SIG=0.0     
C     
C     CALCULATE THE NEW SIGMA OF THE RESIDUALS  
C     
      DO 150 I=1,N
        RESID=TIMING(I) 
        C=0.0     
        IF(IPOLY.LT.NUM)ARG=DMOD(TWOPI*(TIMING(I)-EPLONG)/PLONG,TWOPI)    
        DO 145 J=1,NUM  
          IF(J.LE.IPOLY)THEN  
            IF(DCYC(I).EQ.0.0.AND.J.EQ.1)RESID=RESID-SOLN(1)     
            IF(DCYC(I).NE.0.0)RESID=RESID-SOLN(J)*(DCYC(I)**(J-1))    
            IF(J.LT.NUM)OMINC(I)=RESID/SOLN(2)  
            IF(DCYC(I).EQ.0.0)C=SOLN(1)  
            IF(DCYC(I).NE.0.0)C=C+SOLN(J)*(DCYC(I)**(J-1))
          ELSE    
            IF(J.EQ.IPOLY+1)THEN    
              C=C+SOLN(J)*DCOS(ARG) 
              RESID=RESID-SOLN(J)*DCOS(ARG)     
            END IF
            IF(J.EQ.IPOLY+2)THEN    
              C=C+SOLN(J)*DSIN(ARG) 
              RESID=RESID-SOLN(J)*DSIN(ARG)     
            END IF
          END IF  
  145   CONTINUE  
        OMINC(I)=RESID/SOLN(2)
        PRINT 146,DCYC(I),TIMING(I),C,RESID    
  146   FORMAT(12X,F10.0,2(5X,F15.6),5X,F10.6)  
        SIG=SIG+RESID**2
  150 CONTINUE    
      SIG=SIG/(N-NUM)   
      SIGMA=SQRT(SIG)   
C     
C     CALCULATE THE CORRELATION MATRIX FOR THE ERRORS 
C     
      DO 200 I=1,NUM    
        DO 200 J=1,NUM  
  200     CM(I,J)=COVINV(I,J)*SIG   
      PRINT 250,SIGMA,(SOLN(I),SQRT(CM(I,I)),I=1,IPOLY)     
  250 FORMAT(/12X,'SIGMA OF RESIDUALS ',1PE9.1//
     *        12X,33X,'SIGMA'/12X,38('-')/
     *        12X,'EPOCH',4X,0PF15.6,5X,F9.6/   
     *        12X,'PERIOD',1X,1PE17.10,2X,E12.1/ 
     *        12X,'SECOND',3X,1PE15.3,5X,E9.1/  
     *        12X,'THIRD',4X,1PE14.3,5X,E9.1)   
      IF(NUM.GT.IPOLY)THEN    
        TCOS=SQRT(SOLN(IPOLY+1)**2+SOLN(IPOLY+2)**2) 
        PHICOS=DATAN2(-SOLN(IPOLY+2),SOLN(IPOLY+1))   
        SIGTCS=SQRT(SOLN(IPOLY+1)**2*CM(IPOLY+1,IPOLY+1) + 
     *               SOLN(IPOLY+2)**2*CM(IPOLY+2,IPOLY+2))/TCOS   
        SIGPHC=SQRT(SOLN(IPOLY+2)**2*CM(IPOLY+1,IPOLY+1) + 
     *               SOLN(IPOLY+1)**2*CM(IPOLY+2,IPOLY+2))/(TCOS**2)    
        PRINT 255,(SOLN(I),SQRT(CM(I,I)),I=IPOLY+1,NUM),    
     *             TCOS,SIGTCS,PHICOS,SIGPHC    
  255   FORMAT(12X,'COS TRM',2X,1PE15.3,5X,E9.1/
     *         12X,'SIN TRM',2X,1PE15.3,5X,E9.1/
     *         12X,'T SHIFT',2X,1PE15.3,5X,E9.1/
     *         12X,'PHASE  ',2X,0PF15.5,5X,F9.5)
      END IF
      PRINT 260   
  260 FORMAT(12X,53('-'))     
      EPOCH=SOLN(1)     
      EPER=SOLN(2)
      ESEC=SOLN(3)
      ETHIRD=SOLN(4)    
C
C     WRITE OUT THE RESIDUALS
C
      OPEN(12,FILE='resid',STATUS='UNKNOWN')
      NUM=IPOLY
      DO 300 I=1,N
C
C       SPECIAL JOB FOR AEAQR: WRITE OUT THE RESIDUALS WITHOUT
C       THE ORBITAL VARIATION INCLUDED
C
C       RESID=TIMING(I) 
C       DO 290 J=1,NUM  
C         IF(DCYC(I).EQ.0.0.AND.J.EQ.1)RESID=RESID-SOLN(1)     
C         IF(DCYC(I).NE.0.0)RESID=RESID-SOLN(J)*(DCYC(I)**(J-1))    
C 290   CONTINUE  
C       RESID=RESID/SOLN(2)  
C       WRITE(12,295)TIMING(I),RESID

        WRITE(12,295)TIMING(I),OMINC(I)
  295   FORMAT(F15.6,F10.4)
        CYCLE(I)=DCYC(I)
        OMINC(I)=TIMING(I)-EPOCH-DCYC(I)*EPER
        OMINC(I)=OMINC(I)/EPER               
  300 CONTINUE    
      CLOSE(12)
      DO I=1,90
        CYCCYC(I)=NINT(CYCLE(1)+(CYCLE(N)-CYCLE(1))*(I-3.0)/85.0)
        CALCAL(I)=0.0
        IF(IPOLY.EQ.3)CALCAL(I)=CALCAL(I)+(CYCCYC(I)**2)*ESEC
        IF(IPOLY.LT.NUM)THEN
          ARG=EPOCH+EPER*CYCCYC(I)
          IF(IPOLY.EQ.3)ARG=ARG+(CYCCYC(I)**2)*ESEC
          ARG=DMOD(TWOPI*(ARG-EPLONG)/PLONG,TWOPI)    
          CALCAL(I)=CALCAL(I)+DTLONG*DCOS(ARG)
          CALCAL(I)=CALCAL(I)+SOLN(IPOLY+1)*DCOS(ARG) 
     *                       +SOLN(IPOLY+2)*DSIN(ARG) 
        END IF
        CALCAL(I)=CALCAL(I)/EPER                     
      END DO
      NN=90
      CALL PLOTEPH(CYCLE,OMINC,N,CYCCYC,CALCAL,NN,TITLE)
      RETURN

  400 CONTINUE
C
C     SECTION TO CALCULATE THE EPHEMERIS USING AN ADAPTED VERSION
C     OF THE NONLIN ROUTINE (CALLED EPNLLSQ)
C
C     SET THINGS UP
C
      READ *,IPOLY
      IF(IPOLY.EQ.2)READ *,EPOCH,EPER
      IF(IPOLY.EQ.3)READ *,EPOCH,EPER,ESEC
      READ *,PLONG,EPLONG,DTLONG
      NUM=IPOLY+3
      NFREQ=1
      PRINT 10,FNAME,IPOLY-1    
      PRINT 411,PLONG,EPLONG,DTLONG
  411 FORMAT(12X,'A COSINE OF PERIOD ',F15.7,' D IS INCLUDED IN FIT'/
     *       12X,'WITH PHASING RELATIVE TO EPOCH ',F15.6/
     *       12X,'AND AMPLITUDE                  ',F15.6)   
      PRINT 13    
      CLOSE(11)   
      OPEN(11,FILE=FNAME,STATUS='OLD')     
      READ(11,*)N 
      IF(FORMAT.EQ.' ')THEN
        READ(11,*)(TIMING(J),DCYC(J),J=1,N)
      ELSE
        READ(11,FORMAT)(TIMING(J),DCYC(J),J=1,N)
      END IF
      DO J=1,N
        DCYC(J)=NINT(DCYC(J))
      END DO

      REWIND 11   
      CLOSE(11)   
      
      MD=0
      IPRINT=1
      FREQ=EPER/PLONG
      AMP=DTLONG
      PHI=TWOPI*(EPOCH-EPLONG)/PLONG
      CALL EPNLLSQ(TIMING,DCYC,1,N,CALC,IPOLY,
     *             EPOCH,EPER,ESEC,SEPOCH,SEPER,SESEC,
     *             NFREQ,FREQ,AMP,PHI,SIGMAF,SIGMAA,SIGMAP,
     *             IPRINT,MD)    
      PLONG=EPER/FREQ
      EPLONG=EPOCH-PLONG*PHI/TWOPI
      PRINT 420,EPOCH,SEPOCH,EPER,SEPER,ESEC,SESEC,
     *          PLONG,PLONG*SIGMAF/FREQ,AMP,SIGMAA,EPLONG
  420 FORMAT(12X,'FINAL VALUES : '//12X,35X,'SIGMA'/12X,40('-')/
     *       12X,'EPOCH ',5X,F15.6,5X,F9.6/     
     *       12X,'PERIOD',5X,F15.9,F14.9/       
     *       12X,'SECOND',5X,1PE15.3,E14.3/       
     *       12X,'PLONG ',5X,0PF15.1,5X,F9.1/   
     *       12X,'DT    ',5X,1PE15.8,5X,E9.1/   
     *       12X,'EPLONG',5X,0PF15.1)              
      PRINT 13
      PRINT 14
      DO I=1,N
        PRINT 146,DCYC(I),TIMING(I),CALC(I),TIMING(I)-CALC(I)
      END DO
      OPEN(12,FILE='resid',STATUS='UNKNOWN')
      DO I=1,N
        WRITE(12,295)TIMING(I),(TIMING(I)-CALC(I))/EPER
        CYCLE(I)=DCYC(I)
        OMINC(I)=TIMING(I)-EPOCH-DCYC(I)*EPER
        OMINC(I)=OMINC(I)/EPER
      END DO      
      CLOSE(12)
      DO I=1,90
        CYCCYC(I)=NINT(CYCLE(1)+(CYCLE(N)-CYCLE(1))*(I-5.0)/80.0)
        CALCAL(I)=0.0
        IF(IPOLY.EQ.3)CALCAL(I)=CALCAL(I)+(CYCCYC(I)**2)*ESEC
        ARG=DMOD(TWOPI*FREQ*CYCCYC(I)+PHI,TWOPI)
        CALCAL(I)=CALCAL(I)+AMP*DCOS(ARG)
        CALCAL(I)=CALCAL(I)/EPER          
      END DO
      NN=90
      CALL PLOTEPH(CYCLE,OMINC,N,CYCCYC,CALCAL,NN,TITLE)
      RETURN
      END   
