# CAPACITE 1 : MARCHE AU HASARD A UNE DIMENSION
###################################################################################################################################
# Cellule 1 : Importation des bibliothèques et initialisation des données

# N particules sont placées en x = O à t = 0. 
# La diffusion est modélisée par une marche au hasard avec un libre parcours moyen l et une durée tau entre chaque pas.
# A chaque nouveau pas, la particule se déplace de -l ou +l selon l'axe Ox avec la même probabilité.

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

N = 1000                    # Nombre de particules
l = 0.75                   # Libre parcours moyen (distance parcourue à chaque pas)
tau = 1                    # Durée entre deux collisions (durée d'un pas)
P = 1000                   # Nombre de pas

X = np.zeros(N)             # Tableau des N abscisses initiales nulles
Y = np.linspace(0,1,N)      # Tableau des N ordonnées constantes des N particules, réparties régulièrement entre Y = 0 et Y = 1

# Après chaque pas, on calcule l'abscisse quadratique moyenne : xqm = sqrt(<x²>) (< > moyenne sur les N particules)
# xqm va mesurer la distance caractéristique L de la diffusion des particules. 
# A l'instant t, la théorie de la marche au hasard à une dimensions montre que xqm = sqrt(2Dt) où D est le coefficent de diffusion.
# Le modèle probabiliste développé en cours donne : D = l²/2tau 
# Donc : L² = xqm² = 2Dt = l²t/tau 
# On se propose de vérifier cette loi et d'en déduire D.

Xmoy = np.zeros(P+1)          # Initialisation du tableau des P+1 abscisses moyennes   (1 valeur initiale et P valeurs après chaque pas)
Xqm = np.zeros(P+1)           # Initialisation du tableau des P+1 abscisses quadratiques moyennes (1 valeur initiale et P valeurs après chaque pas)
###################################################################################################################################
# Cellule 2 : Animation de la diffusion des N particules 

start = time.time()          # Lancement d'un chronomètre pour mesurer le temps de calcul

for i in range(1,P+1):       # Pour chaque nouveau pas i (P valeurs de i allant de 1 à P)

    for k in range(N):       # Pour chaque particule k   (N valeurs de k allant de 0 à N-1)
    
        # Lors du pas numéro i, l'abscisse de la particule k varie aléatoirement de +l ou -l avec la même probabilité 1/2
        X[k] = X[k] + random.choice(np.array([-l,l]))  
        
    Xmoy[i] = np.mean(X)               # Moyenne des abscisses des N particules après le pas i 
    Xqm[i] = np.sqrt(np.mean(X**2))    # Abscisse quadratique moyenne au centre O des N particules après le pas i    
    
    plt.clf()                          # Effacement de l'image précédente pour obtenir une animation
    plt.grid()
    plt.axis([-120,120,0,1])
    plt.xlabel('x')
    plt.plot(X,Y,'o',markersize=1)     # Tracé du nuage de points après le pas i
    plt.pause(0.001)                   # Durée entre deux images
    plt.show()
    
end = time.time()         # Fin du chronomètre     
Tcalcul = end - start     # Temps de calcul
print('Le temps de calcul est :',Tcalcul)
####################################################################################################################################
# Cellule 3 : Histogrammes des N abscisses x des N particules 

plt.title("Histogramme de l'abscisse x")
plt.xlabel('x')
plt.hist(X,bins='rice')      # Histogramme des abscisses après les P pas

plt.show()
####################################################################################################################################
 # Cellule 4 : Moyenne des abscisses et abscisse quadratique moyenne en fonction du temps

t = np.arange(0,(P+1)*tau,tau)     # Tableau des temps [0,tau,2tau,3tau,....,Ptau]

plt.subplot(1,2,1)
plt.title("Moyenne des abscisses au cours du temps")
plt.xlabel('temps')
plt.ylabel('x moyen')
plt.plot(t,Xmoy)        
plt.grid()

plt.subplot(1,2,2)
plt.title("Abscisse quadratique moyenne au cours du temps")
plt.xlabel('temps')
plt.ylabel('abscisse quadratique moyenne')
plt.plot(t,Xqm)
plt.grid()

plt.show()
#####################################################################################################################################
# Cellule 5 : Vérification de la loi L² = 2Dt pour un phénomène de diffusion à une dimension

L = Xqm           # Définition de la distance caractéristique L de diffusion, donc d'étalement du nuage de points
plt.plot(t,L**2,'+r',markersize=1)  

p = np.polyfit(t,L**2,1)                                         # Regression linéaire de L² = f(t)
plt.plot(t,p[0]*t+p[1],'k--',linewidth=1,label='Modèle affine')  # Tracé de la droite modèle L² = at + b

plt.title('L² en fonction du temps')
plt.xlabel('temps')
plt.ylabel('L²')
plt.legend()
plt.grid()

print('Le coefficient de diffusion est :',p[0]/2)  
D = l**2/2*tau                                                    # Valeur théorique du coefficient de diffusion
print('La valeur théorique du coefficient de diffusion est :',D) 
   
plt.show()
######################################################################################################################################