C @(#)ftintr.for 17.1.1.1 (ES0-DMD) 01/25/02 17:10:46 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C SUBROUTINE FTINTR(YY,W1,N1,N2,N3,NITER,RELAX,PREC,V1,V2,NACT, + ACTCH,ISTAT) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 17:25 - 13 JAN 1988 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.IDENTIFICATION C C FITLIB.FOR VERSION 1.0 27 MAR 1984 C C.PURPOSE C C INTERFACE ROUTINES FOR THE FITTING STRUCTURES C ITERATE ON THE FITTING PROCESS FOT IMAGE DATA C THE ITERATION IS FINISHED IF C - NITER IS REACHED OR C - CHISQ IS REACHED C - CONVERGENCE IS ACHIEVED C C.ALGORITHM C C USE MIDAS I/O INTERFACES TO FRAMES AND TABLES C C.KEYWORDS C C NON LINEAR FITTING C C C---------------------------------------------------------------- C C C INPUT PARAMETERS C YY REAL MAPPED VALUES OF THE IMAGE C W1 REAL OPTIONAL WEIGHTS C NITER INTG NO. OF ITERATIONS TO BE DONE C RELAX REAL RELAXATION FACTOR. C AUTOMATIC COMPUTATION IF NEGATIVE C PREC REAL PRECISION C V1,V2 REAL OPTIONAL FITTING INTERVAL (IF V1.NE.V2) C C OUTPUT PARAMETERS C NACT INTG ACTUAL NUMBER OF ITERATIONS C ACTCH REAL ACTUAL CHI SQ. C ISTAT INTG STATUS RETURN C INTEGER N1,N2,N3,NITER,NACT,ISTAT,IPRINT,J,K INTEGER NITER1,I,I1,I2,KZ6901,ICOUNT INTEGER IX1,IX2,NOFFSET,IFUN,IP,NN,NP,NNP,IX3,NILAST CHARACTER*80 LINE DOUBLE PRECISION CHI,PP,Q,Y,Y1,Y2,YOUT,CORRNM DOUBLE PRECISION AL,B(50),A(50,50) REAL X(3),YY(N1,N2,N3),W1(N1,N2,N3) REAL RELAX,PREC,V1,V2,ACTCH,V11,V22,CHIOLD,FSUMSQ REAL AC1,AC2,AC3,DELTA,C1,C2,W,CHIPERC LOGICAL AUTO,VCMATR INCLUDE 'MID_INCLUDE:FITI.INC/NOLIST' DOUBLE PRECISION FZDERIV(FZPARMAX) INCLUDE 'MID_INCLUDE:FITC.INC/NOLIST' C CALL STTPUT(' ',ISTAT) CALL STTPUT('Non-linear Least Square Fitting :'// + ' Newton-Raphson Method (defaulted)',ISTAT) CALL STTPUT('_______'// + '___________________________________________________________' + ,ISTAT) CALL STTPUT(' ',ISTAT) CALL STTPUT(' ',ISTAT) CALL STTPUT( + ' Print Max. Fct. eval. Prec. on Param. Relax. Fact.' + ,ISTAT) WRITE (LINE,9000) NINT(FZMETPAR(1)),NINT(FZMETPAR(2)), + FZMETPAR(3),FZMETPAR(4) CALL STTPUT(LINE,ISTAT) CALL STTPUT(' ',ISTAT) VCMATR = FZMETPAR(1) .LE. 0.1 IPRINT = NINT(ABS(FZMETPAR(1))) IF (IPRINT.EQ.0) IPRINT = ABS(NITER) AUTO = RELAX .LT. 0. RELAX = ABS(RELAX) IF (NITER.GT.0) THEN DO 10 J = 1,FZNPTOT FZVALUE(J) = FZGUESS(J) 10 CONTINUE FZNITER = 0 NITER1 = NITER ELSE NITER1 = -NITER END IF IF (V1.NE.V2) THEN V11 = MAX(V1,FZSTART(1)) V22 = MIN(V2,FZSTART(1)+ (FZNPIX(1)-1)*FZSTEP(1)) I1 = (V11-FZSTART(1))/FZSTEP(1) + 1 I2 = (V22-FZSTART(1))/FZSTEP(1) + 1 ELSE I1 = 1 I2 = FZNPIX(1) END IF FZCCHIS = 9999. FZCHISQ = 1.D6 FZRELAX = RELAX CHI = 0.D0 ICOUNT = 0 CORRNM = 1.D6 KZ6901 = 0 99 KZ6901 = KZ6901 + 1 IF ( .NOT. (ICOUNT.LT.NITER1 .AND. CORRNM.GT.PREC)) GO TO 166 DO 30 I = 1,50 DO 20 J = 1,50 A(J,I) = 0.D0 20 CONTINUE B(I) = 0.D0 30 CONTINUE CHIOLD = CHI CHI = 0.D0 C C ITERATION ON IMAGE DATA C DO 100 IX3 = 1,FZNPIX(3) X(3) = FZSTART(3) + (IX3-1)*FZSTEP(3) DO 90 IX2 = 1,FZNPIX(2) X(2) = FZSTART(2) + (IX2-1)*FZSTEP(2) DO 80 IX1 = I1,I2 X(1) = FZSTART(1) + (IX1-1)*FZSTEP(1) IF (FZWEIGHT.EQ.0) THEN W = 1. ELSE W = W1(IX1,IX2,IX3) W = 1./ (W*W) END IF Y = YY(IX1,IX2,IX3) C C ... COMPUTE FUNCTION VALUES AND DERIVATIVES (0 IF FIXED PARAM) C Y1 = 0.D0 NOFFSET = 1 DO 40 IFUN = 1,FZNFUN CALL FTFUNC(FZFCODE(IFUN),FZNIND,X, + FZACTPAR(IFUN),FZVALUE(NOFFSET), + YOUT,FZDERIV(NOFFSET)) Y1 = Y1 + YOUT NOFFSET = NOFFSET + FZACTPAR(IFUN) 40 CONTINUE DO 50 J = 1,FZNPTOT IP = FZFIXED(J) IF (IP.EQ.0) FZDERIV(J) = 0. IF (IP.GT.0) THEN FZDERIV(IP) = FZDERIV(IP) + + FZDERIV(J)/FZPFAC(J) FZDERIV(J) = 0. END IF 50 CONTINUE C C ... COMPUTE RESIDUALS, WEIGHTED RESIDUALS AND CHI**2 C Y2 = Y - Y1 PP = DBLE(W) IF (FZFLAG.EQ.1) PP = 1.D0/AMAX1(1.0,SNGL(Y1)) Q = Y2 CHI = CHI + Q*Q*PP C C ... BUILD NORMAL EQUATION MATRIX C DO 70 J = 1,FZNPTOT AL = FZDERIV(J) DO 60 K = 1,FZNPTOT A(J,K) = A(J,K) + AL*PP*FZDERIV(K) 60 CONTINUE B(J) = B(J) + PP*Q*AL 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE C C ... ADD (PROBABILISTIC) PARAMETER CONSTRAINTS C NN = N1*N2*N3 NP = 0 DO 110 I = 1,FZNPTOT IF (FZFIXED(I).GE.0) THEN A(I,I) = 1.0D0 ELSE NP = NP + 1 IF (FZUNCER(I).GT.0.) THEN NN = NN + 1 Q = 1.0D0/DBLE(FZUNCER(I)*FZUNCER(I)) A(I,I) = A(I,I) + Q B(I) = B(I) + Q* (FZGUESS(I)-FZVALUE(I)) END IF END IF 110 CONTINUE C C ... INVERT NORMAL EQUATION MATRIX TO GET THE COVARIANCE MATRIX C NNP = FZNPTOT CALL DMATIN(A,FZNPTOT,50,FZFIXED,ISTAT) IF (ISTAT.NE.0) THEN CALL STTPUT('*** ERR-NR : Problems inverting matrix ***', + ISTAT) RETURN END IF C C ... COMPUTE REDUCED CHI**2 AND CORRECTIONS FOR THE UNKNOWNS C FSUMSQ = CHI CHI = DMAX1(CHI/DBLE(NN-NP),0.D0) CORRNM = 0.D0 DO 130 I = 1,FZNPTOT Q = 0.0D0 DO 120 J = 1,FZNPTOT Q = Q + A(I,J)*B(J) 120 CONTINUE C C ... COMPUTE ERRORS OF THE CORRECTIONS (I.E. THE UNKNOWNS) C ... AND APPLY CORRECTIONS (DO THE ADJUSTMENT!) C IF (FZFIXED(I).GE.0) THEN FZERROR(I) = 0.D0 ELSE Q = RELAX*Q CORRNM = CORRNM + Q*Q FZERROR(I) = DSQRT(A(I,I)*CHI) FZVALUE(I) = FZVALUE(I) + Q END IF 130 CONTINUE CORRNM = DSQRT(CORRNM) C C ... LINEAR CONSTR. BETWEEN TWO PARAMS C DO 140 I = 1,FZNPTOT IP = FZFIXED(I) IF (IP.GT.0) THEN FZERROR(I) = FZPFAC(I)*FZERROR(IP) FZVALUE(I) = FZPFAC(I)*FZVALUE(IP) END IF 140 CONTINUE ICOUNT = ICOUNT + 1 FZNITER = FZNITER + 1 FZCCHIS = SNGL(CHI) IF (MOD(ICOUNT,IPRINT).EQ.0) THEN CALL STTPUT(' ',ISTAT) CALL STTPUT( + ' Iter Fct. Eval. Sum of Squares Red. Chisq. % Decr.' + ,ISTAT) IF (ICOUNT.GT.1) THEN CHIPERC = 100.* (CHIOLD-CHI)/CHIOLD WRITE (LINE,9010) FZNITER,FZNITER,FSUMSQ,FZCCHIS, + CHIPERC ELSE WRITE (LINE,9010) FZNITER,FZNITER,FSUMSQ,FZCCHIS END IF CALL STTPUT(LINE,ISTAT) CALL STTPUT(' ',ISTAT) CALL STTPUT(' Parameters ',ISTAT) DO 150 I = 1,FZNPTOT WRITE (LINE,9020) FZVALUE(I) CALL STTPUT(LINE,ISTAT) 150 CONTINUE CALL STTPUT(' ',ISTAT) END IF C C ... STORE THE CHISQ C IF (FZNITER.LE.100) THEN FZCHI(FZNITER) = CHI NILAST = FZNITER ELSE DO 160 I = 1,99 FZCHI(I) = FZCHI(I+1) 160 CONTINUE FZCHI(100) = CHI NILAST = 100 END IF C IF (AUTO) THEN IF (NILAST.GT.10) THEN C C ... CHECK FOR CONVERGENCE C IF (FZCHI(NILAST-1).EQ.FZCHI(NILAST)) GO TO 170 AC1 = ALOG(AMAX1(FZCHI(NILAST-2),0.0001)) AC2 = ALOG(AMAX1(FZCHI(NILAST-1),0.0001)) AC3 = ALOG(AMAX1(FZCHI(NILAST),0.0001)) C1 = AC1 - AC2 C2 = AC2 - AC3 IF (ABS(C2).LE.1.E-5) GO TO 170 DELTA = ABS((C1-C2)/C2) IF (C2.GT.0.0 .AND. DELTA.LE.0.25) RELAX = (1.0-DELTA)**4 END IF END IF GOTO 99 166 CONTINUE 170 CALL STTPUT(' ',ISTAT) NACT = FZNITER ACTCH = FZCCHIS FZNDAT = N1*N2*N3 IF (V1.NE.V2) THEN FZNPIX(1) = I2 - I1 + 1 FZSTART(1) = FZSTART(1) + (I1-1)*FZSTEP(1) END IF IF (CORRNM.LE.PREC) THEN CALL STTPUT(' --> NR : Convergence achieved <--',ISTAT) ELSE CALL STTPUT('*** WARN-NR : No convergence reached ***',ISTAT) END IF C C ... PRINT CORREL. MATRIX C IF (VCMATR) CALL DCPRIN(A,NNP,50) RETURN 9000 FORMAT (I6,7X,I7,7X,1PE9.1,8X,0PF4.2) 9010 FORMAT (I5,6X,I5,3X,1PE12.4,2X,1PE12.4,5X,0PF5.2) 9020 FORMAT (2X,1PE15.7) END