# -*- coding: utf-8 -*-
"""TP_3_cinétique_complété (1).ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/15YuceoRF0klIdNdFJS7xaTq0jzXDC7cy

#Présentation expérimentale de la cinétique de décoloration du bleu de bromophénol

## Suivi cinétique

Le bleu de bromophénol est fortement coloré, alors que le composé ayant réagi avec les ions hydroxyde est incolore. Cette propriété va être mise à profit pour suivre la cinétique de la réaction grâce au spectrophotomètre.

La réaction qui a lieu peut être symbolisée par :
				$ BBPH^{-}   +   HO^{–}    \to  BBP^{2-} + H_2O$

On propose d’écrire la loi de vitesse s’écrit selon $v = k [BBPH^{-}]^{\alpha} \cdot [HO^{–}]^{\beta}$

## Rappel de la loi de Beer - Lambert

La spectrophotométrie est l'étude de l'interaction entre la matière et le rayonnement.
Nous profitons donc du fait que le bleu de bromophénol est une molécule colorée et absorbe dans le visible selon la loi de Beer - Lambert.

- L'absorbance d'une solution caractérise son aptitude à absorber une radiation de longueur d'onde donnée. Elle se note $A$ et n'a pas d'unité.

- Soit $I_0$ l'intensité lumineuse en entrée de la cuve contenant la solution, soit $I_{sortant}$ celle en sortie de cuve, on écrit : $\mathrm {A = - log_{10} \dfrac {I_0}{I_{sortant}}}$

- Ici, nous n'avons qu'une seule espèce colorée : le bleu de bromophénol : $\mathrm {A = \varepsilon_i(\lambda) \ell C_i }$
  - la concentration $C_i$ est la concentration molaire de l'espèce colorée $i$,
  - $\ell$ est la longueur de la cuve en $cm$ (souvent $\ell = 1,0 \ cm$)
  - $\varepsilon_i( \lambda)$ est le coefficient d'extinction molaire en $mol^{-1}.L.cm^{-1}$ et dépend de la longueur d'onde de travail $\lambda$

# Etude mathématique

## Ordre partiel $0$ supposé par rapport au $BBPH^{-}$ renommé B dans la suite

- Equation cinétique : $\mathrm { - \dfrac {d[\mathsf {\text {$B$}}]}{\mathsf {\text {$dt$}}} = \mathsf {\text {$k$}}^{'}}$ qui donne après intégration entre l'instant initial $t = 0$ et l'instant $t$ : $\mathrm { [\mathsf {\text {$B$}}] = [\mathsf {\text {$B$}}]_0 - \mathsf {\text {$k^{'}$}}\mathsf {\text {$t$}} }$.
Comme l'absorbance est proportionnelle à la concentration, on peut tracer également $ A = A_0 - \varepsilon (\lambda) \ell k^{'}t$ ou plus simplement $A = f(t)$

## Ordre partiel $1$ supposé par rapport au $BBPH^{-}$

- Equation cinétique : $\mathrm { - \dfrac {d[\mathsf {\text {$B$}}]}{\mathsf {\text {$dt$}}} = \mathsf {\text {$k$}}^{'}[\mathsf {\text {$B$}}]}$ qui donne après intégration entre l'instant initial $t = 0$ et l'instant $t$ : $\mathrm { ln ( \dfrac {[\mathsf {\text {$B$}}]}{[\mathsf {\text {$B$}}]_0}) =  - \mathsf {\text {$k^{'}$}}\mathsf {\text {$t$}} }$
Comme l'absorbance est proportionnelle à la concentration, on peut tracer également $\mathrm { ln (\dfrac {[\mathsf {\text {$A$}}]}{[\mathsf {\text {$A$}}]_0}) =  - \mathsf {\text {$k^{'}$}}\mathsf {\text {$t$}} }$ ou plus simplement $lnA = f(t)$

## Ordre partiel $2$ supposé par rapport au $BBPH^{-}$

- Equation cinétique : $\mathrm { - \dfrac {d[\mathsf {\text {$B$}}]}{\mathsf {\text {$dt$}}} = \mathsf {\text {$k$}}^{'}[\mathsf {\text {$B$}}]^{2}}$ qui donne après intégration entre l'instant initial $t = 0$ et l'instant $t$ : $\mathrm { ( \dfrac {1}{[\mathsf {\text {$B$}}]}- \dfrac {1}{[\mathsf {\text {$B$}}]_0}) =   \mathsf {\text {$k^{'}$}}\mathsf {\text {$t$}} }$.

Comme l'absorbance est proportionnelle à la concentration, on peut tracer également  simplement $1/A = f(t)$ en faisant attention au coefficient de proportionnalité $\varepsilon (\lambda) \ell $

# Suivi cinétique de la solution
"""

