% This is to illustrate diffusion of particles (i.e., Brownian motion) with
% restoring force (Ornstein-Ulenbeck)
clear

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;


K = 4000; %number of particles
N = 200; %number of time steps 
D = 1.0; % diffusion coefficient
v = 1;  % if v=0 this is normal diffusion.  If v ne 0, this is an OU process
x = zeros(K,1); % initial positions

X = x; % store the trajectories of the particles here
dt = .025;  % time step
t = 0; %initial time
var(1) = 0; 
T(1) = 0;
 
diff = sqrt(2*D*dt);
 
for j = 2:N % each time step
     
    x = x -v*x*dt+ diff*randn(K,1);  % move the particles 
    % notice that v ne 0 gives a restoring force
    X = [X,x];  % save the answer
   
      t = t +dt;
      
    T(j) = t; 
    var(j) = sum(x.^2);
end

%Make a histogram of the position of the particles
    figure(1)
    [NN,Y]=hist(x,50);
    dY = Y(2)-Y(1);
      
     if (v~=0)
         sigsq = (D/v)*(1 -exp(-2*t*v ));
     else sigsq = 2*D*t
     end
     
          p = exp(-Y.^2/(2* sigsq))/(sqrt(2* pi*sigsq));
   
         
    plot(Y,NN/(K*dY),'*',Y,p,'--')
    title('Histogram of particle position','fontsize',20)
    xlabel('x','fontsize',20)
    ylabel('Particle count','fontsize',20)
    
     % find slope of msd
    slope =sum(var.^2)/(sum(T.*var)*K)
    
    figure(2)
    if(v==0)
    plot(T,var/K,T,slope*T,'--','linewidth',2)
    else
       
    plot(T,var/K,T,(D/v)*(1 -exp(-2*T*v )),'--','linewidth',2)
    end
     xlabel('t','fontsize',20)
    ylabel('Variance','fontsize',20)
   if(cmykflg==1)
       if(v==0)
print('../../figs_c/chapt_4/particle_var','-deps','-cmyk')
       else
         print('../../figs_c/chapt_10/particle_var','-deps','-cmyk')
       end
    end
    figure(3)
    for j = 1:10
        plot(T,X(j,:),'linewidth',2)
        hold on
    end
    hold off
     xlabel('t','fontsize',20)
    ylabel('x','fontsize',20)
    if(cmykflg==1)
        if(v==0)
           
print('../../figs_c/chapt_4/particle_diff','-deps','-cmyk')
        else
            print('../../figs_c/chapt_10/particle_path','-deps','-cmyk')
        end
    end
 
    
    
        
 