C*******************************************************************************
C
C GOA_LST_SQR generates the Least Squares Best Fit line through a series of
C (x,y) pairs read from a file.
C
C Technique
C =========
C The program reads data pairs (x,y) from a file specified by the user,
C computes the gradient and intercept of the least squares best-fit line
C through the dataset from various statistical quantities described below.
C
C The line is described by it's gradient and intercept as follows:
C y=mx+c; Where: m = gradient; c = intercept;
C x = independent variable; y = dependent variable
C
C It then outputs this least best fit line in form of a plot in the x-y plane
C using the plotting program 'gnuplot'
C
C It also outputs various statistical data such as:
C + Mean of X data, Mean of Y Data
C + Variance of X data, Variance of Y Data
C + Covariance and Correlation Coefficient for the dataset
C + Mean Square Deviation
C + Gradient and Intercept of Least Square Best-Fit Line
C + Uncertainties in the Gradient and Intercept of the Best-Fit Line
C + Equation of the Least-Squares Best Fit line
C
C ALL results are displayed on the screen and saved to file.
C
C
C Options
C =======
C Other options available in the program include:
C + Re-plotting of the Least-Squares Plot
C + Saving the results/plot to a file
C + Printing the results and plot
C + Restarting the program (running more than one analysis without
C having to exit the program)
C
C Various Tests were conducted to ensure that the right file, indices and data
C were used in the analyses.
C
C Error checks were inserted to transfer control to specific locations whenever
C errors are encountered.
C
C
C NOTES
C =====
C - The GNUPLOT program must be available for the graphical output to be displayed.
C - Data entries are in the form x: Column 1
C y: Column 2
C - Columns are separated by a tab or a space.
C
C***********************************************************************************
C
C User Defined Parameters:
C InputFile: Name of file containing the (x,y) data
C X: X Data stored in the 1st column of the file
C Y: Y Dtat stored in the 2nd column of the file
C X_Label: Unit/Label of the X Data (Displayed on plot)
C Y_Label: Unit/Label of the Y Data (Displayed on plot)
C Title: Title of the Dataset (Displayed on plot)
C N_Lower: Index of 1st pair used in the least squares calculation
C N_Upper: Index of last pair used in the least squares calculation
C
C Computed Variables:
C N: Number of Data pairs used in the least squares calculation
C X_Mean: Mean of the X Data
C Y_Mean: Mean of the Y Data
C X_Variance: Variance in the X-Data
C Y_Variance: Variance in the Y-Data
C Cov: Covariance of the dataset
C CorCoeff: Correlation coefficient between the X & Y-Data
C MnSqDev: Mean Square Deviation of the Dataset
C Gradient: Gradient of the computed least squares best-fit line
C Intercept: Intercept of the computed least squares best-fit line
C Uncert_Grad: The uncertainty in the Gradient
C Uncert_Intrcpt: The uncertainty in the Intercept
C
C Other Variables:
C N_Max: Maximum Number of dataset (1000)
C Err1: Variable to determine if Error exists in the Input file
C NDATA: Actual Number of data pairs in the file
C OutputFile: Name of the output File
C RUN: Variable to check if user wants to run another analysis
C
C*******************************************************************************
C
C Last Modified: 28 January 2006
C
C Author: AYENI GBOYEGA
C
C***********************************************************************************
C
C Disclaimer:
C ==========
C This program may be used for simple scientific operations involving least squares.
C But the author would take no liability whatsoever for damages incured from the
C use of the program
C
C Copyright:
C ==========
C All parts of this program may be copied and modified in part or whole for academic
C and research purpose.
C Use of the program for any commercial purpose must be done with a written permission
C from the author.
C
C
C Comments/Enquiries: [email protected]
C
C**********************************************************************************
C

C Program Definition
PROGRAM GOA_LST_SQR

C Implicit Statement
IMPLICIT NONE

C Declaration of variables
CHARACTER InputFile*30,X_Label*30,Y_Label*30,Title*30,
* OutputFile*30,RUN*1
INTEGER N_Max,N,N_Lower,N_Upper,Err1,NDATA
PARAMETER(N_Max=1000)
REAL Gradient,Intercept,MnSqDev,Uncert_Grad,Uncert_Intrcpt,Cov,
* CorCoeff,X(N_Max),Y(N_Max),X_Mean,Y_Mean,X_Variance,Y_Variance

C Common block indicating variables accessible to other subroutines
COMMON/VAL1/N_Lower,N_Upper,Title,X_Label,Y_Label
COMMON/VAL3/N,Gradient,Intercept
COMMON/VAL4/Err1,NDATA
COMMON/VAL5/MnSqDev,Uncert_Grad,Uncert_Intrcpt,Cov,CorCoeff,X,Y,
* X_Mean,Y_Mean,X_Variance,Y_Variance,InputFile,OutputFile

