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
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
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.'
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
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