program fitmain
      implicit real*8 (a-h,o-z)
      parameter(maxpts=1000,maxcol=20,maxex=10)
      character ans,ans1,ans2,ans3,ans4,ans5,yes
      character*14 fn
      dimension x(maxpts),y(maxpts),xx(maxpts),yy(maxpts)
      dimension v(maxcol)
      dimension jxex(maxex),jyex(maxex),jex(maxex)
      yes='y'
      write(6,8)
    8 format(1x,'This is a least-squares fit program.')
  200 write(6,100)
  100 format(1x,'File name ?')
      call fflush(6)
      read(5,101) fn
  101 format(a14)
      open(1,file=fn,status='OLD')
      rewind(1)
      write(6,19)
   19 format(1x,'Set skipping parameters ?')
      call fflush(6)
      read(5,5) ans
      if(ans.eq.yes) then
         write(6,21) 
   21    format(1x,'icount, ibeg, and iskip ?')
         read(5,*) ict,ibeg,iskip
	 ict=min0(ict,maxpts)
	 if(ibeg.lt.0) ibeg=-ibeg
      else
	 ict=maxpts
	 ibeg=0
	 iskip=0
      endif
      write(6,167)
  167 format(1x,'set xcol, ycol ?')
      read(5,5) ans
      if(ans.eq.yes) then
  165    write(6,164)
  164    format(1x,'xcol, ycol ?')
         read(5,*) mx,my
         ncol=max0(mx,my)
         if(ncol.gt.maxcol) go to 165    
      else
         ncol=2
         mx=1
         my=2
      endif
