SUBROUTINE	DRGSUM(X,NORDER,Y,WEIGHT,XSUM,YSUM,Y2SUM)
C
C	SUBROUTINE TO ACCUMULATE SUMS FOR A LEAST SQUARES FIT TO:
C
C		Y= A(1)*X(1) + A(2)*X(2) . . . + A(NORDER)*X(NORDER)
C			WHERE X,Y ARE GIVEN A'S ARE TO BE FOUND
C
C	X	IS THE ARRAY OF LINEAR FUNCIONS.  THE FIRST VALUE OF X MUST BE
C		1.0.  X(1)	= 1.0
C
C	Y	IS THE ACTUAL Y VALUE CORRESPONDING TO THE X'S
C
C	WEIGHT	IS THE WEIGHTING FACTOR = 1./SIGY**2
C		WHERE SIGY	= THE Y STANDARD DEVIATION
C
C	XSUM,YSUM	ARE ARRAYS ACCUMULATED FOR THE LATER FIT.
C		XSUM IS A 2 DIM. ARRAY, AND YSUM 1 DIMENSION.
C
C	THE FINAL ANSWER IS EVALUATED BY SUBROUTINE DRGVAL
C
	IMPLICIT	REAL*8(A-H,O-Z)
	DIMENSION	X(NORDER),YSUM(NORDER),XSUM(NORDER,NORDER)
	Y2SUM	= Y2SUM + WEIGHT*Y*Y
	DO 100 I	= 1,NORDER
	XA	= WEIGHT*X(I)
	YSUM(I)	= YSUM(I)  + Y*XA
	DO 100 J	= 1,I
	XSUM(I,J)	= XSUM(I,J) + XA*X(J)
100	CONTINUE
	END

	SUBROUTINE	DRGVAL(NORDER,PTS,XSUM,YSUM,Y2SUM,ITERM,
	1			A,SIGA,CHISQR,FTEST,B)
C
C	THIS SUBROUTINE TAKES THE ARRAYS XSUM,YSUM GENERATED BY REGSUM
C	AND THE NUMBER OF DATA POINTS (PTS) AND FINISHES THE LEAST
C	SQUARE FIT TO A LINEAR COMBINATION OF FUNCTIONS. IT CALCULATES THE
C	COEFFICIENTS (ARRAY A(NORDER) ), THE ERRORS
C	ON THESE COEFFICIENTS(SIGA), THE REDUCED CHI SQUARED/DEGREE-OF-FREEDOM
C	AND FTEST (SEE BEVINGTON).
C		THE SCRATCH 2 DIMENSIONAL ARRAY B IS USED BY THIS PROGRAM.
C
	IMPLICIT	REAL*8(A-H,O-Z)
	DIMENSION	XSUM(NORDER,NORDER),YSUM(NORDER),A(NORDER)
	DIMENSION	SIGA(NORDER),ITERM(NORDER),B(NORDER,NORDER)
C
C	SYMMETRIZE THE ARRAY A
C
	CHISQR	= -1.
	FTEST	= 0.
	IF(NORDER .LT. 1) RETURN
	IF(XSUM(1,1) .EQ. 0.) RETURN
	IF(NORDER .LT. 2) GO TO 11
	DO 10 I	= 1,NORDER-1
	DO 10 J	= I+1,NORDER
10	XSUM(I,J)	= XSUM(J,I)
11	DO 12 I	= 1,NORDER
	DO 12 J	= 1,NORDER
12	B(I,J)	= XSUM(I,J)
	DO 13 I	= 1,NORDER
	DO 13 J	= 1,NORDER
	IF(ITERM(I) * ITERM(J) .NE. 0) GO TO 13
	B(I,J)	= 0.
	IF(I .EQ. J) B(I,J)	= 1
13	CONTINUE
	CALL	DMATNV(B,NORDER,DET)
	NTERMS	= 0
	DO 14 I	= 1,NORDER
	IF( ITERM(I) .NE. 0 ) GO TO 14
	NTERMS	= NTERMS + 1
	B(I,I)	= 0.
14	CONTINUE
	NTERMS	= NORDER - NTERMS
	IF(DET .EQ. 0.) RETURN			!BAD MATRIX INVERSION
	CALL	DMATML(NORDER,B,NORDER,YSUM,1,A)	!THE RESULT

C
C	CALCULATE THE REDUCED CHI SQUARED/POINT AND THE ERROR
C
	FREEN	= PTS - NTERMS - 1
	SUM	= Y2SUM
	DO 30 I	= 1,NORDER
	SIGA(I)	= SQRT( ABS(B(I,I)) )		!THE ERROR
	IF(ITERM(I) .EQ. 0) GO TO 30
	SUM	= SUM - 2.*A(I)*YSUM(I)
	DO 30 J	= 1,NORDER
	IF(ITERM(J) .EQ. 0) GO TO 30
	SUM	= SUM + A(I)*A(J)*XSUM(I,J)
30	CONTINUE
	IF(FREEN .LT. .5) RETURN
	CHISQR	= SUM/FREEN
	R	= Y2SUM - YSUM(1)**2/SUM
	IF(R .EQ. 0.) RETURN
	R	= SUM/R
	FTEST	= (1.-R) * FREEN/(R * NTERMS)
	END

	SUBROUTINE	DMATML(IX,A,IY,B,IZ,C)
C
C	THIS SUBROUTINE MULTIPLIES 2 MATRICES A,B TOGETHER TO FORM C
C
	IMPLICIT	REAL*8(A-H,O-Z)
	DIMENSION	A(IX,IY),B(IY,IZ),C(IX,IZ)
	DO 10	I	= 1,IX
	DO 10 J	= 1,IZ
	C(I,J)	= 0.
	DO 10 K	= 1,IY
10	C(I,J)	= C(I,J) + A(I,K)*B(K,J)
	END