import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


T=100 #Nombre pas de temps
N=250 #Nombre particule
r=[[] for i in range(N)]
for n in range(N):
    x=[0,]
    y=[0,]
    for i in range(T):
        angle=np.random.uniform(0,2*np.pi) 
        x=np.append(x,x[-1]+1*np.cos(angle))
        y=np.append(y,y[-1]+1*np.sin(angle))
    r[n]=np.append(r[n],x**2+y**2)
    plt.plot(x,y)
    plt.axis('equal')
    teta=np.linspace(0,2*np.pi,100)
X_c=np.sqrt(np.mean(r,axis=0)[-1])*np.cos(teta)
Y_c=np.sqrt(np.mean(r,axis=0)[-1])*np.sin(teta)
plt.plot(X_c,Y_c,'r',linewidth=4)
plt.plot([0,X_c[-10]],[0,Y_c[-10]],'r',linewidth=4)
plt.text(40,53,r'$\propto t$',fontsize=14)
plt.title('Diffusion active de '+str(N)+' algues pendant '+str(T)+' pas de temps')
plt.xlabel('x (en unité arbitraire)')
plt.ylabel('y (en unité arbitraire)')

plt.figure()    
t=np.arange(0,T+1)
def f(x,a,b):
    return a*x+b
popt, pcov = curve_fit(f, t, np.mean(r,axis=0))
plt.plot(t,np.mean(r,axis=0))
plt.plot(t,f(t,popt[0],popt[1]),linewidth=3)
plt.text(40,53,r'$\propto t$',fontsize=14)
plt.xlabel('Temps (en unité arbitraire)')
plt.ylabel('MSD(t) (en unité arbitraire)')
plt.title('MSD de 1000 algues au cours du temps - Diffusion passive')