PROGRAM prog5A
!Initialise the program and give it a name
	IMPLICIT NONE
	!Assume nothing about the variables used within this program
	INTEGER :: npts = 100
	!Declare and assign an integer variable for the number of points to calculate for the plot
	INTEGER, EXTERNAL :: pgbeg
	!Declare an external variable for the pgplot library
	REAL, DIMENSION(100) :: xarray, yarray
	!Declare two 100 element vectors to store the x and f(x) values
	REAL :: a, b, c, h, xmin, xmax, ymin, ymax
	!Declare real variables for the coefficients of the quadratic, the subinterval width and the min/max
	!x and y values
	INTEGER :: i, ierr
	!Declare an integer counting variable and an integer variable for the error return of pgbeg()
	
	PRINT*, 'Input your coefficients, a, b, c, of your quadratic:'
	!Prompt the user to input the coefficients of their quadratic
	READ*, a, b, c
	!Read the coefficients into their respective a, b and c variables
	
	IF (ABS(a) < TINY(a)) THEN
	!Check if a = 0 as this implies the function is not quadratic and so does not need to be dealt with
	!by this program
		PRINT*, 'The function is not quadratic.'
		!Print a comment as such
		GOTO 1
		!Go to the end of the program as we can exit early
	END IF
	!End the if statement
	
	!My xmin and xmax are the points in x at which f'(x) = +/- 1 This gives a nice looking graph
	xmin = (-1 - b)/(2*a)
	!Assign the usual minimum value for x for a > 0
	IF (xmin < 0) THEN
	!Check this value is less than zero
		xmax = (1 - b)/(2*a)
		!If so, calculate the appropriate maximum x value
	ELSE
		xmax = xmin
		xmin = (1 - b)/(2*a)
		!Otherwise calculate min/max x values such that xmax > xmin
	END IF	
	
	IF (a > 0) THEN
		ymin = f(a, b, c, (-b/(2*a)))
		ymax = f(a, b, c, xmin)
		ymin = ymin - 0.25*(ymax - ymin)
		!For a > 0 let ymax be f(xmin) = f(xmax) and ymin be 25% of ymax - f(x) below f(x) with x
		!such that f'(x) = 0
	ELSE
		ymin = f(a, b, c, xmin)
		ymax = f(a, b, c, (-b/(2*a)))
		ymax = ymax + 0.25*(ymax - ymin)
		!For a < 0 let ymin be f(xmin) = f(xmax) and ymax be 25% of ymax - f(x) above f(x) with x
		!such that f'(x) = 0
	END IF
	!End the if statement
	
	h = (xmax - xmin)/(npts - 1)
	!Calculate the size of the subintervals for the array generation below
	
	DO i = 1, npts
	!Start a definite DO loop
		xarray(i) = xmin + (i - 1)*h
		!Calculate the x value in question running from xmin to xmax in steps of h
		yarray(i) = f(a, b, c, xarray(i))
		!Calculate the corresponding y values
	END DO
	!End the DO loop

	ierr = pgbeg(0, 'quadratic_plot.ps/VPS', 1, 1)
	!Call the function pgbeg within the pgplot library to check that the postscript type output file can
	!be created correctly

	IF (ierr == 1) THEN
	!If the output can be correctly created then we can plot the graph
		CALL pgenv(xmin, xmax, ymin, ymax, 0, 1)
		!Call the pgplot subroutine to set the min/max x/y values and how we want the axes and graph to
		!be presented
		CALL pglab('x', 'y = f(x)', 'Graph of a quadratic')
		!Call the pgplot subroutine to set the labels of the graph and axes
		CALL pgline(100, xarray, yarray)
		!Call the pgplot subroutine to plot the graph
		CALL pgend
		!Call the pgplot subroutine to end the pgplot program
	ELSE
		PRINT*, 'Error: Plotting did not begin correctly'
		!If the output file cannot be created then print that an error has occurred
	END IF
	!End the if statement

CONTAINS
!State that the following functions/subroutines will be used by the above program block but we wish to contain
!them within the code

	REAL FUNCTION f(av, bv, cv, xv) RESULT(yv)
	!Declare a function with real result that reads the a, b, c coefficients and the x for which f(x)
	!is to be calculated
		REAL, INTENT(IN) :: av, bv, cv, xv
		!Declare dummy variables for a, b, c and x
		
		yv = (av*(xv**2)) + (bv*xv) + cv
		!Calculate f(x)
	END FUNCTION f
	!End the function, f
1 END PROGRAM
!End the program with the label 1 for the early exit case where a = 0
