# -*- coding: utf-8 -*-
from __future__ import division # Trata como real o resultado de uma divisão entre inteiros
from visual import *            # Importa o módulo de visualização (VPython)
from visual.graph import *      # Importa a parte do VPython que traça gráficos 2D
from random import random       # Importa o gerador de números pseudoaleatórios

# Definição da função a ser integrada
def f(x):
	return sqrt(1-x*x)
	
# Número de pontos a sortear para a estimativa da integral
N = 100

# Traçando o gráfico da função f(x)
a = 0        # Valor mínimo de x
b = 1        # Valor máximo de x
dx = 0.001   # Intervalo de x para os pontos consecutivos
imax = (b-a)/dx+1 # Quantos pontos iremos traçar
#    Criamos uma janela gráfica 2D e indicamos que traçaremos o gráfico de f(x) em amarelo
fig1 = gdisplay(width=500, height=500, xtitle='x', ytitle='f(x)', ymin=0, ymax=1, xmin=a, xmax=b)
grafico = gcurve(color=color.yellow)
#    Varremos os pontos a traçar
i=1
while i<=imax:
    x = a+dx*(i-1)
    grafico.plot(pos=(x,f(x)))
    i = i+1

# Inicializando as variáveis que contam os pontos sorteados e aqueles abaixo da curva da função
n = 0
conta = 0
pontos_acima = gdots(color=color.red, size=4)     # Vamos traçar em vermelho os pontos acima da curva
pontos_abaixo = gdots(color=color.green, size=4)  # Vamos traçar em azul os pontos acima da curva
fig2 = gdisplay(xtitle='n', ytitle='|erro|', ymin=1e-4, ymax=10, logy=True) # Cria uma outra janela para traçar o erro
graf_erro = gcurve(gdisplay=fig2, color=color.cyan) # Vamos traçar o módulo do erro relativo em azul claro
pi_exato = 4*atan(1.)  # Valor exato de pi para comparação
while n<N:         # Vamos sortear ao todo N pontos 
	rate(1000)     # Deixamos a exibição mais lenta
	x = random()   # Sorteamos a coordenada x do ponto
	y = random()   # Sorteamos a coordenada y do ponto
	if y < f(x):   # O ponto sorteado está abaixo da curva da função?
		conta = conta+1		           # Atualizamos o contador desses pontos
		pontos_abaixo.plot(pos=(x,y))  # Registramos o ponto no gráfico
	else:          # O ponto sorteado está acima da curva da função?
		pontos_acima.plot(pos=(x,y))
	n = n+1            # Atualizamos o número de pontos sorteados
	frac = conta/n     # Calculamos a fração abaixo da curva
	pi_est = 4*frac    # Estimamos pi
	erro_rel = (pi_est-pi_exato)/pi_exato # Calculamos o erro relativo
	graf_erro.plot(pos=(n,abs(erro_rel))) # Traçamos a evolução do módulo do erro relativo
	print ('Pontos = ',n,' Estimativa = ',pi_est,' Erro relativo = ',erro_rel)