program ft
      implicit none
      character*80 s80
      real*8,allocatable::x(:),y(:),p(:),q(:,:),e(:,:),qi(:,:),a(:)
      integer*4 i,ic,j,k,N,np,ij,IERR,npa,iargc
      real*8 t,s,tol,u,t2,s2,f
      write(6,*)'polynomial fit'
      npa=iargc()
      if(npa.lt.1)then
       write(6,600)
600    format('usage: ft N_degree input_file [spectrum_orig. corrected]')
       stop
      else
       call getarg(1,s80)
       read(s80,*)N
       call getarg(2,s80)
      endif
      allocate(p(N),q(N,N),qi(N,N),e(N,N+N),a(N))
      
      np=0
      do 9 ic=1,2
      if(ic.eq.2)allocate(x(np),y(np))
      open(9,file=s80)
      np=0
1     read(9,*,err=9,end=9)t,u
      if(abs(u).gt.0.01d0)then
       np=np+1
       if(ic.eq.2)then
        x(np)=t       
        y(np)=u-t       
       endif
      endif
      goto 1
9     close(9)
      write(6,*)np,' points'
      
      do 3 k=1,N
      p(k)=0.0d0
      do 3 j=1,np
3     p(k)=p(k)+y(j)*x(j)**(k-1)
      do 4 k=1,N
      do 4 j=1,N
      q(j,k)=0.0d0
      do 4 ij=1,np
4     q(j,k)=q(j,k)+x(ij)**(j-1)*x(ij)**(k-1)
      IERR=0
      tol=1.0d-90
      call INV(q,qi,N,tol,e,IERR)
      do 5 i=1,N
      a(i)=0.0d0
      do 5 j=1,N
5     a(i)=a(i)+qi(i,j)*p(j)
      do 6 i=1,N
6     write(6,601)a(i),i-1
601   format(f10.2,'  x^',i4)
      open(9,file='T.PRN')
      s2=0.0d0
      do 7 i=1,np
      s=f(x(i),a,N)
      t2=y(i)-s
      s2=s2+t2**2
7     write(9,602)x(i),y(i),s,t2
602   format(4f15.6)
      s2=dsqrt(s2/dble(np))
      write(6,610)s2
610   format(/,'RMS:',f15.6)
      close(9)

      if(npa.gt.2)then
       call getarg(3,s80)
       open(9,file=s80)
       call getarg(4,s80)
       open(91,file=s80)
       np=0
10     read(9,*,end=88,err=88)u,t
       np=np+1
       write(91,990)u-f(u,a,N),t
990    format(2g12.4)
       goto 10
88     close(9)
       close(91)
       write(6,*)np,' points in spectrum'
      endif


      stop
      end

      SUBROUTINE INV(a,ai,n,TOL,e,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension ai(n,n),a(n,n),e(n,2*n)
C
10000 IERR=0
      DO 1  ii=1,n
      DO 1  jj=1,n
      e(ii,jj)=a(ii,jj)
1     e(ii,jj+n)=0.0D0
      do 13 ii=1,n
13    e(ii,ii+n)=1.0D0
c
      DO 2 ii=1,n-1
      if (ABS(e(ii,ii)).LE.TOL) then
       DO 3 io=ii+1,n
3      if (ABS(e(io,ii)).GT.TOL) goto 11
       IERR=1
       write(6,*)ii
       write(6,*)'tol = ',tol
       tol=tol*0.50d0
       if(tol.gt.1.0d-20)goto 10000
       RETURN
c
11     CONTINUE
       DO 4 kk=1,2*n
       w=e(ii,kk)
       e(ii,kk)=e(io,kk)
4      e(io,kk)=w
      ENDIF
      eii=e(ii,ii)
      DO 5 jj=ii+1,n
      e1=e(jj,ii)/eii
      DO 6 kk=ii+1, 2*n
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e1
5     e(jj,ii)=0.0D0
2     CONTINUE
c
      DO 7  i2=n,2,-1
      eii=e(i2,i2)
      DO 7 j2=i2-1,1,-1
      e1=e(j2,i2)/eii
      DO 9 kk=1, n
9     e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1
7     e(j2,i2)=0.0d0
c
      DO 10 ii=1,n
      ei=1.0d0/e(ii,ii)
      DO 12  jj=1,n
12    ai(ii,jj)=e(ii,jj+n)*ei
10    CONTINUE
c
      RETURN
      END

      function f(x,a,N)
      implicit none
      real*8 x,a(*),s,f
      integer*4 N,j
      s=0.0d0
      do 8 j=1,N
8     s=s+a(j)*x**(j-1)
      f=s
      return
      end