!Testing the Monte Carlo routines...

program MC_quasi_random
use MY_STUFF
Implicit none

     integer, parameter :: d=3, N=100000
     integer :: i
     real(KIND=dbl) :: V, f_av1, f_av2, val, val1, val2, err, N_real
     real(KIND=dbl), dimension(d) :: a, b, alfa, beta

alfa(1) = SQRT(3.0) ; alfa(2) = SQRT(5.0) ; alfa(3) = SQRT(6.0)
beta(1) = SQRT(2.0) ; beta(2) = SQRT(10.0) ; beta(3) = SQRT(12.0)

a(1) = 0.0 ; b(1) = pi
a(2) = 0.0 ; b(2) = pi
a(3) = 0.0 ; b(3) = pi

!Integration volume:
V = 1
do i = 1, d
   V = V*(b(i) - a(i))
enddo

call MC_quasi(a, b, N, d, alfa, f_av1)
call MC_quasi(a, b, N, d, beta, f_av2)

val1 = V * f_av1
val2 = V * f_av2

N_real = REAL(N,dbl)

val = (val1 + val2)/2.0

!Error:
err = ABS(val1 - val2)

print '("Number of distributed points:", i8)', N
print '(/)'
print '("Value of Monte Carlo integration (quasi-random sampling):", f10.6)', val
print '(/)'
print '("Estimated error:", f10.4)', err

end program MC_quasi_random