MODULE integrate_library

CONTAINS
	FUNCTION integrate(integrand, integrandParameters, xRange, xStep, yRange, yStep, zRange, zStep)
	!Takes in a function "integrand" that returns a single real value
	!Output will have (xRange(2)-xRange(1))/xStep + 1 elements
		REAL(8), EXTERNAL :: integrand
		REAL, DIMENSION(:) :: integrandParameters
		REAL, DIMENSION(2), INTENT(IN) :: xRange, yRange, zRange
		REAL(8), DIMENSION(:, :, :), ALLOCATABLE :: integral
		REAL(8), DIMENSION(:), ALLOCATABLE :: integrate
		INTEGER :: xMaxIndex, yMaxIndex, zMaxIndex, x, y, z, OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS
		REAL :: xVal, yVal, zVal, xStep, yStep, zStep

		xMaxIndex = getMaxIndex(xRange, xStep)
		yMaxIndex = getMaxIndex(yRange, yStep)
		zMaxIndex = getMaxIndex(zRange, zStep)

		ALLOCATE(integral(xMaxIndex, yMaxIndex, zMaxIndex))
		ALLOCATE(integrate(xMaxIndex))

		!$OMP PARALLEL DO PRIVATE(x,y,z,xVal,yVal,zVal)
		DO x=0,xMaxIndex-1			
			xVal = x * xStep + xRange(1)		
		
			DO y=0,yMaxIndex-1
				yVal = y * yStep + yRange(1)

				DO z=0,zMaxIndex-1
					zVal = z * zStep + zRange(1)		

					integral(x+1,y+1,z+1) = integrand(integrandParameters, xVal, yVal, zVal)
				ENDDO
			ENDDO
		ENDDO
		!$OMP END PARALLEL DO

		integral(:,:,1) = SUM(integral, DIM=3)

		integrate = SUM(integral(:,:,1), DIM=2)

		DEALLOCATE(integral)

		integrate = integrate * yStep * zStep !/ integrate(3) !MAXVAL(integrate) !

		RETURN
	END FUNCTION integrate

	FUNCTION getFrequencyList(varRange, varStep)
		REAL, DIMENSION(2) :: varRange
		REAL :: varStep
		REAL, DIMENSION(:), ALLOCATABLE :: getFrequencyList
		INTEGER :: maxVarIndex

		maxVarIndex = getMaxIndex(varRange, varStep)

		ALLOCATE(getFrequencyList(maxVarIndex))

		DO var=0,maxVarIndex-1 
			getFrequencyList(var+1) = var * varStep + varRange(1)
		END DO

		RETURN
	END FUNCTION getFrequencyList

	FUNCTION getMaxIndex(varRange, varStep)
		REAL, DIMENSION(2) :: varRange
		REAL :: varStep
		INTEGER :: getMaxIndex

		getMaxIndex = NINT(ABS(varRange(2)-varRange(1))/varStep) + 1
	
		RETURN
	END FUNCTION getMaxIndex

	SUBROUTINE writeData(filename, leftColumn, rightColumn, rowNum)
		INTEGER, INTENT(IN) :: rowNum
		REAL, DIMENSION(:), INTENT(IN) :: leftColumn
		REAL(8), DIMENSION(:), INTENT(IN) :: rightColumn
		CHARACTER(LEN=*), INTENT(IN) :: filename

		OPEN(UNIT=1, FILE=filename, ACCESS='SEQUENTIAL', ACTION='WRITE')

		do row=1, rowNum
			WRITE(1, *), leftColumn(row), rightColumn(row)
		enddo

		CLOSE(1)
	END SUBROUTINE writeData

	SUBROUTINE getArgumentsLine(outputFile, alpha, AR, elongation, nu, theta, v)	
	!Provide command line arguments
	!Numbers must be given with a decimal part.
	IMPLICIT NONE
		CHARACTER(LEN=32) :: outputFile, charAlpha, charAR, charElongation, charNu, charTheta, charV
		REAL :: alpha, AR, elongation, nu, theta, v

		CALL GET_COMMAND_ARGUMENT(1, outputFile)
		CALL GET_COMMAND_ARGUMENT(2, charAlpha)
		CALL GET_COMMAND_ARGUMENT(3, charAR)
		CALL GET_COMMAND_ARGUMENT(4, charElongation)
		CALL GET_COMMAND_ARGUMENT(5, charNu)
		CALL GET_COMMAND_ARGUMENT(6, charTheta)
		CALL GET_COMMAND_ARGUMENT(7, charV)

		READ(charAlpha, '(F)') alpha
		READ(charAR, '(F)') AR
		READ(charElongation, '(F)') elongation 
		READ(charNu, '(F)') nu
		READ(charTheta, '(F)') theta
		READ(charV, '(F)') v
	END SUBROUTINE getArgumentsLine

END MODULE integrate_library