
%% 2D diffusion in gravity

c=50*(10^3)^2; %concentration in particles per meter^2
R=10*10^-7; %meters
density_bead=1.04*1000; %kg/m^3
density_water=1.00*1000; %kg/m^3
n=10^-3; %viscocity, pascal seconds
L=2*10^-3; %width in meters
H=5*10^-3; %width in meters
KT=1.38*10^-23*300; %joules
g=9.8; %m^2/s
D=KT/(6*pi*n*R); %meter^s/s

Tsteps=100000;
Nparticles=c*L*H;
deltat=1; %second

%random starting position
X=L/2-L*rand(Nparticles,1);
Y=H*rand(Nparticles,1);

Fnet=g*(density_bead-density_water)*4/3*pi*R^3; %newtons

%make figure to see positions
h=figure;

for tt=1:Tsteps %update positions
    
    %Diffusion step
    delX_Dif=sqrt(2*deltat*D)*randn(Nparticles,1);
    delY_Dif=sqrt(2*deltat*D)*randn(Nparticles,1);
    
    %Drift step
    delY_Drift=-(D*deltat/KT)*Fnet;
    
    %update
    X=X+delX_Dif;
    Y=Y+delY_Dif+delY_Drift;
    
    %apply boundary conditions
    X(X<(-L/2)) =X(X<(-L/2))+L; %periodic
    X(X>(L/2))  =X(X>(L/2)) -L; %periodic
    Y(Y<0)      =-Y(Y<0);       %reflecting
    Y(Y>H)      =2*H-Y(Y>H);    %reflecting
    
    Nwait=3600/deltat; %plot every Nwait updates (here 1 second)
    if (tt/Nwait)==floor(tt/Nwait)
        figure(h);
        plot(X,Y,'o')
        s=sprintf('Time (hours): %g',tt*deltat/3600);
        title(s,'fontsize',18)
        axis equal
        axis([-L/2 L/2 0 H])
        xlabel('X (meters)')
        ylabel('Y (meters)')
        pause(.001)
    end
    
end




