      SUBROUTINE ECLFT1(Q,TQ,NI,NF,IPRINT,MD,TITLE)
C     
C     THIS ROUTINE CALCULATES THE BEST FIT OF AN ECLIPSE MODEL 
C     TO A TOTAL ECLIPSE. 
C     
      common /ecl/modee,param(30),nparam,ai
      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
      CHARACTER*80 TITLE
      REAL*8 FAC,FAC1,FAC2
      N=NF-NI+1   
      NUM=nparam
      IDEBUG=0    
      IF(MD.LT.0)IDEBUG=1     
      IPRINT=1    
      LOOP=0
      IF(AI.LT.0)THEN
C
C       LOOP OVER A RANGE OF INCLINATIONS
C
        AI=-AI   
        LOOP=1
        IPRINT=-1
        DAI=1.0
      END IF

    1 SIG=0.0     
c
c     initialize the parameters
c

c     For 13471-1258
      delt14=0.07244
      delt23=0.06407
      deltw=0.06832
      pi=3.141592653
      x1=(pi*delt14)**2
      x2=(pi*delt23)**2
      x3=(pi*deltw)**2
      r2=sqrt((cos(pi/180.0*ai))**2+x3)
      rho=(x1-x2)/(4*r2*r2)

      param(1)=0.7500 ! This is the eclipse depth
      param(2)=rho
      param(3)=r2
      param(4)=0.0    ! This is phi0
      param(5)=0.250  ! This is the mid-eclipse brightness
      param(6)=0.0900 ! The slope through eclipse
      param(7)=0.0    ! The quadratic term

c     For RR Cae
      delt14=0.03087
      delt23=0.02347
      deltw=0.02717
      pi=3.141592653
      x1=(pi*delt14)**2
      x2=(pi*delt23)**2
      x3=(pi*deltw)**2
      r2=sqrt((cos(pi/180.0*ai))**2+x3)
      rho=(x1-x2)/(4*r2*r2)

      param(1)=0.9000 ! This is the eclipse depth
      param(2)=rho
      param(3)=r2
      param(4)=0.0    ! This is phi0
      param(5)=0.100  ! This is the mid-eclipse brightness
      param(6)=0.0000 ! The slope through eclipse
      param(7)=0.0    ! The quadratic term

