#!/usr/bin/env python           
# -*- coding: iso-8859-1 -*-    

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab

#-----------------------------------------------------------------
def buffon(N, needle_length):
  """Returns an estimate of pi."""
  counter_hit    = 0.

  for i in range(N):
    needle_pos = np.random.uniform(0,1)             #x-position of the center of the needle
    needle_angle = np.random.uniform(0, np.pi/2.)   #angle the needle takes with x-axis
                                                    #We always assume the end of the needle with
                                                    # the smaller angle.
    dist_x = min([needle_pos, 1. -needle_pos])      #The needle may touch more than one line,
                                                    # but the line with the smaller distnance to the
                                                    # center should be the first one to be touched.

    if np.cos(needle_angle) * needle_length/2. > dist_x:
                                                    #Only half the needle's length reach from the center.
      counter_hit += 1.

  return (N/counter_hit) * 2. * needle_length

#-----------------------------------------------------------------
def alternative(N):
  """Make use of alternative method"""
  counter_hit    = 0.

  for i in range(N):
    #A random dot has two independent components x and y.
    #If the x^2 + y^2 <= 1, the dot is in the circle.
    if np.random.uniform(0,1)**2 + np.random.uniform(0,1)**2 < 1:
      counter_hit += 1.

  return (counter_hit/N) * 4.
#----------------------------------------------------------------------------------

N              = 1000000   #Number of throws
needle_length  = 1.       #Length of the needle, needs to be <= 1.

print """
      Proof of Principle.
      Pi is: 
      """
print np.pi
print """
      Buffon with 1 million iterations and a needle length of 1 is:
      """
print buffon(N, needle_length)
print
print """
      Alternative Method with 1 million iterations:
      """
print alternative(N)
print


length_steps  = 5      #step size is 0.2; it starts from 1 and then goes down.
orders_of_mag = 60     #actually an order of magnitude is now only 1.2

collection= []

for ii in range(length_steps):
  list_fixed_length_v = []
  list_fixed_length_N = []
  for jj in range(orders_of_mag):
    list_fixed_length_v.append(abs(buffon(int(10*1.2**(jj+1)), (1-ii/5.))-np.pi))
    list_fixed_length_N.append(int(10*1.2**(jj+1)))

  collection.append([list_fixed_length_N, list_fixed_length_v])

print collection

color_string = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']

for ii in range(length_steps):
  plt.plot(collection[ii][0], collection[ii][1], color = color_string[ii], label=("Needle length ") + str((1.-ii/5.)), linewidth=2.5)
  plt.yscale('log')
  plt.xscale('log')


list_v = []
list_N = []
for jj in range(orders_of_mag):
  list_v.append(abs(alternative(int(10*1.2**(jj+1))) - np.pi))
  list_N.append(int(10*1.2**(jj+1)))

print list_v
print list_N

plt.plot(list_N, list_v, color = color_string[-2], label="Alternative", linewidth=2.5)
plt.yscale('log')
plt.xscale('log')

import pylab
pylab.xlabel("Number of Throws")
pylab.ylabel("Deviation from Pi")
plt.legend(loc='upper right')
pylab.savefig('Deviation_From_Pi')
plt.show()