c c c c
c
      program curvefit
c
c       Take points from a file  and fits the requested
c       order polynomials to them
c
c   First argument on the command line is the name of the input file
c   containing one x,y pair of data values per line
c   If this argument is missing the program assumes data is in a file
c   called "curve.data"
c
c   Second argument on the command line is 'plot' or 'noplot' to
c   request or suppress graphics.  If this argument is absent plots
c   are generated
c
c   The program will request a list of polynomial fits. Type a blank
c   delimited list of integers giving the orders of polynomials you
c   want, or 'help' for more information.
c
      parameter (nmax=20,ndata=1000)
      character*32 infile, arg2
      common/curf/ xd(ndata),yd(ndata),ipvt(nmax),aa(nmax,nmax),
     $ bb(nmax),nl(30),nu(30),isym(30),nd
      character line*80,sym(30)*16
      narg=iargc()
      if(narg.eq.0) then
         infile='curve.data'
         arg2='plot'
c
      else if (narg.eq.1) then
c
         call getarg(1,infile)
         if(infile.eq.'no'.or.infile.eq.'n'.or.infile.eq.'noplot') then
            infile='curve.data'
            arg2='noplot'
         else
            arg2='plot'
            endif
      else
         call getarg(1,infile)
         call getarg(2,arg2)
         if(arg2.eq.'no'.or.arg2.eq.'n'.or.arg2.eq.'noplot') then
            arg2='noplot'
         else
            arg2='plot'
            endif
         endif
      if (infile.eq.' ') infile='curve.data'
      open (11,file=infile)
      open (12,file='curve.fit')
      i=1
  10  read (11,*,end=20) xd(i),yd(i)
      i=i+1
      go to 10
   20 nd=i-1
      write(6,2000) nd
 2000 format(1x,i4,' data points read from input file')
      if(nd.le.0) then
         write(6,*) ' No curve data found.'
         stop
         endif
      if(nd.eq.1) then
         write(6,*) 'More than one data data point required for a fit.'
         stop
         endif
   30 write(6,*) 'Provide a list of polynomial orders for fits.'
      read(5,'(a80)',end=500) line
      call parsl(line,sym,30,nsym,isym,ierr)
      if(sym(1).eq.'end'.or.sym(1).eq.'END'.or.sym(1).eq.'hx'.or.
     $   sym(1).eq.'HX'.or.sym(1).eq.'EXIT'.or.sym(1).eq.'exit'.or.
     $   sym(1).eq.'quit'.or.sym(1).eq.'QUIT'.or.sym(1).eq.' ')go to 500
      if(ierr.ne.0) go to 30
      if(sym(1).eq.'HELP'.or.sym(1).eq.'help'.or.sym(1).eq.'?') then
         write(6,*)
     $ ' This program produces simple least squares polynoial fits to'
         write(6,*)
     $ ' data pairs in the file curve.data.  Each line in this file'
         write(6,*)
     $ ' contains a data (x,y) pair in list directed (free) format.'
         write(6,*)
     $ ' More than one polynomial fit may be requested on a given'
         write(6,*)
     $ ' input prompt by providing a list of blank delimited integers'
         write(6,*)
     $ ' for the desired orders of the polynomials.  A range of fits'
         write(6,*)
     $ ' may be requested by entering two integers separated by a'
         write(6,*)
     $ ' colon.'
         go to 20
         endif
      do 90 i=1,nsym
         call irange(sym(i),jl,ju,ierr)
         if(ierr.ne.0) then
            write(6,*)'Looking for integer in ',sym(i)
            go to 90
            endif
         if(ju.ge.nmax) then
      	    write(6,2028) nmax-1
 2028       format('The largest power this program handles is: ',i3)
            go to 90
            endif
         jl=min0(nd,jl+1)
         ju=min0(nd,ju+1)
         do 80 j=jl,ju
            do 55 ir=1,j
               sumb=0.
               do 45 ij=1,nd
  45              sumb=sumb+yd(ij)*xd(ij)**(ir-1)
               do 50 is=1,j
                  suma=0.
                  do 52 ij=1,nd
  52                 suma=suma+xd(ij)**(is-1)*xd(ij)**(ir-1)
  50              aa(ir,is)=suma
  55           bb(ir)=sumb
            call sgefa(aa,nmax,j,ipvt,info)
            call sgesl(aa,nmax,j,ipvt,bb,0)
            rsid=0.
            do 65 ir=1,nd
               fs=0.
               do 60 is=1,j
   60             fs=fs+bb(is)*xd(ir)**(is-1)
   65          rsid=rsid+(yd(ir)-fs)**2
            rsid=sqrt(rsid/float(nd-1))
            write(6,2001) (j-1),rsid
            write(12,2001) (j-1),rsid
 2001       format(' Fit to polynomial of order ',i3,' with mean',
     $              ' error of',1p,e12.5)
      call wpoly(bb(1),j-1)
      ii=index(infile,' ')-1
      if(arg2.eq.'plot') call plotit(bb(1),j,infile(1:ii))
  80        continue
  90     continue
      go to 30
 500  stop
      end
      SUBROUTINE PARSL (LINE,SYM,NFND,NSYM,ISYM,IERR)
