PROGRAM FITIT
C
C	PROGRAM TO FIT DATA POINTS
C	FIT ROUTINE LEGFIT IS TAKEN FROM BEVINGTON'S TEXT (PAGE 155).
C
	IMPLICIT REAL*16(A-H,O-Z)
	PARAMETER NMAX=20
	PARAMETER IMAX=1000
	BYTE	IB(132)
	CHARACTER IA*80,IT*80, ITYPIN*1,ITYP(10)*1,IFIT(4)*10
	BYTE	TITLE(40)
	REAL*4	X(IMAX),Y(IMAX),S(IMAX)
	REAL*4	XFIT(IMAX),YFIT(IMAX),SCALE,OFFSET
	DIMENSION	A(NMAX),SIGA(NMAX)
	DIMENSION XSUM(NMAX,NMAX),B(NMAX,NMAX),YSUM(NMAX),XF(NMAX),ITERM(NMAX)
	COMMON	/SCRAT/ B
	DATA	ITYP /'P','S','C','L','p','s','c','l',2*'Q'/
	DATA	IFIT/'Polynomial','Sine','Cosine','Legendre'/
	CALL	ERASE
	TYPE 1050
1	CALL	ASK('End of fitting',ITEST,0)
	IF(ITEST .ne. 0) CALL	EXIT
	CALL	ASKSD('Name of file with data to fit ',IA,ISIZE,0,80,
	1 'FITIT.DAT')
	OPEN(UNIT=1,FILE=IA,STATUS='OLD',ERR=1,READONLY)
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
2	READ(1,1007,end=1,err=1) TITLE
	IF(TITLE(1) .eq. '.') go to 2		! Scale adjustment ?
	READ(1,1007,end=1,err=1)
	READ(1,1007,end=1,err=1)
	READ(1,1000,ERR=1,END=1) NPTS,ICHAR,MODE
	READ(1,1007,end=1,err=1)
	IF(NPTS .GT. IMAX) STOP '*** ERROR *** TOO MANY DATA POINTS'
	IF (MODE .GT. 0) THEN			! Statistics included ?
	  READ(1,*,end=1,err=1) (X(J),Y(J),S(J), J=   1,NPTS)
	ELSE
	  READ(1,*,end=1,err=1) (X(J),Y(J), J=  1,NPTS)
	ENDIF
	CLOSE(unit=1)
	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.
	IF(XSTEP .eq. 0.) THEN
	  TYPE *,'BAD input data'
	  GO TO 1
	ENDIF
	NTERMS	= 0
11	WRITE(5,1008)
	CALL	ASKS('Type of fit',ITYPIN,ISIZE,0,1)
	IF(ISIZE .eq. 0) go to 1
	DO 15 J	= 1,10
	IF( ITYPIN .NE. ITYP(J) ) GO TO 15
	ITYPE	= J
	GO TO 16
15	CONTINUE
	GO TO 1					! NOW NEW INPUT FILE
16	CONTINUE
	IF(ITYPE .GT. 4) ITYPE = ITYPE - 4	! HANDLE LOWER CASE
	SCALE	= 1.
	OFFSET	= 0.
	IF(ITYPE .NE. 1) GO TO 22		! Do not rescale non polynomial
	TYPE *,'Unscaled polynomial fits can blow up'
	CALL	ASK('Rescale data ',ITEST,0)
	IF(ITEST .EQ. 0) GO TO 22
	OFFSET	= (XMAX + XMIN) / 2.
	SCALE	= 5./(XMAX-OFFSET)
	CALL	FTRUNC(2,SCALE)
	CALL	FTRUNC(2,OFFSET)
	GO TO 22

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.
	I	= NMAX
	IF(NPTS .lt. I) I=NPTS
	CALL	ASKN('Number of terms to fit ',NTERMS,0,I,0)
	IF(NTERMS .LE. 0) GO TO 11
	XF(1)	= 1.0
	WRITE (5,1013)
	CALL	ASKN('Terms to fit ',NEVEN,-1,1,0)
d	TYPE *,(X(J),Y(J),J=1,NPTS),NPTS,NTERMS,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)

d	type *,'XF'
d	type 1011,(XF(J),J=1,nterms)
d	type *,'YF'
d	type 1011,YF
d1011	FORMAT(1x,1p,8G15.6)

