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
&   'PROVIDE INPUT FILE NAME :'
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'
open (11,file=infile)
open (12,file=outfile)
open (14,file=gnuin, position = 'append')
open (16,file=sampfile)
i=1
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
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
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