C+
C
C  SUBROUTINE POLFT1
C
C  PURPOSE
C	MAKE A LEAST SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE
C	Y = A(1) + A(2)*X + A(3)*X**2 + ...
C
C  USAGE
C	CALL POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR,ARR,IER)
C
C  DESCRIPTION OF PARAMETERS
C	Y	- ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C	SIGMAY	- ARRAY OF STANDARD DEVIATIONS FOR Y-DATA POINTS
C		  (OR - IN CASE MODE = 4 - WEIGHTS FOR POINTS)
C	NPTS	- NO. OF PAIRS OF DATA POINTS
C	NTERMS	- NO. OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
C	MODE	- DETERMINES METHOD OF WEIGHTING LEAST SQUARES FIT
C		  4 (USER DEFINED) WEIGHT(I) = SIGMAY(I)
C		  3 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
C		  2 (NO WEIGHTING) WEIGHT(I) = 1.
C		  1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
C	A	- ARRAY OF COEFFICIENTS OF POLYNOMIAL
C	CHISQR	- REDUCED CHISQUARE FOR FIT
C	ARR	- DOUBLE PRECISION WORK BUFFER; MUST BE AT LEAST
C		  400 WORDS LONG IN THE CALLING ROUTINE
C	IER	- ERROR PARAMETER
C		  -1 DET=0
C		   0 O.K.
C		  +1 NOT ENOUGH POINTS
C
C  SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	DET(ARR,NORD) - EVALUATES THE DETERMINANT OF A SYMMETRIC
C			TWO DIMENSIONAL MATRIX OF ORDER NORD
C
C  COMMENTS
C	DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C	FOR DETAILS OF ALGORITHM SEE "DATA REDUCTION AND ERROR
C	ANALYSIS FOR THE PHYSICAL SCIENCES" - PH. R. BEVINGTON
C	MC GRAW-HILL BOOK COMPANY
C
C	IN THIS SPECIAL VERSION THE ELEMENTS OF COEFFICIENT MATRIX
C	ARR ARE NORMALIZED WITH RESPECT TO A VALUE COMPUTED VIA
C	LOGARITHMIC INTERPOLATION BETWEEN THE SMALLEST AND LARGEST
C	MATRIX ELEMENT
C
C  MODIFICATIONS BY A. SCHERMANN - INST. FOR ASTRONOMY,
C	UNIV. OF VIENNA, AUSTRIA
C
C Modified to assume x-variable is index of Y array
C-
	SUBROUTINE POLFT1(Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR,ARR,IER)
	DOUBLE PRECISION SUMX(19),SUMY(10),ARR(10,10),XTERM,YTERM,
     1	CHISQ
	DIMENSION Y(1),SIGMAY(1),A(1)
	IER=0
C
C  CHECK DEGREE OF FREEDOM
	FREE=NPTS-NTERMS
	IF(FREE.GT.0.) GOTO 11
	IER=1
	GOTO 80
C
C  ACCUMULATE WEIGHTED SUMS
11	NMAX=2*NTERMS-1
	DO 13 N=1,NMAX
13	SUMX(N)=0.
	DO 15 J=1,NTERMS
15	SUMY(J)=0.
	CHISQ=0.
21	DO 50 I=1,NPTS
	XI=I
	YI=Y(I)
31	GOTO(32,37,38,39),MODE
32	IF(YI)35,37,33
33	WEIGHT=1./YI
	GOTO 41
35	WEIGHT=-1./YI
	GOTO 41
37	WEIGHT=1.
	GOTO 41
38	WEIGHT=1./SIGMAY(I)**2
	GOTO 41
39	WEIGHT=SIGMAY(I)
41	XTERM=WEIGHT
	DO 44 N=1,NMAX
	SUMX(N)=SUMX(N)+XTERM
44	XTERM=XTERM*XI
45	YTERM=WEIGHT*YI
	DO 48 N=1,NTERMS
	SUMY(N)=SUMY(N)+YTERM
