% BD_1D.m
% This MATLAB program performs Brownian dynamics
% simulation of a single particle in a potential
% field U(x). A simple explicit Euler method is
% used to propagate the Langevin SDE.
% K. Beers. MIT ChE. 10/27/03
function iflag_main = BD_1D(t_end,dt);
iflag_main = 0;

% First, set the parameters
D = 1;  % diffusivity
zeta = 1;  % drag coefficient
k_sp = 1;  % spring constant for potential energy

% Next, set the vector of time values at which
% we compute the position of the particle.
t_vect = [0:dt:t_end];  num_t_val = length(t_vect);
x_vect = zeros(size(t_vect));

% perform BD iterations to compute trajectory,
% and periodically report result
x = 0;
var_sqrt_2D = sqrt(2*D);  var_sqrt_dt = sqrt(dt);
for iter=1:num_t_val
    x_vect(iter) = x;
    % compute deterministic force from Hookian spring
    F_spring = -k_sp*x;    
    % compute increment of Wiener process
    dWt = var_sqrt_dt*randn;
    % compute BD update
    x = x + F_spring/zeta*dt + var_sqrt_2D*dWt;
end

% Now, plot random trajectory
figure;  plot(t_vect,x_vect);
xlabel('t');  ylabel('x(t)');
title_phrase = ['Brownian dynamics of particle in 1-D quadratic well', ...
        ', \Deltat = ', num2str(dt)];
title(title_phrase);

% Now, we compute the histogram of the trajectory, and plot this in
% comparison to the Boltzmann distribution with kT = 1
figure; num_bins = 50;
[n_hist,x_hist] = hist(x_vect,num_bins);
var1 = trapz(x_hist,n_hist);  n_hist = n_hist ./ var1;
bar(x_hist,n_hist);
x_plot = linspace(-3,3,500);  V_plot = 0.5*k_sp.*(x_plot.^2);
kb_T = 1;  p_Boltzmann = exp(-V_plot/kb_T);
var1 = trapz(x_plot,p_Boltzmann);  p_Boltzmann = p_Boltzmann ./ var1;
hold on;  plot(x_plot,p_Boltzmann);
xlabel('x');  ylabel('p(x)');
title('Probability distribution function of x for 1-D Brownian motion');

iflag_main = 1;
return;