!	call  POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)  !POLYNOM FIT
	SUBROUTINE POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)  
	implicit none
        integer:: i,j,k,l,n,nmax,npts,nterms,mode
	real*8::  XTERM,YTERM,CHISQ,TERM
	real*8:: SUMX(19),SUMY(19),ARRAY(10,10)
	real::  X(NPTS),Y(NPTS),SIGMAY(NPTS),A(NPTS),xi,yi,weight
        real:: determ,delta,chisqr,free
!
!	ACCUMULATE WEIGHTED SUMS
	NMAX=2*NTERMS-1
        sumx(1:nmax)=0.0
        SUMY(1:nterms)=0.0
	CHISQ=0.
	DO I=1,NPTS
	 XI=X(I)
	 YI=Y(I)
	 IF(MODE.lt.0) then !32,37,39
	  IF(YI.gt.0) then !35,37,33
	   WEIGHT=1./YI
	  elseif(yi.lt.0) then
	   WEIGHT=1./(-YI)
	  else
	   WEIGHT=1.
	  endif
         elseif(mode.eq.0) then
          weight=1.
         else 
	  WEIGHT=1./SIGMAY(I)**2
         endif
	 XTERM=WEIGHT

	 DO N=1,NMAX
	  SUMX(N)=SUMX(N)+XTERM
	  XTERM=XTERM*XI
	 enddo
	 YTERM=WEIGHT*YI
	 DO N=1,NTERMS
	  SUMY(N)=SUMY(N)+YTERM
	  YTERM=YTERM*XI
	 enddo
	 CHISQ=CHISQ+WEIGHT*YI**2
	enddo
!
!	CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS
!
51	DO J=1,NTERMS
	DO K=1,NTERMS
	N=J+K-1
54	ARRAY(J,K)=SUMX(N)
	enddo
	enddo
	DELTA=DETERM(ARRAY,NTERMS)
	IF(DELTA) 61,57,61
57	CHISQR=0.
	DO J=1,NTERMS
59	A(J)=0.
	enddo
	GO TO 80
61	DO L=1,NTERMS
62	DO J=1,NTERMS
	DO K=1,NTERMS
	N=J+K-1
65	ARRAY(J,K)=SUMX(N)
	enddo
66	ARRAY(J,L)=SUMY(J)
	enddo
70	A(L)=DETERM(ARRAY,NTERMS)/DELTA
	enddo
!
!	CALCULATE CHI SQUARE
!
	IF(MODE.EQ.0) GO TO 79
71	DO J=1,NTERMS
	CHISQ=CHISQ-2.*A(J)*SUMY(J)
	DO K=1,NTERMS
	N=J+K-1
75	CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
	enddo
	enddo
76	FREE=max(NPTS-NTERMS,1) !avoid div/0
77	CHISQR=CHISQ/FREE
79	CHISQ=0.
	DO I=1,NPTS
	TERM=Y(I)
	DO J=1,NTERMS
	IF (J .EQ. 1) THEN
	TERM = TERM - A(J)
	ELSE
	TERM=TERM-A(J)*X(I)**(J-1)
	ENDIF
200	CONTINUE
	enddo
	CHISQ=CHISQ+TERM**2
100	CONTINUE
	enddo
	FREE=NPTS-NTERMS
	CHISQR=CHISQ/FREE
80	RETURN
	END




	REAL FUNCTION DETERM(ARRAY,NORDER)
	real*8 ARRAY(10,10),SAVE
        integer:: norder,i,j,k,k1
10	DETERM=1.
11	DO K=1,NORDER
!
!	INTERCHANGE COLUMNS IF DIAGONAL ELEMENT IS ZERO
!
	 IF(ARRAY(K,K)) 41,21,41
21	 DO J=K,NORDER
	  IF(ARRAY(K,J)) 31,23,31
23	  CONTINUE
	 enddo
	 DETERM=0.
	 return
31	 DO I=K,NORDER
	  SAVE=ARRAY(I,J)
	  ARRAY(I,J)=ARRAY(I,K)
	  ARRAY(I,K)=SAVE
	 enddo
!
!	SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX
!
41	 DETERM=DETERM*ARRAY(K,K)
	 IF(K.lt.NORDER) then
	  K1=K+1
	  DO I=K1,NORDER
	   DO J=K1,NORDER
 	    ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
	   enddo
	  enddo
	 endif
	enddo
	RETURN
	END