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