Bradley Monk at 09:00, 23 June 2024

2024-06-23T09:00:12Z

← Older revision		Revision as of 02:00, 23 June 2024
Line 1:		Line 1:
	~~==Brownian Motion MATLAB Code==~~

	<syntaxhighlight lang="matlab" line start="1">		<syntaxhighlight lang="matlab" line start="1">
	function [tracks] = BrownianMotion()		function [tracks] = BrownianMotion()

Bradley Monk: Created page with "==Brownian Motion MATLAB Code== function [tracks] = BrownianMotion() % STARTING PARAMETERS %----------------------------------------------------------- D = .3; % Diffusion Rate Ds = .1; % Diffusion Scalar (Ds = Dn) Dn = D/(D/Ds); % new D after scaling L d = 2; % dimensions dT = 1; % time step k = sqrt(d*D); % stdev of..."

2024-06-23T08:59:58Z

Created page with "==Brownian Motion MATLAB Code== <syntaxhighlight lang="matlab" line start="1"> function [tracks] = BrownianMotion() % STARTING PARAMETERS %----------------------------------------------------------- D = .3; % Diffusion Rate Ds = .1; % Diffusion Scalar (Ds = Dn) Dn = D/(D/Ds); % new D after scaling L d = 2; % dimensions dT = 1; % time step k = sqrt(d*D); % stdev of..."

New page

==Brownian Motion MATLAB Code==

<syntaxhighlight lang="matlab" line start="1">
function [tracks] = BrownianMotion()

% STARTING PARAMETERS
%-----------------------------------------------------------
D = .3; % Diffusion Rate
Ds = .1; % Diffusion Scalar (Ds = Dn)
Dn = D/(D/Ds); % new D after scaling L
d = 2; % dimensions
dT = 1; % time step
k = sqrt(d*D); % stdev of D's step size distribution
MSD = 2*d*D; % mean squared displacement
L = sqrt(2*d*D); % average diagonal (2D) step size
Lx = L/sqrt(2); % average linear (1D) step size
Ls = 1/sqrt(D/Ds); % scales Lx values for Dn

MSDtest = [1 0 0]; % test: D, Dn, or L
Scale = 1/10; % scale of model
Ndots = 100;
Nsteps = Ndots;

xyl = ones(2,Ndots);
xyds = ones(2,Ndots);
lims = ((D+1)^2)*10;

% LIVE PARTICLE DIFFUSION
%-----------------------------------------------------------
for t = 1:Nsteps

xyds = STEPxyds(Ndots, k);
[xyl] = AMPARSTEP(Ndots, xyds, xyl);
MAINPLOT(xyl, lims);

end

% MSD RANDOM STEPS ANALYSIS
%-----------------------------------------------------------
tracks = cell(Ndots, 1);

stepN = 1;
for t = 1:Nsteps

xyds = STEPxyds(Ndots, k);
[xyl] = AMPARSTEP(Ndots, xyds, xyl);
[tracks] = MSDfun(stepN, Nsteps, tracks, xyds);

stepN = stepN+1;
end
MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);

% MSD UNIFORM STEPS ANALYSIS
%-----------------------------------------------------------
stepN = 1;
for t = 1:Nsteps

xyds = stepsize(Ndots, Lx);

[xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls);
[tracks] = MSDfun(stepN, Nsteps, tracks, xyds);

stepN = stepN+1;
end
MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);

end

%%

% STEP SIZE GENERATOR
%-----------------------------------------------------------
function xyds = STEPxyds(Ndots, k)

xyds = (k * randn(2,Ndots));

end

% MOVE PARTICLES MAIN FUNCTION
%-----------------------------------------------------------
function [xyl] = AMPARSTEP(Ndots, xyds, xyl)

for j = 1:Ndots
xyl(:,j) = xyl(:,j)+xyds(:,j);
end

end

% LIVE DIFFUSION PLOT
%-----------------------------------------------------------
function [] = MAINPLOT(xyl, lims)

xlim = [-lims lims];
ylim = [-lims lims];
zlim = [-5 5];

figure(1)
subplot(2,1,1),
AMPARPlot = gscatter(xyl(1,:),xyl(2,:));
axis([xlim, ylim]);
set(AMPARPlot,'marker','.','markersize',[6],'color',[1 0 0])

figure(1);
subplot(2,1,2),
gscatter(xyl(1,:),xyl(2,:)); view(20, 30);
axis normal;
grid off
axis([xlim, ylim, zlim]);
set(gca, 'Box', 'on');

end

% MANUAL STEP SIZE FUNCTION
%-----------------------------------------------------------
function xyds = stepsize(Ndots, Lx)

