PROGRAM MAIN
C     LINEFIT.F
C     A PROGRAM TO COMPUTE THE BEST LINE TO DESCRIBE A SET OF X,Y
C     DATA ACCORDING TO LINEAR LEAST SQUARES
C     REFERENCE: 'DATA ANALYSIS' ED. 2, J.P.G. SHEPHERD

C     ORIGINALLY WRITTEN 1994, WHILE A STUDENT AT UW-RF

C     COPYRIGHT (C) 1994, 2007 JIM HALL, JHALL@FREEDOS.ORG
C Permission is hereby granted, free of charge, to any person obtaining a copy
C of this software and associated documentation files (the "Software"), to deal
C in the Software without restriction, including without limitation the rights
C to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
C copies of the Software, and to permit persons to whom the Software is
C furnished to do so, subject to the following conditions:
C 
C The above copyright notice and this permission notice shall be included in
C all copies or substantial portions of the Software.
C 
C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
C IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
C AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
C OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
C THE SOFTWARE.

      PARAMETER (MAXSIZ=300)
      REAL X(MAXSIZ), Y(MAXSIZ)

      N = 0
      XBAR = 0.0
      YBAR = 0.0
      SUMX2 = 0.0
      SUMY2 = 0.0
      SUMXY = 0.0

C     READ ALL VALUES

      PRINT *, 'ENTER X,Y VALUES (END TO QUIT):'
 10   READ (*,*,END=100) X0, Y0

C     THINGS TO COMPUTE:
C     N
C     AVERAGE OF X
C     AVERAGE OF Y
C     SUM OF X**2
C     SUM OF Y**2
C     SUM OF XY

      N = N + 1
      X(N) = X0
      Y(N) = Y0

C     COMPUTE THE SUMS

      SUMX = SUMX + X0
      SUMX2 = SUMX2 + X0**2
      SUMY = SUMY + Y0
      SUMY2 = SUMY2 + Y0**2
      SUMXY = SUMXY + X0 * Y0

C     EXIT THE LOOP EARLY ONLY IF WE HAVE FILLED THE ARRAY

      IF (N.EQ.MAXSIZ) GOTO 100
      GOTO 10

 100  CONTINUE

C     TO COMPUTE THE COEFFICIENTS:
C     M = (SUM(XY) - N * XBAR * YBAR) / DELTA
C     C = YBAR - M * XBAR
C
C     WHERE:
C     DELTA = SUM(X**2) - N * XBAR**2

      XBAR = SUMX / N
      YBAR = SUMY / N
      DELTA = SUMX2 - XBAR**2 * N

      SLOPE = (SUMXY - XBAR * YBAR * N) / DELTA
      YINTCP = YBAR - SLOPE * XBAR

C     TO COMPUTE THE VARIANCE OF THE COEFFICIENTS:
C     VAR(C) = VARNC * SUM(X**2)  / (N * DELTA)
C     VAR(M) = VARNC  / (N * DELTA) = N * VAR(C) / SUM(X**2)
C
C     WHERE:
C     DELTA = SUM(X**2) - N * XBAR**2       (SAME AS ABOVE)
C
C     AND:
C     VARNC = SUM((Y - M * X - C)**2) / N
C           = (SUM(Y**2) + M**2 * SUM(X**2)) / N + C**2
C                  + 2 * (M * C * XBAR - M * SUM(XY) / N - C * YBAR)

      VARNC = 0.0
      DO 110 I = 1, N, 1
         VARNC = VARNC + (Y(I) - SLOPE * X(I) - YINTCP)**2
 110  CONTINUE
      VARNC = VARNC / N

      VARNC2 = YINTCP**2 + (SUMY2 + SLOPE**2 * SUMX2) / N +
     $     2.0 * (SLOPE * YINTCP * XBAR - SLOPE * SUMXY / N -
     $     YINTCP * YBAR)

      VARYIN = VARNC * SUMX2 / (DELTA * N)
      VARSL = VARNC / (DELTA * N)
      VARSL2 = (VARYIN * N) / SUMX2

C     PRINT RESULTS

      PRINT *, 'NUMBER OF POINTS=', N
C     ** MY FIRST METHOD TO COMPUTE UNC IN SLOPE IS BAD
C      PRINT *, 'SLOPE=', SLOPE, '+/-', SQRT(VARSL), '(', VARSL, ')'
      PRINT *, 'SLOPE=', SLOPE, '+/-', SQRT(VARSL2), '(', VARSL2, ')'
C     $     , '= SECOND METHOD'
      PRINT *, 'YINTCP=', YINTCP, '+/-', SQRT(VARYIN), '(', VARYIN, ')'
      PRINT *, 'DELTA=', DELTA
      PRINT *, 'VARNC=', VARNC
      PRINT *, 'VARNC=', VARNC2, '= SECOND METHOD'
      PRINT *, 'XBAR=', XBAR
      PRINT *, 'YBAR=', YBAR

      END