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 *****