      program polyfit
C
C Fit a polynomila of order n to a set of data x,y,dy
C
      integer*4 ma,nld,nldp
      parameter (ma=9)         ! maximum order of fit
      parameter (nld=1000)     ! max number of data points
      parameter (nldp=20)      ! plot fit at this # of x-values
      integer*4 lista(ma)
      real*4 p(ma),pa(ma),covar(ma,ma)
      real*4 x(nld),y(nld),sig(nld),chisq
      real*4 xp(nldp),yp(nldp)  ! plot variables
      character header*80

      external fpoly

   10 write(6,*)'Enter name of file with x,y,dy values!'
      read(5,1)header
    1 format(a80)
      open(1,file=header,status='old',err=10)
      read(1,*,end=10)x(1),y(1),sig(1)
      xmax=x(1)
      xmin=xmax
      ymax=y(1)
      ymin=ymax
      do i=2,nld
       read(1,*,end=11)x(i),y(i),sig(i)
       xmax=max(xmax,x(i))
       xmin=min(xmin,x(i))
       ymax=max(ymax,y(i))
       ymin=min(ymin,y(i))
      enddo
   11 numld=i-1
      if(numld.eq.0)goto 10
      write(6,*)'Enter order of polynomial!'
      read(5,*)mfit
      mfit=mfit+1
      if(mfit.gt.9)stop
c
c     Fit ld with polynomial of third order
c
      do i=1,ma
       lista(i)=i
       pa(i)=0.D0
      enddo
      call lfit(x,y,sig,numld,pa,ma,lista,mfit,covar,ma,chisq,fpoly)

c ....Plot data
      call pgbegin(0,'?',1,1)
      call pgsch(1.5)
      xrange=xmax-xmin
      yrange=ymax-ymin
      xmax=xmax+xrange/10.
      xmin=xmin-xrange/10.
      ymax=ymax+yrange/10.
      ymin=ymin-yrange/10.
      call pgenv(xmin,xmax,ymin,ymax,0,0)
      write(header,2)(pa(j),j=1,mfit)
    2 format('Coeff.: ',9(1x,e10.4))
      call pglabel('x','y',header)
      do i=1,numld
       call pgpoint(1,x(i),y(i),17)
       ylo=y(i)-sig(i)
       yhi=y(i)+sig(i)
       call pgerry(1,x(i),ylo,yhi,1.0)
      enddo
c ....Plot analytical fit
      xstep=(xmax-xmin)/(nldp-1)
      do k=1,nldp
       xp(k)=xmin+(k-1)*xstep
       call fpoly(xp(k),p,mfit)
       yp(k)=0.
       do j=1,mfit
        yp(k)=yp(k)+p(j)*pa(j)
       enddo
      enddo
      ndat=min(k,nldp)
      call pgline(ndat,xp,yp)
      call pgend
      end


      SUBROUTINE FPOLY(X,P,NP)
      save
      DIMENSION P(NP)
      P(1)=1.D0
      DO 11 J=2,NP
        P(J)=P(J-1)*X
11    CONTINUE
      RETURN
      END

      SUBROUTINE LFIT(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQ,
     *FUNCS)
      save
      PARAMETER (MMAX=50)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA),
     *    COVAR(NCVM,NCVM),BETA(MMAX),AFUNC(MMAX)
      KK=MFIT+1
      DO 12 J=1,MA
        IHIT=0
        DO 11 K=1,MFIT
          IF (LISTA(K).EQ.J) IHIT=IHIT+1
11      CONTINUE
        IF (IHIT.EQ.0) THEN
          LISTA(KK)=J
          KK=KK+1
        ELSE IF (IHIT.GT.1) THEN
          PAUSE 'Improper set in LISTA'
        ENDIF
12    CONTINUE
      IF (KK.NE.(MA+1)) PAUSE 'Improper set in LISTA'
      DO 14 J=1,MFIT
        DO 13 K=1,MFIT
          COVAR(J,K)=0.
13      CONTINUE
        BETA(J)=0.
14    CONTINUE
      DO 18 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        YM=Y(I)
        IF(MFIT.LT.MA) THEN
          DO 15 J=MFIT+1,MA
            YM=YM-A(LISTA(J))*AFUNC(LISTA(J))
