## Monte Carlo simulation to estimate the volume of a sphere.
## Points are generated in a box of side 2, centered at the
## origin. The box has volume 8. The sphere is centered at the
## origin. It has radius 1. If the random point lies within the
## sphere, count as "in". Ratio of ins to total number of points
## gives the ratio of the volume of the sphere to the box. So
## the volume of the sphere is the in/total ratio times 8.

from random import *
from visual import *    ## we are not drawing, but will use mag()

## define a function that returns a random number between -1 and 1
def MCcoord():
    x = random()      ## make a random number between 0 and 1
    x *= 2.0          ## double puts into range 0 to 2
    x -= 1.0          ## subtract 1 to put in range -1 to 1
    return x

print "This program uses a Monte Carlo method to estimate the"
print "volume of an object (initially, a sphere, but will change.)"

num_points = input("How many points do you want to throw? ")

num_in = 0.0
for i in range(num_points):
    MCpoint = vector(MCcoord(),MCcoord(),MCcoord())
    ## here is the MC test - the decision/score
    pointsize=0.04
    if mag(MCpoint) < 1:
        s = sphere(radius=pointsize,pos=MCpoint,color=color.red)
        num_in += 1.0

print "Number within the volume is ", num_in
print "Estimated volume is", 8.0 * num_in / num_points
print "Volume of sphere of radius 1 is", 4.0 * 3.1415926535 / 3.0