SUBROUTINE	ANGAD(SIN1,COS1,SIN2,COS2,SIN3,COS3,TAN3)
C
C         CALL ANGAD(SIN1,COS1,SIN2,COS2,SIN3,COS3,TAN3)
C         PRODUCES SIN,COS,TAN OF ANGLE WHERE ANG3 = ANG1+ANG2
C         NOTE RESULT MAY BE STORED IN SOURCEC
C
	SIN3	= (SIN1*COS2 + SIN2*COS1)
	COS3	= (COS1*COS2 - SIN1*SIN2)
	IF(COS3 .NE. 0) GO TO 10
	TAN3	= 1.E30
	RETURN
10	CONTINUE
	TAN3	= SIN3/COS3
	RETURN
	END
	SUBROUTINE	ASKN(STRING,IVAL,MIN,MAX,IDEF)
	BYTE	IA(9),STRING
	CALL	PARAMS(IPAR)
1	IF(IPAR .LT. 5) GO TO 20		! No default available
	ENCODE(8,1000,IA) IDEF
	IA(9)	= 0
1000	FORMAT(I8)
	CALL	ASKS(STRING,IA,ISIZE,0,8,IA)
	IVAL	= IDEF
	IF(ISIZE .EQ. 0) RETURN
10	IA(ISIZE + 1)	= ','
	DECODE(8,1000,IA) IVAL
	IF(IPAR .LT. 3) RETURN		! No max or min specified
	J	= -1
	I	= NULTST(MIN)			! Non zero if no MIN
	IF(IPAR .LT. 4) GO TO 13 		! No max
	J	= NULTST(MAX)			! Non zero if no MAX
	IF(J .NE. 0) GO TO 13
	IF(I .NE. 0) GO TO 12			! No MIN available
	IF(MIN .GT. MAX) RETURN			! MIN and MAX do not agree
12	IF(IVAL .LE. MAX) GO TO 13		! Value is not too big
	WRITE(5,1002) MAX
1002	FORMAT(' *** ERROR *** Value must not be bigger than' I8)
	GO TO 1
13	IF(I .NE. 0) RETURN			! No MAX
	IF(IVAL .GE. MIN) RETURN		! Value not too small
	WRITE(5,1003) MIN
1003	FORMAT(' *** ERROR *** Value must not be smaller than' I8)
	GO TO 1
20	CALL	ASKS(STRING,IA,ISIZE,0,8)
	GO TO 10
	END

	SUBROUTINE	ASKF(STRING,FVAL,FMIN,FMAX,FDEF)
	BYTE	IA(16),STRING
	CALL	PARAMS(IPAR)
1	IF(IPAR .LT. 5) GO TO 20		! No default available
	ENCODE(15,1000,IA) FDEF
	IA(16)	= 0
1000	FORMAT(1P,G15.7)
1001	FORMAT(G15.0)
	CALL	ASKS(STRING,IA,ISIZE,0,15,IA)
	FVAL	= FDEF
	IF(ISIZE .EQ. 0) RETURN
10	IA(ISIZE + 1)	= ','
	DECODE(15,1001,IA) FVAL
	IF(IPAR .LT. 3) RETURN		! No max or min specified
	J	= -1
	I	= NULTST(FMIN)			! Non zero if no FMIN
	IF(IPAR .LT. 4) GO TO 13 		! No max
	J	= NULTST(FMAX)			! Non zero if no FMAX
	IF(J .NE. 0) GO TO 13
	IF(I .NE. 0) GO TO 12			! No FMIN available
	IF(FMIN .GT. FMAX) RETURN		! FMIN and FMAX do not agree
12	IF(FVAL .LE. FMAX) GO TO 13		! Value is not too big
	WRITE(5,1002) FMAX
1002	FORMAT(' *** ERROR *** Value must not be bigger than' 1P,G15.7)
	GO TO 1
13	IF(I .NE. 0) RETURN			! No FMAX
	IF(FVAL .GE. FMIN) RETURN		! Value not too small
	WRITE(5,1003) FMIN
1003	FORMAT(' *** ERROR *** Value must not be smaller than' 1P,G15.7)
	GO TO 1
20	CALL	ASKS(STRING,IA,ISIZE,0,15)
	GO TO 10
	END

	SUBROUTINE	ASK2(STRING,VAL,MIN,MAX,DEF)
	INTEGER*4 VAL,MIN,MAX,DEF
	BYTE	IA(16),STRING
	CALL	PARAMS(IPAR)
1	IF(IPAR .LT. 5) GO TO 20		! No default available
	ENCODE(15,1000,IA) DEF
	IA(16)	= 0
1000	FORMAT(I15)
	CALL	ASKS(STRING,IA,ISIZE,0,15,IA)
	VAL	= DEF
	IF(ISIZE .EQ. 0) RETURN
10	IA(ISIZE + 1)	= ','
	DECODE(15,1000,IA) VAL
	IF(IPAR .LT. 3) RETURN		! No max or min specified
	J	= -1
	I	= NULTST(MIN)			! Non zero if no MIN
	IF(IPAR .LT. 4) GO TO 13 		! No max
	J	= NULTST(MAX)			! Non zero if no MAX
	IF(J .NE. 0) GO TO 13
	IF(I .NE. 0) GO TO 12			! No MIN available
	IF(MIN .GT. MAX) RETURN			! MIN and MAX do not agree
12	IF(VAL .LE. MAX) GO TO 13		! Value is not too big
	WRITE(5,1002) MAX
1002	FORMAT(' *** ERROR *** Value must not be bigger than' I15)
	GO TO 1
13	IF(I .NE. 0) RETURN			! No MAX
	IF(VAL .GE. MIN) RETURN			! Value not too small
	WRITE(5,1003) MIN
1003	FORMAT(' *** ERROR *** Value must not be smaller than' I15)
	GO TO 1
20	CALL	ASKS(STRING,IA,ISIZE,0,15)
	GO TO 10
	END

	FUNCTION DDETRM(ARRAY,NORDER)
C
C	FUNCTION TO FIND DETERMINANT OF DOUBLE PRECISION ARRAY(NORDER,NORDER)
C
	REAL*8 ARRAY(NORDER,NORDER),SAVE,DETSAV
	DETSAV	= 1.			!INITIAL VALUE
	DO 50 K	= 1,NORDER		!DO ALL ROWS
	SAVE	= 1.E-30
	ISAVE	= 0
	DO 23 J	= K,NORDER		!FIND BIGGEST ELEMENT IN ROW
	IF(ABS(ARRAY(K,J)) .LT. SAVE)GO TO 23
	SAVE	= ABS(ARRAY(K,J))
	ISAVE	= J			!COLUMN NUMBER OF BIGGEST
23	CONTINUE
	IF(ISAVE .GT. 0)GO TO 31	!AOK IF FOUND BIGGEST
	DDETRM	= 0.			!NO LARGE ELEMENT, DETERMINANT=0
	RETURN				!NULL DETERMINANT
31	IF(ISAVE .EQ. K)GO TO 41	!IF ROW# = COLUMN DO NOT SWAP
	CALL SWAP(ARRAY(1,ISAVE),ARRAY(1,K),NORDER*4)
	DETSAV	= - DETSAV
41	DETSAV	= DETSAV * ARRAY(K,K)
	IF(K .GE. NORDER) GO TO 50
	KI	= K 
	KI	= K + 1
	DO 46 I	= KI,NORDER		!DO ELEMENTARY ROW MANIPULATION
	DO 46 J	= KI,NORDER		!AND COLUMN MANIPULATION
C		PRODUCE ALL ZERO ELEMENTS IN ROW + COLUMN ON DIAGONAL
46	ARRAY(I,J)	= ARRAY(I,J) - ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
50	CONTINUE
	DDETRM	= DETSAV
	RETURN				!ANSWER IS NON ZERO
	END
	FUNCTION DETERM(ARRAY,NORDER)
C
C	FUNCTION TO FIND DETERMINANT OF DOUBLE PRECISION ARRAY(NORDER,NORDER)
C
	REAL*4 ARRAY(NORDER,NORDER)
	DETERM	= 1.			!INITIAL VALUE
	DO 50 K	= 1,NORDER		!DO ALL ROWS
	SAVE	= 1.E-30
	ISAVE	= 0
	DO 23 J	= K,NORDER		!FIND BIGGEST ELEMENT IN ROW
	IF(ABS(ARRAY(K,J)) .LT. SAVE)GO TO 23
	SAVE	= ABS(ARRAY(K,J))
	ISAVE	= J			!COLUMN NUMBER OF BIGGEST
