#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Oct 11 17:42:58 2023

@author: victorcolas
"""

import numpy as np
import random as rd
import matplotlib.pyplot as plt


def calculPI(N):
    R=1
    DansCercle=0
    for k in range(N):
        ux=2*R*rd.random()-R
        uy=2*R*rd.random()-R
        r=np.sqrt(ux**2+uy**2)
        if r**2<R:
            DansCercle=DansCercle+1

    EstimPi=4*DansCercle/N
    return EstimPi


## Simple calcul
N=1000
EstimPi=calculPI(N)
print(EstimPi)
    
## Convergence vers pi fonction de N
nPoint=50
listeEstimPi=[]
listeNphot=[]
listeN=np.linspace(1,5,nPoint)
for k in range(nPoint):
    print(k)
    Nphot=int(10**listeN[k])
    EstimPi=calculPI(Nphot)
    listeEstimPi.append(EstimPi)
    listeNphot.append(Nphot)

plt.figure()
plt.xlabel('Nombre de photons N')
plt.ylabel('Estimation de pi')
plt.grid(1)
plt.semilogx(listeNphot,listeEstimPi,label='Estimation(N)')
plt.semilogx([listeNphot[0],listeNphot[-1]],[np.pi, np.pi],color='k',label='Valeur de pi',ls='--')
plt.legend()