## MÉTHODE MONTE-CARLO - Relation Z = f(X,Y)
## Armel MARTIN - Septembre 2021
## distribué sous licence CC-by-nc-sa
import numpy as np
import matplotlib.pyplot as plt
import random

#### MESURE D'UNE CAPACITÉ VIA UN TEMPS CARACTÉRISTIQUE
## Mesures (X=R ; Y=tau)
xm,ym = 1032, 0.0012 # Ohm, s
# Précision de la mesure
DeltaX,DeltaY = 50, 0.0001 # Ohm, s
# Incertitude (distribution uniforme supposée)
uX,uY = DeltaX / np.sqrt(3) , DeltaY / np.sqrt(3)

def f(x,y):
    return y/x # Fonction quotient

## Génération aléatoire d'une série de mesures
N = 100000 # nombre de simulations à exécuter.

## En utilisant directement des vecteurs Numpy
vecX,vecY = np.zeros(N),np.zeros(N)
for i in range(0,N):
    vecX[i] = random.uniform(xm-DeltaX,xm+DeltaX) # spécifier l'intervalle
#    vecX[i] = random.gauss(xm,uX) # spécifier moyenne et écart-type
    vecY[i] = random.uniform(ym-DeltaY,ym+DeltaY) # spécifier l'intervalle
#    vecY[i] = random.gauss(ym,uY) # spécifier moyenne et écart-type
vecZ = f(vecX,vecY) # calcul vectorisé

## Statistique sur Z
Zm = np.mean(vecZ)
uZ = np.std(vecZ,ddof=1) # ddof=1 => non biaisé
print('Zm =',Zm,' ; u(Z) =',uZ)

plt.figure(1,figsize=[5,3])
plt.hist(vecX,bins='rice') #,color='grey')
plt.hist(1e6*vecY,bins='rice') #,color='grey')
plt.xlim(0,1.1*max(10**6*vecY))
plt.xlabel('X et 10^6*Y')# (10^-6 s)')
plt.ylabel('frequence')
plt.grid()
titre = "Distribution de X,Y - Incert. relat. ="+ \
        format(100*uX/xm,'.1f')+'%' + ' et '+ format(100*uY/ym,'.1f')+'%' 
plt.title(titre)
plt.tight_layout()
plt.savefig('hist_X,Y_1.png')
plt.show()
#plt.clf()

plt.figure(2,figsize=[5,3])
plt.hist(1e6*vecZ,bins='rice')
#plt.hist(vecZ,range=(0,100000),bins=100)
plt.xlim(0,1.1*max(1e6*vecZ))
plt.xlabel('10^6*Z')
plt.ylabel('frequence')
plt.grid()
titre = "Distribution de Z=f(X,Y)=X/Y"
plt.title(titre)
resultat_MC = 'MC: Moyenne ='+format(1e6*Zm,'.3f')+'; Incert. relat. ='+format(100*uZ/Zm,'.1f')+'%'
resultat_diff = 'Diff: Moyenne ='+format(1e6*f(xm,ym),'.3f')+'; Incert. relat. ='+ \
               format(100*np.sqrt(uX**2/xm**2+uY**2/ym**2),'.1f')+'%'
plt.annotate(resultat_MC,(0,1300))#(40850,700))
plt.annotate(resultat_diff,(0,900))#(40850,500))
plt.tight_layout()
plt.savefig('hist_f(X,Y)_1.png')
plt.show()
#plt.clf()