* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE FITPNL(X,Y,N,NTERMS,IPAR,A,NITER,
     $NUM,YFIT,SIGMA,SIGMAA,EPS,BAVARD, GOOD)
C----------------------------------------------------------------------
C	Fit non lineaire
C
C	N est le nombre de points
C       NTERMS est le nombre de terms a fitter
C       A(J),J=1,NTERMS sont les estimations initiales des parametres
C       X(K,J),J=1,N K=1,NVAR et Y(J),J=1,N sont les valeurs des X et Y
C       NITER est le nombre maximal d'iterations tolere
C       NUM est le nombre d'iterations effectuees
C       YFIT(J),J=1,N  sont les valeurs fittees
C	SIGMA est le sigma du fit (S(Y(J)-YFIT(J))**2)/N
C	SIGMAA(J),J=1,NTERMS sont les sigmas sur les parametres
C	EPS est la tolerance sur SQRT(CHISQ)
C
C-----------------------------------------------------------------------
      INCLUDE 'size.inc'
      INTEGER N                          ! Number of data points
      INTEGER NTERMS                     ! Number of parameters
      REAL*8 X(3,N)                      ! Data points coordinates
     $,Y(N)                              ! Data points values
     $,YFIT(N)                           ! Fitted values
      REAL*8 SIGMAA(NTERMS)              ! Errors on parameters
     $,A(MVAR)                           ! Parameters
      REAL*8 SIGMA                       ! Rms of residuals
      REAL*8 EPS                         ! Stability of minimum
      INTEGER NITER,NUM                  ! Maximum and current Number of iterations
      LOGICAL BAVARD                     ! Write down progress
      LOGICAL GOOD(N)
      REAL*8 FUNPNL,FCHISQ
      EXTERNAL FUNPNL,FCHISQ
*
      REAL*8 ALPHA(MVAR,MVAR),ARRAY(MVAR,MVAR),BETA(MVAR),
     $DERIV(MVAR),B(MVAR),AA(MVAR)
      INTEGER IPAR(MVAR),NTERMX,I,J,K,JJ,KK
      REAL*8 FLAMDA,RTS,CHISQ1,CHISQR
*
      NTERMX=MVAR
      FLAMDA=0.001
      RTS=57.29577951*3600.
*
      DO NUM=1,NITER
         DO J=1,NTERMS
            BETA(J)=0.
            DO K=1,J
               ALPHA(J,K)=0.
            ENDDO
         ENDDO
         DO I=1,N
            IF (GOOD(I)) THEN
               CALL DPOINL(X(1,I),A,DERIV)
               DO  J=1,NTERMS
                  JJ=IPAR(J)
                  BETA(J)=BETA(J)+(Y(I)-FUNPNL(X(1,I),A))*DERIV(JJ)
                  DO K=1,J
                     KK=IPAR(K)
                     ALPHA(J,K)=ALPHA(J,K)+DERIV(JJ)*DERIV(KK)
                  ENDDO
               ENDDO
            ENDIF
         ENDDO
         DO J=1,NTERMS
            DO K=1,J
               ALPHA(K,J)=ALPHA(J,K)
            ENDDO
         ENDDO
         DO I=1,N
            IF (GOOD(I)) YFIT(I)=FUNPNL(X(1,I),A)
         ENDDO
         CHISQ1=FCHISQ(Y,N,NTERMS,YFIT,GOOD)
*
100      CONTINUE
         DO J=1,NTERMS
            DO K=1,NTERMS
               ARRAY(J,K)=ALPHA(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K))
            ENDDO
            ARRAY(J,J)=1.+FLAMDA
         ENDDO
         CALL AMTINV(ARRAY,NTERMS)
         DO J=1,NTERMS
            JJ=IPAR(J)
            B(J)=A(JJ)
            DO  K=1,NTERMS
               B(J)=B(J)+BETA(K)*ARRAY(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K))
            ENDDO
         ENDDO
         DO J=1,NTERMX
            AA(J)=A(J)
         ENDDO
         DO J=1,NTERMS
            JJ=IPAR(J)
            AA(JJ)=B(J)
         ENDDO
         DO I=1,N
            IF (GOOD(I)) YFIT(I)=FUNPNL(X(1,I),AA)
         ENDDO
         CHISQR=FCHISQ(Y,N,NTERMS,YFIT,GOOD)
         IF (CHISQ1.LT.CHISQR) THEN
            FLAMDA = 10*FLAMDA
            GOTO 100
         ENDIF
*
* Converged ...
         FLAMDA=FLAMDA/10.
         SIGMA=DSQRT(CHISQR)
         DO J=1,NTERMS
            JJ=IPAR(J)
            A(JJ)=B(J)
            SIGMAA(J)=DSQRT(ARRAY(J,J)/ALPHA(J,J))*SIGMA
         ENDDO
         IF (BAVARD) THEN
            WRITE(6,1001) NUM,SIGMA*RTS
            DO K=1,NTERMS
               JJ=IPAR(K)
               WRITE(6,1000) NUM,JJ,A(JJ)*RTS,SIGMAA(K)*RTS
            ENDDO
         ENDIF
         IF ((CHISQ1-CHISQR).LT.EPS*CHISQ1) RETURN
      ENDDO
1000  FORMAT(' ite = ',I3,' a(',I2,') =',F12.2,' err = ',F12.3)
1001  FORMAT(' ite = ',I3,' sigma = ',F12.2)
      END