C
C    SEPARATE THE INDIVIDUAL WORDS IN THE COMMND LINE FOR
C    LATER INTERPRETATION.
C
      DIMENSION ISYM(*)
      CHARACTER LINE*(*),SYM(*)*(*)
      LSYM=LEN(SYM(1))
      LLINE=LEN(LINE)
      NSYM=0
      IERR=0
      DO 10 I=1,NFND
   10 SYM(I)=' '
      I=1
      J=1
C
 20   IF(LINE(I:I).EQ.' ') GO TO 24
      IF(LINE(I:I).EQ.',') GO TO 24
      IF(LINE(I:I).EQ.'!')  GO TO 80
      GO TO 30
 24   I=I+1
      IF (I.GT.LLINE) GO TO 80
C
C   CHANGE ABOVE TO AN ERROR CALL.
C
      J=1
      NEW=1
      GO TO 20
   30 IF (1.LT.J) GO TO 40
      IF(SYM(NSYM).EQ.'!') GO TO 80
      NSYM=NSYM+1
      ISYM(NSYM)=I
      IF(NSYM.GT.NFND) GO TO 80
   40 IF (LSYM.LT.J) GO TO 55
      SYM(NSYM)(J:J)=LINE(I:I)
      I=I+1
      J=J+1
      GO TO 20
 55   IERR=1
      SYM(NSYM)=' '
      WRITE(6,2002) LSYM
 2002 FORMAT(' THIS WORD HAS MORE THAN',I3,' CHARACTERS')
      CALL POINAT(LINE,9,I)
      IL=I+1
      DO 60 IS=IL,LLINE
      I=IS
      IF(LINE(IS:IS).EQ.' ') GO TO 24
   60 CONTINUE
 80   RETURN
C
      END
*DECK POINAT
      SUBROUTINE POINAT(LINE,NW,IC)
C
C   WRITE LINE AND A POINTER UNDER CHARACTER IC.
C
      CHARACTER LINE*80
      CHARACTER F*40
      ICM=IC-1
      WRITE(F,2000) ICM
 2000 FORMAT('(',I3,'X,1H)')
      WRITE(6,2020)LINE
 2020 FORMAT(/(1X,A80))
      WRITE(6,F)
      RETURN
      END
      SUBROUTINE IRANGE (SYM,INT1,INT2,IERR)
C
C     GET TWO INTEGER VALUES FROM A RANGE EXPRESSION (INT1:INT2)
C     RETURN IERR=1 IF EITHER IS NOT A VALID INTEGER
C
      CHARACTER SYM*(*),SYM1*16,SYM2*16
      LSYM=LEN(SYM)
      NSYM=2
      IERR=0
      DO 10 I=1,LSYM
      IF(SYM(I:I).EQ.' ') GO TO 20
      IC=I
      IF(SYM(I:I).EQ.':') GO TO 30
   10 CONTINUE
   20 LM=MIN0(LSYM,16)
      SYM1=SYM(1:LM)
      NSYM=1
      GO TO 40
   30 ICM=IC-1
      ICP=IC+1
      IF(SYM(ICP:ICP).EQ.' ') GO TO 104
      LM=MIN0(LSYM,ICP+15)
      SYM1=SYM(1:ICM)
   40 CALL GETINT(SYM1,INT1,IERR1)
      IF(NSYM.EQ.2) THEN
         SYM2=SYM(ICP:LM)
         CALL GETINT(SYM2,INT2,IERR2)
      ELSE
         IERR2=0
         INT2=INT1
         ENDIF
      IF(IERR1.GT.0.OR.IERR2.GT.0) IERR=1
      RETURN
  104 IERR=1
      RETURN
      END
      SUBROUTINE GETINT(SYM,INTG,IERR)
