# Numerical integration to calculate value of pi
from pylab import *

r = 1.0					# radius of circle
dx = .05				# integration delta x
x = arange(0, r+dx, dx)	# array ox x-coordinates
y = sqrt(r**2 - x**2)   # equation of circle, r^2 = x^2 + y^2
plot(x,y, 'x-')
show()

N = len(y)
A = 0.0
for k in range(N-1):			# N-1 trapezoids 
	at = (y[k] + y[k+1])* dx/2	# area of one trapezoid 
	A += at						# add up the areas

print A*4, pi		# compare area of circle calculated and pi*r^2