SUBROUTINE CGRIDLS
C This subroutine is adapted from GRIDLS by R.S. Simon, which in
C turn is loosely based on the subroutine GRIDLS in
C "Data Reduction and Error Analysis for the Physical Sciences" by
C P.R. Bevington (1969), pg. 212.
C By trial and error, it minimizes the A.F.
C
C J.T. Armstrong      4 Jan 1991
C
C-----------------------------------------------------------------------

      INCLUDE 'BINFIT.INC'

      REAL*8      CSQLIM, CSQ1, CSQ2, OCHISQ
      REAL*8      DELTA, SAVE, PREAF, PREPRM, PREDEL, FN
      INTEGER*4   WHERE
      LOGICAL*1 QUERY
C
      METHOD = 'BRUTE FORCE !!!'
      CSQLIM = 0.003
D     WRITE(OUTC,*) ' Entering CGRIDLS, ready to go to BINVIS'
      CALL CADJMOD
      CALL BINVIS         ! Calculate visibility due to model
D     WRITE(OUTC,*) ' Calculating parameter increments'
      DO IG = 1, NFILT
          DO IC = 1, NCOMP
              IF((IC.EQ.1) .AND. (IG.EQ.1)) THEN
                  DELTAP(IC,IG,1) = 0.05*SEP(IC,IG)
                  DELTAP(IC,IG,2) = 5.0 *DEGRAD
                  DELTAP(IC,IG,3) = 0.05*DIAM(IC,IG)
                  DELTAP(IC,IG,8) = 0.05*FLUXTOT(IC,IG)
              ELSE IF((IC.EQ.1) .AND. (IG.NE.1)) THEN
                  DELTAP(IC,IG,4) = 0.05*DRATIO(IC,IG)
                  DELTAP(IC,IG,8) = 0.05*FLUXTOT(IC,IG)
              ELSE IF((IC.NE.1) .AND. (IG.EQ.1)) THEN
                  DELTAP(IC,IG,1) = 0.05*SEP(IC,IG)
                  DELTAP(IC,IG,2) = 5.0 *DEGRAD
                  DELTAP(IC,IG,3) = 0.05*DIAM(IC,IG)
                  DELTAP(IC,IG,5) = 0.01            ! Mag. diff.
              ELSE IF((IC.NE.1) .AND. (IG.NE.1)) THEN
                  DELTAP(IC,IG,4) = 0.05*DRATIO(IC,IG)
                  DELTAP(IC,IG,6) = 0.01            ! Color
              END IF
          END DO
      END DO
D     WRITE(OUTC,*) ' Increments for parameters:'
D     CALL CWRITEMD(OUTC,' ',DELTAP,0,0)
D     WRITE(OUTC,*) ' Calculating reduced chi squared'
      CALL CAGRFAC
      DO ITER = 1,MAXIT
          JITER = ITER - 1
          OCHISQ = TCHISQ
C         -- Print current model parameters and fit
          WHERE = 1
          CALL CAFOUT(WHERE)
C
          DO IC=1,NCOMP
              DO IG=1,NFILT
                  DO K=1,8
                      IF (VARY(IC,IG,K)) THEN
                          DELTA = DELTAP(IC,IG,K)
                          CSQ1 = TCHISQ
                          PREAF  = TCHISQ
                          PREPRM = PARM(IC,IG,K)
                          PREDEL = DELTA
                          PARM(IC,IG,K) = PREPRM + DELTA
D                         WRITE(OUTC,*) ' Before CADJMOD'
D                         CALL CWRITEMD(OUTC,' ',PARM,IC,IG)
                          CALL CADJMOD
D                         WRITE(OUTC,*) ' After CADJMOD'
D                         CALL CWRITEMD(OUTC,' ',PARM,IC,IG)
                          CALL BINVIS
                          CALL CAGRFAC
D                         WRITE(OUTC,'(A,F7.3)') ' TCHISQ=',TCHISQ
                          CSQ2 = TCHISQ