23	CONTINUE
	IF(ISAVE .GT. 0)GO TO 31	!AOK IF FOUND BIGGEST
	DETERM	= 0.			!NO LARGE ELEMENT, DETERMINANT=0
	RETURN				!NULL DETERMINANT
31	IF(ISAVE .EQ. K)GO TO 41	!IF ROW# = COLUMN DO NOT SWAP
	CALL SWAP(ARRAY(1,ISAVE),ARRAY(1,K),NORDER*2)
	DETERM	= - DETERM
41	DETERM	= DETERM * ARRAY(K,K)
	IF(K .GE. NORDER) GO TO 50
	KI	= K 
	KI	= K + 1
	DO 46 I	= KI,NORDER		!DO ELEMENTARY ROW MANIPULATION
	DO 46 J	= KI,NORDER		!AND COLUMN MANIPULATION
C		PRODUCE ALL ZERO ELEMENTS IN ROW + COLUMN ON DIAGONAL
46	ARRAY(I,J)	= ARRAY(I,J) - ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
50	CONTINUE
	RETURN				!ANSWER IS NON ZERO
	END
	SUBROUTINE DMATNV(ARRAY,N,DET)
C  
C  MATRIX INVERSION ROUTINE--BEVINGTON,P.302
C
C  ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED
C  N IS THE ORDER, DET IS THE CALCULATED DETERMINANT
C  THE INVERTED MATRIX REPLACES THE ORIGINAL ONE
C
	IMPLICIT	REAL*8 (A-H, O-Z)
	REAL	ARRAY(N,N)
	DIMENSION IK(10),JK(10)
	DET	= 1.
	DO 100 K	= 1,N
