C   25/06/78 601302202  MEMBER NAME  CURFIT5  (S)           FORTRAN
      SUBROUTINE CURFIT
     ,(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,SIGMAA,FLAMDA,YFIT,CHISQR,
     , IPRO)
C      SUBROUTINE CURFIT
C
C    COPIED FROM F22YAM.UTILTYS(CURFIT). MODIFIED TO PRINT FOR DEBUGGING
C     NON CONVERGENCE SITUATION CHANGED  7.10.85  J.OLSSON
C     DSQRT ARGUMENT MADE ABSOLUTE...   30.01.86  J.OLSSON
C
C      ZWECK:
C      MAKE A LEAST-SQUAREA FIT TO A NON LINEAR FUNCTION
C      WITH A LINEARIZATION OF THE FITTING FUNCTION
C
C      USAGE:
C      CALL CURFIT(X,Y,SIGMAY,NPYS,NTERMS,MODE,A,DELTAA,
C
C        SIGMAA,FLAMDA,YFIT,CHISQR)
C      DESCRICPTION OF PARAMETERS
C
C      X           ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C      Y           ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C      SIGMAY      ARRAY OF STANDART DEVIATIONS FOR Y DATA POINTS
C      NPTS        NUMBER OF PAIRS OF DATA POINTE
C      NTERMS      NUMBER OF PARAMETERS
C      MODE        DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
C                  +1 (INSTRUMENTAL) WEIGHT(I)=1./SIGMAY(I)**2
C                  0  (NO WEIGHTING) WEIGHT(I)=1.
C                  -1 (STATISTICAL) WEIGHT(I)=1./Y(I)
C      A           ARRAY OF PARAMETERS
C      DELTAA      ARRAY OF INCREMENTS FOR PARAMETER A (DERIVATION)
C      SIGMAA      ARRAY OF STANDART DEVIATIONS FOR PARAMETER A
C      FLAMDA      PROPORTION OF GRADIENT SEARCH INCLUDET
C      YFIT        ARRAY OF CALCULATED VALUES OF Y
C      CHISQR      REDUCED CHI SQUARE FOR FIT
C
C     IPRO:  PRINT STEERING. IF NONZERO, DEBUG PRINT IS PRODUCED
C
C
C      SUBROUTINES AND FUNCTIONS NEEDED
C      FUNCTN(X,I,A)
C         EVALUATES THE FITTING FUNCTION FOR THE I TH TERM
C      FCHISQ(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
C         EVALUATES REDUCED CHI SQUARE FOR FIT DATA
C      FDERIV(X,I,A,DELTAA,NTERMS,DERIV)
C         EVALUATES THE DERIVATIVES OF THE FITTING FUNCTION
C         FOR THE I TH TERM WITH RESECT TO EACH PARAMETER
C         STEPWIDTHS IST DELTAA
C      MATINV(ARRAY,NTERMS,DET)
C         AND CALCULATES ITS DETERNINANT
C         INVERTS A SYMMETRIC 2 DIMENSIONAL MATRIX OF DEGREE NTERMS
C
C      COMMENTS
C      THE SUBROUTINE RETURNS AFTER EACH SUCSESSFUL STEP
C      THE MAIN PROGRAMM HAS TO CHECK CHISQR, WETHER CONVERGENZ IS REACH
C
C      NTERMS <= 10
C      SET FLAMDA = .001 AT THE BEGINNING OF SEARCH
C
C***********************************************************************
C
C
      REAL * 8 ARRAY(10,10),FLAMDA
      REAL X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1),
     ,     YFIT(1), WEIGHT(300),ALPHA(15,15),BETA(15),DERIV(15),B(15)
      REAL YFITSA(300),ASAVE(10),ASISA(10)
C      NPTS <= 300
C
      IF(NTERMS.GT.10.OR.NTERMS.LT.1) GO TO 13
      IF(NPTS.GT.300) GO TO 13
      NFREE=NPTS-NTERMS
      IF(NFREE.GT.0) GO TO 20
13    CHISQR=0.
      GO TO 110
C
C         EVALUATE WEIGHTS
C
20    DO 30 I=1,NPTS
21     IF(MODE)22,27,29
22    IF(Y(I))25,27,23
23    WEIGHT(I)=1./Y(I)
      GO TO 30