C Welcome Statement/Instructions
WRITE(*,*)'WELCOME TO GOA_LST_SQR'
WRITE(*,*)'AUTHOR: AYENI GBOYEGA'
WRITE(*,*)'JANUARY 2006'
WRITE(*,*)
WRITE(*,*)'-----------------------------------------------------'
WRITE(*,*)'THIS PROGRAM COMPUTES THE LEAST SQUARES BEST FIT LINE'
WRITE(*,*)'THROUGH A SERIES OF (X,Y) PAIRS READ FROM A FILE'
WRITE(*,*)
WRITE(*,*)'TO USE THE PROGRAM, FOLLOW INSTRUCTIONS ON THE SCREEN'
WRITE(*,*)
WRITE(*,*)'Comments/Enquiries: [email protected]'
WRITE(*,*)'======================================================'

C Call Subroutine USERINPUT to get all the user defined Variables
1 CALL USERINPUT(InputFile,X_Label,Y_Label,Title,N_Lower,N_Upper)

C Set variable Run to 'N', i.e. No.
C (This would be used later to determine whether to run the program after each operation)
RUN='N'

C Call Subroutine GETDATA to read (x,y) data from the specified Input file
CALL GETDATA(X,Y,N_Lower,N_Upper,N,InputFile)

C If there was an error in the specified Input file, then return to beginning of file
C and start program again (Err1=1 when there was an error in reading the file)
IF(Err1.EQ.1) then
REWIND(1)
GOTO 1
END IF

C Call Subroutine Mean() to compute the mean of the x-data "X_Mean"
CALL MEAN(X,N,X_Mean)

C Call Subroutine Mean() to compute the mean of the y-data "Y_Mean"
CALL MEAN(Y,N,Y_Mean)

C Call Subroutine VAR() to compute the variance of the x-data "X_Variance"
CALL VAR(X,N,X_Mean,X_Variance)

C Call Subroutine VAR() to compute the variance of the y-data "Y_Variance"
CALL VAR(Y,N,Y_Mean,Y_Variance)

C Call Subroutine COVAR() to compute the covariance "Cov" between the x- & y- data
CALL COVAR(X,Y,N,X_Mean,Y_Mean,Cov)

C Call Subroutine CORRCOEF() to compute the correlation coefficient "CorCoeff"
C between the x- & y- data
CALL CORRCOEF(X_Variance,Y_Variance,Cov,CorCoeff)

C Call Subroutine GRAD() to compute the gradient "Gradient" of the least squares
C best-fit line
CALL GRAD(Cov,X_Variance,Gradient)

C Call Subroutine INTRCPT() to compute the intercept "Intercept" of the least squares
C best-fit line
CALL INTRCPT(X_Mean,Y_Mean,Gradient,Intercept)

C Call Subroutine MSQDEV() to compute the Mean Squared Deviation "MnSqDev" of the
C least squares best-fit line
CALL MSQDEV(X,Y,Gradient,Intercept,N,MnSqDev)

C Call Subroutine UNCTGRAD() to compute the uncertainty in the gradient of the
C least squares best-fit line "Uncert_Grad"
CALL UNCTGRAD(MnSqDev,N,X_Variance,Uncert_Grad)

C Call Subroutine UNCTINT() to compute the uncertainty in the intercept of the
C least squares best-fit line "Uncert_Intrcpt"
CALL UNCTINT(Uncert_Grad,X_Variance,N,X_Mean,Uncert_Intrcpt)

C Call Subroutine GNUFILE() to save the x- & y-data used in the analysis (i.e.
C from 1st to last indices specified by user); y-values for the best fit line
C and the equation of the best fit line into appropriate files.
CALL GNUFILE(X,Y)

C Call Subroutine RESULT_OUT() to screen and to file "OutputFile"
CALL RESULT_OUT()

C Call Subroutine GNU_PLOT() to plot the (x.y) pairs and the least squares best
C fit line through the data pairs
CALL GNU_PLOT()

C Close all open files
CLOSE (1)
CLOSE (7)
CLOSE (8)
CLOSE (9)
CLOSE(11)

C Ask the user whether to run the program again. Exit program except when Run='Y'
C Recall default was Run = 'N' at beginning of the program
WRITE(*,*)
WRITE(*,*)'======================================='
WRITE(*,*)'Would You like to run another Analysis?'
WRITE(*,*)'Y=Yes'
WRITE(*,*)'Any other key to exit'
WRITE(*,*)'======================================='
WRITE(*,*)

C Get the user's response
READ(*,'(A1)')RUN

C Transfer control to Statement labelled 1 if Run='Y'
IF(RUN.EQ.'Y'.OR.RUN.EQ.'y') GOTO 1

C Quit program
STOP
END


C END OF MAIN PROGRAM




C*******************************************************************************
C Subroutine USERINPUT
C This subroutine is used to acquire all the necessary User-defined variables
C including InputFile, X_Label, Y_Label, Title, N_Lower, N_Upper (all defined
C at the beginning of main program).
C
C Local Variables
C =================
C EXIST: Logical variable used to check if the Input File Exists