C
C     GET THE INTEGER VALUE FROM SYMBOL
C     RETURN IERR=1 IF IT'S NOT AN INTEGER
C
      CHARACTER SYM*(*)
      LSYM=LEN(SYM)
      I1=1
      IF(SYM(1:1).EQ.'+'.OR.SYM(1:1).EQ.'-') I1=2
      DO 10 I=I1,LSYM
      IF(SYM(I:I).EQ.' ') GO TO 20
      IF(SYM(I:I).LT.'0'.OR.SYM(I:I).GT.'9') GO TO 100
   10 CONTINUE
   20 READ(SYM,*) INTG
      IERR=0
      RETURN
  100 IERR=1
      RETURN
      END
      subroutine wpoly(a,n)
c
c     Write out a polynomial of order n with coefficients a(i).
c
      dimension a(0:n)
      character line*80,num*16,pow*2
      line=' '
      line(2:5)='y = '
      if (a(0).lt.0.) line(6:6)='-'
      ax=abs(a(0))
c 
      write(num,'(1p,e16.10)') ax
c 
      line(8:23)=num
      if(a(1).lt.0.) then
         line(31:31)='-'
      else
         line(31:31)='+'
         endif
      ax=abs(a(1))
      write(num,'(1p,e16.10)') ax
      line(33:48)=num
      line(50:50)='x'
      if(n.lt.2) go to 10
      if(a(2).lt.0.) then
         line(56:56)='-'
      else
         line(56:56)='+'
         endif
      ax=abs(a(2))
      write(num,'(1p,e16.10)') ax
      line(58:73)=num
      line(75:78)='x**2'
   10 write(6,2020) line
      write(12,2020) line
      nl=n/3
      i=2
      do 60 il=1,nl
      line=' '
      do 50 jl=0,2
      i=i+1
      ic=6+25*jl
      if(a(i).lt.0.) then
         line(ic:ic)='-'
      else
         line(ic:ic)='+'
         endif
      ic=ic+2
      icu=ic+15
      ax=abs(a(i))
      write(num,'(1p,e16.10)')  ax
      line(ic:icu)=num
      ic=icu+2
      icu=ic+4
      write(pow,'(i2)') i
      if(pow(1:1).eq.' ' ) then
         pow(1:1)=pow(2:2)
         pow(2:2)=' '
         endif
      line(ic:icu)='x**'//pow
      if (i.eq.n) go to 55
   50 continue
   55 write(6,2020) line
      write(12,2020) line
 2020 format(a80)
   60 continue
      do 70 i=1,80
      line(i:i)='-'
   70 continue
      write(6,2020) line
      write(12,2020)line
      return
      end
      subroutine sgefa(a, lda, n, ipvt, info)
