

c=50; %concentration in particles per micron^2
sigma=.1; %micron
D=0.1; %micron^2/s
L=10; %size of membrane in micron
KT=1.38*10^-23*300; %joules
Tsteps=100000;
Nparticles=c*L^2;
deltat=.01*sigma^2/D; %second


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

B=2*KT;
%make figure to see positions
h=figure;
plot(X,Y,'o')
axis([-L/2 L/2 -L/2 L/2])
%coordinates for a circle
r=2*sigma;
theta=0:.1:2*pi;
xc=r*sin(theta);
yc=r*cos(theta);
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
    U=-B*exp(-X.^2/(2*sigma^2)).*exp(-Y.^2/(2*sigma^2));
    delX_Drift=(D*deltat/KT)*X.*U/sigma^2;
    delY_Drift=(D*deltat/KT)*Y.*U/sigma^2;
    
    X=X+delX_Dif+delX_Drift;
    Y=Y+delY_Dif+delY_Drift;
    
    %apply periodic boundary conditions
    X(X<(-L/2)) =X(X<(-L/2))+L;
    X(X>(L/2))  =X(X>(L/2)) -L;
    Y(Y<(-L/2)) =Y(Y<(-L/2))+L;
    Y(Y>(L/2))  =Y(Y>(L/2)) -L;
    
    Nwait=100; %plot every 100 updates
    if (tt/Nwait)==floor(tt/Nwait)  
        %count the number withing 2 sigma
        N2sig=sum( (X.^2+Y.^2)<((2*sigma)^2));
        plot(X,Y,'o')
        hold on
        plot(xc,yc,'r');hold off
        s=sprintf('Number within 2 sigma: %d',N2sig);
        title(s)
        axis equal
        axis([-L/8 L/8 -L/8 L/8])
        pause(.001)
    end
    
end


%% Compare number with r<2*sigma to that calculated from principle of constant chemical potential
c=50; %concentration in particles per micron^2
sigma=.1; %micron
KT=1.38*10^-23*300; %joules

dr=sigma/100;
r=(0:dr:2*sigma);

B=1*KT;
Nexpected=sum(dr*r.*2*pi.*c.*exp(1/KT*B*exp(-r.^2/(2*sigma^2))))

B=2*KT;
Nexpected=sum(dr*r.*2*pi.*c.*exp(1/KT*B*exp(-r.^2/(2*sigma^2))))

B=5*KT;
Nexpected=sum(dr*r.*2*pi.*c.*exp(1/KT*B*exp(-r.^2/(2*sigma^2))))



