#program to compare grid and Monte Carlo point sampling to est. area
#T.L. Einstein, Dec. 2003
from math import *
from random import *
print "Volume of unit-sphere in N-space is pi**(N/2)/(N/2)! \n"
print "(4/3)pi/8 = pi/6 is %f \n" %(pi/6.) 
print "Integral of r^2 over octant of unit sphere is (1/8)(4pi/5),"\
" so  %f \n" %(0.1*pi)
for n in range(1,9):
    xyzstart=1./2**n #starting position in each direction
    xyzincr=2*xyzstart #increment in each direction
    iestpio4=0  #running total of grid points inside curve
    rsqest=0.
    iranpio4=0  #running total of randomly-poicked points inside curve
    rransqMC=0.
    npts=0
    for ix in range(2**(n-1)):
        xval=xyzstart+ix*xyzincr
        for iy in range(2**(n-1)):
            yval=xyzstart+iy*xyzincr
            for iz in range(2**(n-1)):
                zval=xyzstart+iz*xyzincr
                rvalsq=xval**2 + yval**2 + zval**2
                xran=random()
                yran=random()
                zran=random()
                rransq=xran**2 + yran**2 + zran**2
                if rvalsq<=1: 
                    iestpio4 += 1
                    #test of whether the point is inside the sphere
                    rsqest += rvalsq
                    #sum of integrand r**2 at randomly-selected points
                if rransq <= 1: 
                    iranpio4 += 1
                    rransqMC += rransq
    print "At level %d, estimate of pi/4 by grid is %f & by MC is %f"\
    %(n,float(iestpio4)/(8**(n-1)),float(iranpio4)/(8**(n-1)))
    #2**(n-1) points in each direction
    print "   At level %d, estimate of integral r^2 by grid is %f & by MC is %f"\
    %(n,rsqest/(8**(n-1)),rransqMC/(8**(n-1)))