#!/usr/bin/env python

"""Computes an estimate of pi using Monte Carlo trials

USAGE: (from shell prompt)
    python ./monte_carlo_pi.py

AUTHOR:
    Jonathan Senning <jonathan.senning@gordon.edu>
    Gordon College
    Python Version: March 8, 2008
    Adapted to shell program: October 27, 2010
"""

from pylab import *

def monte_carlo_pi( n ):
    """Estimate pi by Monte Carlo trials

    USAGE:
        p = monte_carlo_pi( n )

    INPUT:
        n       - positive integer; number of trials ("throws at unit square")

    OUTPUT:
        float   - estimate of pi

    NOTES:

        This is based on a description in an introduction to Monte-Carlo
        methods found at http://www.chem.unl.edu/zeng/joy/mclab/mcintro.html.
        The idea is to estimate pi by noting that the area of a quarter of a
        unit circle is pi/4 and that if n darts are thrown randomly into a
        unit square the fraction of them that land in a quarter of a unit
        circle is equal to the fraction of the square area taken up by the
        quarter of the unit circle.

            |******--------------
            |       ****  "miss"|
            |            **     |
            |               *   |
            |                *  |
            |  "hit"          * |
            |                 * |
            |                  *|
            |__________________*|

    AUTHOR:
        Jonathan R. Senning <jonathan.senning@gordon.edu>
        Gordon College
        Written February 15, 2005
        Converted to Python May 21, 2008
    """
    x = rand( n, 1 )
    y = rand( n, 1 )

    # Here is an explanation of the calculation which follows:
    #  a = sqrt( x**2 + y**2 )         Computes the array of distance values
    #  b = 1 - sign( a )	           2 if distance is less than 1,
    #                                         0 if distance is greater than 1
    #  c = sum( 1 - sign( b ) ) / 2    Number of "hits"
    #  4 * ( c ) / n                   Estimate of pi.

    return 4 * ( sum( 1 - sign( sqrt( x**2 + y**2 ) - 1 ) ) / 2 ) / n

if __name__ == "__main__":

    n = 10000000
    p = monte_carlo_pi( n )
    print "Using %d trials, the Monte Carlo estimate of pi is %f" % (n, p)