25    WEIGHT(I)=1./(-Y(I))
      GO TO 30
27    WEIGHT(I)=1.
      GO TO 30
29    IF(SIGMAY(I)) 27,27,290
290   WEIGHT(I)=1./SIGMAY(I)/SIGMAY(I)
30    CONTINUE
C
C  PRINT INPUT AND CALCULATED WEIGHTS
      IF(IPRO.EQ.0) GO TO 7701
      WRITE(6,7702) NPTS,NTERMS,MODE
7702  FORMAT(' CURFIT INPUT:::  NPTS NTERMS MODE ',3I5)
      WRITE(6,7703) (A(K),K=1,NTERMS)
7703  FORMAT(' A: ',10(F10.4,1X))
      WRITE(6,7704) (DELTAA(K),K=1,NTERMS)
7704  FORMAT(' DELTAA: ',10(F10.4,1X))
      WRITE(6,7706) (X(K),K=1,NPTS)
7706  FORMAT(' X: ',10(F10.4,1X))
      WRITE(6,7707) (Y(K),K=1,NPTS)
7707  FORMAT(' Y: ',10(F10.4,1X))
      WRITE(6,7705) (SIGMAY(K),K=1,NPTS)
7705  FORMAT(' SIGMAY: ',10(F10.4,1X))
      WRITE(6,7709) (WEIGHT(K),K=1,NPTS)
7709  FORMAT(' WEIGHT: ',10(F10.4,1X))
7701  CONTINUE
C
C         EVALUATE ALPHA AND BETA MATRICES
C
31    DO 34 J=1,NTERMS
      BETA(J)=0.
      DO 34 K=1,J
34    ALPHA(J,K)=0.
41    DO 50 I=1,NPTS
      IF(IPRO.EQ.1) WRITE(6,1431) I
1431  FORMAT('   CALLING FDERIV,                  I = ',I10)
      CALL FDERIV(X,I,A,DELTAA,NTERMS,DERIV)
      FUNXIA = FUNCTN(X,I,A)
      IF(IPRO.EQ.1) WRITE(6,1461) I,X(I),FUNXIA,(A(II),II=1,NTERMS)
1461  FORMAT(' I X FUN  A.. ',I3,F8.2,E14.6,8F12.4)
      IF(IPRO.EQ.1) WRITE(6,1432) (DERIV(JK),JK=1,NTERMS)
1432  FORMAT('   DERIV: ',10(F10.4,1X))
      DO 46 J=1,NTERMS
      BETA(J)=BETA(J)+WEIGHT(I)*(Y(I)-FUNXIA)*DERIV(J)
      DO 46 K=1,J
46    ALPHA(J,K)=ALPHA(J,K)+WEIGHT(I)*DERIV(J)*DERIV(K)
50    CONTINUE
51    DO 53 J=1,NTERMS
      DO 53 K=1,J
53    ALPHA(K,J)=ALPHA(J,K)
      IF(IPRO.EQ.1) WRITE(6,7801) ((ALPHA(J,K),J=1,NTERMS),K=1,NTERMS)
7801  FORMAT(' ALPHA: ',10(F10.4,1X))
      IF(IPRO.EQ.1) WRITE(6,7802) (BETA(J),J=1,NTERMS)
7802  FORMAT(' BETA: ',10(F10.4,1X))
C
C         EVALUATE CHI SQUARE AT STARTING POINT
C
61    DO 62 I=1,NPTS
      IF(IPRO.EQ.1) WRITE(6,1311) I,X(I),(A(JJ),JJ=1,NTERMS)
