c carray1.f c c
      program curvefit
c
      implicit double precision (a-h,o-z)
c
c
c       Take points from specified data file and develop a Quadratic spline
c       fit. This assumes that the data file has an extension ".txt"
c
      parameter (nmax=20,ndata=1000)
      character*32 infile, arg2, data_name, outfile, sampfile, gnuin
      common/curf/ xd(ndata),yd(ndata),b(ndata), c(ndata), d(ndata),
     $ xs(1001),ys(1001),nd
      integer*4 narg
      integer*2 istatus
c
      WRITE(6,'(a)',advance='no')
     &   'PROVIDE INPUT FILE NAME :'
      READ (*,*) infile
      name_len=index(infile,'.')-1
      infile=infile(1:name_len)//'.txt'
      outfile=infile(1:name_len)//'.fit'
      sampfile=infile(1:name_len)//'.sam'
c     gnuin=infile(1:name_len)//'.plt'
      gnuin = 'quadfit.plt'
      open (11,file=infile)
      open (12,file=outfile)
      open (14,file=gnuin, position = 'append')
      open (16,file=sampfile)
      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
c
c   Generate splines
c
      call quad (nd, xd, yd, b, c)
c
c   Sample the interpolation table for a plot
c
      write(16,*) xd(1),yd(1)
      x=xd(1)
      dx=(xd(nd)-x)*.001
      do 100 i=2,1001
         x=x+dx
         call quadint (x,y,xd,yd,b,c,nd)
         write(16,*) x,y
 100  continue
c
c   Write information for data table
c
      do 200 i=1,nd
         write(12,2020) xd(i),yd(i),b(i),c(i)
 200  continue
 2020 format(1p,5x,'& ',4(e15.8,',') )
         call plotit(infile(1:name_len))
 500  stop
      end
      subroutine plotit(infile)
      implicit real*8 (a-h,o-z)
      character*(*) infile
      character*80 line
      write(14,*) 'set data style lines'
      write(14,*) 'set nokey'
      line = 'set title '//
     $ '''Quadratic Fit to Data for : '//infile//''''
      ii=index(line,':')
      ii=index(line(ii+2:80),' ')+ii
      write(14,'(a)') line(1:ii)
      line = 'plot '''//infile//'.sam'' using 1:2 , \\'
      ii=index(line,'\\')
      write(14,'(a)' ) line(1:ii)
      line = ' '''//infile//'.txt'' using 1:2 with points'
      ii=index(line,'ts')
      write(14,'(a)' ) line(1:ii+1)
      write(14,*) 'pause -1'
c      call system ('gnuplot gnuin-fit')
      return
      end
      subroutine quad (n, x, y, b, c)
      implicit real*8 (a-h,o-z)
      integer n
      real*8 x(n), y(n), b(n), c(n)
c
c  the coefficients b(i), c(i) i=1,2,...,n are computed
c  for a quadratic with continuous derivatives
c
c    s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2
c
c    for  x(i) .le. x .le. x(i+1)
c
c  input..
c
c    n = the number of data points or knots (n.ge.2)
c    x = the abscissas of the knots in strictly increasing order
c    y = the ordinates of the knots
c
c  output..
c
c    b, c     = arrays of quadratic coefficients as defined above.
c
      c(1)=((y(3)-y(1))/(x(3)-x(1))-(y(2)-y(1))/(x(2)-x(1)))/(x(3)-x(2))
      b(1)=(y(2)-y(1))/(x(2)-x(1))-c(1)*(x(2)-x(1))
      do 100 i=2,n-1
         b(i)=b(i-1)+2*c(i-1)*(x(i)-x(i-1))
         c(i)=((y(i+1)-y(i))/(x(i+1)-x(i))-b(i))/(x(i+1)-x(i))
  100 continue
         b(n)=0.
         c(n)=0.
c
c
      return
      end
c
      subroutine quadint(x,y,xtab,ytab,bcoef,ccoef,ntab)
      implicit real*8 (a-h,o-z)
c
c    Perform Quadratic interpolation on a table of data with
c    precomputed coefficients
c
      real*8 xtab(ntab),ytab(ntab),bcoef(ntab),ccoef(ntab)
      save ilast
      data ilast/1/

c    Start the search from the last point of table use index
c
      if (x.le.xtab(ilast+1)) then
c
c    Search down the table from point of last use
c
          do 20 i1=ilast,1,-1
              if(x.ge.xtab(i1)) go to 60
  20          continue
          write(6,*) 'x = ', x, '  is below the table range'  
          i1=1
          go to 60
      else
c
c    Search up the table from point of last use
c
          do 40 i1=ilast+1,ntab-1
              if(x.le.xtab(i1+1)) go to 60
  40          continue
          write(6,*) 'x = ', x, '  is above the table range'
          i1=ntab-1
          go to 60
      endif
c
c   Bounding points found, interpolate
c
  60  dx=(x-xtab(i1))
      y=ytab(i1)+bcoef(i1)*dx+ccoef(i1)*dx**2
      ilast=i1
      return
      end
c
c c