% Numerical solution of the Langevin equation in one dimension
% Marc R. Roussel (roussel@uleth.ca)

% Clear all variables.
clear

% Calculate frictional coefficient and Gamma
D = 5.7e-10; % sucrose, in SI units
m = 5.68e-25;   % sucrose, in kg
kB = 1.38e-23;
T = 293.15;
f = kB*T/D
Gamma = 2*f*kB*T

% Simulation for tmax time units
tmax = 10e-9; % in seconds
dt = m/(10*f); % Time step

% Initial conditions stored as first point of vectors for later plotting.
npts = 1;
t(npts) = 0;
v(npts) = 0;
x(npts) = 0;

% Simulate until tmax
while (t(npts) <= tmax)
    % Calculate Delta_v
    % Note: randn() generates Gaussian-distributed random numbers.
    Delta_v = -f*v(npts)/m*dt + sqrt(Gamma)*randn(1)/m*sqrt(dt);
    % Update velocity, position and time; store all for future plotting.
    npts = npts + 1;
    t(npts) = t(npts-1) + dt;
    v(npts) = v(npts-1) + Delta_v;
    % Note: Delta_x = v*dt
    x(npts) = x(npts-1) + v(npts-1)*dt;
end

% Plot velocity and position vs time
figure('DefaultAxesFontSize',18)
plot(t,x)
xlim([0 tmax])
xlabel('t/s')
ylabel('x/m')

figure('DefaultAxesFontSize',18)
plot(t,v)
xlim([0 tmax])
xlabel('$t$/s','interpreter','latex')
ylabel('$v/$m s$^{-1}$','interpreter','latex')

% Check if <v> is approximately zero.
v_av = mean(v);
check_v_av = ['<v> = ',num2str(v_av),'.'];
disp(check_v_av)

% Check if <v^2> is approximately kB*T/m.
% Note: v.^2 squares each element of v
v2av = mean(v.^2);
compare_vals = ['<v^2> is ',num2str(v2av),' and kB T/m is ',num2str(kB*T/m),'.'];
disp(compare_vals)
