c<html>
c<head><title>array1.f</title></head>
c<body>
c<pre>
program curvefit
c<a name="dp"><font color="FF0000">
implicit double precision (a-h,o-z)
c</a></font>
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</pre>
c</body>
c</html>
Copyright ©2017 - STI - Todos os direitos reservados