C                         -- Reverse direction of search if AF increasing
                          IF (CSQ2 .GT. CSQ1) THEN
                              DELTA = -DELTA
                              PARM(IC,IG,K) = PARM(IC,IG,K) + DELTA
                              SAVE = CSQ1
                              CSQ1 = CSQ2
                              CSQ2 = SAVE
                          END IF
C                         -- Increment PARM(IC,IG,K) until AF increases
                          FN = 0.
   60                     FN = FN + 1.
                          PARM(IC,IG,K) = PARM(IC,IG,K) + DELTA
D                         WRITE(OUTC,*) ' Before CADJMOD'
D                         CALL CWRITEMD(OUTC,' ',PARM,IC,IG)
                          CALL CADJMOD
D                         WRITE(OUTC,*) ' After CADJMOD'
D                         CALL CWRITEMD(OUTC,' ',PARM,IC,IG)
                          CALL BINVIS
                          CALL CAGRFAC
D                         WRITE(OUTC,'(A,F7.3)') ' TCHISQ=',TCHISQ
                          IF (TCHISQ .LT. CSQ2) THEN
                              CSQ1 = CSQ2
                              CSQ2 = TCHISQ
                              GO TO 60
                          ELSE IF (TCHISQ .GT. CSQ2) THEN
C                             Find min. of parabola thru last three points
                              DELTA = DELTA 
     1                          * (1./(1.+(CSQ1-CSQ2)/(TCHISQ-CSQ2))+.5)
                              PARM(IC,IG,K) = PARM(IC,IG,K) - DELTA
                              DELTAP(IC,IG,K) = DELTAP(IC,IG,K) * FN/3.
                          END IF
D                         WRITE(OUTC,*) ' Before CADJMOD'
D                         CALL CWRITEMD(OUTC,' ',PARM,IC,IG)
                          CALL CADJMOD
D                         WRITE(OUTC,*) ' After CADJMOD'
D                         CALL CWRITEMD(OUTC,' ',PARM,IC,IG)
                          CALL BINVIS
C                         -- AF for final set of parameters
                          CALL CAGRFAC
D                         WRITE(OUTC,'(A,F7.3)') ' TCHISQ=',TCHISQ
                          IF (TCHISQ .GE. PREAF) THEN
C                             Model diverged, so reset model and decrease step size.
                              PARM(IC,IG,K) = PREPRM
                              DELTAP(IC,IG,K) = PREDEL/5.
                              CALL CADJMOD
                              CALL BINVIS
                              CALL CAGRFAC
                          END IF
                      END IF
                      IF (CCFLAG) GO TO 225
                  END DO
              END DO
          END DO
          JITER = JITER + 1
C
C If AF doesn't change significantly by 4th iteration, quit here.
C
          IF (ITER.GT.3 .AND. (OCHISQ-TCHISQ)/OCHISQ .LT. CSQLIM .AND.
     +        OCHISQ .GE. TCHISQ) THEN
              WRITE (OUTC,'(A,F6.3,A)',ERR=220)
     +            ' Convergence is slow (',TCHISQ,')'
              IF (QUERY(' Continue? (Y or N): ')) THEN
                  CSQLIM = CSQLIM/2.
              ELSE
                  GO TO 220
              END IF
          END IF
          OCHISQ = TCHISQ
      END DO
C
      CALL PUTOUT(' ')
      CALL PUTOUT(' Iteration halted after maximum no. of iterations.')
      GOTO 230
C
  220 CALL PUTOUT(' ')
      CALL PUTOUT(' Iteration halted;  model seems to have converged.')
      GOTO 230
C
  225 CALL PUTOUT(' ')
      CALL PUTOUT(' Iteration halted after an <interrupt>.')
      CALL CTRLCEN
      GOTO 230
C
  230 CALL CFIXMOD

      RETURN
      END