SUBROUTINE USERINPUT(InputFile,X_Label,Y_Label,Title,
* N_Lower,N_Upper)

C Implicit Statement
IMPLICIT NONE

C Variable Declaration
CHARACTER InputFile*30,X_Label*30,Y_Label*30,Title*30
INTEGER N_Lower,N_Upper
LOGICAL EXIST

C Write instructions
WRITE(*,*)
WRITE(*,*)'Please Enter the following parameters'

C Get the Input Filename
2 WRITE(*,*)'Input Filename (+ extension):'
READ(*,'(A30)')InputFile

C Check if File Exists
INQUIRE(FILE=InputFile, EXIST=EXIST)

C If Input File does not exist in the current directory, output
C an error message and go back to read the input Filename again
IF(.NOT. EXIST) THEN
WRITE(*,*)'Input file does not exist. Please try again.'
WRITE(*,*)'Ensure that File is saved in your current path'
WRITE(*,*)
GOTO 2
ELSE
C Else if the Input File exists, then Open the File
OPEN(1,STATUS='OLD',FILE=InputFile)
ENDIF

C Get the Label for the X-axis
WRITE(*,*)
WRITE(*,*)'X LABEL: (Maximum 30 Characters)'
READ(*,'(A30)')X_Label

C Get the Label for the Y-axis
WRITE(*,*)
WRITE(*,*)'Y LABEL: (Maximum 30 Characters)'
READ(*,'(A30)')Y_Label

C Get the Title of the Analysis
WRITE(*,*)
WRITE(*,*)'TITLE: (Maximum 30 Characters)'
READ(*,'(A30)')Title

C Get the Index of the First (x,y) pair to use in the analysis
C If there is an error while reading the Index, go to statement labelled 3
WRITE(*,*)
4 WRITE(*,*)'Index of First pair to read'
READ(*,'(I3)',ERR=3)N_Lower

C Get the Index of the Last (x,y) pair to use in the analysis
C If there is an error while reading the Index, go to statement labelled 3
WRITE(*,*)
WRITE(*,*)'Index of Last pair to read'
READ(*,'(I3)',ERR=3)N_Upper

C****Check the indices
C If the Index of the first (x,y) pair is greater than or equal
C to the index of the last (x,y) pair then there is an error!!
C Go to Statement 3.
IF(N_Lower.GT.N_Upper.OR.N_Lower.EQ.N_Upper) then
WRITE(*,*)'Fisrt Index:',N_Lower,'Last Index:',N_Upper
WRITE (*,*)
GOTO 3
ELSE

C Else, everything is OK......transfer control to statement labelled 5.
GOTO 5
END IF


C Output reason(s) for error and give probable solutions.
3 WRITE(*,*)
WRITE(*,*)'Error: Invalid Index(ices)!!!'
WRITE(*,*)'Indices MUST be an Integer AND'
WRITE(*,*)'Index of Last pair MUST be greater than Index
*of the 1st pair'
WRITE(*,*)

C then go back to read the indices.
GOTO 4

C Return to main program.
5 RETURN
END




C*******************************************************************************
C Subroutine GETDATA is used is used to acquire all (x,y) data pairs from
C the specified Input file

SUBROUTINE GETDATA(X,Y,N_Lower,N_Upper,N,InputFile)

C Variable Declaration
INTEGER N,NDATA,Err1
CHARACTER InputFile*30
PARAMETER (N_Max=100)
REAL X(N_Max),Y(N_Max)

C Common Block
COMMON/VAL4/Err1,NDATA

C Set Err1 to zero (Err1 remains zero except if an error is encountered
C while reading the file)
Err1=0

C Set N (Number of Data points) to zero
N=0

C Read all the (x,y) data pairs in the input file until the end of file.
C Increment number of data points after reading each (x,y) pair in the file
C Transfer control to Statement 12 If Error is encountered
C Transfer control to Statement 11 If end of file is encountered
DO I=1,N_Max
READ(1,*, END=11, ERR=12)X(I),Y(I)
N=N+1
END DO

C Store actual number of data pairs in the file in the variable NDATA
C This is important because the number to be used for the Least Squares
C computation might differ

11 NDATA=N

C If 1st Index is not equal to 1, or Last Index is not equal to the total
C number of data pairs, then the user would not be using the whole dataset
IF(N_Lower.NE.1.OR.N_Upper.NE.N)THEN
C
C Therefore a new data range must be defined.
C
C Set a new counter to zero
L=0
C Set the data pairs from 1st to last indices as the new set of data.
C Count the new number of data pairs to be used in the analysis
DO J=N_Lower,N_Upper
L=L+1
X(L)=X(J)
Y(L)=Y(J)
END DO

C Replace the old number of data pairs with the new.
N=L

