""" From "A SURVEY OF COMPUTATIONAL PHYSICS", Python eBook Version
   by RH Landau, MJ Paez, and CC Bordeianu
   Copyright Princeton University Press, Princeton, 2012; Book  Copyright R Landau, 
   Oregon State Unv, MJ Paez, Univ Antioquia, C Bordeianu, Univ Bucharest, 2012.
   Support by National Science Foundation , Oregon State Univ, Microsoft Corp"""                             

#  StoneThrowing: Monte-Carlo integration via stone throwing 

import random
from visual.graph import *

N        = 1000   # points to plot the function
graph    = display(width=600,height=600,title='Integration by Stone Throwing')
xsinx    = curve(x=list(range(0,N)), color=color.yellow, radius=0.5)    
pts      = label(pos=(-30, -50), text='x ->', box=0) 
pts      = label(pos=(-60, -60), text='Number Stones = ', box=0)          # labels
pts2     = label(pos=(-30, -60), box=0)
inside   = label(pos=(30,-60), text='Number Splashes = ', box=0)
inside2  = label(pos=(60,-60), box=0)
arealbl  = label(pos=(-65,60), text='Computed A = ', box=0)
arealbl2 = label(pos=(-35,60), box=0)
areanal  = label(pos=(30,60), text='Analytic A = ', box=0)
zero     = label(pos=(-85,-48), text='0', box=0)
fofx     = label(pos=(-85,0), text='f(x)', box=0)
five     = label(pos=(-85,50), text='5', box=0)
twopi    = label(pos=(90,-48), text='2pi',box=0)

def fx (x):  return x*sin(x)*sin(x)      # Function to integrate
		
def plotfunc():             # Plot function & box
    incr = 2.0*pi/N
    for i in range(0,N):
        xx         = i*incr
        xsinx.x[i] = ((80.0/pi)*xx-80)
        xsinx.y[i] = 20*fx(xx)-50
    box            = curve(pos=[(-80,-50), (-80,50), (80,50)
                   ,(80,-50), (-80,-50)], color=color.white)        # box
    
plotfunc()                         # area of box = height * width = 5*2pi
j            = 0
Npts         = 100                    # points generated inside the box
analyt       = pi**2                              # Analytical integral
areanal.text = 'Analytic A = %8.5f'%analyt     # Output analytical integral
genpts       = points(size=2)
for i in range(1,Npts):                     # generates points inside box
    rate(500)                 # slow process down to see point generation
    x = 2.*pi*random.random()                               # 0=< x <= 2pi
    y = 5.*random.random()                                    # 0=< x <= 5
    xp = x*80./pi-80.
    yp = 20.0*y-50. 
    pts2.text = '%4s' %i
    if y <= fx(x):                                               # splash
           j += 1                               # increase numbr splashes
           genpts.append(pos=(xp,yp), color=color.cyan)
           inside2.text='%4s'%j                      # Label for splashes
    else:  genpts.append(pos=(xp,yp), color=color.green)
    boxarea = 2.*pi*5.                          # Area of box propto Npts
    area = boxarea*j/(Npts-1)                    # Area of curve propto j
    arealbl2.text = '%8.5f'%area            # write current computed area