400	CALL	HRGSUM(XF,NTERMS,YF,WEIGHT,XSUM,YSUM,Y2SUM)
	PTS	= NPTS
	CALL	HRGVAL(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		(J,A(J),SIGA(J), J=1, NTERMS)
510	CALL	ASK('Plot result',ITEST,-1)
	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
c	CALL	DISCKP				! Disable checkpointing
	CALL	QPLOT(X,Y,S,NPTS,-1,ICHAR)	!PLOT RAW DATA POINTS
	IF(NTERMS .LE. 0) GO TO 600
	ENCODE(60,1003,IT)IFIT(ITYPE),CHISQR
	CALL	PLCHAR(220,764,%DESCR(TITLE),40)
	CALL	PLCHAR(256,0,%DESCR(IT),60)		!WRITE OUT THE CHISQR
	CALL	QPLOT(XFIT,YFIT,S,200,1,0)	!PLOT THE FIT
c	CALL	TIDET				! Detach terminal
c	CALL	ENACKP

C
C	SECTION TO OUTPUT RESULT
C
600	CALL	ASK('Output ',ITEST,0)
	IF (ITEST .EQ. 0) GO TO 21
	CALL	PLTOLP			!OUTPUT PLOT TO LP:
	IF( NTERMS .LE. 0 ) GO TO 21
	IF( IOUT6 .eq. 0)
	1 OPEN(6,FILE='FITIT.OUT',DISPOSE='print',status='new')
	IOUT6	= IOUT6 + 1
	WRITE (6,1001) CHISQR,FTEST,OFFSET,SCALE,
	1		(J,A(J),SIGA(J), J=1, NTERMS)
620	CALL	ERASE
	CALL	ASK('Append fit to input file',ITEST,0)
	IF(ITEST .EQ. 0) GO TO 22
	J	= 200
	OPEN(UNIT=1,FILE=IA,STATUS='OLD',ERR=901,readonly)
	OPEN(UNIT=2,FILE=IA,STATUS='new',ERR=902,CARRIAGECONTROL='list')
621	READ(1,1004,END=622) ISIZ,(IB(J),J=1,ISIZ)
	WRITE(2,1006) (IB(J),J=1,ISIZ)
	GO TO 621
622	CLOSE(UNIT=1)
	WRITE(2,1016)
	WRITE(2,1003) IFIT(ITYPE),CHISQR
625	WRITE(2,*) (XFIT(J),YFIT(J),J=1,200)
	CLOSE(UNIT=2)
	GO TO 21
C
C	SECTION TO END PROGRAM
C
901	CALL MSG('ERROR - Unit 1')
	go to 900
902	CALL MSG('ERROR - Unit 2')
900	CLOSE(UNIT=6,DISPOSE='PRINT')
c	CALL	DISCKP
	CALL	EXIT

1000	FORMAT(8I10)
1001	FORMAT('  CHI SQUARED/DEG-OF-FREEDOM =	'1P,G10.2
	1	/'  FTEST =			'G10.3
	2	/'  OFFSET ='G10.3'	SCALE ='G10.3
	2	/' TERM COEFFICIENTS:	ERRORS: '
	3	/(0p,I5,1P,G14.6,'	',G12.2) )
1002	FORMAT(1P,G13.6)	
1003	FORMAT(A,' Fit   chi sq/freedom= '1P,G10.2)
1004	FORMAT(Q,132(A1))
1005	FORMAT(' 		Title of dataset:'/ 1X,60A1 //)
1006	FORMAT(132(A1))
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')
1050	FORMAT(//' 	FITIT - A General fitting routine'
	1 /' 	Input data must be in RICEPLOT format'
	2/' 	SIN/COS/LEGENDRE poly fits X is in degrees'
	3/' 	No answer goes back to previous question'///)
	END

	SUBROUTINE	XFUNC (XF,NTERMS,XX,ITYPE)
	IMPLICIT REAL*16(A-H,O-Z)
	DOUBLE PRECISION	XF(NTERMS),XX,RAD,TEMP,FL
	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)	= SIND(XX*(J-1))
125	CONTINUE
	RETURN
130	CONTINUE
	DO 135 J	= 2,NTERMS
135	XF(J)	= COSD(XX*(J-1))
	RETURN
140	CONTINUE
	TEMP	= COSD(XX)		!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