C===========================================================================
C Copyright (C) 1995-2005 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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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                                                                               
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.VERSION: 1.0  ESO-FORTRAN Conversion, AA  17:25 - 13 JAN 1988                 
C 051109	last modif
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                            

      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)                                    
      DOUBLE PRECISION AL,B(FZPARMAX),A(FZPARMAX,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,FZPARMAX                                                            
          DO 20, J = 1,FZPARMAX                                                        
              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,FZPARMAX,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,FZPARMAX)
      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                                                                       