15        CONTINUE
        ENDIF
        SIG2I=1./SIG(I)**2
        DO 17 J=1,MFIT
          WT=AFUNC(LISTA(J))*SIG2I
          DO 16 K=1,J
            COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K))
16        CONTINUE
          BETA(J)=BETA(J)+YM*WT
17      CONTINUE
18    CONTINUE
      IF (MFIT.GT.1) THEN
        DO 21 J=2,MFIT
          DO 19 K=1,J-1
            COVAR(K,J)=COVAR(J,K)
19        CONTINUE
21      CONTINUE
      ENDIF
      CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1)
      DO 22 J=1,MFIT
        A(LISTA(J))=BETA(J)
22    CONTINUE
      CHISQ=0.
      DO 24 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        SUM=0.
        DO 23 J=1,MA
          SUM=SUM+A(J)*AFUNC(J)
23      CONTINUE
        CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
24    CONTINUE
      CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
      RETURN
      END

      SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
      DIMENSION COVAR(NCVM,NCVM),LISTA(MFIT)
      DO 12 J=1,MA-1
        DO 11 I=J+1,MA
          COVAR(I,J)=0.
11      CONTINUE
12    CONTINUE
      DO 14 I=1,MFIT-1
        DO 13 J=I+1,MFIT
          IF(LISTA(J).GT.LISTA(I)) THEN
            COVAR(LISTA(J),LISTA(I))=COVAR(I,J)
          ELSE
            COVAR(LISTA(I),LISTA(J))=COVAR(I,J)
          ENDIF
13      CONTINUE
14    CONTINUE
      SWAP=COVAR(1,1)
      DO 15 J=1,MA
        COVAR(1,J)=COVAR(J,J)
        COVAR(J,J)=0.
15    CONTINUE
      COVAR(LISTA(1),LISTA(1))=SWAP
      DO 16 J=2,MFIT
        COVAR(LISTA(J),LISTA(J))=COVAR(1,J)
16    CONTINUE
      DO 18 J=2,MA
        DO 17 I=1,J-1
          COVAR(I,J)=COVAR(J,I)
17      CONTINUE
18    CONTINUE
      RETURN
      END

      SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
      PARAMETER (NMAX=50)
      DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
      DO 11 J=1,N
        IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
        BIG=0.0
        DO 13 J=1,N
          IF(IPIV(J).NE.1)THEN
            DO 12 K=1,N
              IF (IPIV(K).EQ.0) THEN
                IF (ABS(A(J,K)).GE.BIG)THEN
                  BIG=ABS(A(J,K))
                  IROW=J
                  ICOL=K
                ENDIF
              ELSE IF (IPIV(K).GT.1) THEN
                PAUSE 'Singular matrix'
              ENDIF
12          CONTINUE
          ENDIF
13      CONTINUE
        IPIV(ICOL)=IPIV(ICOL)+1
        IF (IROW.NE.ICOL) THEN
          DO 14 L=1,N
            DUM=A(IROW,L)
            A(IROW,L)=A(ICOL,L)
            A(ICOL,L)=DUM
14        CONTINUE
          DO 15 L=1,M
            DUM=B(IROW,L)
            B(IROW,L)=B(ICOL,L)
            B(ICOL,L)=DUM
15        CONTINUE
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.'
        PIVINV=1.0/A(ICOL,ICOL)
        A(ICOL,ICOL)=1.D0
        DO 16 L=1,N
          A(ICOL,L)=A(ICOL,L)*PIVINV
16      CONTINUE
        DO 17 L=1,M
          B(ICOL,L)=B(ICOL,L)*PIVINV
17      CONTINUE
        DO 21 LL=1,N
          IF(LL.NE.ICOL)THEN
            DUM=A(LL,ICOL)
            A(LL,ICOL)=0.0
            DO 18 L=1,N
              A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18          CONTINUE
            DO 19 L=1,M
              B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19          CONTINUE
          ENDIF
21      CONTINUE
22    CONTINUE
      DO 24 L=N,1,-1
        IF(INDXR(L).NE.INDXC(L))THEN
          DO 23 K=1,N
            DUM=A(K,INDXR(L))
            A(K,INDXR(L))=A(K,INDXC(L))
            A(K,INDXC(L))=DUM
23        CONTINUE
        ENDIF
24    CONTINUE
      RETURN
      END