# importation des biliothèques
import matplotlib.pyplot as plt
import numpy as np
from math import *

A = np.array([0.815, 0.753, .691, .639, .585, .537, .494, .454, .417, .383, .351, .322, .298, .274, .252, .230, .210, .195, .176, .162, .148, .135, .126, .114, .105, .097 ]) # ordre 0

time = np.array([0.17, 0.67, 1.17, 1.67, 2.17, 2.67, 3.17, 3.67, 4.17, 4.67, 5.17, 5.67, 6.17, 6.67, 7.17, 7.67, 8.17, 8.67, 9.17, 9.67, 10.17, 10.67, 11.17, 11.67, 12.17, 12.67]) # temps en min

Z1 = np.log(A) # ordre 1

Z2 = (1/A) # ordre 2

# Affichage graphique : ordre 0 supposé
plt.figure(1)
plt.title(r'Tracé de $A$ en fonction du temps ')
plt.ylabel(r'$A$')
plt.xlabel('temps en min (min)')
plt.plot(time, A,'r+', label='points expérimentaux')
plt.xlim(0,20)
plt.ylim(0,1)
plt.grid(linestyle="-.")
plt.legend()
plt.show()

# Affichage graphique ; ordre 1 supposé
plt.figure(2)
plt.title(r'Tracé de $lnA$ en fonction du temps ')
plt.ylabel(r'$lnA$')
plt.xlabel('temps en minutes (min)')
plt.plot(time, Z1,'r+', label='points expérimentaux')
plt.grid(linestyle="-.")
plt.legend()
plt.show()

# Affichage graphique : ordre 2 supposé
plt.figure(5)
plt.title(r'Tracé de $1/A$ en fonction du temps ')
plt.ylabel(r'$1/A$')
plt.xlabel('temps en minutes (min)')
plt.plot(time, Z2,'r+', label='points expérimentaux')
#plt.xlim(0,15)
#plt.ylim(0,1)
plt.grid(linestyle="-.")
plt.legend()
plt.show()

"""# Régression linéaire sur la courbe sélectionnée"""

p1 = np.polyfit(time, Z1, 1)                      #Régression linéaire de lnA en fonction de t

delta_A =(2/100)*A # demi - étendue du spectrophotomètre donné par le constructeur
u_A = delta_A/np.sqrt(3) # incertitude - type par rapport à l'absorbance
u_Z1=u_A/A # incertitude - type par rapport au calcul du ln : axe des Y

modele1=("ln A = "+ format(p1[0], "#.3e")+"x t +"+ format(p1[1], "#.3g"))   # on fabrique l'équation de la droite de régression
plt.plot(time,np.polyval(p1,time),color='green',label=modele1)

plt.plot(time, np.polyval(p1,time))
plt.errorbar(time, Z1, yerr=u_Z1, fmt='.r')   #On ajoute les barres d'incertitude
plt.xlabel('t/min')
plt.ylabel('ln A' )
plt.title('Évolution de $lnA$ en fonction de t')
plt.legend()
plt.grid()
plt.show()

print("kapp =",-p1[0],'min^-1')                   #Affichage de kapp, l'opposé de la pente de la droite de régression