1311  FORMAT(' LOOP62: I X A..',I3,F10.2,8E12.4)
62    YFIT(I)=FUNCTN(X,I,A)
63    CHISQ1=FCHISQ(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
C
C         INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARAMETERS
C
71    DO 74 J=1,NTERMS
      IF(ABS(ALPHA(J,J)).LE.1.E-30) ALPHA(J,J)=1.E-30
      DO 73 K=1,NTERMS
      IF(ABS(ALPHA(K,K)).LE.1.E-30) ALPHA(K,K)=1.E-30
73    ARRAY(J,K)=ALPHA(J,K)/SQRT(ALPHA(J,J)*ALPHA(K,K))
74    ARRAY(J,J)=1.D0+FLAMDA
80    CALL MATINV(ARRAY,NTERMS,DET)
81    DO 84 J=1,NTERMS
      B(J)=A(J)
      IF(IPRO.EQ.1) WRITE(6,1502) J,B(J)
1502  FORMAT(' LOOP84: J BJ ',I3,E12.4)
      DO 84 K=1,NTERMS
      IF(IPRO.EQ.1) WRITE(6,1501) J,K,BETA(K),ARRAY(J,K),ALPHA(J,J),
     $ ALPHA(K,K),B(J)
1501  FORMAT(' J K,'2I2,' BETAK ARRAYJK ALFAJJ,KK ',4E12.4,' BJ ',E12.4)
84    B(J)=B(J)+BETA(K)*ARRAY(J,K)/SQRT(ALPHA(J,J)*ALPHA(K,K))
C
C         IF CHI SQUARE INCREASED , INCREASE FLAMDA AND TRY AGAIN
C
91    DO 92 I=1,NPTS
      IF(IPRO.EQ.1) WRITE(6,1312) I,X(I),(B(JJ),JJ=1,NTERMS)
1312  FORMAT(' LOOP92: I X B..',I3,F10.2,8E12.4)
92    YFIT(I)=FUNCTN(X,I,B)
      IF(IPRO.EQ.1) WRITE(6,4803) (YFIT(I),I=1,NPTS)
4803  FORMAT(' YFIT: ',10F10.4)
93    CHISQR=FCHISQ(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      IF(IPRO.EQ.1) WRITE(6,7803) CHISQ1,CHISQR,FLAMDA
7803  FORMAT(' CURFIT: CHISQ1,CHISQR,FLAMDA ',2E14.6,D12.4)
      IF(CHISQ1-CHISQR) 95,101,101
95    CONTINUE
C     FLAMDA=10.D0*FLAMDA
C     GO TO 71
C
C INCREASING FLAMDA DOES NOT WORK, RESTORE PREVIOUS RESULT INSTEAD.
C
      DO 1641  J1 = 1,NTERMS
      A(J1) = ASAVE(J1)
1641  SIGMAA(J1) = ASISA(J1)
      DO 1642  J2 = 1,NPTS
1642  YFIT(J2) = YFITSA(J2)
      GO TO 110
C
C         EVALUATE PARAMETERS AND UNCERTAINTIES
C
101   DO 103 J=1,NTERMS
      A(J)=B(J)
103   SIGMAA(J)=DSQRT(DABS(ARRAY(J,J))/ABS(ALPHA(J,J)))
      FLAMDA=FLAMDA/10.D0
      IF(IPRO.EQ.1) WRITE(6,7805) FLAMDA,(SIGMAA(JK),JK=1,8)
7805  FORMAT(' FLAMDA ',D12.4,' SIGAA ',8(F10.4,1X))
      IF(IPRO.EQ.1) WRITE(6,7807) (A(JK),JK=1,NTERMS)
7807  FORMAT(' NEW A: ',8(F10.4,1X))
C
C  SAVE FIT RESULT
C
      DO 1712  J1 = 1,NTERMS
      ASAVE(J1) = A(J1)
1712  ASISA(J1) = SIGMAA(J1)
      CHISQO = CHISQR
      DO 1713  J2 = 1,NPTS
1713  YFITSA(J2) = YFIT(J2)
110   CONTINUE
C
      IF(IPRO.EQ.1) WRITE(6,1714) (ASAVE(JJ),JJ=1,NTERMS)
1714  FORMAT(' ASAVE AT EXIT ',8E12.4)
      IF(IPRO.EQ.1) WRITE(6,1715) (ASISA(JJ),JJ=1,NTERMS)
1715  FORMAT(' ASISA AT EXIT ',8E12.4)
      IF(IPRO.EQ.1) WRITE(6,1716) (YFITSA(JJ),JJ=1,NPTS)
1716  FORMAT(' YFITSA EXIT ',10F10.2)
C
      RETURN
      END