C
C  FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL
C  PARTIAL PIVOTING
C
	AMAX	= 0.0
	DO 30 I	= K,N
	DO 30 J	= K,N
	IF(ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30
24	AMAX	= ARRAY(I,J)
	IK(K)	= I
	JK(K)	= J
30	CONTINUE
	IF(AMAX)41,32,41
32	DET	= 0.
	GO TO 140
41	I	= IK(K)
	IF(I-K) 21,51,43
43	DO 50 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= ARRAY(I,J)
50	ARRAY(I,J)	= -SAVE
51	J	= JK(K)
	IF(J-K)21,61,53
53	DO 60 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= ARRAY(I,J)
60	ARRAY(I,J)	= -SAVE
C
C  INVERT MATRIX
C
61	DO 70 I	= 1,N
	IF(I-K)63,70,63
63	ARRAY(I,K)	= -ARRAY(I,K)/AMAX
70	CONTINUE
71	DO 80 I	= 1,N
	DO 80 J	= 1,N
	IF(I-K)74,80,74
74	IF(J-K)75,80,75
75	ARRAY(I,J)	= ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80	CONTINUE
81	DO 90 J	= 1,N
	IF(J-K)83,90,83
83	ARRAY(K,J)	= ARRAY(K,J)/AMAX
90	CONTINUE
	ARRAY(K,K)	= 1./AMAX
100	DET	= DET*AMAX
C
C  RESTORE ORDERING OF MATRIX
C
101	DO 130 L	= 1,N
	K	= N-L+1
	J	= IK(K)
	IF(J-K)111,111,105
105	DO 110 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= -ARRAY(I,J)
110	ARRAY(I,J)	= SAVE
111	I	= JK(K)
	IF(I-K)130,130,113
113	DO 120 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= -ARRAY(I,J)
120	ARRAY(I,J)	= SAVE
130	CONTINUE
140	RETURN
	END
	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
	PROGRAM FITIT
C
C	PROGRAM TO FIT DATA POINTS
C	FIT ROUTINE LEGFIT IS TAKEN FROM BEVINGTON'S TEXT (PAGE 155).
C
	BYTE	ITYP(10),ITYPIN(10),IA(80),TITLE(40)
	DOUBLE PRECISION	A,SIGA,CHISQR,FTEST,B,XX,OFFSET
	DOUBLE PRECISION	XSUM,YSUM,Y2SUM,PTS,XF,YF,WEIGHT
	DIMENSION	X(200),Y(200),S(200),A(20),SIGA(20)
	DIMENSION	XFIT(200),YFIT(200)
	DIMENSION XSUM(20,20),B(20,20),YSUM(20),XF(20),ITERM(20)
	COMMON	/SCRAT/ B
	DATA	ITYP /'P','S','C','L',6*0/
	DATA	NMAX/ 20/
1	CALL	ASKS('Name of file with data to fit ',IA,ISIZE,0,80,
	1 'SY:FITIT.DAT')
	CALL	ASSIGN(1,IA,ISIZE)
C
C	 READ IN WEIGHTING MODE AND DATA POINTS.
C	IF MODE	=
C		  -1 (STATISTICAL WEIGHTING--WEIGHTING FACTOR	= 1/Y)
C		   0 (NO WEIGHTING--WEIGHTING FACTOR	= 1)
C		  +1 (INSTRUMENTAL WEIGHTING--WEIGHTING FACTOR	= 1/S**2)
C
	READ(1,1007) TITLE
	READ(1,1007)
	READ(1,1007)
	READ(1,1000,ERR=1,END=1) NPTS,ICHAR,MODE
	READ(1,1007)
	IF(NPTS .GT. 200) STOP '*** ERROR *** ONLY 200 DATA POINTS ALLOWED'
	IF (MODE .GT. 0) READ(1,*) (X(J),Y(J),S(J), J=   1,NPTS)
	IF (MODE .LE. 0) READ(1,*) (X(J),Y(J), J=  1,NPTS)
	WRITE(5,1005) TITLE
C
C	SET UP TO GRAPH RESULTS.
C
	XMAX 	= X(1)
	XMIN 	= XMAX
	YMAX	= Y(1)
	YMIN	= YMAX
	DO 10  J	= 1, NPTS
	IF(MODE .EQ. 0) S(J)	= 1.
	IF(MODE .LT. 0) S(J)	= SSQRT( Y(J) )	!STATISTICAL WEIGHT
	IF(S(J) .LE. 0.) STOP '*** ERROR *** BAD ERRORS ON Y'
	IF (XMAX .LT. X(J)) XMAX	= X(J)
	IF (XMIN .GT. X(J)) XMIN	= X(J)
	IF(YMIN .GT. Y(J) ) YMIN	= Y(J)
	IF(YMAX .LT. Y(J) ) YMAX	= Y(J)
10	CONTINUE
	XSTEP 	= (XMAX - XMIN)/200.
	NTERMS	= 0
11	WRITE(5,1008)
	CALL	ASKS('Type of fit',ITYPIN,ISIZE,0,10)
	DO 15 J	= 1,10
	IF( ITYPIN(1) .NE. ITYP(J) ) GO TO 15
	ITYPE	= J
	GO TO 16
15	CONTINUE
	CALL	CLOSE(1)
	GO TO 1					! NOW NEW INPUT FILE
16	CONTINUE
	SCALE	= 1.
	OFFSET	= 0.
	IF(ITYPE .NE. 1) GO TO 21		! Do not rescale non polynomial
	CALL	ASK('Rescale data ',ITEST)
	IF(ITEST .EQ. 0) GO TO 21
	SCALE	= 5./(ABS(XMAX) + ABS(XMIN) )
	OFFSET	= (XMAX + XMIN) / 2.
	GO TO 18
17	CONTINUE
	OFFSET	= XMIN
	SCALE	= 360./(XMAX - XMIN)
	IF(ITYPE .EQ. 4) SCALE	= SCALE/2.
18	CONTINUE
	CALL	FTRUNC(2,SCALE)
	CALL	FTRUNC(2,OFFSET)

C
C	START THE FIT BY INIT VARIABLES
C
21	CONTINUE
 	CALL	ERASE
22	DO 25 I	= 1,NMAX
	YSUM(I)	= 0.
	A(I)	= 0.
	SIGA(I)	= 0.
	DO 25 J	= 1,NMAX
25	XSUM(I,J)	= 0.
	Y2SUM	= 0.

	CALL	ASKN('Number of terms to fit',NTERMS,0,NMAX)
	IF(NTERMS .LE. 0) GO TO 11
	XF(1)	= 1.0
	WRITE (5,1013)
	CALL	ASKN('Terms to fit ',NEVEN)
	DO 30 J	= 1,NTERMS
	ITERM(J)	= 1
30	CONTINUE
	IF(NEVEN .EQ. 0) GO TO 35
	MIN	= 3
	IF(NEVEN .LT. 0) MIN	= 2
	DO 31 J	= MIN,NTERMS,2
31	ITERM(J)	= 0
35	CONTINUE
C 
C	DO FITTING.
C
	DO 400 I	= 1,NPTS
	XX	= X(I)
	XX	= (XX - OFFSET) * SCALE
	CALL	XFUNC(XF,NTERMS,XX,ITYPE)	!GET THE SERIES
	WEIGHT	= 1./S(I)**2
	YF	= Y(I)
400	CALL	DRGSUM(XF,NTERMS,YF,WEIGHT,XSUM,YSUM,Y2SUM)
	PTS	= NPTS
	CALL	DRGVAL(NTERMS,PTS,XSUM,YSUM,Y2SUM,ITERM,
	1		A,SIGA,CHISQR,FTEST,B)
	DO 450 I	= 1,200			!CALCULATE THE FITTED CURVE
	XX	= XMIN + XSTEP*(I-1)	!THE X VALUE OF THE CURVE
	XFIT(I)	= XX
	XX	= (XX - OFFSET) * SCALE
	CALL	XFUNC(XF,NTERMS,XX,ITYPE)
	YFIT(I)	= 0.
	DO 450 J	= 1,NTERMS
	YFIT(I)	= YFIT(I) + A(J) * XF(J)	!THE Y VALUE OF CURVE
450	CONTINUE

C
C	WRITE OUT RESULTS.
C
	WRITE (5,1001) CHISQR,FTEST,OFFSET,SCALE,
	1		(A(J),SIGA(J), J=1, NTERMS)
510	CALL	ASK('Plot result',ITEST)
	IF(ITEST .EQ. 0) GO TO 22
C
C	PLOT THE RESULTS ON THE ADM-3A
C
500	ICHAR  =  -4
	IF (MODE .NE. 0) ICHAR = -104	!PLOT WITH ERROR BARS
	CALL	SSET		! Attach the terminal
	CALL	DISCKP
	CALL	QPLOT(X,Y,S,NPTS,-1,ICHAR)	!PLOT RAW DATA POINTS
	IF(NTERMS .LE. 0) GO TO 600
	ENCODE(60,1003,IA)CHISQR
	CALL	CHAR(0,760,TITLE,40)
	CALL	CHAR(256,0,IA,60)		!WRITE OUT THE CHISQR
	CALL	QPLOT(XFIT,YFIT,S,200,1,0)	!PLOT THE FIT
	CALL	TIDET				! Detach terminal
	CALL	ENACKP

C
C	SECTION TO OUTPUT RESULT
C
600	CALL	ASK('Output ',ITEST)
	IF (ITEST .EQ. 0) GO TO 21
	CALL	PLTOLP			!OUTPUT PLOT TO LP:
	CALL	SSET
	IF( NTERMS .LE. 0 ) GO TO 21
	WRITE (6,1001) CHISQR,FTEST,OFFSET,SCALE,
	1		(A(J),SIGA(J), J=1, NTERMS)
620	CALL	ASK('Append fit to input file',ITEST)
	IF(ITEST .EQ. 0) GO TO 21
	J	= 200
	WRITE(1,1016)
	WRITE(1,1003) CHISQR
625	WRITE(1,*) (XFIT(J),YFIT(J),J=1,200)
	GO TO 21
C
C	SECTION TO END PROGRAM
C
900	CLOSE(UNIT=6,DISPOSE='PRINT')
	CALL	DISCKP
	CALL	SSET
	CALL	ERASE
	CALL	EXIT

1000	FORMAT(8I10)
1001	FORMAT('  CHI SQUARED/DEG-OF-FREEDOM =	'1P,G12.3
	1	/'  FTEST =			'G12.3
	2	/'  OFFSET ='G15.7'	SCALE ='G15.7
	2	/' COEFFICIENTS:	ERRORS: '
	3	/(1P,G14.6,'	',G12.2) )
1002	FORMAT(1P,G13.6)	
1003	FORMAT('POLYNOMIAL FIT   CHI SQ/FREEDOM= '1P,G10.3)
1005	FORMAT(' 		Title of dataset:'/ 1X,60A1 //)
1007	FORMAT(80A1)
1008	FORMAT(' P	= POLYNOMIAL'/' S	= SINE'/' C	= COSINE'
	1	/' L	= LEGENDRE POLY.'//)
1013	FORMAT( ' 	+1 - Fit to EVEN terms ONLY'/
	2 ' 	 0 - Fit to ALL terms'/
	3 ' 	-1 - Fit to ODD terms ONLY'/ )
1016	FORMAT(' 200,0,0')
	END

	SUBROUTINE	XFUNC (XF,NTERMS,XX,ITYPE)
	DOUBLE PRECISION	XF(NTERMS),XX,RAD,TEMP,FL
	DATA	RAD /.017453292519943296/	!CONVERTS DEGREES TO RADIANS
	DATA	PI  /3.141592653589793/
	XF(1)	= 1.0
	IF(NTERMS .LE. 1) RETURN
	GO TO (110,120,130,140,150) ITYPE
110	TEMP	= XX
	DO 115 J	= 2,NTERMS
	XF(J)	= TEMP
115	TEMP	= TEMP*XX
	RETURN
120	CONTINUE
	DO 125 J	= 2,NTERMS
	XF(J)	= SIN(XX*RAD*(J-1))
125	CONTINUE
	RETURN
130	CONTINUE
	DO 135 J	= 2,NTERMS
135	XF(J)	= COS(XX*RAD*(J-1))
	RETURN
140	CONTINUE
	TEMP	= COS(XX*RAD)		!FIND TEMP OF VARIABLE
	XF(2)	= TEMP
	DO 145 J	= 3,NTERMS
	FL	= J-1
145	XF(J)	= ( (2.*FL-1)*TEMP*XF(J-1) - (FL-1)*XF(J-2) )/FL
	RETURN
150	CONTINUE
	END
	SUBROUTINE	FROUND(IDIGIT,F)
	INTEGER*4 ITEST,JTEST
	DOUBLE PRECISION VALUE,SCALE
	IF(IDIGIT .GT. 6) RETURN
	IF( F .EQ. 0) RETURN
	VALUE	= ABS(F)
	SCALE	= LOG10(VALUE)
	IF(SCALE .LT. 0.) SCALE	= SCALE - 1.
	IEXP	= SCALE
	IEXP	= IEXP - IDIGIT
	SCALE	= 10.
	SCALE	= SCALE**IEXP
	ITEST	= VALUE/SCALE + .5
	JTEST	= 10
	JTEST	= JTEST**IDIGIT
	IF(ITEST .LT. JTEST)  GO TO 20
	ITEST	= ITEST/10
	SCALE	= SCALE*10.
20	CONTINUE
	IF( F .LT. 0) ITEST	= -ITEST
	F	= ITEST*SCALE
	END
	SUBROUTINE	FTRUNC(IDIGIT,F)
	INTEGER*4 ITEST,JTEST
	DOUBLE PRECISION VALUE,SCALE
	IF(IDIGIT .GT. 6) RETURN
	IF( F .EQ. 0) RETURN
	VALUE	= ABS(F)
	SCALE	= LOG10(VALUE)
	IF(SCALE .LT. 0.) SCALE	= SCALE - 1.
	IEXP	= SCALE
	IEXP	= IEXP - IDIGIT
	SCALE	= 10.
	SCALE	= SCALE**IEXP
	ITEST	= VALUE/SCALE
	JTEST	= 10
	JTEST	= JTEST**IDIGIT
	IF(ITEST .LT. JTEST)  GO TO 20
	ITEST	= ITEST/10
	SCALE	= SCALE*10.
20	CONTINUE
	IF( F .LT. 0) ITEST	= -ITEST
	F	= ITEST*SCALE
	END
	SUBROUTINE MATINV(ARRAY,N,DET)
C  
C  MATRIX INVERSION ROUTINE--BEVINGTON,P.302
C
C  ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED
C  N IS THE ORDER, DET IS THE CALCULATED DETERMINANT
C  THE INVERTED MATRIX REPLACES THE ORIGINAL ONE
C
	DIMENSION	ARRAY(N,N)
	DIMENSION IK(20),JK(20)
	DET	= 1.
	DO 100 K	= 1,N
C
C  FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL
C  PARTIAL PIVOTING
C
	AMAX	= 0.0
21	DO 30 I	= K,N
	DO 30 J	= K,N
	IF(ABS(AMAX) .GT. ABS(ARRAY(I,J))) GO TO 30
	AMAX	= ARRAY(I,J)
	IK(K)	= I
	JK(K)	= J
30	CONTINUE
	IF(AMAX .NE. 0.) GO TO 41
	DET	= 0.
	GO TO 140
41	I	= IK(K)
	IF(I-K) 21,51,43
43	DO 50 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= ARRAY(I,J)
50	ARRAY(I,J)	= -SAVE
51	J	= JK(K)
	IF(J-K)21,61,53
53	DO 60 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= ARRAY(I,J)
60	ARRAY(I,J)	= -SAVE
C
C  INVERT MATRIX
C
61	DO 70 I	= 1,N
	IF(I .EQ. K) GO TO 70
	ARRAY(I,K)	= -ARRAY(I,K)/AMAX
70	CONTINUE
71	DO 80 I	= 1,N
	DO 80 J	= 1,N
	IF(I .EQ. K) GO TO 80
	IF(J .EQ. K) GO TO 80
	ARRAY(I,J)	= ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80	CONTINUE
81	DO 90 J	= 1,N
	IF(J .EQ. K) GO TO 90
	ARRAY(K,J)	= ARRAY(K,J)/AMAX
90	CONTINUE
	ARRAY(K,K)	= 1./AMAX
100	DET	= DET*AMAX
C
C  RESTORE ORDERING OF MATRIX
C
101	DO 130 L	= 1,N
	K	= N-L+1
	J	= IK(K)
	IF(J .LE. K) GO TO 111
	DO 110 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= -ARRAY(I,J)
110	ARRAY(I,J)	= SAVE
111	I	= JK(K)
	IF(I .LE. K) GO TO 130
	DO 120 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= -ARRAY(I,J)
120	ARRAY(I,J)	= SAVE
130	CONTINUE
140	RETURN
	END
	SUBROUTINE	 MMULT(AMAT,BMAT,CMAT,IROWA,ICOLA,ICOLB)
C
C	CALL MMULT(AMAT,BMAT,CMAT,IROWA,ICOLA,ICOLB)
C
C	CMAT	= AMAT X BMAT
C
C	MATRIX MULTIPLICATION
	REAL	AMAT(IROWA,ICOLA),BMAT(ICOLA,ICOLB),CMAT(IROWA,ICOLB)
	DO 10 I	= 1,IROWA
	DO 10 J	= 1,ICOLB
	CMAT(I,J)	= 0
	DO 10 K = 1,ICOLA
	CMAT(I,J)	= CMAT(I,J) + AMAT(I,K)*BMAT(K,J)
10	CONTINUE
	END
	SUBROUTINE POLFIT(X,Y,S,NPTS,NTERM,MODE,A,CHISQR)
	REAL X(NPTS),Y(NPTS),S(NPTS),A(NTERM)
	REAL*8 SUMX(19),SUMY(10),ARRAY(100)
	REAL*8 XTERM,YTERM,CHISQ
	CALL ZERO(A,NTERM*2)
	CHISQR	= 0.
	NTERMS	= NTERM
	IF(NTERMS .GT. 10)NTERMS	= 10
	IF(NTERMS .GE. NPTS)NTERMS	= NPTS - 1
	NMAX	= 2*NTERMS - 1
	CALL ZERO(SUMX,76)
	CALL ZERO(SUMY,40)
	CHISQ	= 0.
	XSAV	= 0.
	DO 20 I	= 1,NPTS
	IF(XSAV .GT. ABS(X(I)) )GO TO 20
	XSAV	= ABS(X(I))
20	CONTINUE
	XSAV	= XSAV/1.5
	DO 50 I	= 1,NPTS
	XI	= X(I)/XSAV
	YI	= Y(I)
	IF(MODE) 32,37,39
32	WEIGHT	= 1./(.5 + ABS(YI) )
	GO TO 41
37	WEIGHT	= 1.
	GO TO 41
39	WEIGHT	= 1./S(I)**2
41	XTERM	= WEIGHT
	DO 44 N	= 1,NMAX
	SUMX(N)	= SUMX(N) + XTERM
44	XTERM	= XTERM * XI
	YTERM	= WEIGHT* YI
	YT	= WEIGHT* XT
	DO 48 N	= 1,NTERMS
	SUMY(N)	= SUMY(N) + YTERM
48	YTERM	= YTERM * XI
49	CHISQ	= CHISQ + WEIGHT*YI**2
50	CONTINUE
C
C	CONSTRUCT MATRICES + CALCULATE COEFFICIENTS
C
	DO 54 J	= 1,NTERMS
	DO 54 K	= 0,NTERMS-1
	N	= J + K
	M	= J + NTERMS*K
54	ARRAY(M)	= SUMX(N)
	DELTA	= DDETRM(ARRAY,NTERMS)
	IF(ABS(DELTA) .LT. 1.E-30)RETURN
61	DO 70 L	= 1,NTERMS
	DO 66 J	= 1,NTERMS
	DO 65 K	= 0,NTERMS-1
	N	= J + K
	M	= J + NTERMS*K
65	ARRAY(M)	= SUMX(N)
	M	= J + NTERMS*(L-1)
66	ARRAY(M)	= SUMY(J)
70	A(L)	= DDETRM(ARRAY,NTERMS) /( XSAV**(L-1) * DELTA )
C
C	CALCULATE CHI SQUARE
C
	DO 75 J	= 1,NTERMS
	CHISQ	= CHISQ - 2.*A(J)*SUMY(J)*XSAV**(J-1)
	DO 75 K	= 1,NTERMS
	N	= J + K - 1
75	CHISQ	= CHISQ + A(J)*A(K)*SUMX(N)*XSAV**(J+K-2)
	CHISQR	= CHISQ / (NPTS - NTERMS)
	RETURN
	END

	SUBROUTINE	REGSUM(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 REGVAL
C
	DIMENSION	X(NORDER),YSUM(NORDER),XSUM(NORDER,NORDER)
	Y2SUM	= Y2SUM + WEIGHT*Y*Y		!SUM OF SQUARES OF Y
	DO 100 I	= 1,NORDER
	XA	= WEIGHT*X(I)
	YSUM(I)	= YSUM(I)  + Y*XA		!SUM OVER XjY
	DO 100 J	= 1,I
	XSUM(I,J)	= XSUM(I,J) + XA*X(J)	!SUM OVER XiXj
100	CONTINUE
	END

	SUBROUTINE	REGVAL(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
	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.			!DEFAULT VALUES IF BAD INPUT
	FTEST	= 0.
	IF(NORDER .LT. 1) RETURN	!BAD NUMBER OF TERMS
	IF(XSUM(1,1) .EQ. 0.) RETURN	!BAD SUM OVER WEIGHTING FUNCTION
	IF(NORDER .LT. 2) GO TO 11	!ONLY SINGLE CONSTATNT TO DETERMINE
	DO 10 I	= 1,NORDER-1
	DO 10 J	= I+1,NORDER
10	XSUM(I,J)	= XSUM(J,I)	!SYMMETRIZE THE MATRIX
11	DO 12 I	= 1,NORDER
	DO 12 J	= 1,NORDER
12	B(I,J)	= XSUM(I,J)		!MOVE MATRIX TO SCRATCH MATRIX
C
C	GET RID OF UNWANTED TERMS
C
	DO 13 I	= 1,NORDER
	DO 13 J	= 1,NORDER
	IF(ITERM(I) * ITERM(J) .NE. 0) GO TO 13
	B(I,J)	= 0.			!UNWANTED ROW/COLUMN
	IF(I .EQ. J) B(I,J)	= 1.	!SO IS UNIT MATRIX
13	CONTINUE
	CALL	MATINV(B,NORDER,DET)	!INVERT THE SCRATCH MATRIX
	NTERMS	= 0
	DO 14 I	= 1,NORDER
	IF( ITERM(I) .NE. 0 ) GO TO 14
	NTERMS	= NTERMS + 1
	B(I,I)	= 0.			!NULL VALUE FOR UNWANTED TERMS
14	CONTINUE
	NTERMS	= NORDER - NTERMS
	IF(DET .EQ. 0.) RETURN		!BAD MATRIX INVERSION
	CALL	MATMUL(NORDER,B,NORDER,YSUM,1,A)	!THE RESULT

C
C	CALCULATE THE REDUCED CHI SQUARED/POINT AND THE ERROR
C
	FREEN	= PTS - NTERMS - 1		!THE DEGREES-OF-FREEDOM
	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)		!FIND CHI SQUARED
	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	MATMUL(IX,A,IY,B,IZ,C)
C
C	THIS SUBROUTINE MULTIPLIES 2 MATRICES A,B TOGETHER TO FORM C
C
	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
      FUNCTION SSQRT(A) 
      IF(A)  10,20,20 
   10 SSQRT =-SQRT(-A)
      RETURN
   20 SSQRT = SQRT(A) 
      RETURN
      END 
	SUBROUTINE	ANGAD(SIN1,COS1,SIN2,COS2,SIN3,COS3,TAN3)
C
C         CALL ANGAD(SIN1,COS1,SIN2,COS2,SIN3,COS3,TAN3)
C         PRODUCES SIN,COS,TAN OF ANGLE WHERE ANG3 = ANG1+ANG2
C         NOTE RESULT MAY BE STORED IN SOURCEC
C
	SIN3	= (SIN1*COS2 + SIN2*COS1)
	COS3	= (COS1*COS2 - SIN1*SIN2)
	IF(COS3 .NE. 0) GO TO 10
	TAN3	= 1.E30
	RETURN
10	CONTINUE
	TAN3	= SIN3/COS3
	RETURN
	END
	SUBROUTINE	ASKN(STRING,IVAL,MIN,MAX,IDEF)
	BYTE	IA(9),STRING
	CALL	PARAMS(IPAR)
1	IF(IPAR .LT. 5) GO TO 20		! No default available
	ENCODE(8,1000,IA) IDEF
	IA(9)	= 0
1000	FORMAT(I8)
	CALL	ASKS(STRING,IA,ISIZE,0,8,IA)
	IVAL	= IDEF
	IF(ISIZE .EQ. 0) RETURN
10	IA(ISIZE + 1)	= ','
	DECODE(8,1000,IA) IVAL
	IF(IPAR .LT. 3) RETURN		! No max or min specified
	J	= -1
	I	= NULTST(MIN)			! Non zero if no MIN
	IF(IPAR .LT. 4) GO TO 13 		! No max
	J	= NULTST(MAX)			! Non zero if no MAX
	IF(J .NE. 0) GO TO 13
	IF(I .NE. 0) GO TO 12			! No MIN available
	IF(MIN .GT. MAX) RETURN			! MIN and MAX do not agree
12	IF(IVAL .LE. MAX) GO TO 13		! Value is not too big
	WRITE(5,1002) MAX
1002	FORMAT(' *** ERROR *** Value must not be bigger than' I8)
	GO TO 1
13	IF(I .NE. 0) RETURN			! No MAX
	IF(IVAL .GE. MIN) RETURN		! Value not too small
	WRITE(5,1003) MIN
1003	FORMAT(' *** ERROR *** Value must not be smaller than' I8)
	GO TO 1
20	CALL	ASKS(STRING,IA,ISIZE,0,8)
	GO TO 10
	END

	SUBROUTINE	ASKF(STRING,FVAL,FMIN,FMAX,FDEF)
	BYTE	IA(16),STRING
	CALL	PARAMS(IPAR)
1	IF(IPAR .LT. 5) GO TO 20		! No default available
	ENCODE(15,1000,IA) FDEF
	IA(16)	= 0
1000	FORMAT(1P,G15.7)
1001	FORMAT(G15.0)
	CALL	ASKS(STRING,IA,ISIZE,0,15,IA)
	FVAL	= FDEF
	IF(ISIZE .EQ. 0) RETURN
10	IA(ISIZE + 1)	= ','
	DECODE(15,1001,IA) FVAL
	IF(IPAR .LT. 3) RETURN		! No max or min specified
	J	= -1
	I	= NULTST(FMIN)			! Non zero if no FMIN
	IF(IPAR .LT. 4) GO TO 13 		! No max
	J	= NULTST(FMAX)			! Non zero if no FMAX
	IF(J .NE. 0) GO TO 13
	IF(I .NE. 0) GO TO 12			! No FMIN available
	IF(FMIN .GT. FMAX) RETURN		! FMIN and FMAX do not agree
12	IF(FVAL .LE. FMAX) GO TO 13		! Value is not too big
	WRITE(5,1002) FMAX
1002	FORMAT(' *** ERROR *** Value must not be bigger than' 1P,G15.7)
	GO TO 1
13	IF(I .NE. 0) RETURN			! No FMAX
	IF(FVAL .GE. FMIN) RETURN		! Value not too small
	WRITE(5,1003) FMIN
1003	FORMAT(' *** ERROR *** Value must not be smaller than' 1P,G15.7)
	GO TO 1
20	CALL	ASKS(STRING,IA,ISIZE,0,15)
	GO TO 10
	END

	SUBROUTINE	ASK2(STRING,VAL,MIN,MAX,DEF)
	INTEGER*4 VAL,MIN,MAX,DEF
	BYTE	IA(16),STRING
	CALL	PARAMS(IPAR)
1	IF(IPAR .LT. 5) GO TO 20		! No default available
	ENCODE(15,1000,IA) DEF
	IA(16)	= 0
1000	FORMAT(I15)
	CALL	ASKS(STRING,IA,ISIZE,0,15,IA)
	VAL	= DEF
	IF(ISIZE .EQ. 0) RETURN
10	IA(ISIZE + 1)	= ','
	DECODE(15,1000,IA) VAL
	IF(IPAR .LT. 3) RETURN		! No max or min specified
	J	= -1
	I	= NULTST(MIN)			! Non zero if no MIN
	IF(IPAR .LT. 4) GO TO 13 		! No max
	J	= NULTST(MAX)			! Non zero if no MAX
	IF(J .NE. 0) GO TO 13
	IF(I .NE. 0) GO TO 12			! No MIN available
	IF(MIN .GT. MAX) RETURN			! MIN and MAX do not agree
12	IF(VAL .LE. MAX) GO TO 13		! Value is not too big
	WRITE(5,1002) MAX
1002	FORMAT(' *** ERROR *** Value must not be bigger than' I15)
	GO TO 1
13	IF(I .NE. 0) RETURN			! No MAX
	IF(VAL .GE. MIN) RETURN			! Value not too small
	WRITE(5,1003) MIN
1003	FORMAT(' *** ERROR *** Value must not be smaller than' I15)
	GO TO 1
20	CALL	ASKS(STRING,IA,ISIZE,0,15)
	GO TO 10
	END

	FUNCTION DDETRM(ARRAY,NORDER)
C
C	FUNCTION TO FIND DETERMINANT OF DOUBLE PRECISION ARRAY(NORDER,NORDER)
C
	REAL*8 ARRAY(NORDER,NORDER),SAVE,DETSAV
	DETSAV	= 1.			!INITIAL VALUE
	DO 50 K	= 1,NORDER		!DO ALL ROWS
	SAVE	= 1.E-30
	ISAVE	= 0
	DO 23 J	= K,NORDER		!FIND BIGGEST ELEMENT IN ROW
	IF(ABS(ARRAY(K,J)) .LT. SAVE)GO TO 23
	SAVE	= ABS(ARRAY(K,J))
	ISAVE	= J			!COLUMN NUMBER OF BIGGEST
23	CONTINUE
	IF(ISAVE .GT. 0)GO TO 31	!AOK IF FOUND BIGGEST
	DDETRM	= 0.			!NO LARGE ELEMENT, DETERMINANT=0
	RETURN				!NULL DETERMINANT
31	IF(ISAVE .EQ. K)GO TO 41	!IF ROW# = COLUMN DO NOT SWAP
	CALL SWAP(ARRAY(1,ISAVE),ARRAY(1,K),NORDER*4)
	DETSAV	= - DETSAV
41	DETSAV	= DETSAV * ARRAY(K,K)
	IF(K .GE. NORDER) GO TO 50
	KI	= K 
	KI	= K + 1
	DO 46 I	= KI,NORDER		!DO ELEMENTARY ROW MANIPULATION
	DO 46 J	= KI,NORDER		!AND COLUMN MANIPULATION
C		PRODUCE ALL ZERO ELEMENTS IN ROW + COLUMN ON DIAGONAL
46	ARRAY(I,J)	= ARRAY(I,J) - ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
50	CONTINUE
	DDETRM	= DETSAV
	RETURN				!ANSWER IS NON ZERO
	END
	FUNCTION DETERM(ARRAY,NORDER)
C
C	FUNCTION TO FIND DETERMINANT OF DOUBLE PRECISION ARRAY(NORDER,NORDER)
C
	REAL*4 ARRAY(NORDER,NORDER)
	DETERM	= 1.			!INITIAL VALUE
	DO 50 K	= 1,NORDER		!DO ALL ROWS
	SAVE	= 1.E-30
	ISAVE	= 0
	DO 23 J	= K,NORDER		!FIND BIGGEST ELEMENT IN ROW
	IF(ABS(ARRAY(K,J)) .LT. SAVE)GO TO 23
	SAVE	= ABS(ARRAY(K,J))
	ISAVE	= J			!COLUMN NUMBER OF BIGGEST
23	CONTINUE
	IF(ISAVE .GT. 0)GO TO 31	!AOK IF FOUND BIGGEST
	DETERM	= 0.			!NO LARGE ELEMENT, DETERMINANT=0
	RETURN				!NULL DETERMINANT
31	IF(ISAVE .EQ. K)GO TO 41	!IF ROW# = COLUMN DO NOT SWAP
	CALL SWAP(ARRAY(1,ISAVE),ARRAY(1,K),NORDER*2)
	DETERM	= - DETERM
41	DETERM	= DETERM * ARRAY(K,K)
	IF(K .GE. NORDER) GO TO 50
	KI	= K 
	KI	= K + 1
	DO 46 I	= KI,NORDER		!DO ELEMENTARY ROW MANIPULATION
	DO 46 J	= KI,NORDER		!AND COLUMN MANIPULATION
C		PRODUCE ALL ZERO ELEMENTS IN ROW + COLUMN ON DIAGONAL
46	ARRAY(I,J)	= ARRAY(I,J) - ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
50	CONTINUE
	RETURN				!ANSWER IS NON ZERO
	END
	SUBROUTINE DMATNV(ARRAY,N,DET)
C  
C  MATRIX INVERSION ROUTINE--BEVINGTON,P.302
C
C  ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED
C  N IS THE ORDER, DET IS THE CALCULATED DETERMINANT
C  THE INVERTED MATRIX REPLACES THE ORIGINAL ONE
C
	IMPLICIT	REAL*8 (A-H, O-Z)
	REAL	ARRAY(N,N)
	DIMENSION IK(10),JK(10)
	DET	= 1.
	DO 100 K	= 1,N
C
C  FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL
C  PARTIAL PIVOTING
C
	AMAX	= 0.0
	DO 30 I	= K,N
	DO 30 J	= K,N
	IF(ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30
24	AMAX	= ARRAY(I,J)
	IK(K)	= I
	JK(K)	= J
30	CONTINUE
	IF(AMAX)41,32,41
32	DET	= 0.
	GO TO 140
41	I	= IK(K)
	IF(I-K) 21,51,43
43	DO 50 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= ARRAY(I,J)
50	ARRAY(I,J)	= -SAVE
51	J	= JK(K)
	IF(J-K)21,61,53
53	DO 60 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= ARRAY(I,J)
60	ARRAY(I,J)	= -SAVE
C
C  INVERT MATRIX
C
61	DO 70 I	= 1,N
	IF(I-K)63,70,63
63	ARRAY(I,K)	= -ARRAY(I,K)/AMAX
70	CONTINUE
71	DO 80 I	= 1,N
	DO 80 J	= 1,N
	IF(I-K)74,80,74
74	IF(J-K)75,80,75
75	ARRAY(I,J)	= ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80	CONTINUE
81	DO 90 J	= 1,N
	IF(J-K)83,90,83
83	ARRAY(K,J)	= ARRAY(K,J)/AMAX
90	CONTINUE
	ARRAY(K,K)	= 1./AMAX
100	DET	= DET*AMAX
C
C  RESTORE ORDERING OF MATRIX
C
101	DO 130 L	= 1,N
	K	= N-L+1
	J	= IK(K)
	IF(J-K)111,111,105
105	DO 110 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= -ARRAY(I,J)
110	ARRAY(I,J)	= SAVE
111	I	= JK(K)
	IF(I-K)130,130,113
113	DO 120 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= -ARRAY(I,J)
120	ARRAY(I,J)	= SAVE
130	CONTINUE
140	RETURN
	END
	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
	PROGRAM FITIT
C
C	PROGRAM TO FIT DATA POINTS
C	FIT ROUTINE LEGFIT IS TAKEN FROM BEVINGTON'S TEXT (PAGE 155).
C
	BYTE	ITYP(10),ITYPIN(10),IA(80),TITLE(40)
	DOUBLE PRECISION	A,SIGA,CHISQR,FTEST,B,XX,OFFSET
	DOUBLE PRECISION	XSUM,YSUM,Y2SUM,PTS,XF,YF,WEIGHT
	DIMENSION	X(200),Y(200),S(200),A(20),SIGA(20)
	DIMENSION	XFIT(200),YFIT(200)
	DIMENSION XSUM(20,20),B(20,20),YSUM(20),XF(20),ITERM(20)
	COMMON	/SCRAT/ B
	DATA	ITYP /'P','S','C','L',6*0/
	DATA	NMAX/ 20/
1	CALL	ASKS('Name of file with data to fit ',IA,ISIZE,0,80,
	1 'SY:FITIT.DAT')
	CALL	ASSIGN(1,IA,ISIZE)
C
C	 READ IN WEIGHTING MODE AND DATA POINTS.
C	IF MODE	=
C		  -1 (STATISTICAL WEIGHTING--WEIGHTING FACTOR	= 1/Y)
C		   0 (NO WEIGHTING--WEIGHTING FACTOR	= 1)
C		  +1 (INSTRUMENTAL WEIGHTING--WEIGHTING FACTOR	= 1/S**2)
C
	READ(1,1007) TITLE
	READ(1,1007)
	READ(1,1007)
	READ(1,1000,ERR=1,END=1) NPTS,ICHAR,MODE
	READ(1,1007)
	IF(NPTS .GT. 200) STOP '*** ERROR *** ONLY 200 DATA POINTS ALLOWED'
	IF (MODE .GT. 0) READ(1,*) (X(J),Y(J),S(J), J=   1,NPTS)
	IF (MODE .LE. 0) READ(1,*) (X(J),Y(J), J=  1,NPTS)
	WRITE(5,1005) TITLE
C
C	SET UP TO GRAPH RESULTS.
C
	XMAX 	= X(1)
	XMIN 	= XMAX
	YMAX	= Y(1)
	YMIN	= YMAX
	DO 10  J	= 1, NPTS
	IF(MODE .EQ. 0) S(J)	= 1.
	IF(MODE .LT. 0) S(J)	= SSQRT( Y(J) )	!STATISTICAL WEIGHT
	IF(S(J) .LE. 0.) STOP '*** ERROR *** BAD ERRORS ON Y'
	IF (XMAX .LT. X(J)) XMAX	= X(J)
	IF (XMIN .GT. X(J)) XMIN	= X(J)
	IF(YMIN .GT. Y(J) ) YMIN	= Y(J)
	IF(YMAX .LT. Y(J) ) YMAX	= Y(J)
10	CONTINUE
	XSTEP 	= (XMAX - XMIN)/200.
	NTERMS	= 0
11	WRITE(5,1008)
	CALL	ASKS('Type of fit',ITYPIN,ISIZE,0,10)
	DO 15 J	= 1,10
	IF( ITYPIN(1) .NE. ITYP(J) ) GO TO 15
	ITYPE	= J
	GO TO 16
15	CONTINUE
	CALL	CLOSE(1)
	GO TO 1					! NOW NEW INPUT FILE
16	CONTINUE
	SCALE	= 1.
	OFFSET	= 0.
	IF(ITYPE .NE. 1) GO TO 21		! Do not rescale non polynomial
	CALL	ASK('Rescale data ',ITEST)
	IF(ITEST .EQ. 0) GO TO 21
	SCALE	= 5./(ABS(XMAX) + ABS(XMIN) )
	OFFSET	= (XMAX + XMIN) / 2.
	GO TO 18
17	CONTINUE
	OFFSET	= XMIN
	SCALE	= 360./(XMAX - XMIN)
	IF(ITYPE .EQ. 4) SCALE	= SCALE/2.
18	CONTINUE
	CALL	FTRUNC(2,SCALE)
	CALL	FTRUNC(2,OFFSET)

C
C	START THE FIT BY INIT VARIABLES
C
21	CONTINUE
 	CALL	ERASE
22	DO 25 I	= 1,NMAX
	YSUM(I)	= 0.
	A(I)	= 0.
	SIGA(I)	= 0.
	DO 25 J	= 1,NMAX
25	XSUM(I,J)	= 0.
	Y2SUM	= 0.

	CALL	ASKN('Number of terms to fit',NTERMS,0,NMAX)
	IF(NTERMS .LE. 0) GO TO 11
	XF(1)	= 1.0
	WRITE (5,1013)
	CALL	ASKN('Terms to fit ',NEVEN)
	DO 30 J	= 1,NTERMS
	ITERM(J)	= 1
30	CONTINUE
	IF(NEVEN .EQ. 0) GO TO 35
	MIN	= 3
	IF(NEVEN .LT. 0) MIN	= 2
	DO 31 J	= MIN,NTERMS,2
31	ITERM(J)	= 0
35	CONTINUE
C 
C	DO FITTING.
C
	DO 400 I	= 1,NPTS
	XX	= X(I)
	XX	= (XX - OFFSET) * SCALE
	CALL	XFUNC(XF,NTERMS,XX,ITYPE)	!GET THE SERIES
	WEIGHT	= 1./S(I)**2
	YF	= Y(I)
400	CALL	DRGSUM(XF,NTERMS,YF,WEIGHT,XSUM,YSUM,Y2SUM)
	PTS	= NPTS
	CALL	DRGVAL(NTERMS,PTS,XSUM,YSUM,Y2SUM,ITERM,
	1		A,SIGA,CHISQR,FTEST,B)
	DO 450 I	= 1,200			!CALCULATE THE FITTED CURVE
	XX	= XMIN + XSTEP*(I-1)	!THE X VALUE OF THE CURVE
	XFIT(I)	= XX
	XX	= (XX - OFFSET) * SCALE
	CALL	XFUNC(XF,NTERMS,XX,ITYPE)
	YFIT(I)	= 0.
	DO 450 J	= 1,NTERMS
	YFIT(I)	= YFIT(I) + A(J) * XF(J)	!THE Y VALUE OF CURVE
450	CONTINUE

C
C	WRITE OUT RESULTS.
C
	WRITE (5,1001) CHISQR,FTEST,OFFSET,SCALE,
	1		(A(J),SIGA(J), J=1, NTERMS)
510	CALL	ASK('Plot result',ITEST)
	IF(ITEST .EQ. 0) GO TO 22
C
C	PLOT THE RESULTS ON THE ADM-3A
C
500	ICHAR  =  -4
	IF (MODE .NE. 0) ICHAR = -104	!PLOT WITH ERROR BARS
	CALL	SSET		! Attach the terminal
	CALL	DISCKP
	CALL	QPLOT(X,Y,S,NPTS,-1,ICHAR)	!PLOT RAW DATA POINTS
	IF(NTERMS .LE. 0) GO TO 600
	ENCODE(60,1003,IA)CHISQR
	CALL	CHAR(0,760,TITLE,40)
	CALL	CHAR(256,0,IA,60)		!WRITE OUT THE CHISQR
	CALL	QPLOT(XFIT,YFIT,S,200,1,0)	!PLOT THE FIT
	CALL	TIDET				! Detach terminal
	CALL	ENACKP

C
C	SECTION TO OUTPUT RESULT
C
600	CALL	ASK('Output ',ITEST)
	IF (ITEST .EQ. 0) GO TO 21
	CALL	PLTOLP			!OUTPUT PLOT TO LP:
	CALL	SSET
	IF( NTERMS .LE. 0 ) GO TO 21
	WRITE (6,1001) CHISQR,FTEST,OFFSET,SCALE,
	1		(A(J),SIGA(J), J=1, NTERMS)
620	CALL	ASK('Append fit to input file',ITEST)
	IF(ITEST .EQ. 0) GO TO 21
	J	= 200
	WRITE(1,1016)
	WRITE(1,1003) CHISQR
625	WRITE(1,*) (XFIT(J),YFIT(J),J=1,200)
	GO TO 21
C
C	SECTION TO END PROGRAM
C
900	CLOSE(UNIT=6,DISPOSE='PRINT')
	CALL	DISCKP
	CALL	SSET
	CALL	ERASE
	CALL	EXIT

1000	FORMAT(8I10)
1001	FORMAT('  CHI SQUARED/DEG-OF-FREEDOM =	'1P,G12.3
	1	/'  FTEST =			'G12.3
	2	/'  OFFSET ='G15.7'	SCALE ='G15.7
	2	/' COEFFICIENTS:	ERRORS: '
	3	/(1P,G14.6,'	',G12.2) )
1002	FORMAT(1P,G13.6)	
1003	FORMAT('POLYNOMIAL FIT   CHI SQ/FREEDOM= '1P,G10.3)
1005	FORMAT(' 		Title of dataset:'/ 1X,60A1 //)
1007	FORMAT(80A1)
1008	FORMAT(' P	= POLYNOMIAL'/' S	= SINE'/' C	= COSINE'
	1	/' L	= LEGENDRE POLY.'//)
1013	FORMAT( ' 	+1 - Fit to EVEN terms ONLY'/
	2 ' 	 0 - Fit to ALL terms'/
	3 ' 	-1 - Fit to ODD terms ONLY'/ )
1016	FORMAT(' 200,0,0')
	END

	SUBROUTINE	XFUNC (XF,NTERMS,XX,ITYPE)
	DOUBLE PRECISION	XF(NTERMS),XX,RAD,TEMP,FL
	DATA	RAD /.017453292519943296/	!CONVERTS DEGREES TO RADIANS
	DATA	PI  /3.141592653589793/
	XF(1)	= 1.0
	IF(NTERMS .LE. 1) RETURN
	GO TO (110,120,130,140,150) ITYPE
110	TEMP	= XX
	DO 115 J	= 2,NTERMS
	XF(J)	= TEMP
115	TEMP	= TEMP*XX
	RETURN
120	CONTINUE
	DO 125 J	= 2,NTERMS
	XF(J)	= SIN(XX*RAD*(J-1))
125	CONTINUE
	RETURN
130	CONTINUE
	DO 135 J	= 2,NTERMS
135	XF(J)	= COS(XX*RAD*(J-1))
	RETURN
140	CONTINUE
	TEMP	= COS(XX*RAD)		!FIND TEMP OF VARIABLE
	XF(2)	= TEMP
	DO 145 J	= 3,NTERMS
	FL	= J-1
145	XF(J)	= ( (2.*FL-1)*TEMP*XF(J-1) - (FL-1)*XF(J-2) )/FL
	RETURN
150	CONTINUE
	END
	SUBROUTINE	FROUND(IDIGIT,F)
	INTEGER*4 ITEST,JTEST
	DOUBLE PRECISION VALUE,SCALE
	IF(IDIGIT .GT. 6) RETURN
	IF( F .EQ. 0) RETURN
	VALUE	= ABS(F)
	SCALE	= LOG10(VALUE)
	IF(SCALE .LT. 0.) SCALE	= SCALE - 1.
	IEXP	= SCALE
	IEXP	= IEXP - IDIGIT
	SCALE	= 10.
	SCALE	= SCALE**IEXP
	ITEST	= VALUE/SCALE + .5
	JTEST	= 10
	JTEST	= JTEST**IDIGIT
	IF(ITEST .LT. JTEST)  GO TO 20
	ITEST	= ITEST/10
	SCALE	= SCALE*10.
20	CONTINUE
	IF( F .LT. 0) ITEST	= -ITEST
	F	= ITEST*SCALE
	END
	SUBROUTINE	FTRUNC(IDIGIT,F)
	INTEGER*4 ITEST,JTEST
	DOUBLE PRECISION VALUE,SCALE
	IF(IDIGIT .GT. 6) RETURN
	IF( F .EQ. 0) RETURN
	VALUE	= ABS(F)
	SCALE	= LOG10(VALUE)
	IF(SCALE .LT. 0.) SCALE	= SCALE - 1.
	IEXP	= SCALE
	IEXP	= IEXP - IDIGIT
	SCALE	= 10.
	SCALE	= SCALE**IEXP
	ITEST	= VALUE/SCALE
	JTEST	= 10
	JTEST	= JTEST**IDIGIT
	IF(ITEST .LT. JTEST)  GO TO 20
	ITEST	= ITEST/10
	SCALE	= SCALE*10.
20	CONTINUE
	IF( F .LT. 0) ITEST	= -ITEST
	F	= ITEST*SCALE
	END
	SUBROUTINE MATINV(ARRAY,N,DET)
C  
C  MATRIX INVERSION ROUTINE--BEVINGTON,P.302
C
C  ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED
C  N IS THE ORDER, DET IS THE CALCULATED DETERMINANT
C  THE INVERTED MATRIX REPLACES THE ORIGINAL ONE
C
	DIMENSION	ARRAY(N,N)
	DIMENSION IK(20),JK(20)
	DET	= 1.
	DO 100 K	= 1,N
C
C  FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL
C  PARTIAL PIVOTING
C
	AMAX	= 0.0
21	DO 30 I	= K,N
	DO 30 J	= K,N
	IF(ABS(AMAX) .GT. ABS(ARRAY(I,J))) GO TO 30
	AMAX	= ARRAY(I,J)
	IK(K)	= I
	JK(K)	= J
30	CONTINUE
	IF(AMAX .NE. 0.) GO TO 41
	DET	= 0.
	GO TO 140
41	I	= IK(K)
	IF(I-K) 21,51,43
43	DO 50 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= ARRAY(I,J)
50	ARRAY(I,J)	= -SAVE
51	J	= JK(K)
	IF(J-K)21,61,53
53	DO 60 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= ARRAY(I,J)
60	ARRAY(I,J)	= -SAVE
C
C  INVERT MATRIX
C
61	DO 70 I	= 1,N
	IF(I .EQ. K) GO TO 70
	ARRAY(I,K)	= -ARRAY(I,K)/AMAX
70	CONTINUE
71	DO 80 I	= 1,N
	DO 80 J	= 1,N
	IF(I .EQ. K) GO TO 80
	IF(J .EQ. K) GO TO 80
	ARRAY(I,J)	= ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80	CONTINUE
81	DO 90 J	= 1,N
	IF(J .EQ. K) GO TO 90
	ARRAY(K,J)	= ARRAY(K,J)/AMAX
90	CONTINUE
	ARRAY(K,K)	= 1./AMAX
100	DET	= DET*AMAX
C
C  RESTORE ORDERING OF MATRIX
C
101	DO 130 L	= 1,N
	K	= N-L+1
	J	= IK(K)
	IF(J .LE. K) GO TO 111
	DO 110 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= -ARRAY(I,J)
110	ARRAY(I,J)	= SAVE
111	I	= JK(K)
	IF(I .LE. K) GO TO 130
	DO 120 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= -ARRAY(I,J)
120	ARRAY(I,J)	= SAVE
130	CONTINUE
140	RETURN
	END
	SUBROUTINE	 MMULT(AMAT,BMAT,CMAT,IROWA,ICOLA,ICOLB)
C
C	CALL MMULT(AMAT,BMAT,CMAT,IROWA,ICOLA,ICOLB)
C
C	CMAT	= AMAT X BMAT
C
C	MATRIX MULTIPLICATION
	REAL	AMAT(IROWA,ICOLA),BMAT(ICOLA,ICOLB),CMAT(IROWA,ICOLB)
	DO 10 I	= 1,IROWA
	DO 10 J	= 1,ICOLB
	CMAT(I,J)	= 0
	DO 10 K = 1,ICOLA
	CMAT(I,J)	= CMAT(I,J) + AMAT(I,K)*BMAT(K,J)
10	CONTINUE
	END
	SUBROUTINE POLFIT(X,Y,S,NPTS,NTERM,MODE,A,CHISQR)
	REAL X(NPTS),Y(NPTS),S(NPTS),A(NTERM)
	REAL*8 SUMX(19),SUMY(10),ARRAY(100)
	REAL*8 XTERM,YTERM,CHISQ
	CALL ZERO(A,NTERM*2)
	CHISQR	= 0.
	NTERMS	= NTERM
	IF(NTERMS .GT. 10)NTERMS	= 10
	IF(NTERMS .GE. NPTS)NTERMS	= NPTS - 1
	NMAX	= 2*NTERMS - 1
	CALL ZERO(SUMX,76)
	CALL ZERO(SUMY,40)
	CHISQ	= 0.
	XSAV	= 0.
	DO 20 I	= 1,NPTS
	IF(XSAV .GT. ABS(X(I)) )GO TO 20
	XSAV	= ABS(X(I))
20	CONTINUE
	XSAV	= XSAV/1.5
	DO 50 I	= 1,NPTS
	XI	= X(I)/XSAV
	YI	= Y(I)
	IF(MODE) 32,37,39
32	WEIGHT	= 1./(.5 + ABS(YI) )
	GO TO 41
37	WEIGHT	= 1.
	GO TO 41
39	WEIGHT	= 1./S(I)**2
41	XTERM	= WEIGHT
	DO 44 N	= 1,NMAX
	SUMX(N)	= SUMX(N) + XTERM
44	XTERM	= XTERM * XI
	YTERM	= WEIGHT* YI
	YT	= WEIGHT* XT
	DO 48 N	= 1,NTERMS
	SUMY(N)	= SUMY(N) + YTERM
48	YTERM	= YTERM * XI
49	CHISQ	= CHISQ + WEIGHT*YI**2
50	CONTINUE
C
C	CONSTRUCT MATRICES + CALCULATE COEFFICIENTS
C
	DO 54 J	= 1,NTERMS
	DO 54 K	= 0,NTERMS-1
	N	= J + K
	M	= J + NTERMS*K
54	ARRAY(M)	= SUMX(N)
	DELTA	= DDETRM(ARRAY,NTERMS)
	IF(ABS(DELTA) .LT. 1.E-30)RETURN
61	DO 70 L	= 1,NTERMS
	DO 66 J	= 1,NTERMS
	DO 65 K	= 0,NTERMS-1
	N	= J + K
	M	= J + NTERMS*K
65	ARRAY(M)	= SUMX(N)
	M	= J + NTERMS*(L-1)
66	ARRAY(M)	= SUMY(J)
70	A(L)	= DDETRM(ARRAY,NTERMS) /( XSAV**(L-1) * DELTA )
C
C	CALCULATE CHI SQUARE
C
	DO 75 J	= 1,NTERMS
	CHISQ	= CHISQ - 2.*A(J)*SUMY(J)*XSAV**(J-1)
	DO 75 K	= 1,NTERMS
	N	= J + K - 1
75	CHISQ	= CHISQ + A(J)*A(K)*SUMX(N)*XSAV**(J+K-2)
	CHISQR	= CHISQ / (NPTS - NTERMS)
	RETURN
	END

	SUBROUTINE	REGSUM(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 REGVAL
C
	DIMENSION	X(NORDER),YSUM(NORDER),XSUM(NORDER,NORDER)
	Y2SUM	= Y2SUM + WEIGHT*Y*Y		!SUM OF SQUARES OF Y
	DO 100 I	= 1,NORDER
	XA	= WEIGHT*X(I)
	YSUM(I)	= YSUM(I)  + Y*XA		!SUM OVER XjY
	DO 100 J	= 1,I
	XSUM(I,J)	= XSUM(I,J) + XA*X(J)	!SUM OVER XiXj
100	CONTINUE
	END

	SUBROUTINE	REGVAL(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
	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.			!DEFAULT VALUES IF BAD INPUT
	FTEST	= 0.
	IF(NORDER .LT. 1) RETURN	!BAD NUMBER OF TERMS
	IF(XSUM(1,1) .EQ. 0.) RETURN	!BAD SUM OVER WEIGHTING FUNCTION
	IF(NORDER .LT. 2) GO TO 11	!ONLY SINGLE CONSTATNT TO DETERMINE
	DO 10 I	= 1,NORDER-1
	DO 10 J	= I+1,NORDER
10	XSUM(I,J)	= XSUM(J,I)	!SYMMETRIZE THE MATRIX
11	DO 12 I	= 1,NORDER
	DO 12 J	= 1,NORDER
12	B(I,J)	= XSUM(I,J)		!MOVE MATRIX TO SCRATCH MATRIX
C
C	GET RID OF UNWANTED TERMS
C
	DO 13 I	= 1,NORDER
	DO 13 J	= 1,NORDER
	IF(ITERM(I) * ITERM(J) .NE. 0) GO TO 13
	B(I,J)	= 0.			!UNWANTED ROW/COLUMN
	IF(I .EQ. J) B(I,J)	= 1.	!SO IS UNIT MATRIX
13	CONTINUE
	CALL	MATINV(B,NORDER,DET)	!INVERT THE SCRATCH MATRIX
	NTERMS	= 0
	DO 14 I	= 1,NORDER
	IF( ITERM(I) .NE. 0 ) GO TO 14
	NTERMS	= NTERMS + 1
	B(I,I)	= 0.			!NULL VALUE FOR UNWANTED TERMS
14	CONTINUE
	NTERMS	= NORDER - NTERMS
	IF(DET .EQ. 0.) RETURN		!BAD MATRIX INVERSION
	CALL	MATMUL(NORDER,B,NORDER,YSUM,1,A)	!THE RESULT

C
C	CALCULATE THE REDUCED CHI SQUARED/POINT AND THE ERROR
C
	FREEN	= PTS - NTERMS - 1		!THE DEGREES-OF-FREEDOM
	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)		!FIND CHI SQUARED
	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	MATMUL(IX,A,IY,B,IZ,C)
C
C	THIS SUBROUTINE MULTIPLIES 2 MATRICES A,B TOGETHER TO FORM C
C
	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
      FUNCTION SSQRT(A) 
      IF(A)  10,20,20 
   10 SSQRT =-SQRT(-A)
      RETURN
   20 SSQRT = SQRT(A) 
      RETURN
      END