c
c     sgefa factors a real matrix by gaussian elimination.
c
c     sgefa is usually called by sgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for sgeco) = (1 + 9/n)*(time for sgefa) .
c
c     on entry
c
c        a       real(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that sgesl or sgedi will divide by zero
c                     if called.  use  rcond  in sgeco for a reliable
c                     indication of singularity.
c
c     this version is written in cdc compass for the ftn compiler
c     communication convention.  results may vary slightly from those
c     obtained with a fortran version, and such a version is available
c     upon request.  this program simulates the fortran sequence..
c
c     subroutine sgefa (a, lda, n, ipvt, info)
      dimension a(lda,n), ipvt(n)
      info = 0
      if(n.gt.1) go to 10
      ipvt(1) = 1
      if(a(1,1).eq.0.) info = 1
      return
 10   nm1 = n - 1
      do 80 i = 1,nm1
        m = i
        ip1 = i + 1
        do 20 k = ip1,n
 20       if(abs(a(k,i)).gt.abs(a(m,i))) m = k
        ipvt(i) = m
        if (a(m,i).ne.0.) go to 25
        info = i
        go to 80
 25     if(m.eq.i) go to 40
        do 30 k = i,n
          t = a(i,k)
          a(i,k) = a(m,k)
 30       a(m,k) = t
 40     rp = -1./a(i,i)
        do 50 k = ip1,n
 50       a(k,i) = a(k,i)*rp
        do 70 j = ip1,n
          c = a(i,j)
          if(c.eq.0.) go to 70
          do 60 k = ip1,n
 60         a(k,j) = a(k,j) + a(k,i)*c
 70     continue
 80   continue
      ipvt(n) = n
      if(a(n,n).eq.0.) info = n
      return
      end
      subroutine sgesl(a,lda,n,ipvt,b,job)
c
c     sgesl solves the real system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by sgeco or sgefa.
c
c     on entry
c
c        a       real(lda, n)
c                the output from sgeco or sgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from sgeco or sgefa.
c
c        b       real(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if sgeco has set rcond .gt. 0.0
c        or sgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call sgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call sgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     this version is written in cdc assembly language for the ftn
c     compiler communication convention.  results may vary slightly
c     from those obtained with a fortran version, and such a version
c     is available upon request.  this program simulates the
c     following fortran sequence..
c
c     subroutine sgesl (a, lda, n, ipvt, b, job)
      dimension a(lda,n), b(n), ipvt(n)
      nm1 = n - 1
      if(n.gt.1) go to 10
      b(1) = b(1)/a(1,1)
      return
 10   if (job.ne.0) go to 70
      do 40 i = 1,nm1
        m = ipvt(i)
        c = b(m)
        if (m.eq.i) go to 20
        b(m) = b(i)
        b(i) = c
 20     if(c.eq.0.) go to 40
        ip1 = i + 1
        do 30 k = ip1,n
 30       b(k) = b(k) + a(k,i)*c
 40     continue
      do 60 i = 1,nm1
        j = n - i + 1
        b(j) = b(j)/a(j,j)
        c = b(j)
        if(c.eq.0.) go to 60
        jm1 = j - 1
        do 50 k = 1,jm1
 50       b(k) = b(k) - a(k,j)*c
 60     continue
      b(1) = b(1)/a(1,1)
      return
 70   do 90 k = 1,n
        t = 0.
        if (k.eq.1) go to 90
        km1 = k - 1
        do 80 i = 1,km1
 80       t = t + a(i,k)*b(i)
 90     b(k) = (b(k) - t)/a(k,k)
      do 110 kb = 1,nm1
        k = n - kb
        t = 0.
        kp1 = k + 1
        do 100 i = kp1,n
 100      t = t + a(i,k)*b(i)
        b(k) = b(k) + t
        l = ipvt(k)
        if (l.eq.k) go to 110
        t = b(l)
        b(l) = b(k)
        b(k) = t
 110    continue
      return
      end
      subroutine plotit(bb,n,infile)
      dimension bb(n)
      character*(*) infile
      character*80 line
      character*20 term
      character*4 ord
      open (14,file='gnuin-fit')
      if (n.eq.2) then
         ord=' 1st'
      else if (n.eq.3) then
         ord=' 2nd'
      else if (n.eq.4) then
         ord=' 3rd'
      else if (n.le.10) then
         write(ord,'(1x,i1,''th'')') n-1
      else
         write(ord,'(i2,''th'')') n-1
         endif
c
c   Uncomment the next statement if you are running with
c   a Tektronix 40xx emulation for graphics
c
c     write(14,*) 'set terminal tek40xx'
      line = 'set title '''//ord//' Order Fit to Data in file :'//infile//''''
      write(14,*) line
      line = 'plot '''//infile//''' using 1:2 t ''Data'',\\'
      ii=index(line,'\\')
      write(14,*) line(1:ii)
      line=' '
      write(line,'(1p,e14.6)') bb(1)
      iline=15
      do 30 j=2,n
         write(term,'(1p,e14.7)') bb(j)
         if(bb(j).ge.0.0) term(1:1)='+'
         term(15:18)='*x**'
         if (j.le.10) then
            write(term(19:20),'(i1,1x)') j-1
         else
            write(term(19:20),'(i2)') j-1
         endif
         line(iline:(iline+19))=term
         iline=iline+20
         if(iline.gt.50) then
            write(line(iline:(iline+1)),'(a2)') ' \\'
            ii=index(line,'\\')
            write(14,*) line(1:ii)
            line=' '
            iline=1
            endif
c
   30    continue
c
      write (line(iline:(iline+7)),'(a8)') ' t ''Fit'''
      write(14,*) line
      write(14,*) 'pause -1'
c
      close(14)
c
      call system ('gnuplot gnuin-fit')
      return
      end
c 
c c c