!+ Subroutine to derive an ESFT.
!
! Method:
!	An array of geiven transmissivities is fitted using
!	an algorithm originally devised by W. J. Wiscombe and
!	modified by J.-F. Geleyn.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE ESFT_FIT(IERR, MAX_N_TERM, TOL_FIT
     &   , N_DATA_IN, U_IN, TRANS_IN
     &   , N_TERM, K, A
     &   , TRANS_FIT
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/esft_parameters.acf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/precision_esft.tcf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS
      INTEGER 	!        , INTENT(OUT)
     &     IERR                               
!             ERROR FLAG
     &   , N_TERM                             
!             NUMBER OF TERMS USED
      INTEGER 	!        , INTENT(IN)
     &     MAX_N_TERM                         
!             MAXIMUM NUMBER OF TERMS
     &                                        
!         TO BE USED IN FIT
     &   , N_DATA_IN                          
!             NUMBER OF INPUT DATA POINTS
      REAL*8	!        , INTENT(IN)
     &     TOL_FIT                            
!             TOLERANCE FOR ESFT FIT
     &   , U_IN(NPD_AMOUNT)               
!             AMOUNTS OF ABSORBER
     &   , TRANS_IN(NPD_AMOUNT)           
!             ACTUAL TRANSMISSIVITIES
      REAL*8	!        , INTENT(OUT)
     &     K(NPD_ESFT_TERM)                     
!             ESFT EXPONENTS
     &   , A(NPD_ESFT_TERM)                     
!             ESFT WEIGHTS
     &   , TRANS_FIT(NPD_AMOUNT+1)        
!             FITTED TRANSMISSIONS
!
!     LOCAL VARIABLES
      INTEGER
     &     I                                  
!             LOOP VARIABLE
     &   , N_DATA                             
!             NUMBER OF DATA POINTS USED
     &   , ITERATION                          
!             ITERATION COUNT
      LOGICAL
     &     L_END                              
!             TRUE TO TERMINATE
     &   , L_CONVERGE                         
!             TRUE IF CONVERGED
     &   , L_RESOLVE                          
!             TRUE OF SVD TO BE RESOLVED
      REAL*8
     &     WEIGHT(NPD_AMOUNT+1)           
!             WEIGHTS APPLIED TO RESIDUALS
     &   , EXP_KU(NPD_ESFT_TERM, NPD_AMOUNT+1)
!             PRECALCULATED EXPONENTIALS
     &   , RESIDUAL(NPD_AMOUNT+1)         
!             RESIDUALS IN FIT
     &   , RELATIVE_ERROR(NPD_AMOUNT+1)   
!             RELATIVE ERRORS
     &   , U(NPD_AMOUNT+1)                
!             AUGMENTED ARRAY OF AMOUNTS
     &   , TRANS(NPD_AMOUNT+1)            
!             AUGMENTED ARRAY OF TRANS.
     &   , A_OLD(NPD_ESFT_TERM)                 
!             OLD ESFT WEIGHTS
     &   , K_MINIMUM
     &   , R_MINIMUM
     &   , SUM_WEIGHT
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     EXPONENT_KU, LEAST_SQUARES, REMOVE_NEGATIVE
     &   , RESIDUAL_FIT, LOCATE_MINIMUM, TERMINATION
     &   , REORDER, ADJUST_FIT, CALC_FIT
!
!
!
!     THE INPUT DATA ARE AUGMENTED BY AN EXTRA POINT FOR WHICH THERE IS
!     NO ABSORBER: THIS ALLOWS US TO FORCE THE TRANSMISSIONS TO SUM TO 1
!     AS CLOSELY AS REQUIRED.
      N_DATA=N_DATA_IN+1
      SUM_WEIGHT=0.0D+00
      DO I=1, N_DATA_IN
         U(I)=U_IN(I)
         TRANS(I)=TRANS_IN(I)
         WEIGHT(I)=1.0D+00/(1.0D+00-WEIGHT_ABSORB*TRANS(I)
     &      -WEIGHT_TRANS*(1-TRANS(I)))**2
         SUM_WEIGHT=SUM_WEIGHT+WEIGHT(I)
      ENDDO
      DO I=1, N_DATA_IN
         WEIGHT(I)=WEIGHT_ESFT_DATA*WEIGHT(I)/SUM_WEIGHT
      ENDDO
      U(N_DATA)=0.0D+00
      TRANS(N_DATA)=1.0D+00
      WEIGHT(N_DATA)=1.0D+00-WEIGHT_ESFT_DATA
!
!     INITIALIZE FOR THE FIT.
      N_TERM=1
      K(1)=0.0D+00
      ITERATION=0
!
!     THIS IS THE MAIN LOOP
100   CONTINUE
!     ADVANCE THE ITERATION COUNT.
      ITERATION=ITERATION+1
!
!     PRINT OUT DETAILS OF THE ITERATION CYCLE.
      WRITE(IU_MONITOR, '(/A, I5)') 'ITERATION = ', ITERATION
!
200   CONTINUE
!        CALCULATE THE VALUES OF THE EXPONENTIALS FOR THE K'S TO HAND.
         CALL EXPONENT_KU(N_TERM, N_DATA
     &      , K, U
     &      , EXP_KU
     &      )
!
!        SOLVE THE LINEAR LEAST SQUARES EQUATIONS TO GET THE FIT.
         CALL LEAST_SQUARES(IERR, N_TERM, N_DATA
     &      , TRANS, WEIGHT, EXP_KU
     &      , A
     &      )
         IF (IERR.NE.I_NORMAL) RETURN
!
!
!        IF THE FIT PRODUCES NEGATIVE VALUES OF A THESE MUST BE REMOVED
!        AND THE LEAST SQUARES PROBLEM RESOLVED.
         CALL REMOVE_NEGATIVE(N_DATA, N_TERM
     &      , A, K, A_OLD
     &      , L_RESOLVE
     &      , ITERATION
     &      )
      IF (L_RESOLVE) GOTO 200
!
!     FIND THE RESIDUALS AT THE DATA POINTS RESULTING FROM THIS FIT.
      CALL RESIDUAL_FIT(N_TERM, N_DATA
     &   , A, EXP_KU, TRANS
     &   , RESIDUAL
     &   )
!
      CALL LOCATE_MINIMUM(N_TERM, N_DATA
     &   , K, WEIGHT, RESIDUAL, U
     &   , K_MINIMUM, R_MINIMUM
     &   )
!
!     ASSESS THE CRITERIA FOR TERMINATION.
      CALL TERMINATION(N_DATA, ITERATION, N_TERM
     &      , WEIGHT, RESIDUAL, R_MINIMUM
     &      , TOL_FIT, L_END, L_CONVERGE
     &      )
!
!     IF A TERMINATING CONDITION HAS NOT BEEN REACHED EXECUTE ANOTHER
!     CYCLE.
      IF (.NOT.L_END) THEN
!        ADD A NEW VALUE OF K TO THE LIST AT THE MINIMUM OF THE
!        RESIDUAL POLYNOMIAL.
!        IF WE HAVE ALREADY USED THE MAXIMUM NUMBER OF TERMS, THE
!        LAST TERM IS ARBITRARILY REPLACED.
         IF ( (N_TERM.LT.MAX_N_TERM).AND.(N_TERM.LT.N_DATA) )
     &         N_TERM=N_TERM+1
         K(N_TERM)=K_MINIMUM
         A(N_TERM)=0.0D+00
!        REORDER THE ROOTS SO THAT K(N_TERM) IS THE LARGEST.
         CALL REORDER(N_TERM, K, A)
!        RECORD THE OLD VALUES OF A FOR LATER COMPARISON.
         DO I=1, N_TERM
            A_OLD(I)=A(I)
         ENDDO
         GOTO 100
      ENDIF
!
!
!     ADJUST THE FIT TO FORCE THE SUM OF THE WEIGHTS TO UNITY.
      CALL ADJUST_FIT(N_TERM, A, K)
!
!     PRINT OUT THE TERMINATING INFORMATION
      IF (L_CONVERGE) THEN
         WRITE(IU_MONITOR, '(/A/)')
     &      'CONVERGENCE HAS BEEN ACHIEVED: '
     &      //'THE FINAL FIT IS AS FOLLOWS.'
      ELSE
         WRITE(IU_MONITOR, '(/A/)')
     &      'CONVERGENCE HAS NOT BEEN ACHIEVED: '
     &      //'THE LAST FIT OBTAINED WAS:'
      ENDIF
      CALL CALC_FIT(N_DATA, N_TERM
     &   , A, K, U, TRANS
     &   , TRANS_FIT, RELATIVE_ERROR
     &   )
      WRITE(IU_MONITOR, '(3X, A6, 13X, A10, 9X, A6)')
     &   'AMOUNT', 'CALCULATED', 'ACTUAL'
      DO I=1, N_DATA-1
         WRITE(IU_MONITOR, '(3(3X,1PE16.9))')
     &      U(I), TRANS_FIT(I), TRANS(I)
      ENDDO
      WRITE(IU_MONITOR, '(/, 3(3X,1PE16.9), /)')
     &   U(N_DATA), TRANS_FIT(N_DATA), TRANS(N_DATA)
      WRITE(IU_MONITOR, '(/A,I4)') 'NUMBER OF TERMS IN FIT = ',N_TERM
      WRITE(IU_MONITOR, '(''NO. OF TERM'', 3X
     &, ''COEFFICIENT OF EXPONENTIAL'' ,3X, ''VALUE OF K'')')
      DO I=1, N_TERM
         WRITE(IU_MONITOR, '(I5, 9X,1PE16.9, 13X, 1PE16.9)')
     &      I, A(I), K(I)
      ENDDO
!
!
      RETURN
      END