"""Pour s'assurer efficacement que la droite de régression passe par les barres d'incertitude, on peut tracer :
- les résidus avec les barres d'incertitude ;
- les résidus normalisés directement.

Notre critère de décision est le dépassement, en valeur absolue, de la valeur 2 pour l'écart normalisé.
"""

#Tracé des résidus

plt.figure(3)
plt.title("Résidus")
plt.ylabel('lnA - modèle affine')
plt.xlabel('temps en (s)')
plt.errorbar(time,Z1-np.polyval(p1,time), yerr = u_Z1, fmt = 'ro',label="Résidus")
plt.plot((np.min(time), np.max(time)), (0, 0) )
plt.grid(linestyle="-.")
plt.legend()
plt.show()

# Tracé des résidus normalisés z
plt.figure(4)
plt.title('Résidus normalisés')
z = (Z1-np.polyval(p1,time))/(u_Z1)
plt.plot(time, z,'bo')
plt.xlabel('temps en (s)')
plt.ylabel('Résidus normalisés de $lnA$' )
plt.plot((np.min(time), np.max(time)), (0, 0) )
plt.grid(linestyle="-.")
plt.show()

"""On valide le résultat si, sur le graphe des résidus normalisés, les points représentés sont inférieurs à 2. Cela signifie qu'il y a moins de deux incertitudes types entre le modèle affine et les points expérimentaux.
Si ce n'est pas le cas, il est possible que :
* le modèle d'ordre choisi ne permette pas de décrire correctement ces résultats expérimentaux et la loi est autre ;
* certaines incertitudes aient été sous-estimées ;
* une erreur de manipulation ait été commise.

# Méthode de MC pour la régression linéaire
# Méthode de Monte Carlo
Une autre méthode (appelée Monte-Carlo) consiste à réaliser un grand nombre d'expériences virtuelles, en faisant varier les paramètres qu'on sait ne pas être connus de manière certaine. Pour simuler un résultat de mesure (l’ensemble des valeurs raisonnablement attribuables à une grandeur), on génère un ensemble de valeurs aléatoires. Selon les cas, pour modéliser les grandeurs mal connues, on choisit une distribution de probabilité  ou uniforme :
- normale (gaussienne) si on a des raisons de connaître l'écart-type de la distribution ;
- uniforme si on ne pense que connaitre l'étendue de la distribution.

### Nombres aléatoires

On utilise le sous-module `random` de la bibliothèque Numpy :
1. dans sa version la plus simple, l' expression  `rd.uniform(x-a, x+a, N)` permet de créer un array de $N$ valeurs tirées aléatoirement entre $x−a$ et $x+a$ et issues d'une distribution uniforme (autrement dit rectangulaire) ; ici `x`et `a` sont des scalaires ;
2. l' expression  `X + rd.uniform(-a, +a, N)` permet de rajouter à `X`, un scalaire __ou un `array`__, un array de $N$ valeurs tirées aléatoirement entre $−a$ et $+a$, issues d'une distribution uniforme ; cette expression permet de gérer le cas où on a un ensemble de différentes valeurs mesurées, auxquelles on rajoute des variables aléatoires uniformes ayant toute la même étendue `a` ;
3. l' expression  `X + A * rd.uniform(-1, +1, N)`  permet de gérer le cas où on a un ensemble de différentes valeurs mesurées, auxquelles on rajoute des variables aléatoires uniforme de __différentes__ étendues, contenues dans un array `A` ; ainsi peut-on prendre en compte le fait que la « précision » varie en fonction de la valeur de la lecture à l'appareil.

4. l'expression `rd.normal(loc=0.0, scale=1.0, size=None)` crée un array de taille `size` qu’on peut préciser (si on ne le précise pas on obtient un scalaire), contenant des valeurs tirées aléatoirement issues d’une distribution gaussienne d’espérance `loc` et d’écart-type `scale`.
    - C’est ce qu’on utlise habituellement si on a un processus qu’on peut raisonnablement penser comme Gaussien, par exemple en application du théorème de la limite centrale (un processus additif de nombreuses variables aléatoires).
    - On fait souvent **l’hypothèse gaussienne lorsqu’on souhaite modéliser un résultat de mesurage dont on ne connait pas la distribution exacte, mais dont on connait une estimation de l’espérance et de l’écart-type** (la valeur mesurée et l’incertitude-type, en termes métrologiques).
    
### Détermination de l'incertitude - type des paramètres de la droite de régression

* Le résultat d'une mesure est l'ensemble des valeurs raisonnablement attribuables à la grandeur mesurée. On caractérise cet ensemble par la valeur mesurée (parfois une moyenne), et l'écart-type expérimental ; ce dernier est par définition l'incertitude-type associée à la valeur mesurée.

* L'idée pour estimer l'incertitude-type associée à la pente de la droite   consiste à simuler *in silico* le résultat de $N$ expériences virtuelles. Pour cela on  ajoute aux observations de la grandeur portée en $Y$ des valeurs aléatoires dont l'écart-type est égal à l'incertitude-type qu'on a évalué lors de l'expérience. On fait la même chose pour la grandeur portée en $X$.

* Pour chaque graphe de $Y$ en fonction de $X$, on fait une régression linéaire et on obtient une valeur la pente. L'ensemble des valeurs de la pente  représente l'ensemble de des valeurs raisonnablement attribuables à cette pente, donc le résultat de la mesure.

* On la caractérise par :
    * sa valeur moyenne, autrement dit la valeur mesurée ;
    * son écart-type, autrement dit l'incertitude-type.

* Bien entendu, pour chaque graphe on peut récupérer également l'ordonnée à l'origine (oo) dont on peut calculer la moyenne et l'écart-type.
"""

