#On essaie ici de modeliser la marche au hasard vue en cours ; le processus est ici bidimensionnel (les courbes seront plus jolies) ; la seule différence est l'expression du coefficient qui sera D=a**2/(4*tau)


import numpy as np
from scipy.integrate import odeint
from scipy import fft
import matplotlib.pyplot as plt

#Variables d'une trajectoire
a = 1 # distance parcourue à chaque saut
nd = 1001 # nb de sauts et durée d'un saut
Tmax = 1 # durée totale de la marche au hasard

t=np.linspace(0.0, Tmax, nd) # instants de mesure

Delta_t = Tmax/(nd-1) # pas de temps numérique

D_th=a**2/(4*Delta_t) # Diffusion théorique attendue




##IMPLEMENTATION NUMERIQUE

#Pour y parvenir on va commencer par créer et visualiser une première trajectoire utilisant la fonction random.choice() (doc.) qui permet de tirer au hasard, avec une distribution uniforme, une valeur comprise entre 4 choix possible (vers le haut, vers le bas, à droite ou à gauche.)
#Les coordonnées (x,y) de chaque étape seront stockés dans un np.array et la trajectoire visualisée grâce à matplotlib.
#Pour que la modélisation ait un sens il faut multiplier les particules suivies et en tirer des valeurs statistiques moyennes. On crée donc un tableau position identique à pos_1 mais contenant Num particules au lieu d'une seule.


pos_1=np.zeros((nd,2)) # tableau contenant nd couples de deux valeurs (x,y)
for i in range(1,nd):
    step=[[a,0],[-a,0],[0,a],[0,-a]]
    choix=np.array([0,1,2,3])
    ind = np.random.choice(choix,1)[0]#Un pas en avant ou un pas en arrière?
    pos_1[i,:]=pos_1[i-1,:] + step[ind]

plt.figure()
plt.plot(0,0,'o', label='départ')
plt.plot(pos_1[nd-1,0],pos_1[nd-1,1],'x', label='arrivée')
plt.legend()
plt.plot(pos_1[:,0],pos_1[:,1])
plt.grid(True)
plt.show()


#Nombre de trajectoires
Num = 1000 #Nb de trajectoires calculées
position = np.zeros((Num,nd,2)) #Tableau contenant toutes les positions 2D au cours

# Simulation des trajectoires
for j in range(Num):
    for i in range(1,nd) :
        step=[[a,0],[-a,0],[0,a],[0,-a]]
        choix=np.array([0,1,2,3])
        ind = np.random.choice(choix,1)[0]#Un pas en avant ou un pas en arrière?
        position[j,i,:] = position[j,i-1,:] + step[ind]


#On trace les trajectoires individuelles, en veillant à conserver le rapport d'aspect des deux axes du graphique grâce à la commande plt.axis('equal').
#Attention le tracé peut-être long si Num est élevé.


# Tracé de la trajectoire de chaque particule suivie
plt.figure('Evolution spatiale individuelle')
plt.axis('equal')
plt.grid(True)
j = 0
for j in range(Num):
    x=position[j,:,0]
    y=position[j,:,1]
    plt.plot(y,x)
plt.xlim(-80,80)
plt.ylim(-80,80)
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.show()



#On peut alors faire une évaluation statistique du coefficient de diffusion pour en donner une valeur expérimentale, accompagnée de son incertitude-type.


# Tracé de l'évolution temporelle de la moyenne instantannée et quadratique du coefficient de diffusion

#C'est le rayon moyen que l'on va calculer en mesurant, à chaque instant, la moyenne des distances l**2=x**2+y**2 parcourue par les diverses trajectoires. Cette moyenne sera stockée à chaque instant dans la liste xyquad, puis tracée dans le graphe ci-dessous.

xyquad=np.zeros(nd)
for i in range(nd):
    xyquad[i]=np.mean(position[:,i,0]**2+position[:,i,1]**2)

Diff=xyquad[1:]/t[1:]

D_mean = np.mean(Diff/4)
D_u = np.std(Diff/4, ddof=1)

plt.figure()
plt.hist(Diff/4,bins='rice')
plt.grid(True)
plt.xlabel('Valeur de D')
plt.ylabel('Occurrences')
plt.show()

print("La valeur moyenne mesurée du coefficient de diffusion est ",format(D_mean,"#.3f"),"+/-",format(D_u,"#.3f"),"m^2/s")
print("La valeur théorique du coefficient de diffusion est ",format(D_th,"#.3f"),"m^/s")