% The aim of this matlab is to simulate the displacement of n particules
% for x steps of 1 second each. This simulation is here to help understand
% the principle of diffusion, and the applicability of mean displacement to
% find the coefficient of diffusion of a system. All equations are coming
% from professor Kolinski course, Experimental Method in Engineering
% Mechanics, given at EPFL in the Mechanical Engineering section.

clc
clear all
close all
%%
%Those are the variable you have to choose.PLEASE ONLY CHANGE THOSE
%VARIABLES
number = 300; %The number of particules to simulate
total_step = 40; %The number of steps to simulate
displ = 1; %The displacement of one step 

%Please note that for those simulations, t = 1 seconde. So, D = displ^2/2 !
%%
for i =1:number
    x_2(i,1) = 0;
    y_2(i,1) = 0;
end

for step = 2:(total_step+1)
    for i = 1:number
        %Each step, the particule goes -1 or +1 both in the x and y direction
        x_2(i,step) = x_2(i,step-1) + (randi([1,2],1,1)*2-3)*displ;
        y_2(i,step) = y_2(i,step-1) + (randi([1,2],1,1)*2-3)*displ;
    end
end

figure(1);
hold on
%A color is randomly picked for each particule
for i=1:number
h(i) = animatedline('color',[rand rand rand],'LineStyle',':');
end

%We draw the displacement of each particule
for k = 1:(total_step+1)
    tic;
    for i=1:number
        addpoints(h(i),x_2(i,k),y_2(i,k));
        drawnow
    end
    toc;
end
%We draw the final position of the particules
plot(x_2(:,total_step+1),y_2(:,total_step+1),'d','color','red')
hold off

%We calculate the mean displacement for each particule at each step !
for step = 1:(total_step+1)
    for i = 1:number
        r(i,step) = x_2(i,step)^2 + y_2(i,step)^2;
    end
    mean_x(step) = 1/number * sum(x_2(:,step).^2);
    mean_y(step) = 1/number * sum(y_2(:,step).^2);
    mean_r(step) = mean_x(step)+mean_y(step);
end

figure(2);
hold on
plot(0:total_step,mean_r)
ylabel('Mean displacement squarred, <r^2>')
xlabel('Lag time, T')

%This is the theoretical D you should obtain
D = displ^2/2

%This is the D obtained by doing the slope of the curve, only taking into
%account first and last point
D_simu = (mean_r(total_step)-mean_r(1))/total_step/4

p=polyfit(0:total_step,mean_r,1);
plot(0:total_step,polyval(p,0:total_step));
%This is the D obtained by linearizing all the points
D_simu_lin = p(1,1)/4