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 FTTNTR(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                                                
C THE ITERATION IS FINISHED IF                                                  
C - NITER IS REACHED OR                                                         
C - CHISQ IS REACHED                                                            
C - CONVERGENCE IS ACHIEVED                                                     
C                                                                               
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 INPUT PARAMETERS                                                              
C NITER INTG NO. OF ITERATIONS TO BE DONE                                       
C   <0 USE PREVIOUS COMPUTED VALUES                                             
C   >0 USE GUESS                                                                
C RELAX REAL RELAXATION FACTOR                                                  
C CHISQ REAL CHI SQUARE                                                         
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 NITER,NACT,ISTAT,IPRINT,I,J,NITER1,ICOUNT,NN,NP,NNP
      INTEGER ISTAR,NCOL,NS,KZ6941,N,IDAT,NOFFSET,IFUN,K,IP,NILAST
      INTEGER ICOL(10),ACOL,AROW,NC,NROW

      REAL RELAX,PREC,V1,V2,ACTCH,VS,VE,CHIOLD,CHIPERC,AC1,AC2,AC3
      REAL DELTA,C1,C2
      REAL XVAL(10),X(8),W,FSUMSQ                                               

      CHARACTER*80 LINE                                                         
      CHARACTER*6 TYPE                                                          

      DOUBLE PRECISION CHI,PP,Q,Y,Y1,Y2,YOUT,CORRNRM                            

      LOGICAL NULL(10),VALID,NEXT,AUTO,ISEL,VCMATPR                             

      INCLUDE 'MID_INCLUDE:FITI.INC/NOLIST'                                     
      DOUBLE PRECISION FZDERIV(FZPARMAX)                                    
      DOUBLE PRECISION AL,B(FZPARMAX),A(FZPARMAX,FZPARMAX)
      INCLUDE 'MID_INCLUDE:FITC.INC/NOLIST'                                     

      EQUIVALENCE (XVAL(1),W)                                                   
C      EQUIVALENCE (XVAL(2),Y)
      EQUIVALENCE (XVAL(3),X(1))     

      DATA TYPE/'R*4'/                
                                                                                
      DO 200, I=1,10
         NULL(I) = .FALSE.
200   CONTINUE

      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)                                                    
      VCMATPR = FZMETPAR(1) .LE. 0.1                                          
      IPRINT = NINT(ABS(FZMETPAR(1)))                                         
      IF (IPRINT.EQ.0) IPRINT = ABS(NITER)                                      
      AUTO   = RELAX .LT. 0.                                                    
      RELAX  = ABS(RELAX)                                                       
      DO 10, I = 1,8  
          NULL(I) = .FALSE.                                                     
   10 CONTINUE                                                                  
      NEXT   = .TRUE.                                                           
      IF (NITER.GT.0) THEN                                                      
          DO 20, J = 1,FZNPTOT                                                 
              FZVALUE(J) = FZGUESS(J)                                       
   20     CONTINUE                                                              
          FZNITER = 0                                                         
          NITER1 = NITER 
      ELSE                                                                      
          NITER1 = -NITER                                                       
      END IF                                                                    
                                                                                
      VS     = MIN(V1,V2)                                                       
      VE     = MAX(V1,V2)                                                       
      FZCCHIS = 9999.                                                         
      FZCHISQ = 1.D-6                                                         
      FZRELAX = RELAX                                                         
      CORRNRM = 1.D6                                                            
      CHI    = 0.D0                                                             
      ICOUNT = 0                                                                
      IF (FZWEIGHT.EQ.0) THEN                                                 
          ISTAR  = 2                                                            
          W      = 1.                                                           
          NCOL   = FZNIND + 1                                                 
          ICOL(1) = 0                                                           
                                                                                
      ELSE                                                                      
          ISTAR  = 1                                                            
          NCOL   = FZNIND + 2                                                 
          ICOL(1) = FZWEIGHT                                                  
      END IF                                                                    
                                                                                
      ICOL(2) = FZDVAR                                                        
      NULL(1) = .FALSE.                                                         
      DO 30, I = 1,FZNIND                                                      
          ICOL(I+2) = FZIVAR(I)                                               
   30 CONTINUE                                                                  

      CALL TBIGET(FZIDEN,NC,NROW,NS,ACOL,AROW,ISTAT)                        
      KZ6941 = 0
99    KZ6941 = KZ6941 + 1
          IF ( .NOT. (NEXT)) GO TO 180                                          
          DO 50, I = 1,FZPARMAX 
              DO 40, J = 1,FZPARMAX
                  A(J,I) = 0.D0                                                 
   40         CONTINUE                                                          
              B(I)   = 0.D0                                                     
   50     CONTINUE                                                              
          CHIOLD = CHI                                                          
          CHI    = 0.D0                                                         
C                                                                               
C ITERATION FOR TABLES                                                          
C                                                                               
          N      = 0                                                            
          DO 100 IDAT = 1,NROW                                                  
              CALL TBSGET(FZIDEN,IDAT,ISEL,ISTAT)                             
