# -*- coding: utf-8 -*-
"""
Propagation des incertitudes par methode Monte-Carlo :
exemple de y = x1x2 /(x1+x2)
"""

# Bibliotheques
import random as rd                # Pour tirage aleatoire
import numpy as np                 # Pour sqrt, log, cos...
import matplotlib.pyplot as plt    # graphes

#valeur des donnees et leurs incertitudes dans la meme unite !!
x1 = 10 
u_x1 = 0.5
x2 = 5
u_x2 = 0.3



# Calcul de l'incertitude sur y 
N = 1000        # nb d'iterations Monte Carlo
L_y = []        # initialisation liste valeurs de y calculees
for i in range(N):
    #tirage aléatoire uniforme de x1_MC de valeur moyenne x1 et dans l'intervalle +/- sqrt(3)*u_x1
    x1_MC = x1 + rd.uniform(-1,1)*np.sqrt(3)*u_x1 
    #tirage aléatoire uniforme de x2_MC de valeur moyenne x2 et dans l'intervalle +/- sqrt(3)*u_x2
    x2_MC = x2 + rd.uniform(-1,1)*np.sqrt(3)*u_x2 
    y_MC = x1_MC * x2_MC / (x1_MC + x2_MC)             # calcul de y
    L_y.append(y_MC)                      # stocker cette valeur


ymoy = np.mean(L_y)                    # calcul valeur moyenne de L_y
u_y = np.std(L_y,ddof=1)               # calcul ecart-type de L_y

# Afficher resultat
print('y= {} +/- {} unité'.format(ymoy,u_y))


# Histogramme des valeurs de y 
plt.hist(L_y,bins='rice') 
plt.xlabel('valeur de y')
plt.show()