END IF

C If the IF statement is FALSE, the block above is not executed because
C the whole dataset in file would be used in the computations
C (so return to to main program)
C
C Jump to Statement Label 13 (ie return to main program)
GOTO 13

C Error block (Only Executed if there is an error while reading the input file
C Give advice on probable cause(s) of the error
12 WRITE(*,*)'---------------------------------------------------'
WRITE(*,*)'An Error occurred while trying to read the data'
WRITE(*,*)'HINTS:'
WRITE(*,*)'Make sure Input file exists and data properly saved'
WRITE(*,*)'Save the x and y-values in columns'
WRITE(*,*)'Separate each x-y pair with a TAB'
WRITE(*,*)'No Heading should be present at the top'
WRITE(*,*)'---------------------------------------------------'

C Set Err1 to 1 (this sends a message to the main program that an error was
C encountered while reading the input file). Only executed if there is an error.
Err1=1

C Return to Main program
13 RETURN

END




C*******************************************************************************
C Subroutine MEAN() returns the means for each of the two data columns x- & y-
C The subroutine computes the mean of the dataset 'A' ("A" corresponds to either
C the x- or y- data depending on the entry into the call statement)

SUBROUTINE MEAN(A,NA,A_MEAN)

C Variable declaration
INTEGER NA
REAL A(NA),A_MEAN, A_SUM

C Set sum of entries "A_SUM" to zero
A_SUM=0.

C Add up all the entries in the data column
DO I=1,NA
A_SUM=A_SUM+A(I)
END DO

C Mean of the data is equal to the Sum of the data divided by number of data points
A_MEAN=A_SUM/NA

C Return to the main program
RETURN
END




C*******************************************************************************
C Subroutine VAR() returns the variance for each of the two data columns x- & y-
C The subroutine computes the variance of the dataset 'A' ("A" corresponds to either
C the x- or y- data depending on the entry into the call statement)

SUBROUTINE VAR(A,NA,A_MEAN,A_VAR)

C Variable declaration
INTEGER NA
REAL A(NA),A_MEAN,A_VAR,VAR_SUM

C Set the Numerator (Summation) in the variance equation to zero
VAR_SUM=0.

C Set the Variance to zero (Not important...used as a control while
C building the program)
A_VAR=0.

C Compute the Numerator (Summation) in the variance equation
DO I=1,NA
VAR_SUM=(VAR_SUM+((A(I)-A_MEAN)**2))
END DO

C Compute the variance of the data column
A_VAR=(VAR_SUM/(NA-1))

C Return to main program
RETURN
END


C******************************************************************************
C Subroutine COVAR() returns the covariance between the two data columns x- & y-
C The subroutine computes the covariance of the dataset 'A' and 'B' ("A" & "B"
C correspond to x- & y- data resp.)

SUBROUTINE COVAR(A,B,N,A_MEAN,B_MEAN,Cov)

C Variable declaration
INTEGER N
REAL A(N),B(N),A_MEAN,A_VAR,B_MEAN,Cov,COV_SUM

C Set the Numerator (Summation) in the covariance equation to zero
COV_SUM=0.

C Compute the Numerator (Summation) in the covariance equation
DO I=1,N
COV_SUM=COV_SUM+((A(I)-A_MEAN)*(B(I)-B_MEAN))
END DO

C Compute the covariance between the two data columns

Cov=COV_SUM/(N-1)

C Return to the main program
RETURN
END


C*******************************************************************************
C Subroutine CORRCOEF() returns the correlation coefficient between the two data
C columns x- & y- .

SUBROUTINE CORRCOEF(A_VAR,B_VAR,Cov,CorCoeff)

C Variable declaration
REAL A_VAR,B_VAR,Cov,CorCoeff

C Compute the correlation coefficient between the two data columns
C A_VAR & B_VAR correspond to the variances in A nad B respectively
CorCoeff=Cov/((A_VAR*B_VAR)**0.5)

C Return to the main program
RETURN
END


C*******************************************************************************
C Subroutine GRAD() returns the gradient of the least squares best fit line through
C the (x,y) pairs

SUBROUTINE GRAD(Cov,A_VAR,Gradient)

C Variable declaration
REAL Cov,A_VAR,Gradient

C Compute the gradient of the least squares best fit line
Gradient=Cov/A_VAR

C Return to the main program
RETURN
END

C*******************************************************************************
C Subroutine INTRCPT() returns the intercept of the least squares best fit line through
C the (x,y) pairs
SUBROUTINE INTRCPT(A_MEAN,B_MEAN,Gradient,Intercept)
REAL A_MEAN,B_MEAN,Gradient,Intercept

C Compute the Intercept of the least squares best fit line
Intercept=B_MEAN-Gradient*A_MEAN

C Return to the main program
RETURN
END


C*******************************************************************************
C Subroutine MSQDEV() returns the mean squared deviation of the the best fit line through
C the (x,y) pairs
SUBROUTINE MSQDEV(A,B,Gradient,Intercept,N,MnSqDev)

C Variable declaration
INTEGER N
REAL A(N),B(N),Gradient,Intercept,MnSqDev,MnSqDev_Sum

C Set the Numerator (Summation) in the mean squared deviation equation to zero
MnSqDev_Sum=0.

C Compute the Numerator (Summation) in the mean squared deviation equation
DO I=1,N
MnSqDev_Sum=MnSqDev_Sum+((B(I)-Gradient*A(I)-Intercept)**2)
END DO

C Compute the mean squared deviation of the best fit line
MnSqDev=MnSqDev_Sum/(N-2)

C Return to the main program
RETURN
END

C*******************************************************************************
C Subroutine UNCTGRAD() returns the uncertainty in the gradient "Uncert_Grad" of
C the best fit line through the (x,y) pairs

SUBROUTINE UNCTGRAD(MnSqDev,N,A_VAR,Uncert_Grad)

C Variable declaration
INTEGER N
REAL Uncert_Grad,MnSqDev,A_VAR,UASQ

C Compute the uncertainty in the gradient of the best fit line
Uncert_Grad=(MnSqDev/((N-1)*A_VAR))**0.5

C Return to the main program
RETURN
END

C*******************************************************************************
C Subroutine UNCTINT() returns the uncertainty in the intercept "Uncert_Intrcpt" of
C the best fit line through the (x,y) pairs

SUBROUTINE UNCTINT(Uncert_Grad,A_VAR,N,A_MEAN,Uncert_Intrcpt)

C Variable declaration
INTEGER N
REAL Uncert_Intrcpt,MnSqDev,A_VAR,UBSQ

C Compute the uncertainty in the intercept of the best fit line
Uncert_Intrcpt=(Uncert_Grad*(A_VAR*(1-(1/N))+A_MEAN**2))**0.5

C Return to the main program
RETURN
END


C*******************************************************************************
C Subroutine UNCTINT() writes the (x,y) pairs used in the analysis, the (x,y2) values
C for the best fit line, and the equation of the best fit line into appropriate files
C
C
C Variable(s) not yet defined
C =================
C EQ: character variable used to store the equation of the best fit line
C Y2: data column used to store the y- values for the best fit line
C
SUBROUTINE GNUFILE(X,Y)

C Variable declaration
INTEGER N
CHARACTER EQ*17
REAL Gradient, Intercept, Y2(N), X(N), Y(N)


C Common Block
COMMON/VAL2/EQ
COMMON/VAL3/N,Gradient,Intercept

C Open output files "RawData", "BestFit", "EQN", & "DATA"
OPEN (UNIT=7, FILE='RawData', STATUS='UNKNOWN')
OPEN (UNIT=8, FILE='BestFit', STATUS='UNKNOWN')
OPEN (UNIT=9, FILE='EQN', STATUS='UNKNOWN')
OPEN (UNIT=11, FILE='DATA', STATUS='UNKNOWN')

DO I=1,N

C Compute the Y2 values for the best fit line using the same X values
C y=mx+c (already defined in main program)
Y2(I)=Gradient*X(I)+Intercept

C Write the (x,y) pairs to file: "RawData"
WRITE(7,*)X(I),Y(I)

C Write the (x,y2) pairs to file: "BestFit"
WRITE(8,*)X(I),Y2(I)

C Write the y2 values to file: "DATA"
WRITE(11,*)Y2(I)
END DO


C Write the equation of the best fit line to file: "EQN"
WRITE(9,'(A2,F6.2,A1,A1,F6.2)')'Y=',Gradient,'X','+',Intercept
REWIND (9)

C Read equation from file "EQN" and store in variable "EQ"
READ(9,'(A17)')EQ

C Return to the main program
RETURN
END

C*******************************************************************************
C Subroutine RESULT_OUT() outputs the results and variables to the screen and to a file.
C
SUBROUTINE RESULT_OUT()

C Variable declaration
CHARACTER File1*30,XLab*30,YLab*30,Title*30,OutputFile*30,EQ*17
INTEGER NMax,N,NLow,NUpp
PARAMETER(NMax=1000)
REAL Grad,Inter,MSDev,UnGrad,UnInter,Cov,CorCoeff,X(NMax),
* Y(NMax),Y2(NMax),XMean,YMean,XVar,YVar

C Common Block
COMMON/VAL1/NLow,NUpp
COMMON/VAL2/EQ
COMMON/VAL3/N,Grad,Inter
COMMON/VAL4/Err1,NDATA
COMMON/VAL5/MSDev,UnGrad,UnInter,Cov,CorCoeff,X,Y,
* XMean,YMean,XVar,YVar,File1,OutputFile


C Go to beginning of file "DATA" and read off data column "Y2"
REWIND (11)
DO I=1,N
READ(11,*,END=41)Y2(I)
END DO

C Open the Output file
OPEN(UNIT=12, FILE='OutputFile',Status='UNKNOWN')

C Call the Subroutine DTIME() to output the date and time to screen and output file
CALL DTIME()

C Output all the results to the screen
WRITE(*,*)
41 WRITE(*,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(*,*) 'SUMMARY OF RESULTS'
WRITE(*,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(*,*) 'Input Filename: ',File1
WRITE(*,*) 'Total Number of data in File: ', NDATA
WRITE(*,*) 'Index of 1st pair of x-y used: ', NLow
WRITE(*,*) 'Index of Last pair of x-y used:', NUpp
WRITE(*,*) 'Number of data pairs used : ', N
WRITE(*,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(*,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(*,*)
WRITE(*,*) 'S/N INDEX X-DATA Y-DATA Y-LSQR'

NL=NLow

DO I=1,N
WRITE(*,'(I3,3x,I3,2x,3(F7.2,3x))')I,NL,X(I),Y(I),Y2(I)
NL=NL+1
END DO

WRITE(*,*)
WRITE(*,*) '++++++++++++++++++++++++++++++++'
WRITE(*,'(A13,F9.2)') ' Mean X: ',XMean
WRITE(*,'(A13,F9.2)') ' Mean Y: ',YMean
WRITE(*,'(A13,F9.2)') ' Variance X: ',XVar
WRITE(*,'(A13,F9.2)') ' Variance Y: ',YVar
WRITE(*,'(A13,F9.2)') ' Covariance: ',Cov
WRITE(*,'(A13,F9.2)') ' Corr. Coeff: ',CorCoeff
WRITE(*,'(A13,F9.2)') ' Gradient: ',Grad
WRITE(*,'(A13,F9.2)') ' Intercept: ',Inter
WRITE(*,'(A13,F9.2)') ' Mean Sq. Dev:',MSDev
WRITE(*,'(A25,F6.2)') ' Uncertainty in Gradient :',UnGrad
WRITE(*,'(A25,F6.2)') ' Uncertainty in Intercept:',UnInter
WRITE(*,*) '++++++++++++++++++++++++++++++++'
WRITE(*,*) 'Equation of Least Square Line:'
WRITE(*,'(2x,A17)')EQ
WRITE(*,*) '==============================='
WRITE(*,*)


C Save a summary of the results in the file OutputFile

WRITE(12,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(12,*) 'SUMMARY OF RESULTS'
WRITE(12,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(12,*) 'Input Filename: ',File1
WRITE(12,*) 'Total Number of data in File: ', NDATA
WRITE(12,*) 'Index of 1st pair of x-y used: ', NLow
WRITE(12,*) 'Index of Last pair of x-y used:', NUpp
WRITE(12,*) 'Number of data pairs used : ', N
WRITE(12,*) '+++++++++++++++++++++++++++++++++++++++++++'
WRITE(12,*)
WRITE(12,*) 'S/N INDEX X-DATA Y-DATA Y-LSQR'

NL=NLow

DO I=1,N
WRITE(12,'(I3,3x,I3,2x,3(F7.2,3x))')I,NL,X(I),Y(I),Y2(I)
NL=NL+1
END DO

WRITE(12,*)
WRITE(12,*) '===================================='
WRITE(12,'(A13,F9.2)') ' Mean X: ',XMean
WRITE(12,'(A13,F9.2)') ' Mean Y: ',YMean
WRITE(12,'(A13,F9.2)') ' Variance X: ',XVar
WRITE(12,'(A13,F9.2)') ' Variance Y: ',YVar
WRITE(12,'(A13,F9.2)') ' Covariance: ',Cov
WRITE(12,'(A13,F9.2)') ' Corr. Coeff: ',CorCoeff
WRITE(12,'(A13,F9.2)') ' Gradient: ',Grad
WRITE(12,'(A13,F9.2)') ' Intercept: ',Inter
WRITE(12,'(A13,F9.2)') ' Mean Sq. Dev:',MSDev
WRITE(12,'(A25,F6.2)') ' Uncertainty in Gradient :',UnGrad
WRITE(12,'(A25,F6.2)') ' Uncertainty in Intercept:',UnInter
WRITE(12,*) '==============================='
WRITE(12,*) 'Equation of Least Square Line:'
WRITE(12,'(2x,A17)')EQ
WRITE(12,*) '==============================='

WRITE(12,*)
WRITE(*,*) 'NOTE: Summary of Results saved in File "OutputFile"'
WRITE(12,*)
CLOSE(12)

C Return to the main program
RETURN
END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE GNU_PLOT()
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C This Subroutine outputs a plot of the original data and the computed Least
C Squares Best Fit line.
C It also displays the equation of the line, indices of the data used in the
C Analyses (i.e. the part of the data used) and other useful parameters

C Discussion:
C This subroutine reads and plots the data from files "RawData" and
C "BestFit" for the original data and the Least Squares Best Fit
C Line respectively.

C The subroutine write the necessary commands for the graph into the file "CMD"
C and the invokes GNUPLOT to display your picture.

C Once this routine invokes GNUPLOT, a graphics window opens up,
C and the FORTRAN program pauses. Hitting RETURN allows the
C FORTRAN routine to proceed.
C
C The subroutine also gives an option to save a postscript version of the plot,
C re-plot the last graph, or to print a hardcopy of the graph. Control is
C transferred to the appropriate point depending on the user's choice
C

C Modified: 2 January 2006
C Author: Ayeni Gboyega
C Parameters:
C Title,X_Label,Y_Label,N_Lower,N_Upper (Same meanings as described in Main Program)
C EQ (Same meaning as described in subroutine GNUFILE)
C OPN: OPTION (P-PRINT,R-REPEAT,S-SAVE OR X-EXIT)
C STAT: Integer used to check if Gnuplot returns error when called with "SYSTEM"
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

IMPLICIT NONE

C Variable Declaration
CHARACTER X_Label*30,Y_Label*30,Title*30
CHARACTER EQ*17,OPN*1
INTEGER N_Lower,N_Upper,SYSTEM,STAT

C Common block indicating variables common to other subroutines
COMMON/VAL1/N_Lower,N_Upper,Title,X_Label,Y_Label
COMMON/VAL2/EQ

C Open the File "CMD" to write in the gnuplot commands for plotting the graph
21 OPEN(UNIT=10,FILE='CMD',STATUS='UNKNOWN')

C Write the gnuplot commands to set the title, labels for x- & y-axes, equation
C of the least squares fit line and indices of the data used in the analyses
WRITE(10,*)'set title ','"',Title,'"'
WRITE(10,*)'set xlabel ','"',X_Label,'"'
WRITE(10,*)'set ylabel ','"',Y_Label,'"'
WRITE(10,*)'set label 1 ','"EQN OF BESTFIT LINE: ',EQ,'"',' at
* graph .01,.97 left'
WRITE(10,*)'set label 2 ','"NLOWER:',N_Lower,'"',' at
* graph .01,.93 left'
WRITE(10,*)'set label 3 ','"NUPPER:',N_Upper,'"',' at
* graph .01,.89 left'

C Write the gnuplot command to plot the Raw Data and Best-Fit line
WRITE ( 10,'(A)')'plot "RawData" using 1:2,"BestFit"
* using 1:2 with lines'

C Write the gnuplot command to pause while the graph is displayed
WRITE ( 10, '(A)' ) 'pause -1 "Hit return to continue"'

C Write the gnuplot command to quit after return is hit by the user
WRITE ( 10, '(A)' ) 'q'

C Call program GNUPLOT to display plot using the commands in file "CMD"
STAT = SYSTEM ('gnuplot CMD')

C If an Error was generated while the gnuplot command is issued,
C the output an Error message and exit. If No error occurs here
C then there would not be need to check in subsequent sections"
IF(STAT /= 0 ) THEN
WRITE(*,*)'********************NOTE*************************'
WRITE(*,*)'FATAL ERROR!'
WRITE(*,*)' AN ERROR CODE WAS GENERATED WHEN THE GNUPLOT '
WRITE(*,*)' COMMAND WAS ISSUED. MAYBE "GNUPLOT" IS NOT INSTALLED!! '
WRITE(*,*)' CHECK TO ENSURE THAT THE PROGRAM IS IN YOUR PATH.'
WRITE(*,*)'*************************************************'
WRITE(*,*)' '
STOP
END IF

CLOSE (10)

C Ask User for what to do next: (P)rint, (S)ave, (R)epeat Plot or E(x)it
30 WRITE(*,*)
WRITE(*,*)'OPTIONS:'
WRITE(*,*)'---------'
WRITE(*,*)'R: Repeat Plot'
WRITE(*,*)'S: Save Plot'
WRITE(*,*)'P: Print Plot'
WRITE(*,*)'X: Exit'

C Read User's response. If Error occurs while reading, show options again.
READ(*,'(A1)',ERR=30)OPN

C If User would like to see plot again, transfer control to statement label 21.
IF (OPN.EQ.'R'.OR.OPN.EQ.'r') GOTO 21


C If User would like to save the plot then follow the same procedure
C used to display the plot but with the terminal set to postscript and
C output set to a the file "LSQFile.ps"
IF (OPN.EQ.'S'.OR.OPN.EQ.'s') then
OPEN(UNIT=10,FILE='CMD',STATUS='OLD')
WRITE(10,*)'set output "LSQFile.ps"'
WRITE(10,*)'set term postscript'
WRITE(10,*)'set title ','"',Title,'"'
WRITE(10,*)'set xlabel ','"',X_Label,'"'
WRITE(10,*)'set ylabel ','"',Y_Label,'"'
WRITE(10,*)'set label 1 ','"EQN OF BESTFIT LINE: ',EQ,'"',' at
* graph .01,.97 left'
WRITE(10,*)'set label 2 ','"NLOWER:',N_Lower,'"',' at
* graph .01,.93 left'
WRITE(10,*)'set label 3 ','"NUPPER:',N_Upper,'"',' at
* graph .01,.89 left'

WRITE ( 10,'(A)')'plot "RawData" using 1:2,"BestFit"
* using 1:2 with lines'


C Write Message for Fortran to pause and inform user that the file was saved
WRITE ( 10, '(A)' ) 'pause -1 "File Saved"'
WRITE ( 10, '(A)' ) 'q'

C Call GNUPLOT to output the plot to file "LSQFile.ps"
STAT = SYSTEM ('gnuplot CMD')
CLOSE (10)

C Inform User that the file was saved
WRITE(*,*)'The graph has been saved as "LSQFile.ps" in your
* present working directory'

WRITE(*,*) 'Summary of Results saved in File "OutputFile"'

WRITE(*,*)

C Return to Options Statement
GOTO 30

ENDIF


C If User would like to print the plot then follow the same procedure
C used to save the plot and inform user to open the file
C LSQFile.ps and print
IF (OPN.EQ.'P'.OR.OPN.EQ.'p') then


C Open the file "CMD" to write the gnuplot commands
OPEN(UNIT=10,FILE='CMD',STATUS='OLD')

WRITE(10,*)'set output "LSQFile.ps"'
WRITE(10,*)'set term postscript'
WRITE(10,*)'set title ','"',Title,'"'
WRITE(10,*)'set xlabel ','"',X_Label,'"'
WRITE(10,*)'set ylabel ','"',Y_Label,'"'
WRITE(10,*)'set label 1 ','"EQN OF BESTFIT LINE: ',EQ,'"',' at
* graph .01,.97 left'
WRITE(10,*)'set label 2 ','"NLOWER:',N_Lower,'"',' at
* graph .01,.93 left'
WRITE(10,*)'set label 3 ','"NUPPER:',N_Upper,'"',' at
* graph .01,.89 left'
WRITE ( 10,'(A)')'plot "RawData" using 1:2,"BestFit"
* using 1:2 with lines'
WRITE(10,'(A)')'!lpr -Pps49 LSQFile.ps'
WRITE(10,'(A)')'pause -1 "To Print, open the file LSQFile.ps"'
WRITE(10,'(A)')'q'
STAT = SYSTEM ('gnuplot CMD')
CLOSE (10)
WRITE(*,*)
WRITE(*,*)'To Print, open the file LSQFile.ps and print'
WRITE(*,*)

C Return to Options Statement
GOTO 30

ENDIF

C If User wishes to exit, the return to main program
IF (OPN.EQ.'X'.OR.OPN.EQ.'x') GOTO 24

C If User's Option is invalid, then Inform User and return to Options
IF(OPN.NE.'P'.AND.OPN.NE.'p') then
IF(OPN.NE.'S'.AND.OPN.NE.'s') then
IF(OPN.NE.'R'.AND.OPN.NE.'r') then
IF(OPN.NE.'X'.AND.OPN.NE.'x') then
WRITE(*,*)
WRITE(*,*) 'INVALID ENTRY'
WRITE(*,*) 'SELECT P,S,R OR X'
WRITE(*,*)
GOTO 30
END IF
END IF
END IF
END IF

C Return to the main program
24 RETURN
END


C*******************************************************************************
C Subroutine DTIME() outputs the current time in the Screen and Output File
SUBROUTINE DTIME()

C Implicit Statement
IMPLICIT NONE

C Variable declaration
CHARACTER DATE*8,TIME*10,ZONE*5
INTEGER D,H,M,MM,N,S,Y,VALUES(8)
CHARACTER(4), PARAMETER, DIMENSION(12) :: MONTH = (/
* 'JAN. ','FEB. ','MAR. ','APR. ','MAY ','JUN. ',
* 'JUL. ','AUG. ','SEP. ','OCT. ','NOV. ','DEC. ' /)

C Call function DATE_AND_TIME() to get the date, time ,etc
CALL DATE_AND_TIME (DATE,TIME,ZONE,VALUES)
Y=VALUES(1)
M=VALUES(2)
D=VALUES(3)
H=VALUES(5)
N=VALUES(6)
S=VALUES(7)
C Output the date and time to the screen and file
WRITE(*,*)
WRITE(*,'(1X,A,I2,1X,I4,I2,A1,I2,A1,I2)' )
* MONTH(M), D, Y, H,':',N,':',S
WRITE(12,*)
WRITE(12,'(1X,A,I2,1X,I4,I2,A1,I2,A1,I2)' )
* MONTH(M), D, Y, H,':',N,':',S

C Return to the main program
RETURN
END