C              XVAL(2) = Y     ! to avoid problems with the equivalence
              CALL TBRRDR(FZIDEN,IDAT,NCOL,ICOL(ISTAR),                  
     +                    XVAL(ISTAR),NULL(ISTAR),ISTAT)                        
              Y=XVAL(2)
              VALID  = ISEL .AND. ( .NOT. NULL(1)) .AND.                        
     +                 ( .NOT. NULL(2)) .AND. ( .NOT. NULL(3))                  
C                                                                               
C ...       COMPUTE FUNCTION VALUES AND DERIVATIVES (0 IF FIXED PARAM)          
C                                                                               
              IF (V1.NE.V2) VALID  = VALID .AND. XVAL(1) .GE. VS .AND.          
     +                               XVAL(1) .LE. VE                            
              IF (VALID) THEN                                                   
                  N      = N + 1                                                
                  Y1     = 0.                                                   
                  NOFFSET = 1                                                   
                  DO 60 IFUN = 1,FZNFUN                                       
                      CALL FTFUNC(FZFCODE(IFUN),FZNIND,X,                   
     +                            FZACTPAR(IFUN),FZVALUE(NOFFSET),          
     +                            YOUT,FZDERIV(NOFFSET))                      
                      Y1     = Y1 + YOUT                                        
                      NOFFSET = NOFFSET + FZACTPAR(IFUN)                      
   60             CONTINUE                                                      
                  DO 70 J = 1,FZNPTOT                                         
                      IP     = FZFIXED(J)                                     
                      IF (IP.EQ.0) FZDERIV(J) = 0.D0                          
                      IF (IP.GT.0) THEN                                         
                          FZDERIV(IP) = FZDERIV(IP) +                       
     +                                    FZDERIV(J)/FZPFAC(J)              
                          FZDERIV(J) = 0.D0                                   
                      END IF                                                    
                                                                                
   70             CONTINUE                                                      
C                                                                               
C ...         COMPUTE RESIDUALS, WEIGHTED RESIDUALS AND CHI**2                  
C                                                                               
                  W      = 1./ (W*W)                                            
                  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 90 J = 1,FZNPTOT                                         
                      AL     = FZDERIV(J)                                     
                      DO 80 K = 1,FZNPTOT                                     
                          A(J,K) = A(J,K) + AL*PP*FZDERIV(K)                  
   80                 CONTINUE                                                  
                      B(J)   = B(J) + PP*Q*AL                                   
   90             CONTINUE                                                      
              END IF                                                            
                                                                                
  100     CONTINUE                                                              
C                                                                               
C ...     ADD (PROBABILISTIC) PARAMETER CONSTRAINTS                             
C                                                                               
          NN     = N                                                            
          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))            
C     ^                                                                         
C   To be modified for soft constraints                                         
                  END IF                                                        
                                                                                
              END IF                                                            
                                                                                
  110     CONTINUE                                                              
C                                                                               
C ...     INVERT NORMAL EQUATION MATRIX TO GET THE COVARIANCE MATRIX            
C                                                                               
          CALL DMATIN(A,FZNPTOT,FZPARMAX,FZFIXED,ISTAT) 
          NNP    = FZNPTOT                                                    
          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)                                
          CORRNRM = 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).EQ.0) THEN                                       
                  FZERROR(I) = 0.D0                                           
                                                                                
              ELSE                                                              
                  Q      = RELAX*Q                                              
                  CORRNRM = CORRNRM + Q*Q                                       
                  FZERROR(I) = DSQRT(A(I,I)*CHI)                              
                  FZVALUE(I) = FZVALUE(I) + Q                               
              END IF                                                            
                                                                                
  130     CONTINUE                                                              
          CORRNRM = DSQRT(CORRNRM)                                              
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)                                                 
          NEXT   = ICOUNT .LT. NITER1 .AND. CORRNRM .GT. PREC                   
          IF ( .NOT. NEXT .OR. 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                                                                               
C ...                                                                           
C                                                                               
          IF (AUTO) THEN                                                        
              IF (NILAST.GT.10) THEN                                            
C                                                                               
C ... CHECK FOR CONVERGENCE                                                     
C                                                                               
                  IF (FZCHI(NILAST-1).EQ.FZCHI(NILAST)) GO TO 190           
                  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 190                               
                  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
  180 CONTINUE                                                                  
  190 CALL STTPUT(' ',ISTAT)                                                    
      NACT   = FZNITER                                                        
      ACTCH  = FZCCHIS                                                        
      FZNDAT = N                                                              
      IF (CORRNRM.LE.PREC) THEN                                                 
          CALL STTPUT('  --> NR  : Convergence achieved <--',ISTAT)             
                                                                                
      ELSE                                                                      
          CALL STTPUT('*** WARN-NR  : No convergence reached ***',ISTAT)        
      END IF                                                                    
                                                                                
      CALL STTPUT(' ',ISTAT)                                                    
C                                                                               
C ... PRINT CORREL. MATRIX                                                      
C                                                                               
      IF (VCMATPR) 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,4X,0PF5.2)                         
 9020 FORMAT (2X,1PD15.7)                                                       
      END                                                                       
