      PROGRAM LEASTQ
C
C     LEAST SQUARES STRIGHT LINE FIT  (MISN-0-359)
C     -----------------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U    1973)
C     Revised: Peter Signell (Mich. St. U.      3/86)
C     Revised: Gary Herzenstiel (Mich. St. U.   3/86)
C     Revised: E.J.D.Kales (Mich. St. U.       11/87)
C     Revised: Peter Signell (Mich. St. U.      3/02)
C
C     Find a two-parameter fit to a straight line
C
C     SUBROUTINES:
C          GETINF  - get values for N,X(1),Y(1),DY(1)
C          RESIDU  - calculate residuals A1, A2, and
C                    CHI-SQUARE
C          GRAPHI  - display the graph of X1 vs Y1
C          QUERY   - ask if user wants to try again
C
C                            L E A S T S Q
C                                  :
C            -----------------------------------------------
C            :                :             :             :
C         GETINF           RESIDU        GRAPHI       QUERY
C       (N,X,Y,DY)                                   (AGAIN)
C
C     ***** DECLARATIONS *****
      REAL N,X(100),Y(100),DY(100)
      REAL SXX,S1,CHI,XMAX,YMAX,YMIN,A1,A2,SX,XMIN
      CHARACTER AGAIN
      OPEN(UNIT=3,FILE='(T:81)m359p1f.out')
C                    
C     ***** M A I N   B O D Y *****
C
 10   WRITE(3,'(A1)') 12
      CALL GETINF  (                                        N,X,Y,DY)
      CALL RESIDU  (N,X,Y,DY,SXX,S1,SX,CHI,XMAX,XMIN,YMAX,YMIN,A1,A2)
      CALL GRAPHI  (N,X,Y,DY,SXX,S1,SX,CHI,XMAX,XMIN,YMAX,YMIN,A1,A2)
      CALL QUERY   (                                           AGAIN)
      IF(AGAIN.EQ.'Y') GOTO 10
C     else
      STOP
      END
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(    N,X,Y,DY)
C
C     gets values for N, the first set of coordinates,
C     and an error for the first point.
C
C     ***** DECLARATIONS *****
      REAL N,X(100),Y(100),DY(100)
      INTEGER M
      CHARACTER READFI*12
C
C     N   = number of data points
C     X,Y = coordinates of first point
C     DY  = error for the first point
C
      PRINT *,'The data will be read from your file'
      PRINT *,'It must be in the following order:'
      PRINT *,'First line has N, number of data, in I5 format.'
      PRINT *,'Succeeding lines have X(n),Y(n),DY(n) in 3F5.2'

      READFI='m359-p1f.dat'
      OPEN(UNIT=2,FILE=READFI)
      READ( 2,20) M
 20   FORMAT(I5)
      N=M
      READ( 2,22) (X(I),Y(I),DY(I),I=1,M)
 22   FORMAT(3F5.2)
      WRITE(*,24) N,X(1),Y(1),DY(1)
 24   FORMAT('m359p1f.for:    Least Squares Straight Line Fit:'/
     +'   N=',F7.3,' X(1)=',F7.3,' Y(1)=',F7.3,' DY(1)=',F7.3/)
      RETURN
      END
C
C     ***** R E S I D U *****
C
      SUBROUTINE RESIDU
     +           (N,X,Y,DY,SXX,S1,SX,CHI,XMAX,XMIN,YMAX,YMIN,A1,A2)
C
C     calculate residual and least squares parameters A1, 
C     A2, and CHI
C
C     ***** DECLARATIONS *****
      REAL N,X(100),Y(100),DY(100)
      REAL SXX,S1,CHI,XMAX,YMAX,YMIN,A1,A2,SX,XMIN
      REAL SIG,SY,SXY,Y1,R
      INTEGER M,J
C
C     S??? = summing variables
C     XMAX = maximum of all the values of X
C     XMIN = minimum of all the values of X
C     R    = residual
C     CHI  = chi-square
C
      IF(DY(1).GT.0) THEN
        SIG=DY(1)**2
      ELSE
        SIG=1.0
      ENDIF
      S1=1.0/SIG
      SX=X(1)/SIG
      SXX=X(1)**2/SIG
      SY=Y(1)/SIG
      SXY=X(1)*Y(1)/SIG
      XMAX=X(1)
      XMIN=X(1)
      YMAX=Y(1)
      YMIN=Y(1)
      M=N
C
C     loop to compute sums of all data points and errors
C
      DO 36 J=2,M
C
C       compute sums
C
        IF(DY(J).GT.0) THEN
          SIG=DY(J)**2
        ELSE
          SIG=1.0
        ENDIF
        S1=S1+1.0/SIG
        SX=SX+X(J)/SIG
        SXX=SXX+X(J)**2/SIG
        SY=SY+Y(J)/SIG
        SXY=SXY+X(J)*Y(J)/SIG
