% 2 d diffusion

clear all
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

N = 1000; % number of particles
% initial position
c = zeros(N,2);  % the particle positions are stored in c
 
dt = 0.01;  %time step
D = 1; % specify the diffusion coefficient  
sc = sqrt(2*D*dt);
 
n_time_steps = 1000;  %specify the number of time steps
mn = [];
T = [];


for j = 1:n_time_steps % number of time steps
displ = sc*randn(N,2);  % the displacement is normally distributed
%take a forward time step
 
cnew = c +   displ;
c = cnew;
rsq = c(:,1).^2+c(:,2).^2;  % find the squared displacement
mn(j) = sum(rsq)/N;  % mean squared displacement
T(j) = j*dt;
end

figure(1)
plot(c(:,1),c(:,2),'*')
xlabel('x','fontsize', 20)
ylabel('y','fontsize',20)

 if(cmykflg==1)
print('../../figs_c/chapt_4/two_d_diff','-deps','-cmyk')
 end

  
% do statistics on the radial distibution  
rsq = c(:,1).^2+c(:,2).^2;
rstat = sort(rsq);

[p,s] = polyfit(rstat(1:N-1),log(1-[1:N-1]'/N),1)
y=polyval(p,rstat);
figure(2)
plot(rstat,[1:N]/N,'linewidth',2)
xlabel('R^2','fontsize',20)
ylabel('P(R) = fraction of points with r<R','fontsize',20)

 if(cmykflg==1)
print('../../figs_c/chapt_4/two_d_dist','-deps','-cmyk')
 end

  
figure(3)
plot(rstat(1:N-1),log(1-[1:N-1]/N),rstat(1:N-2),-0.025*rstat(1:N-2),'--','linewidth',2)
xlabel('R^2','fontsize',20)
ylabel('ln(1-P(R))','fontsize',20)

 if(cmykflg==1)
print('../../figs_c/chapt_4/two_d_log_dist','-deps','-cmyk')
 end

  

figure(4)
plot(T,mn,'linewidth',2)
xlabel('t','fontsize',20)
ylabel('Mean Squared Displacement','fontsize',20)

 if(cmykflg==1)
print('../../figs_c/chapt_4/two_d_msd','-deps','-cmyk')
  end

