Qual: Difference between revisions
Jump to navigation
Jump to search
Bradley Monk (talk | contribs) No edit summary |
Bradley Monk (talk | contribs) No edit summary |
||
| Line 1: | Line 1: | ||
{{Box|font= | {{Box|font=150%|width=60%|float=left|text=14px|boarder=solid #aaa 1px|Qual Project Source Code|All source code used in my qualifying exam project can be found [https://github.com/subroutines/djhbrm here in my github repo for the project]]}} | ||
{{Clear}} | |||
==SI Fig 1== | |||
== | <html><iframe width="560" height="315" src="//www.youtube.com/embed/JH-hGjzhEFQ" frameborder="0" allowfullscreen></iframe></html> | ||
(hover→ {{Popup|{{#widget:Html5media | |||
|url=http://youtu.be/JH-hGjzhEFQ | |||
|width=200 | |||
|height=200 | |||
}}}} ). Instead, '''D''' is often calculated from the ''mean squared diffusion'' (MSD) path of the particle, defined below. | |||
Revision as of 11:10, 5 January 2015
Qual Project Source Code
SI Fig 1
<html><iframe width="560" height="315" src="//www.youtube.com/embed/JH-hGjzhEFQ" frameborder="0" allowfullscreen></iframe></html>
(hover→ {{#info: {{#widget:Html5media
|url=http://youtu.be/JH-hGjzhEFQ
|width=200
|height=200
}} }} ). Instead, D is often calculated from the mean squared diffusion (MSD) path of the particle, defined below.
Brownian Motion Matlab Code
function [boom] = BradsBrownianMotionScript()
format compact;
format short;
close all;
% Note: This matlab script requires 'msdanalyzer' toolbox
%-------------############################------------------%
% 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);
boom = 'done';
%-------------############################------------------%
end % ## END MAIN FUNCTION ##
%-------------############################------------------%
%%
%%
%-------------############################------------------%
% ## SUBFUNCTIONS ##
%-------------############################------------------%
%-------------------------------------------%
% 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];
%=================================%
% MAIN 2D PLOT
%---------------------------------%
figure(1)
subplot(2,1,1),
AMPARPlot = gscatter(xyl(1,:),xyl(2,:));
axis([xlim, ylim]);
set(AMPARPlot,'marker','.','markersize',[6],'color',[1 0 0])
%=================================%
% 3D PLOT
%---------------------------------%
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