C
C       find maximum and minimum values of X and Y
C
        IF(X(J).LT.XMIN) XMIN=X(J)
        IF(X(J).GT.XMAX) XMAX=X(J)
        IF(Y(J).LT.YMIN) YMIN=Y(J)
        IF(Y(J).GT.YMAX) YMAX=Y(J)
 36   CONTINUE
      WRITE(3,38) M
 38   FORMAT(//'m359p1f.for: Least Squares Straight Line Fit:'//
     +'Number of data points =',I5//
     +12X,'X',12X,'Y',11X,'DY',10X,'Residuals'/)
C
C     COMPUTE PARAMETERS A1 AND A2
C
      A1=(SXX*SY-SX*SXY)/(S1*SXX-SX*SX)
      A2=(SX*SY-S1*SXY)/(SX*SX-SXX*S1)
C
C     compute R (residuals) and CHI (chi-square)
C
      CHI=0.0
      DO 44,J=1,M
        Y1=A1+A2*X(J)
        R=Y(J)-Y1
        IF(DY(J).GT.0) THEN
          SIG=DY(J)**2
        ELSE
          SIG=1.0
        ENDIF
        CHI=CHI+R**2/SIG
        WRITE(3,42) X(J),Y(J),DY(J),R
 42     FORMAT(1X,4F14.5)
 44   CONTINUE
      RETURN
      END
C
C     ***** G R A P H I *****
C
      SUBROUTINE GRAPHI
     +       (N,X,Y,DY,SXX,S1,SX,CHI,XMAX,XMIN,YMAX,YMIN,A1,A2)
C
C     calculates and displays uncertainties DA1 and DA2
C     and plots Y1 against X1 using fitted X values
C
C     ***** DECLARATIONS *****
      REAL N,X(100),Y(100),DY(100)
      REAL SXX,S1,CHI,XMAX,YMAX,YMIN,A1,A2,SX,XMIN
      REAL DA1,DA2,RAV,DY1,DX,Y1,X1,II,XSCLE
      INTEGER J,K,M,M1,SETM2,M2,RMAR
      CHARACTER LINE*101,LINCHR*1
C
      M=N
C
C     compute uncertainties DA1 and DA2
C
      DA1=SQRT(SXX/(S1*SXX-SX*SX))
      DA2=SQRT(S1/(S1*SXX-SX*SX))
C
C     if DY(1) = ZERO, assume errors are unknown
C
      IF(DY(1).EQ.0)THEN
        RAV=SQRT(CHI)/N
        DA1=DA1*RAV
        DA2=DA2*RAV
      ENDIF
      LINCHR=' '  
      XSCLE=6.0
      RMAR=61
C
C     DX and DY1 are step sizes used to make the plot
C
      DX=1.2*(XMAX-XMIN)/50.0
      DY1=1.2*(YMAX-YMIN)/(RMAR-1.)
C
C     Within the dashed lines is a Kludge to attempt to make
C     the *'S form a straight line.  Remove if troublesome:
C     ----------------------------------------------------
      II=INT(A2*DX/DY1+0.49999)
      DX=II*DY1/A2
C     ----------------------------------------------------
C
C     set first X value for plot
C
      X1=XMIN-5.0*DX
C
      LINE='
     +                                                                 '
C
      WRITE(3,54) A1,DA1,A2,DA2,CHI,12
 54   FORMAT(//'Least Squares Fit Parameters:'//
     +       '    A1 =',F10.5,'  +/-',F10.5/
     +       '    A2 =',F10.5,'  +/-',F10.5/
     +       ' CHISQ =',F10.5/
     +A1///'      X       Y     ')
      WRITE(3,58)
 58   FORMAT(19X,50H..................................................,
     +11H...........)
C
C     loop to display graph
C
      DO 70,J=1,51
        LINE(1:1)='.'
        LINE(RMAR:RMAR)='.'
C
C       Y1 is the fitted Y value
C
        Y1=A1+A2*X1
        M1=(Y1-YMIN)/DY1+XSCLE
C
C       if M1 is between 1 AND RMAR, store it for the
C       fitted point
C
        IF(M1.GE.1.AND.M1.LE.RMAR) LINE(M1:M1)='*'
C
        SETM2=0
        DO 60,K=1,M
C
C         if this X1 is within DX of any X(K),K=1,N then
C         store "O" 
C
          IF(ABS(X1-X(K)).LE.DX/2.0) THEN
            M2=(Y(K)-YMIN)/DY1+XSCLE
            LINE(M2:M2)='O'
            SETM2=1
          ENDIF
 60     CONTINUE
C
        WRITE(3,64) X1,Y1,LINE(1:61)
 64     FORMAT(1X,2F8.3,2X,A61)
C
C       Restore LINE to blanks for next plot
C
        IF(M1.GE.1.AND.M1.LE.RMAR) LINE(M1:M1)=' '
        IF(SETM2.EQ.1) LINE(M2:M2)=' '
C
C       do next fitted X value
C
        X1=X1+DX
 70   CONTINUE
      WRITE(3,58)
      RETURN
      END
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C
C     asks if user would like to do another
C
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      PRINT *
      WRITE(*,92)
 92   FORMAT(
     +'Would you like to do this again using new values? (Y/N)? ')
      READ( *,'(A1)')AGAIN
      RETURN
      END
C     ***** P R O G R A M    E N D *****