from numpy import random as rd
# importation des biliothèques
import matplotlib.pyplot as plt
import numpy as np
from math import *

A = np.array([0.815, 0.753, .691, .639, .585, .537, .494, .454, .417, .383, .351, .322, .298, .274, .252, .230, .210, .195, .176, .162, .148, .135, .126, .114, .105, .097 ])

time = np.array([0.17, 0.67, 1.17, 1.67, 2.17, 2.67, 3.17, 3.67, 4.17, 4.67, 5.17, 5.67, 6.17, 6.67, 7.17, 7.67, 8.17, 8.67, 9.17, 9.67, 10.17, 10.67, 11.17, 11.67, 12.17, 12.67])

#Proposition des valeurs des précisions :
p_time = 0.2/(2*60) # demi - étendue du déclenchement et arrêt du chronomètre
delta_A=(2/100)*A # demi - étendue de l'absorbance donnée par le constructeur


# Paramètres
N_sim = 20000     # Nombre d'itérations souhaitées
N = len(time)        # Nombre de données de notre expérience
p1_sim=np.zeros(N_sim)


for i in range(0, N_sim):   # Création d'une boucle for permettant de générer les N_sim expériences virtuelles

    time_sim = time + rd.uniform(-p_time, p_time, N)                 # Génère les N valeurs aléatoires de pH selon une loi rectangulaire

    A1_sim = A + delta_A*rd.uniform(-1, 1, N)        # Si on connaissait une précision, on génererait les N valeurs aléatoires de A selon une loi rectangulaire

    Z1_sim = np.log(A1_sim)

    p1_sim [i]= np.polyfit(time_sim, Z1_sim, 1)[0]


plt.figure(2)
plt.title('Pour 20 000 iterations')
plt.hist(-p1_sim, bins='rice', color='brown', label = "$k'_{1}$") # tracé d'un histogramme
plt.xlabel(r"$k'_{1}$ ")
plt.ylabel('effectif')
plt.legend()
plt.show()

print("k'1 = ", format(np.average(-p1_sim), "#.5g"), 'min^-1')
print("u(k'1) = ", format(np.std(-p1_sim, ddof=1), "#.3g"), 'min^-1')

print("k'1 = ", format(np.average(-p1_sim), "#.4e"), 'min^-1')
print("u(k'1) = ", format(np.std(-p1_sim, ddof=1), "#.3g"), 'min^-1')