c-----------------------------------------------------------------------
c Find the least squares linear fit
c   y = a + b*x
c
c Read (x,y) data from a file in a two-column format or from the keyboard
c
c Variables ...
c   x(i)   ... independent variable
c   y(i)   ... dependent variable
c   ipt    ... number of independent measurements
c   a      ... y-intercept
c   b      ... slope
c   nmax   ... maximum number of data points
c-----------------------------------------------------------------------
c Program Note:
c   Call a subroutine slope to find the slope and intercept.
c-----------------------------------------------------------------------
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

c Change the following parameter to handle more data points ------------
      parameter (nmax=1000)
      dimension x(nmax), y(nmax)
      character fname*40, resp*1, mresp*1

c Print program header -------------------------------------------------
      write(*,*)'************************************'
      write(*,*)'Find the least squares linear fit***'
      write(*,*)'  y = a + b*x                       '
      write(*,*)'************************************'
      write(*,*)

      write(*,600)' Read data from a file (Y/n)? '
        read(*,500)resp

c Read (x,y) manually from a terminal ----------------------------------
      if(resp .eq. 'n' .or. resp .eq. 'N')then
        write(*,*)'Enter x,y:'
        do 10 i=1, nmax
          write(*,608)' x(', i, ')= '
            read(*,*)x(i)
          write(*,608)' y(', i, ')= '
            read(*,*)y(i)
          write(*,600)' More data (Y/n)? '
            read(*,500)mresp
            if(mresp .eq. 'n' .or. mresp .eq. 'N')goto 11
10      continue
11      ipt=i
c Read data x, y from a file (may include a header) --------------------
      else
        write(*,600)' Filename containing x-y data in 2-column format: '
        read(*,500)fname
        open(1, file=fname, status='old')
        do 20 i=1, nmax
19        read(1,*,err=19,end=21)x(i), y(i)
20      continue
21      close(1)
        ipt=i-1
      endif

c Report the number of points read -------------------------------------
      write(*,*)' Number of points read = ',ipt

c Call a routine to find the slope and intercept
      call slope(x,y,ipt,a,b)

c Print the intercept and slope ----------------------------------------
      write(*,*)' intercept = ', a
      write(*,*)' slope =     ', b

c Some formats ---------------------------------------------------------
500   format(a)
600   format(a\)
608   format(a,i2,a\)
      end


c-----------------------------------------------------------------------
      subroutine slope(x,y,ipt,a,b)
c-----------------------------------------------------------------------
c Find the least squares linear fit
c   y = a + b*x
c
c Variables ...
c   x(i)   ... independent variable (input)
c   y(i)   ... dependent variable (input)
c   ipt    ... number of independent measurements (input)
c   a      ... y-intercept (output)
c   b      ... slope (output)
c   sumx   ... sum of x(i)
c   sumy   ... sum of y(i)
c   sumxx  ... sum of x(i)*x(i)
c   sumxy  ... sum of x(i)*y(i)
c   nmax   ... maximum number of data points
c-----------------------------------------------------------------------
c Program Note:
c   Calculate slope and intercept from the sum of x, y, xy, etc.
c-----------------------------------------------------------------------

c Declare arrays -------------------------------------------------------
      dimension x(1), y(1)

c Find the various sums needed for calculating slope and intercept -----

      sumx = 0.
      sumy = 0.

      do 10 i=1, ipt
        sumx  = sumx  + x(i)
        sumy  = sumy  + y(i)
        sumxx = sumxx + x(i)*x(i)
        sumxy = sumxy + x(i)*y(i)
10    continue

c Find the average of x and y ------------------------------------------
      xave = sumx / float(ipt)
      yave = sumy / float(ipt)

c Find the slope (b) and y-intercept (a)
      b = (sumxy - sumx*yave) / (sumxx - sumx*xave)
      a = yave - b*xave

      end