48	YTERM=YTERM*XI
49	CHISQ=CHISQ+WEIGHT*YI**2
50	CONTINUE
C
C  GET LARGEST AND SMALLEST MATRIX ELEMENT (FOR NORMALIZATION)
	XTERM=SUMX(1)
	YTERM=XTERM
	DO 100 I=2,NMAX
	IF(SUMX(I).GT.XTERM) XTERM=SUMX(I)
	IF(SUMX(I).LT.YTERM) YTERM=SUMX(I)
100	CONTINUE
	DO 110 I=1,NTERMS
	IF(SUMY(I).GT.XTERM) XTERM=SUMY(I)
	IF(SUMY(I).LT.YTERM) YTERM=SUMY(I)
110	CONTINUE
	IF(YTERM.LE.0.) YTERM=1.D0
C
C  LOGARITHMIC INTERPOLATION OF NORMALIZATION VALUE
	XTERM=1.D1**((DLOG10(XTERM)+DLOG10(YTERM))/2.)
C
C  CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS
51	DO 54 J=1,NTERMS
	DO 54 K=1,NTERMS
	N=J+K-1
54	ARR(J,K)=SUMX(N)/XTERM
	DELTA=DET(ARR,NTERMS)
	IF(DELTA.NE.0.) GOTO 61
57	CHISQR=0.
	DO 59 J=1,NTERMS
59	A(J)=0.
	IER=-1
	GOTO 80
61	DO 70 L=1,NTERMS
62	DO 66 J=1,NTERMS
	DO 65 K=1,NTERMS
	N=J+K-1
65	ARR(J,K)=SUMX(N)/XTERM
66	ARR(J,L)=SUMY(J)/XTERM
70	A(L)=DET(ARR,NTERMS)/DELTA
C
C  CALCULATE CHISQUARE
71	DO 75 J=1,NTERMS
	CHISQ=CHISQ-2.*A(J)*SUMY(J)
	DO 75 K=1,NTERMS
	N=J+K-1
75	CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
77	CHISQR=CHISQ/FREE
80	RETURN
	END
c-----------------------------------------------------------
	FUNCTION POLYNO (COE,NPOL,IX)
C
C EVALUATE A POLYNOMIAL OF ORDER NPOL WITH COEFFICIENTS COE(1),...,
C COE (NPOL+1)  FOR AN INDEX IX
C
	DIMENSION COE(*)
C
	IF(NPOL.GT.0) GO TO 10
	POLYNO=COE(1)
	RETURN
C
10	POLYNO=COE(NPOL+1)
	DO 20 I=NPOL,1,-1
20	    POLYNO=POLYNO*IX+COE(I)
C
	RETURN
	END
c-----------------------------------------------------------
C+
C
C  FUNCTION DET
C
C  PURPOSE
C	CALCULATE DETERMINANT OF A SQUARE MATRIX
C
C  USAGE
C	DET = DET (ARR,NORDER)
C
C  DESCRIPTION OF PARAMETERS
C	ARRAY	- MATRIX
C	NORDER	- DEGREE OF MATRIX
C
C  SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C		NONE
C
C  COMMENTS
C	THIS SUBPROGRAM DESTROYS THE INPUT MATRIX ARRAY
C
	FUNCTION DET(ARRAY,NORDER)
	DOUBLE PRECISION ARRAY(10,10),SAVE
	DET=1.
	DO 90 K=1,NORDER
C  INTERCHANGE COLUMNS, IF DIAGONAL ELEMENT IS ZERO
	IF(ARRAY(K,K).NE.0.) GOTO 40
	DO 10 J=K,NORDER
	IF(ARRAY(K,J).NE.0.) GOTO 20
10	CONTINUE
	DET=0.
	GOTO 100
20	DO 30 I=K,NORDER
	SAVE=ARRAY(I,J)
	ARRAY(I,J)=ARRAY(I,K)
30	ARRAY(I,K)=SAVE
	DET=-DET
C  SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX
40	DET=DET*ARRAY(K,K)
	IF(K-NORDER.GE.0) GOTO 90
	K1=K+1
	DO 50 I=K1,NORDER
	DO 50 J=K1,NORDER
50	ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
90	CONTINUE
100	RETURN
	END