import numpy as np
import matplotlib.pyplot as plt

### Pour 1 seule particule

def pos1(K):
    # Entrée : K : nombre de pas de temps
    # Sortie : x : liste de taille K des positions de la particule
    
    x=np.zeros(K)
    for i in range(1,K):
        x[i]=x[i-1]+np.random.choice([-1,1])
    return x

### Pour N particules

def posN(K,N):
    # Entrée : K : nombre de pas de temps
    # Entrée : N : nombre de particules
    # Sortie : tab_x : tableau à 2 dimensions des positions des N particules
    # (1 particule par ligne)
    
    tab_x=np.zeros((N,K))
    for i in range(0,N):
        tab_x[i]=pos1(K)
    return tab_x

### Histogrammes

K=500
N=1000  # On choisit N de sorte que la position moyenne soit environ nulle.

tabx=posN(K,N)
xf=tabx[:,K-1]

plt.figure()
plt.hist(xf)
plt.xlabel('Numero du site occupe')
plt.ylabel('Nombre de particules')
plt.title('Histogramme des positions finales')
plt.show()


### Caracterisation de l'etalement des particules

k_list=np.linspace(0,K-1,100)   # Liste des pas de temps pour le trace graphique
# Calcul de la liste des moyennes et des ecarts-types
moy_list=[]
ecart_type_list=[]
for k in k_list:
    moy_list.append(np.mean(tabx[:,int(k)]))
    ecart_type_list.append(np.std(tabx[:,int(k)]))

plt.figure()
plt.subplot(2,1,1)
plt.plot(k_list,moy_list)
plt.title('Position moyenne pour '+str(N)+' particules')
plt.xlabel('Numero du pas de temps')
plt.ylabel('Position')
plt.grid()
plt.subplot(2,1,2)
plt.plot(k_list,ecart_type_list)
plt.title('Ecart-type pour '+str(N)+' particules')
plt.xlabel('Numero du pas de temps')
plt.ylabel('Ecart-type')
plt.grid()
plt.tight_layout()
plt.show()

# Tracé de graphe pour valider le modèle ecart_type=sqrt(k)
plt.figure()
plt.plot(k_list,np.array(ecart_type_list)**2)
plt.title('Ecart-type au carré pour '+str(N)+' particules')
plt.xlabel('Numéro du pas de temps')
plt.ylabel('Ecart-type au carré')
plt.grid()