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

' Change the following parameter to handle more data points ------------
      const nmax=1000
      dim x(nmax), y(nmax)

' Print program header -------------------------------------------------
      print "************************************"
      print "Find the least squares linear fit***"
      print "  y = a + b*x                       "
      print "************************************"
      print

      input "Read data from a file (Y/n)? ", resp$

' Read (x,y) manually from a terminal ----------------------------------
      if resp$="n" or resp$="N" then
        print "Enter x,y:"
        for i=1 to nmax
          i$ = ltrim$(str$(i))            'no extra spaces
          print "x("; i$; ") = ";
          input "", x(i)
          print "y("; i$; ") = ";
          input "", y(i)
          input "More data (Y/n)? ", mresp$
          if mresp$="n" or mresp$="N" then exit for
        next
        ipt=i

' Read data x, y from a file (may include a header) --------------------
      else
        input "Filename containing x-y data in 2-column format: ", filename$
        open filename$ for input as #1
        for i=1 to nmax
          if eof(1) then exit for
          input #1, x(i), y(i)
        next
        close #1
        ipt=i-1
      endif

' Report the number of points read -------------------------------------
      print "Number of points read = "; ipt

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

' Print the intercept and slope ----------------------------------------
      print "intercept = "; a
      print "slope =     "; b


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

' Declare arrays -------------------------------------------------------
'     dim x(1), y(1)

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

      sumx = 0.
      sumy = 0.

      for i=1 to ipt
        sumx  = sumx  + x(i)
        sumy  = sumy  + y(i)
        sumxx = sumxx + x(i)*x(i)
        sumxy = sumxy + x(i)*y(i)
      next

' Find the average of x and y ------------------------------------------
      xave = sumx / ipt
      yave = sumy / ipt

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

      end sub