Lx(1:2,1:Ndots) = Lx;
xyd = randi([0 1],Ndots,2)';
xyd(xyd == 0) = -1;
xyds = (Lx.*xyd);

end

% MSD SCALED STEPS FUNCTION
%-----------------------------------------------------------
function [xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls)

for j = 1:Ndots
xyds(:,j) = xyds(:,j)*Ls;
xyl(:,j) = xyl(:,j)+xyds(:,j);
end

end

% MSD TRACKS GENERATOR
%-----------------------------------------------------------
function [tracks] = MSDfun(stepN, Nsteps, tracks, xyds)

time = (0:Nsteps-1)';
xymsd = xyds';
xymsd = cumsum(xymsd,1);
tracks{stepN} = [time xymsd];

end

% MSD TRACKS ANALYSIS
%-----------------------------------------------------------
function [] = MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest)

SPACE_UNITS = 'µm';
TIME_UNITS = 's';
N_PARTICLES = Ndots;
N_TIME_STEPS = Nsteps;
N_DIM = 2;

oD = D; % raw µm^2/s
D = D*Scale; % to-scale µm^2/s

oDn = Dn; % raw µm^2/s
Dn = Dn*Scale; % to-scale µm^2/s

oL = L; % raw µm
L = L*Scale; % to-scale µm

dTbase = dT; % raw time-step
dT = dT*Scale; % to-scale time-step
k = k; % stdv of step distribution

ma = msdanalyzer(2, SPACE_UNITS, TIME_UNITS);
ma = ma.addAll(tracks);
disp(ma)

figure
ma.plotTracks;
ma.labelPlotTracks;

ma = ma.computeMSD;
ma.msd;

t = (0 : N_TIME_STEPS)' * dT;
[T1, T2] = meshgrid(t, t);
all_delays = unique( abs(T1 - T2) );

figure
ma.plotMSD;

cla
ma.plotMeanMSD(gca, true)

mmsd = ma.getMeanMSD;
t = mmsd(:,1);
x = mmsd(:,2);
dx = mmsd(:,3) ./ sqrt(mmsd(:,4));
errorbar(t, x, dx, 'k')

[fo, gof] = ma.fitMeanMSD;
plot(fo)
ma.labelPlotMSD;
legend off

ma = ma.fitMSD;

good_enough_fit = ma.lfit.r2fit > 0.8;
Dmean = mean( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;
Dstd = std( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;

Dheader1 = ['Raw Unscaled Values'];
Dhead1 = [' D Dn L'];
Ddat1 = [oD oDn oL];
disp(' ')
disp(Dheader1)
disp(Dhead1)
disp(Ddat1)

yourtesthead = ['YOU ARE TESTING DIFFUSION FOR:'];
if MSDtest(1)
yourtest = [' D: original diffusion rate'];
elseif MSDtest(2)
yourtest = [' Dn: new diffusion rate'];
elseif MSDtest(3)
yourtest = [' L: step length'];
else
yourtest = [' generic diffusion rate'];
end
disp(yourtesthead)
disp(yourtest)

disp(' ')
fprintf('Estimation of raw D coefficient from MSD:\n')
fprintf('D = %.3g ± %.3g (mean ± std, N = %d)\n', ...
Dmean, Dstd, sum(good_enough_fit));

% Retrieve instantaneous velocities, per track
trackV = ma.getVelocities;

% Pool track data together
TV = vertcat( trackV{:} );

% Velocities are returned in a N x (nDim+1) array: [ T Vx Vy ...]. So the
% velocity vector in 2D is:
V = TV(:, 2:3);

% Compute diffusion coefficient
varV = var(V);
mVarV = mean(varV); % Take the mean of the two estimates
Dest = mVarV / 2 * dT;

Dheader2 = ['Scaling to model...'];
Dhead2 = [' D Dn L'];
Ddat2 = [D Dn L];

disp(' ')
disp(Dheader2)
disp(Dhead2)
disp(Ddat2)
fprintf('Estimation from velocities histogram:\n')
fprintf('Tested D = %.3g %s, compare to scaled Des value of %.3g %s\n', ...
Dest, [SPACE_UNITS '²/' TIME_UNITS], D, [SPACE_UNITS '²/' TIME_UNITS]);

% printf('D.psd target value was %.3g %s\n', ...
% Dest, msdDpsd, [SPACE_UNITS '²/' TIME_UNITS]);

end

</syntaxhighlight>

<br /><br />
{{Clear}}

Brownian Motion MATLAB Code - Revision history

Bradley Monk at 09:00, 23 June 2024