SUBROUTINE CGRADLS C This subroutine is adapted from GRADLS by R.S. Simon, which in C turn is loosely based on the subroutine GRADLS in C "Data Reduction and Error Analysis for the Physical Sciences" by C P. R. Bevington (1969), pg. 221. C By moving along the gradient of the "agreement factor", minimizes C the A.F. C C J.T. Armstrong 29 Dec 1990 C C----------------------------------------------------------------------- INCLUDE 'BINFIT.INC' REAL*8 CSQLIM, CSQ1, CSQ2, OCHISQ, GRAD(MXMOD,MXFILT,8) REAL*8 SUM, DELTA, SAVE INTEGER*4 NBACK, WHERE, IGG LOGICAL*1 QUERY METHOD = 'Gradient search' CSQLIM = 0.01 D WRITE(OUTC,*) ' Entering CGRADLS, 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 ! Find reduced CHISQ DO 210 ITER = 1, MAXIT JITER = ITER - 1 OCHISQ = TCHISQ CSQ1 = TCHISQ SUM = 0. WHERE = 1 CALL CAFOUT(WHERE) ! Print model parameters and fit to terminal C Find components of gradient: DO IC=1,NCOMP DO IG=1,NFILT DO K=1,8 IF (VARY(IC,IG,K)) THEN D IF(ITER.EQ.1) THEN D WRITE(OUTC,'(A,I1,A,I4,A,I2,A)') D 1 ' VARY(',K,',',IC,',',IG,') is TRUE' D END IF DELTA = 0.1 * DELTAP(IC,IG,K) PARM(IC,IG,K) = PARM(IC,IG,K) + DELTA IF( (K.EQ.1 .OR. K.EQ.2 .OR. K.EQ.3) 1 .AND. IG.EQ.1 ) THEN DO IGG=2,NFILT PARM(IC,IGG,K) = PARM(IC,IG,K) END DO END IF CALL CADJMOD CALL BINVIS CALL CAGRFAC PARM(IC,IG,K) = PARM(IC,IG,K) - DELTA IF( (K.EQ.1 .OR. K.EQ.2 .OR. K.EQ.3 ) 1 .AND. IG.EQ.1 ) THEN DO IGG=2,NFILT PARM(IC,IGG,K) = PARM(IC,IG,K) END DO END IF GRAD(IC,IG,K) = CSQ1 - TCHISQ SUM = SUM + GRAD(IC,IG,K)*GRAD(IC,IG,K) END IF END DO END DO END DO D WRITE(OUTC,*) ' Components of the gradient:' D CALL CWRITEMD(OUTC,' ',GRAD,0,0) IF (SQRT(SUM) .EQ. 0.) GO TO 220 C Normalize the gradient components and change each variable parameter. DO IC = 1,NCOMP DO IG = 1,NFILT DO K = 1,8 IF (VARY(IC,IG,K)) THEN GRAD(IC,IG,K) = DELTAP(IC,IG,K) 1 *GRAD(IC,IG,K)/SQRT(SUM) C Move to a new point. PARM(IC,IG,K) = PARM(IC,IG,K) + GRAD(IC,IG,K) END IF END DO END DO END DO D WRITE(OUTC,*) ' Components of the normalized gradient:' D CALL CWRITEMD(OUTC,' ',GRAD,0,0) CALL CADJMOD CALL BINVIS CALL CAGRFAC C Make sure AF decreases. IF (CSQ1 .LE. TCHISQ) THEN DO NBACK = 1, 15 C Moved too far, so decrease step size. C Moved too far, so backup in steps of 1/2, until AF has decreased C to less than its original value, then rearrange parameters so that C it looks like step size was correct. CSQ2 = TCHISQ DO IC=1,NCOMP DO IG=1,NFILT DO K=1,8 IF (VARY(IC,IG,K)) THEN GRAD(IC,IG,K) = GRAD(IC,IG,K) / 2. PARM(IC,IG,K) = PARM(IC,IG,K)-GRAD(IC,IG,K) DELTAP(IC,IG,K) = DELTAP(IC,IG,K) / 2. END IF END DO END DO END DO CALL CADJMOD CALL BINVIS CALL CAGRFAC IF (TCHISQ .LT. CSQ1) THEN SAVE = CSQ2 CSQ2 = TCHISQ TCHISQ = SAVE DO IC = 1,NCOMP DO IG = 1,NFILT DO K = 1, 8 IF (VARY(IC,IG,K)) THEN PARM(IC,IG,K) 1 = PARM(IC,IG,K)+GRAD(IC,IG,K) END IF END DO END DO END DO GO TO 80 END IF END DO GO TO 200 END IF C Increment parameters until AF starts to increase. CSQ2 = TCHISQ 70 DO IC = 1, NCOMP DO IG = 1, NFILT DO K = 1,8 IF (VARY(IC,IG,K)) THEN PARM(IC,IG,K) = PARM(IC,IG,K) + GRAD(IC,IG,K) END IF END DO END DO END DO CALL CADJMOD CALL BINVIS CALL CAGRFAC IF (TCHISQ .LT. CSQ2) THEN CSQ1 = CSQ2 CSQ2 = TCHISQ GO TO 70 END IF C Find minimum of parabola defined by last three points 80 IF (TCHISQ .GT. CSQ2) THEN DELTA = 1./ (1.+(CSQ1-CSQ2)/(TCHISQ-CSQ2)) + 0.5 ELSE C Backup if AF hasn't decreased DELTA = 1. END IF DO IC = 1,NCOMP DO IG = 1,NFILT DO K = 1,8 IF (VARY(IC,IG,K)) THEN PARM(IC,IG,K) 1 = PARM(IC,IG,K)-GRAD(IC,IG,K)*DELTA END IF END DO END DO END DO CALL CADJMOD CALL BINVIS CALL CAGRFAC IF (TCHISQ .GT. CSQ2) THEN C Model diverged, so backup to lowest minimum. DO IC = 1, NCOMP DO IG = 1, NFILT DO K = 1,8 IF (VARY(IC,IG,K)) THEN PARM(IC,IG,K) 1 = PARM(IC,IG,K)+GRAD(IC,IG,K)*(DELTA-1.) END IF END DO END DO END DO CALL CADJMOD CALL BINVIS CALL CAGRFAC END IF 200 CONTINUE JITER = JITER + 1 C If AF doesn't change significantly by 2nd iteration, quit here IF (ITER.GT.1 .AND. (OCHISQ-TCHISQ)/OCHISQ .LT. CSQLIM .AND. + OCHISQ .GT. TCHISQ) THEN WRITE (OUTC,'(A,F7.3,A)',ERR=220) + ' Convergence is slow (',TCHISQ,')' IF (QUERY(' Continue? (Y or N): ')) THEN CSQLIM = CSQLIM/2. ELSE GOTO 220 END IF END IF IF (CCFLAG) GO TO 225 210 CONTINUE CALL PUTOUT(' ') CALL PUTOUT(' Iteration halted after maximum no. '// 1 'of iterations along gradient.') GO TO 230 220 CALL PUTOUT(' ') CALL PUTOUT(' Iteration halted; model seems to have '// 1 'converged along gradient.') GO TO 230 225 CALL PUTOUT(' ') CALL PUTOUT(' Iteration halted by user.') CALL CTRLCEN 230 CALL CFIXMOD RETURN END