C SUBROUTINE CHIFIT.F 
C
C SOURCE
C   BEVINGTON, PAGES 228-231.
C
C PURPOSE
C   MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION
C      WITH A PARABOLIC EXPANSION OF CHI SQUARE 
C
C USAGE 
C   CALL CHIFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA,
C      SIGMAA, YFIT, CHISQR)
C
C DESCRIPTION OF PARAMETERS
C   X      - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C   Y      - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C   SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS
C   NPTS   - NUMBER OF PAIRS OF DATA POINTS
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 PARAMETERS A
C   SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A
C   YFIT   - ARRAY OF CALCULATED VALUES OF Y
C   CHISQR - REDUCED CHI SQUARE FOR FIT 
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C   FUNCTN (X, I, A)
C      EVALUATES THE FITTING FUNCTION FOR THE ITH TERM
C   FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
C      EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA
C   MATINV (ARRAY, NTERMS, DET) 
C      INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE NTERMS
C      AND CALCULATES ITS DETERMINANT
C
C COMMENTS
C   DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C
      SUBROUTINE CHIFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,
     *SIGMAA,YFIT,CHISQR)
      DOUBLE PRECISION ALPHA
      DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1), 
     *YFIT(1)
      DIMENSION ALPHA(10,10),BETA(10),DA(10)
      REAL FUNCTN
      EXTERNAL FUNCTN
C 
11    NFREE=NPTS-NTERMS
      FREE=NFREE
      IF (NFREE) 14,14,16
14    CHISQR=0.
      GOTO 120
16    DO 17 I=1,NPTS
17    YFIT(I)=FUNCTN (X,I,A)
      CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
C 
C EVALUATE ALPHA AND BETA MATRICES
C 
20    DO 60 J=1,NTERMS
C 
C A(J) + DELTAA(J)
C 
21    AJ=A(J) 
      A(J)=AJ+DELTAA(J)
      DO 24 I=1,NPTS
24    YFIT(I)=FUNCTN (X,I,A)
      CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      ALPHA(J,J)=CHISQ2-2.*CHISQ1
      BETA(J)=-CHISQ2 
31    DO 50 K=1,NTERMS
      IF (K-J) 33,50,36
33    ALPHA(K,J)=(ALPHA(K,J)-CHISQ2)/2.
      ALPHA(J,K)=ALPHA(K,J)
      GOTO 50 
36    ALPHA(J,K)=CHISQ1-CHISQ2
C 
C A(J) + DELTAA(J) AND A(K) + DELTAA(K) 
C 
41    AK=A(K) 
      A(K)=AK+DELTAA(K)
      DO 44 I=1,NPTS
44    YFIT(I)=FUNCTN (X,I,A)
      CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      ALPHA(J,K)=ALPHA(J,K)+CHISQ3
      A(K)=AK 
50    CONTINUE
C 
C A(J) - DELTAA(J)
C 
51    A(J)=AJ-DELTAA(J)
      DO 53 I=1,NPTS
53    YFIT(I)=FUNCTN (X,I,A)
      CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      A(J)=AJ 
      ALPHA(J,J)=(ALPHA(J,J)+CHISQ3)/2.
      BETA(J)=(BETA(J)+CHISQ3)/4.
60    CONTINUE
C 
C ELIMINATE NEGATIVE CURVATURE
C 
61    DO 70 J=1,NTERMS
      IF (ALPHA(J,J)) 63,65,70
63    ALPHA(J,J)=-ALPHA(J,J)
      GOTO 66 
65    ALPHA(J,J)=0.01 
C Changed from DO 70
66    DO 700 K=1,NTERMS		
C Changed from 68, 70, 68
      IF (K-J) 68,700,68
68    ALPHA(J,K)=0.
      ALPHA(K,J)=0.
C New continuation statement
700   CONTINUE
70    CONTINUE
C 
C INVERT MATRIX AND EVALUATE PARAMETER INCREMENTS
C 
71    CALL MATINV (ALPHA,NTERMS,DET)
      DO 76 J=1,NTERMS
      DA(J)=0.
74    DO 75 K=1,NTERMS
75    DA(J)=DA(J)+BETA(K)*ALPHA(J,K)
76    DA(J)=0.2*DA(J)*DELTAA(J)
C 
C MAKE SURE CHI SQUARE DECREASES
C 
81    DO 82 J=1,NTERMS
82    A(J)=A(J)+DA(J) 
83    DO 84 I=1,NPTS
84    YFIT(I)=FUNCTN (X,I,A)
      CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      IF (CHISQ1-CHISQ2) 87,91,91
87    DO 89 J=1,NTERMS
      DA(J)=DA(J)/2.
89    A(J)=A(J)-DA(J) 
      GOTO 83 
C 
C INCREMENT PARAMETERS UNTIL CHI SQUARE STARTS TO INCREASE
C 
91    DO 92 J=1,NTERMS
92    A(J)=A(J)+DA(J) 
      DO 94 I=1,NPTS
94    YFIT(I)=FUNCTN (X,I,A)
      CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      IF (CHISQ3-CHISQ2) 97,101,101
97    CHISQ1=CHISQ2
      CHISQ2=CHISQ3
99    GOTO 91 
C 
C FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS 
C 
101   DELTA=1./(1.+(CHISQ1-CHISQ2)/(CHISQ3-CHISQ2))+0.5
      DO 104 J=1,NTERMS
      A(J)=A(J)-DELTA*DA(J)
104   SIGMAA(J)=DELTAA(J)*SQRT(FREE*ALPHA(J,J))
      DO 106 I=1,NPTS 
106   YFIT(I)=FUNCTN (X,I,A)
      CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
111   IF (CHISQ2-CHISQR) 112,120,120
112   DO 113 J=1,NTERMS
113   A(J)=A(J)+(DELTA-1.)*DA(J)
      DO 115 I=1,NPTS 
115   YFIT(I)=FUNCTN (X,I,A)
      CHISQR=CHISQ2
120   RETURN
      END