      SUBROUTINE ECLFIT(Q,TQ,NI,NF,NUM,A,B,C,SIGA,SIGB,SIGC,
     *                  IPRINT,MD,MODE)    
C     
C     THIS ROUTINE CALCULATES THE BEST FIT OF AN ECLIPSE MODEL 
C     TO A TOTAL ECLIPSE. MODE=1 FOR INGRESSES, =-1 FOR EGRESSES
C     
      dimension iwork2(33)
      DOUBLE PRECISION COV(33,33),COVINV(33,33),Y(33),DEL(33),CM(33,33) ! DS001
      DOUBLE PRECISION DERIV(33),UNIT(33,33),PARMS(33)
      DIMENSION JUNK(33),Q(NF) 
      DOUBLE PRECISION TQ(NF)
      DOUBLE PRECISION ARG,TI,QI,SI,CI,Q0,tmid
      REAL*8 FAC,FAC1,FAC2
      N=NF-NI+1   
      IDEBUG=0    
      IF(MD.LT.0)IDEBUG=1     
      IPRINT=1    
      SIG=0.0     
      PARMS(1)=A
      PARMS(2)=B
      PARMS(3)=C   
      DO 2 J=NUM+1,33     
    2 PARMS(J)=0.0
      ICONV=0     
      ITER=0
      SIG=0.0     
      DO 5 I=NI,NF
        TI=TQ(I)  
        QI=Q(I)   
        Q0=func(ti,num,parms,mode)
        SIG=SIG+(QI-Q0)**2    
    5 CONTINUE    
      SIG=SIG/(N-NUM)   
      IF(IPRINT.EQ.1)PRINT 8,NI,NF  
    8 FORMAT(12X,'NON-LINEAR LEAST SQUARES ECLIPSE FIT OF POINTS '/
     *       12x,I6,' TO ',I6)
      IF(IPRINT.EQ.1.AND.NUM.GT.0)PRINT 10,NUM,(PARMS(J),J=1,NUM)      
   10 FORMAT(/12X,I2,' PARAMETERS INCLUDED IN THE FIT WITH'/ 
     *       12X,'STARTING VALUES GIVEN BY :'// 
     *       12X,55('-')/(12X,3(1PE15.8,5X)))   
      PRINT 12    
   12 FORMAT(12X,55('-')////12X,'ITER',5X,'AVE SIG**2'/12X,109('-'))    
      PRINT 225,ITER,SIG,(PARMS(J),J=1,NUM)     
   15 ITER=ITER+1 
      DO 20 I=1,NUM     
        DEL(I)=0.0
        Y(I)=0.0  
        DO 20 J=1,NUM   
          CM(I,J)=0.0   
          COV(I,J)=0.0  
   20 CONTINUE    
C     
C     CALCULATE THE NORMAL EQUATION MATRIX AND Q VECTOR     
C     
      DO 50 I=NI,NF     
        TI=TQ(I)  
        QI=Q(I)   
        Q0=func(ti,num,parms,mode)
        DO 30 J=1,NUM 
          DERIV(J)=dfunc(ti,num,j,parms,mode)
   30   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)+(QI-Q0)*DERIV(J)
   40   CONTINUE  
   50 CONTINUE    
C     
C     INVERT THE NORMAL MATRIX 
C     
      DO 60 J=1,NUM     
      DO 60 K=1,NUM     
   60 COVINV(J,K)=COV(J,K)    
C     
C     INVERT THE NORMAL MATRIX USING DINV.FOR
C     
      call dinv(num,covinv,33,iwork2,ifail)
      if (ifail.ne.0) then
        print *,'--------------------------------------------------'
        print *,'MATRIX INVERSION PROBLEMS'
        print *,'ifail = ',ifail
        idebug=1
        ierror=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   
        print *,'--------------------------------------------------'
      END IF
C     
C     CALCULATE THE CORRECTIONS TO PARMS AND STORE IN DEL   
C     
      DO 140 I=1,NUM    
        DO 140 J=1,NUM  
          DEL(I)=DEL(I)+COVINV(I,J)*Y(J)  
  140 CONTINUE    
C     
C     APPLY CORRECTIONS 
C     
      DO 145 I=1,NUM    
  145 PARMS(I)=PARMS(I)+DEL(I)
C     
C     CALCULATE THE NEW SIGMA**2 OF THE RESIDUALS     
C     
      SIGNEW=0.0  
      DO 150 I=NI,NF    
        TI=TQ(I)  
        QI=Q(I)   
        Q0=func(ti,num,parms,mode)
        SIGNEW=SIGNEW+(QI-Q0)**2    
  150 CONTINUE    
      SIGNEW=SIGNEW/(N-NUM)   
