Qual: Difference between revisions

Revision as of 12:05, 5 January 2015

Qual Project Source Code
All source code used in my qualifying exam project can be found here in my github repo for the project]
SI Fig 1

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
Qual: Difference between revisions

Revision as of 12:05, 5 January 2015

SI Fig 1

Brownian Motion Matlab Code

Navigation menu

Search
