PROGRAM MAIN

*     A FORTRAN PROGRAM FOR THE LINEAR LEAST SQUARES FIT.
*     ACCEPTS DATA ON STDIN AND WRITES RESULTS ON STDOUT.
*     REFERENCE: JOHN P.G. SHEPHERD'S "DATA ANALYSIS", ED. 1.

*     ORIGINALLY WRITTEN APRIL 16, 1994 WHILE A STUDENT AT UW-RF.

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

*     INITIALIZE VARIABLES.

      PARAMETER( SMALL = 1E-6 )

      SUMX = 0.
      SUMY = 0.
      SUMXY = 0.
      SUMXSQ = 0.
      SUMYSQ = 0.
      NCOUNT = 0

*     READ THE DATA FROM STDIN.

      PRINT 900
 900  FORMAT( ' 1994 JAMES HALL' / )
      PRINT 901
 901  FORMAT( ' A FORTRAN PROGRAM FOR THE LINEAR LEAST SQUARES FIT.' )
      PRINT 902
 902  FORMAT( ' ACCEPTS DATA ON STDIN AND WRITES RESULTS ON STDOUT.' )

      PRINT 1000
 1000 FORMAT( / ' ENTER DATA:' / )
 10   READ (*,*,IOSTAT=IE) X,Y

*     UPDATE THE SUMS AND THE COUNTER.

      IF( IE.EQ.0 ) THEN
         SUMX = SUMX + X
         SUMY = SUMY + Y
         SUMXY = SUMXY + X * Y
         SUMXSQ = SUMXSQ + X**2
         SUMYSQ = SUMYSQ + Y**2
         NCOUNT = NCOUNT + 1
         GOTO 10
      ENDIF

*     CHECK FOR ENOUGH DATA PAIRS.

      IF( NCOUNT.LE.1 ) THEN
         PRINT *, 'NOT ENOUGH DATA.'
         STOP
      ENDIF

*     COMPUTE THE AVERAGE X AND Y.  DETERMINE DELTA AND CHECK THAT IT
*     IS NON-ZERO.

      XBAR = SUMX / NCOUNT
      YBAR = SUMY / NCOUNT
      DELTA = SUMXSQ - XBAR**2 * NCOUNT

      IF( DELTA.LT.SMALL ) THEN
         PRINT *, 'NOT ENOUGH VARIATION.'
         STOP
      ENDIF

*     PRINT SOME RESULTS AND GENERAL INFORMATION.

      PRINT 1010
 1010 FORMAT( / ' RESULTS:' / )

      PRINT 2000, NCOUNT, 'NUMBER OF PAIRS'
 2000 FORMAT( ' ', I4, T20, A )
      PRINT 2010, SUMX, 'SUM OF X'
      PRINT 2010, SUMY, 'SUM OF Y'
      PRINT 2010, SUMXY, 'SUM OF X*Y'
      PRINT 2010, SUMXSQ, 'SUM OF X**2'
      PRINT 2010, SUMYSQ, 'SUM OF Y**2'
 2010 FORMAT( ' ', G13.6, T20, A )

*     CALCULATE SLOPE AND INTERCEPT.  PRINT RESULTS.

      SLOPE = ( SUMXY - XBAR * YBAR * NCOUNT ) / DELTA
      YINT = YBAR - SLOPE * XBAR
      PRINT 2010, SLOPE, 'SLOPE'
      PRINT 2010, YINT, 'INTERCEPT'

*     CALCULATE THE VARIANCE IN THE COEFFICIENTS AND PRINT RESULTS.
*     THE VARIANCES SHOULD ALL BE POSITIVE, BUT WE'LL APPLY THE ABS()
*     FUNCTION ANYWAY TO AVOID PROBLEMS.

      VAR = ( SUMYSQ + SUMXSQ * SLOPE**2 - 2. * SLOPE * SUMXY ) /
     $     NCOUNT + YINT**2 + 2. * ( SLOPE * YINT * XBAR - YINT * YBAR )
      VARSLO = VAR / ( DELTA * NCOUNT )
      VARINT = VARSLO * SUMXSQ
      PRINT 2010, VARINT, 'VARIANCE OF INTERCEPT'
      PRINT 2010, SQRT( ABS( VARINT ) ), 'UNCERTAINTY OF INTERCEPT'
      PRINT 2010, VARSLO, 'VARIANCE OF SLOPE'
      PRINT 2010, SQRT( ABS( VARSLO ) ), 'UNCERTAINTY OF SLOPE'
      PRINT 2010, VAR, 'VARIANCE OF DATA'

      STOP
      END