C     
C     CALCULATE THE CORRELATION MATRIX FOR THE ERRORS 
C     
      DO 180 I=1,NUM    
        CM(I,I)=COVINV(I,I)*SIGNEW  
  180 CONTINUE    
      DO 200 I=1,NUM    
      DO 200 J=1,I-1    
        FAC1=COVINV(I,I)
        FAC2=COVINV(J,J)
        FAC=SQRT(FAC1*FAC2)
        IF(FAC.EQ.0.0)PRINT 190,FAC1,FAC2,FAC,I,J
  190   FORMAT(' FAC1= ',1PE15.6,' FAC2= ',1PE15.6,
     *         ' FAC= ',1PE15.6,' I= ',I5,' J= ',I5)
        CM(I,J)=COVINV(I,J)/FAC
  200 CONTINUE
C     
C     CHECK FOR CONVERGENCE   
C     
      ERROR=1.0E-7
      DELSIG=ABS(SIGNEW-SIG)/SIG    
      IF(DELSIG.LE.ERROR)ICONV=1    
      SIG=SIGNEW  
C     
C     PRINT RESULTS SECTION   
C     
      IF(IPRINT.EQ.1)THEN     
        IF(ICONV.EQ.1)PRINT 210,ITER
  210   FORMAT(/12X,'*** CONVERGENCE AFTER ',I3,' ITERATIONS ***'/)     
        IF(ITER.GT.9)PRINT 220
  220   FORMAT(/12X,'*** NO CONVERGENCE AFTER 10 ITERATIONS ***'/)
        PRINT 225,ITER,SIG,(PARMS(J),J=1,NUM)   
  225   FORMAT(12X,I4,1P7E15.6/(31X,6E15.6))    
        IF(ICONV.EQ.1.OR.ITER.GT.9)PRINT 230    
  230   FORMAT(12X,109('-')//)
      END IF
      IF(ICONV.EQ.1.OR.ITER.GT.9)THEN     
        IF(MD.EQ.1)THEN 
C     
C         PREWHITEN THE DATA VIA THE FINAL SOLUTION   
C     
          DO 235 I=NI,NF
            TI=TQ(I)    
            QI=Q(I)     
            Q0=func(ti,num,parms,mode)
            Q(I)=Q(I)-Q0
  235     CONTINUE
        END IF    
C     
C       CALCULATE ERRORS ON PARMS   
C     
        SIGA=SQRT(CM(1,1))  
        SIGB=SQRT(CM(2,2))  
        SIGC=SQRT(CM(3,3))  
        A=PARMS(1)  
        B=PARMS(2)  
        C=PARMS(3)  
        IF(IPRINT.EQ.1)PRINT 245    
  245   FORMAT(12X,'FINAL VALUES : '//12X,33X,'SIGMA'/12X,38('-'))
        IF(IPRINT.EQ.1.AND.NUM.GT.0)    
     *    PRINT 252,(PARMS(J),SQRT(CM(J,J)),J=1,NUM)
  252   FORMAT(12X,'   A',5X,1PE15.8,5X,E9.1/   
     *         12X,'   B',5X,E15.8,5X,E9.1/     
     *         12X,'   C',5X,E15.8,5X,E9.1)     
        PRINT 254,SQRT(SIG)   
  254   FORMAT(12X,38('-')//12X,'STANDARD DEVIATION OF RESIDUALS IS ',  
     *         1PE14.6/)
        IF(IPRINT.EQ.1)THEN   
          PRINT 255     
  255     FORMAT(12X,'CORRELATION COEFFICIENTS ARE :'/)     
          DO 270 I=2,NUM
            PRINT 260,(CM(I,J),J=1,I-1)   
  260       FORMAT(12X,15F8.3)
  270     CONTINUE
        END IF    
        RETURN
      END IF
C     
C     RETURN FOR THE NEXT ITERATION 
C     
      GO TO 15    
      END   


      REAL FUNCTION func(ti,num,parms,mode)
      DOUBLE PRECISION ti,tmid,parms(*)

      width=50.0
      tmid=parms(3)
      if(ti.le.tmid-width/2.0)then
        func=parms(1)
        return
      else if(ti.gt.tmid-width/2.0.and.
     *        ti.lt.tmid+width/2.0)then
        func=parms(1)-mode*parms(2)*
     *                (ti-tmid+width/2.0)/width
        return
      else if(ti.ge.tmid+width/2.0)then
        func=parms(1)-mode*parms(2)
        return
      end if
      end


      REAL FUNCTION dfunc(ti,num,j,parms,mode)
      DOUBLE PRECISION ti,tmid,parms(*)

      width=50.0
      tmid=parms(3)
      if(ti.le.tmid-width/2.0)then
        if(j.eq.1)dfunc=1.0
        if(j.eq.2)dfunc=0.0
        if(j.eq.3)dfunc=0.0
        return
      else if(ti.gt.tmid-width/2.0.and.
     *        ti.lt.tmid+width/2.0)then
        if(j.eq.1)dfunc=1.0
        if(j.eq.2)dfunc=mode*(tmid-ti-width/2.0)/width
        if(j.eq.3)dfunc=mode*parms(2)/width
        return
      else if(ti.ge.tmid+width/2.0)then
        if(j.eq.1)dfunc=1.0
        if(j.eq.2)dfunc=-1.0*mode
        if(j.eq.3)dfunc=0.0
        return
      end if
      end