c
c
c
      npts=0
      mpts=0
      jct=(ict-1)*(iskip+1)+ibeg+1
      do 30 i=1,jct
      if(i.le.ibeg.or.mod(i-ibeg-1,iskip+1).ne.0) then
	 read(1,5,end=31) ans
      else
         mpts=mpts+1
         read(1,*,end=31) (v(k),k=1,ncol)
         x(mpts)=v(mx)
         y(mpts)=v(my)
         npts=npts+1
      endif
   30 continue
   31 write(6,18)
   18 format(1x,'Exhibit data points ?')
      call fflush(6)
      read(5,5) ans
    5 format(a1)
      if(ans.eq.yes) then
         write(6,14)
   14    format(/' Data read in :'/)
         do 10 i=1,npts
         write(6,7) i,x(i),y(i)
    7    format(1x,i3,3x,1p,d12.5,3x,d12.5)
   10    continue
    4    write(6,50)
   50    format(///)
      endif
      write(6,61)
   61 format(1x,'Subtract offsets ?')
      read(5,5) ans5
   60 if(ans5.eq.yes) then
         write(6,62)
   62    format(1x,'xoff, yoff ?')
         read(5,*) xoff,yoff
      else
         xoff=0.d0
         yoff=0.d0
      endif
   70 write(6,71)
   71 format(1x,'Take log(x) ?')
      read(5,5) ans1
      if(ans1.ne.yes) then
         write(6,81)
   81    format(1x,'Take exp(x) ?')
         read(5,5) ans3
      endif
      jx=0
      if(ans1.eq.yes) then
         do 65 i=1,npts
         xt=dabs(x(i)-xoff)
         if(xt.ne.0.d0) then
            xx(i)=dlog10(xt)
         else
            jx=jx+1
            if(jx.gt.maxex) go to 60
            jxex(jx)=i
         endif
   65    continue
      elseif(ans3.eq.yes) then
         do 85 i=1,npts
         xx(i)=dexp(x(i)-xoff)
   85    continue
      else
         do 66 i=1,npts
         xx(i)=x(i)-xoff
   66    continue
      endif
      write(6,72)
   72 format(1x,'Take log(y) ?')
      read(5,5) ans2
      if(ans2.ne.yes) then
         write(6,83)
   83    format(1x,'Take exp(y) ?')
         read(5,5) ans4
      endif
      jy=0
      if(ans2.eq.yes) then
         do 67 i=1,npts
         yt=dabs(y(i)-yoff)
         if(yt.ne.0.d0) then
            yy(i)=dlog10(yt)
         else
            jy=jy+1
            if(jy.gt.maxex) go to 60
            jyex(jy)=i
         endif
   67    continue
      elseif(ans4.eq.yes) then
         do 87 i=1,npts
         yy(i)=dexp(y(i)-yoff)
   87    continue
      else
         do 68 i=1,npts
         yy(i)=y(i)-yoff
   68    continue
      endif
c
c
c
      jxy=0
      do 90 i=1,jx
      jxy=jxy+1
      jex(jxy)=jxex(i)
   90 continue
      do 95 i=1,jy
         do 96 j=1,jxy
         if(jyex(i).eq.jex(j)) go to 95
   96    continue
         jxy=jxy+1
         if(jxy.gt.maxex) go to 60
         jex(jxy)=jyex(i)
   95 continue
c
c
c
    3 write(6,40)
   40 format(1x,'Beginning and ending indices.')
      call fflush(6)
      read(5,*) mm,nn
      if(mm.le.0.or.nn.gt.npts) go to 3
      l=0
      do 97 i=mm,nn
      do 98 j=1,jxy
         if(i.eq.jex(j)) go to 97
   98 continue
      l=l+1
   97 continue
      if(l.lt.3) go to 3
      write(6,53) mm,nn,l,xoff,yoff
   53 format(1x,'Fit: #',i3,'-- #',i3,',',i4,' pts. '
     c             ,'offsets: ',1p,2(d12.5,1x))
c
c
c
      call fit(mm,nn,xx,yy,a,b,r,sa,sb,jex,jxy)
      write(6,15) b,sb,a,sa
   15 format(' fit: (',1p,e12.5,' +_',e12.5,
     c       ') x + (',e12.5,' +_',e12.5,').')
      write(6,16) r
   16 format(' *** correlation coefficient: ',e12.5,' ***')
      a0=10.d0**a
      a1=10.d0**(a+sa)
      a2=10.d0**(a-sa)
      write(6,17) a0,a2,a1
   17 format(' amplitude: ',1p,e12.5,' (',e12.5,'< A <',e12.5,').')
c
c
c
      write(6,52)
   52 format(/' Change range ?')
      call fflush(6)
      read(5,5) ans
      if(ans.eq.yes) go to 3
      write(6,63)
   63 format(' Change offsets ?')
      read(5,5) ans5
      if(ans5.eq.yes) go to 60
      write(6,69)
   69 format(' Change log/no log ?')
      read(5,5) ans
      if(ans.eq.yes) go to 70
      write(6,201)
  201 format(' Another file ?')
      call fflush(6)
      read(5,5) ans
      if(ans.ne.yes) stop
      close(1)
      go to 200
      end
c
c
c
      subroutine fit(m,n,x,y,a,b,r,sa,sb,jex,jxy)
      parameter(maxex=10)
      implicit real*8 (a-h,o-z)
      dimension x(1),y(1)
      dimension jex(maxex)
      l=0
      sumx=0.
      sumy=0.
      sumxy=0.
      sumx2=0.
      sumy2=0.
      do 1 i=m,n
      do 10 j=1,jxy
         if(i.eq.jex(j)) go to 1
   10 continue
      l=l+1
      sumx=sumx+x(i)
      sumy=sumy+y(i)
      sumxy=sumxy+x(i)*y(i)
      sumx2=sumx2+x(i)**2
      sumy2=sumy2+y(i)**2
    1 continue
      b=(dfloat(l)*sumxy-sumx*sumy)/(dfloat(l)*sumx2-sumx**2)
      a=(sumy-b*sumx)/dfloat(l)
      error=sumy2-sumy**2/dfloat(l)-b**2*(sumx2-sumx**2/dfloat(l))
      sa=dsqrt(error/dfloat(l-2))
      sb=sa*dsqrt(dfloat(l))/dsqrt(dfloat(l)*sumx2-sumx**2)
      r=dsign(dsqrt(1.-error/(sumy2-sumy**2/dfloat(l))),b)
      return
      end