program main
	use montecarlo
	implicit none
	real(8), external :: f
	real(8), parameter :: pi = 3.141592635
	real(8), parameter :: trueint = 1.3932039296856768591842462603255
	integer, parameter :: N = 1000000
	real(8), parameter :: acc = 0.001
	real(8) :: integral, err
	real(8) :: a(3), b(3)
	
	! Integration limit:
	a(1:3) = 0.0
	b(1:3) = pi
	
	! Run Plain Monte Carlo:
	print *, " Finding the integral of f(x,y,z) = 1/(1-cos(x)cos(y)cos(z))*(1/pi^3)"
	print *, " N = ", N
	print *, " True value: ", trueint
	print *, "=================================="
	print *, " Plain Monte Carlo: "
	
	call plainmc(f, a, b, N, integral, err)
	
	print *, "Integral = ", integral
	print *, "Error = ", err
	print *, "Integral fejl = ", abs(integral - trueint)


	print *, "=================================="
	print *, " Quasi-Random: "

	! Run Quasi-Random Monte Carlo:
	call quasirandom(f, a, b, N, integral, err)
	
	print *, "Integral = ", integral
	print *, "Error = ", err
	print *, "Integral fejl = ", abs(integral - trueint)

	print *, "=================================="
	print *, " Stratified: "

	! Run Stratified Monte Carlo:
	call strata(f, a, b, N, acc, integral, err)
	
	print *, "Integral = ", integral
	print *, "Error = ", err
	print *, "Integral fejl = ", abs(integral - trueint)
	print *, "=================================="
end program main

real(8) function f(x)
	implicit none
	real(8), intent(in) :: x(3)
	real(8), parameter :: pi = 3.141592635

	f = (1.0/pi**3) * 1.0/(1-cos(x(1))*cos(x(2))*cos(x(3)))
end function