c     For SDSS-0155 
      delt14=0.061000
      delt23=0.059000
      deltw=0.060000
      pi=3.141592653
      x1=(pi*delt14)**2
      x2=(pi*delt23)**2
      x3=(pi*deltw)**2
      r2=sqrt((cos(pi/180.0*ai))**2+x3)
      rho=(x1-x2)/(4*r2*r2)
      print *,'rho= ',rho

      param(1)=0.9500 ! This is the eclipse depth
      param(2)=rho
      param(3)=r2
      param(4)=-0.001 ! This is phi0
      param(5)=0.050  ! This is the mid-eclipse brightness
      param(6)=0.0000 ! The slope through eclipse
      param(7)=0.0    ! The quadratic term


      PARMS(1)=param(1)
      PARMS(2)=param(2)
      PARMS(3)=param(3)
      PARMS(4)=param(4)
      PARMS(5)=param(5)
      PARMS(6)=param(6)
      PARMS(7)=param(7)
      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=func1(ti,num,parms,ai)
        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)))   
      IF(IPRINT.EQ.1)PRINT 12    
   12 FORMAT(12X,55('-')////12X,'ITER',5X,'AVE SIG**2'/12X,109('-'))    
      IF(IPRINT.EQ.1)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=func1(ti,num,parms,ai)
        DO 30 J=1,NUM 
          DERIV(J)=dfunc1(ti,num,j,parms,ai)
   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   
      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    
      IF(IDEBUG.EQ.1)THEN
        PRINT 142 
  142   FORMAT(///' THE CORRECTION VECTOR IS : '/) 
        PRINT 80,(DEL(J),J=1,NUM)     
        PRINT 80  
        PRINT *,'--------------------------------------------------'
      END IF
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=func1(ti,num,parms,ai)
        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=func1(ti,num,parms,ai)
            Q(I)=Q(I)-Q0
  235     CONTINUE
        END IF    
C     
C       CALCULATE ERRORS ON PARMS   
C     
        IF(IPRINT.NE.0)THEN   
          open(11,file='results',status='unknown',access='append')
          write(11,240)(parms(j),sqrt(cm(j,j)),j=1,num)
  240     format(2f6.3,2f6.3,2f7.4,1p2e8.1,0p2f6.3,2f6.3,2f6.3)
          write(11,241)sqrt(sig),title
  241     format(1x,f7.4,1x,a18)
          close(11)
C         
C         recalculate the eclipse widths
C
          rho=parms(2)
          r2=parms(3)
          cosi2=cos(pi/180.0*ai)**2
          x3=r2*r2-cosi2
          x1=(4.0*r2*r2*rho+2.0*r2*r2*rho*rho+2.0*x3)/2.0
          x2=(2.0*r2*r2*rho*rho+2.0*x3-4.0*r2*r2*rho)/2.0
          r1=rho*r2
          delt14=sqrt(x1)/pi
          delt23=sqrt(x2)/pi
          delt12=(delt14-delt23)/2.0
          deltw=sqrt(x3)/pi
          IF(IPRINT.EQ.1)PRINT 245,TITLE(1:60),ai
  245     FORMAT(12X,'FINAL VALUES : '/
     *           12X,A60/12X,'INCLINATION: ',F6.3,13X,'SIGMA'/
     *           12X,39('-'))
          IF(IPRINT.EQ.1)
     *    PRINT 252,(PARMS(J),SQRT(CM(J,J)),J=1,NUM)
  252     FORMAT(12X,'DEPTH',5X,1PE15.8,5X,E9.1/   
     *           12X,'RHO  ',5X,E15.8,5X,E9.1/     
     *           12X,'R2   ',5X,E15.8,5X,E9.1/     
     *           12X,'PHI0 ',5X,E15.8,5X,E9.1/     
     *           12X,'A0   ',5X,E15.8,5X,E9.1/     
     *           12X,'A1   ',5X,E15.8,5X,E9.1/     
     *           12X,'A2   ',5X,E15.8,5X,E9.1)     
          IF(IPRINT.EQ.1)
     *    PRINT 253,X1,X2,X3,R1,DELT14,DELT23,DELTW,DELT12
  253     FORMAT(12X,39('-')/
     *           12X,'X1    ',4X,F15.5/12X,'X2    ',4X,F15.5/
     *           12X,'X3    ',4X,F15.5/12X,'R1    ',4X,F15.5/
     *           12X,'DELT14',4X,F15.5/12X,'DELT23',4X,F15.5/
     *           12X,'DELTW ',4X,F15.5/12X,'DELT1/2',3X,F15.5)
          IF(IPRINT.EQ.1)
     *    PRINT 254,SQRT(SIG)   
  254     FORMAT(12X,39('-')//12X,'STANDARD DEVIATION OF RESIDUALS IS ',  
     *           1PE14.6/)
          IF(IPRINT.EQ.1)
     *    PRINT 255     
  255     FORMAT(12X,'CORRELATION COEFFICIENTS ARE :'/)     
          DO 270 I=2,NUM
            IF(IPRINT.EQ.1)PRINT 260,(CM(I,J),J=1,I-1)   
  260       FORMAT(12X,15F8.3)
  270     CONTINUE
          IF(IPRINT.EQ.-1)THEN
            IF(LOOP.EQ.1)PRINT 280
  280       FORMAT(12X,'  i  ','   Depth ','   Rho  ','   r2   ',
     *             '  Phi0  ','   a0   ','   a1   ','   r1   ',
     *             '  DelT14 ',' DelT23 ','  DelTW ',' Delt12 ')
            PRINT 290,AI,(PARMS(J),J=1,NUM),R1,
     *                DELT14,DELT23,DELTW,DELT12
  290       FORMAT(12X,F5.2,F8.3,F8.3,F8.4,F8.4,F8.3,F8.3,
     *             F8.5,F8.5,F8.5,F8.5,F8.5)
          END IF    
        END IF    
        do j=1,7
          param(j)=PARMS(J)
        end do
        IF(LOOP.GE.1)THEN
C
C         LOOP OVER THE NEXT INCLINATION ANGLE
C
          LOOP=LOOP+1
          AI=AI+DAI
          IF(AI.LE.90.0)GO TO 1
        END IF
        RETURN
      END IF
C     
C     RETURN FOR THE NEXT ITERATION 
C     
      GO TO 15    
      END   


      REAL FUNCTION func1(ti,num,parms,ai)
      DOUBLE PRECISION ti,parms(*)
      double precision pi,f,val,rho,r2,phi0,x,s
      double precision rho2,s2
      rho=parms(2)
      rho2=rho*rho
      r2=parms(3)
      phi0=parms(4)
c
c     ti is actually phase
c
      val=parms(5)
      if(num.ge.6)val=val+parms(6)*(ti-phi0)
      if(num.eq.7)val=val+parms(7)*(ti-phi0)**2
      f=parms(1)
      pi=3.141592653
      cosi=cos(ai*pi/180.0)
      x=1.-sqrt((2.*pi*(ti-phi0))**2+(cosi*cosi))/r2
      s=(2.*x-x*x-rho2)/(2.*(1.-x))
      s2=s*s
      if(x.ge.-rho.and.x.le.rho)
     *   f=f*(0.5
     *        -(rho2*asin(s/rho) +
     *          s*sqrt(rho2-s2) -
     *          sqrt((rho2-s2)*(1.-rho2+s2)) +
     *          asin(sqrt(rho2-s2)))/(pi*rho2)) 
      if(x.gt.rho)f=0.0
      val=val+f
      func1=val
      return
      end


      REAL FUNCTION dfunc1(ti,num,j,parms,ai)
      DOUBLE PRECISION ti,tmid,parms(*)
      double precision pi,f,val,rho,r2,phi0,x,s
      double precision rho2,s2,dfdrho,dfdr2,dfdp0
      rho=parms(2)
      rho2=rho*rho
      r2=parms(3)
      phi0=parms(4)
c
c     ti is actually phase
c
      pi=3.141592653
      cosi=cos(ai*pi/180.0)
      x=1.-sqrt((2.*pi*(ti-phi0))**2+(cosi*cosi))/r2
      s=(2.*x-x*x-rho2)/(2.*(1.-x))
      s2=s*s
      dfunc1=0.0
      if(j.eq.5)then
        dfunc1=1.0
        return
      end if
      if(j.eq.6)then
        dfunc1=ti-phi0
        return
      end if
      if(j.eq.7)then
        dfunc1=(ti-phi0)*(ti-phi0)
        return
      end if
      if(j.eq.1.or.j.eq.2)then
        f=parms(1)
        if(x.ge.-rho.and.x.le.rho)
     *     f=f*(0.5
     *          -(rho2*asin(s/rho) +
     *            s*sqrt(rho2-s2) -
     *            sqrt((rho2-s2)*(1.-rho2+s2)) +
     *            asin(sqrt(rho2-s2)))/(pi*rho2)) 
        if(x.gt.rho)f=0.0
      end if
      if(j.eq.1)then
        dfunc1=f
        return
      end if
      if(x.lt.-rho.or.x.gt.rho)then
        dfunc1=0.0
        return
      end if
      if(j.eq.2)then
        dfdrho=2.*sqrt(1-s2/rho2)/(pi*(1.-x))*
     *         (1.-s/sqrt(1.-rho2+s2)) -
     *         2.*(f-0.5)/rho -
     *         2./(pi*rho)*(asin(s/rho)+
     *                      sqrt((rho2-s2)/(1.-rho2+s2)))
        dfunc1=parms(1)*dfdrho
        return
      end if
      if(j.eq.3)then
        dfdr2=-2.*sqrt(1-s2/rho2)/(pi*rho)*
     *         (1.-s/sqrt(1.-rho2+s2))*
     *         (1.+s/(1-x))*
     *         (sqrt((2.*pi*(ti-phi0))**2+cosi*cosi)/(r2*r2))
        dfunc1=parms(1)*dfdr2
        return
      end if
      if(j.eq.4)then
        dfdp0=-2.*sqrt(1-s2/rho2)/(pi*rho)*
     *         (1.-s/sqrt(1.-rho2+s2))*
     *         (1.+s/(1-x))*
     *         (4*pi*pi*(ti-phi0)/(r2*r2*(1.-x))) -
     *         parms(6)
        dfunc1=parms(1)*dfdp0
        return
      end if
      print *,'dfunc1 error: j=',j
      stop
      end

