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