# Qual

Qual Paper SI Material

Qual Project Source Code

All source code used in my qualifying exam project can be found here in my github repo for the project. An online version of my qual with enhanced features like citation data popup graphics and animations can be found on this page

## SI Fig. 2 - Brownian Motion Matlab Code

Brownian Motion Matlab Code

```  1 function [boom] = BradsBrownianMotionScript()
2 format compact;
3 format short;
4 close all;
5
6 % Note: This matlab script requires 'msdanalyzer' toolbox
7 %-------------############################------------------%
8 %                 STARTING PARAMETERS
9 %-------------############################------------------%
10 D = .3;                     % Diffusion Rate
11 Ds = .1;                    % Diffusion Scalar (Ds = Dn)
12 Dn = D/(D/Ds);              % new D after scaling L
13 d = 2;                      % dimensions
14 dT = 1;                     % time step
15 k = sqrt(d*D);	            % stdev of D's step size distribution
16 MSD = 2*d*D;                % mean squared displacement
17 L = sqrt(2*d*D);            % average diagonal (2D) step size
18 Lx = L/sqrt(2);             % average linear (1D) step size
19 Ls = 1/sqrt(D/Ds);          % scales Lx values for Dn
20
21
22 MSDtest = [1 0 0];			% test: D, Dn, or L
23 Scale = 1/10;				% scale of model
24 Ndots = 100;
25 Nsteps = Ndots;
26
27
28 xyl = ones(2,Ndots);
29 xyds = ones(2,Ndots);
30 lims = ((D+1)^2)*10;
31
32
33 %===========================================================%
34 %               LIVE PARTICLE DIFFUSION
35 %-----------------------------------------------------------%
36 for t = 1:Nsteps
37
38     xyds = STEPxyds(Ndots, k);
39 	[xyl] = AMPARSTEP(Ndots, xyds, xyl);
40 	MAINPLOT(xyl, lims);
41
42 end
43 %===========================================================%
44
45
46 %===========================================================%
47 %               MSD RANDOM STEPS ANALYSIS
48 %-----------------------------------------------------------%
49 tracks = cell(Ndots, 1);
50
51 stepN = 1;
52 for t = 1:Nsteps
53
54 	xyds = STEPxyds(Ndots, k);
55 	[xyl] = AMPARSTEP(Ndots, xyds, xyl);
56     [tracks] = MSDfun(stepN, Nsteps, tracks, xyds);
57
58 stepN = stepN+1;
59 end
60 MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);
61
62
63
64
65
66 %===========================================================%
67 %               MSD UNIFORM STEPS ANALYSIS
68 %-----------------------------------------------------------%
69 stepN = 1;
70 for t = 1:Nsteps
71
72 xyds = stepsize(Ndots, Lx);
73
74 	[xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls);
75     [tracks] = MSDfun(stepN, Nsteps, tracks, xyds);
76
77 stepN = stepN+1;
78 end
79 MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest);
80
81
82 boom = 'done';
83 %-------------############################------------------%
84 end %         ##    END MAIN FUNCTION   ##
85 %-------------############################------------------%
86 %%
87
88
89
90
91
92 %%
93 %-------------############################------------------%
94 %             ##      SUBFUNCTIONS      ##
95 %-------------############################------------------%
96
97
98 %-------------------------------------------%
99 % STEP SIZE GENERATOR
100 %-------------------------------------------%
101 function xyds = STEPxyds(Ndots, k)
102
103     xyds = (k * randn(2,Ndots));
104
105 end
106
107
108 %-------------------------------------------%
109 % MOVE PARTICLES MAIN FUNCTION
110 %-------------------------------------------%
111 function [xyl] = AMPARSTEP(Ndots, xyds, xyl)
112
113 	for j = 1:Ndots
114         xyl(:,j) = xyl(:,j)+xyds(:,j);
115 	end
116
117 end
118
119
120 %-------------------------------------------%
121 % LIVE DIFFUSION PLOT
122 %-------------------------------------------%
123 function [] = MAINPLOT(xyl, lims)
124 %-------------------------------------------%
125 xlim = [-lims lims];
126 ylim = [-lims lims];
127 zlim = [-5 5];
128
129 %=================================%
130 %       MAIN 2D PLOT
131 %---------------------------------%
132 figure(1)
133 subplot(2,1,1),
134 AMPARPlot = gscatter(xyl(1,:),xyl(2,:));
135 axis([xlim, ylim]);
136 set(AMPARPlot,'marker','.','markersize',[6],'color',[1 0 0])
137
138
139 %=================================%
140 %           3D PLOT
141 %---------------------------------%
142 figure(1);
143 subplot(2,1,2),
144 gscatter(xyl(1,:),xyl(2,:)); view(20, 30);
145 axis normal;
146 grid off
147 axis([xlim, ylim, zlim]);
148 set(gca, 'Box', 'on');
149
150 end
151
152
153 %-------------------------------------------%
154 % MANUAL STEP SIZE FUNCTION
155 %-------------------------------------------%
156 function xyds = stepsize(Ndots, Lx)
157 %-------------------------------------------%
158
159    Lx(1:2,1:Ndots) = Lx;
160    xyd = randi([0 1],Ndots,2)';
161    xyd(xyd == 0) = -1;
162    xyds = (Lx.*xyd);
163
164 end
165
166
167 %-------------------------------------------%
168 % MSD SCALED STEPS FUNCTION
169 %-------------------------------------------%
170 function [xyl xyds] = MSDAMPARSTEP(Ndots, xyds, xyl, Ls)
171
172 	for j = 1:Ndots
173         xyds(:,j) = xyds(:,j)*Ls;
174         xyl(:,j) = xyl(:,j)+xyds(:,j);
175 	end
176
177 end
178
179
180 %-------------------------------------------%
181 % MSD TRACKS GENERATOR
182 %-------------------------------------------%
183 function [tracks] = MSDfun(stepN, Nsteps, tracks, xyds)
184     time = (0:Nsteps-1)';
185     xymsd = xyds';
186     xymsd = cumsum(xymsd,1);
187     tracks{stepN} = [time xymsd];
188
189 end
190
191
192 %-------------------------------------------%
193 % MSD TRACKS ANALYSIS
194 %-------------------------------------------%
195 function [] = MSDfunction(tracks,Ndots,Nsteps,D,Dn,L,dT,k,Scale,MSDtest)
196
197 SPACE_UNITS = 'µm';
198 TIME_UNITS = 's';
199 N_PARTICLES = Ndots;
200 N_TIME_STEPS = Nsteps;
201 N_DIM = 2;
202
203 oD = D;				% raw		µm^2/s
204 D  = D*Scale;       % to-scale	µm^2/s
205
206 oDn = Dn;			% raw		µm^2/s
207 Dn = Dn*Scale;		% to-scale	µm^2/s
208
209 oL = L;				% raw		µm
210 L = L*Scale;		% to-scale	µm
211
212 dTbase = dT;		% raw time-step
213 dT = dT*Scale;		% to-scale time-step
214 k = k;				% stdv of step distribution
215
216 ma = msdanalyzer(2, SPACE_UNITS, TIME_UNITS);
218 disp(ma)
219
220 figure
221 ma.plotTracks;
222 ma.labelPlotTracks;
223
224 ma = ma.computeMSD;
225 ma.msd;
226
227 t = (0 : N_TIME_STEPS)' * dT;
228 [T1, T2] = meshgrid(t, t);
229 all_delays = unique( abs(T1 - T2) );
230
231 figure
232 ma.plotMSD;
233
234
235 cla
236 ma.plotMeanMSD(gca, true)
237
238 mmsd = ma.getMeanMSD;
239 t = mmsd(:,1);
240 x = mmsd(:,2);
241 dx = mmsd(:,3) ./ sqrt(mmsd(:,4));
242 errorbar(t, x, dx, 'k')
243
244 [fo, gof] = ma.fitMeanMSD;
245 plot(fo)
246 ma.labelPlotMSD;
247 legend off
248
249
250 ma = ma.fitMSD;
251
252 good_enough_fit = ma.lfit.r2fit > 0.8;
253 Dmean = mean( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;
254 Dstd  =  std( ma.lfit.a(good_enough_fit) ) / 2 / ma.n_dim;
255
256 Dheader1 = ['Raw Unscaled Values'];
257 Dhead1 = ['    D        Dn        L'];
258 Ddat1 = [oD oDn oL];
259 disp(' ')
262 disp(Ddat1)
263
264
265 yourtesthead = ['YOU ARE TESTING DIFFUSION FOR:'];
266 if MSDtest(1)
267 	yourtest = ['   D:   original diffusion rate'];
268 elseif MSDtest(2)
269 	yourtest = ['   Dn:  new diffusion rate'];
270 elseif MSDtest(3)
271 	yourtest = ['   L:  step length'];
272 else
273 	yourtest = ['   generic diffusion rate'];
274 end
276 disp(yourtest)
277
278 disp(' ')
279 fprintf('Estimation of raw D coefficient from MSD:\n')
280 fprintf('D = %.3g ± %.3g (mean ± std, N = %d)\n', ...
281     Dmean, Dstd, sum(good_enough_fit));
282
283
284
285
286 % Retrieve instantaneous velocities, per track
287  trackV = ma.getVelocities;
288
289  % Pool track data together
290  TV = vertcat( trackV{:} );
291
292  % Velocities are returned in a N x (nDim+1) array: [ T Vx Vy ...]. So the
293  % velocity vector in 2D is:
294  V = TV(:, 2:3);
295
296  % Compute diffusion coefficient
297 varV = var(V);
298 mVarV = mean(varV); % Take the mean of the two estimates
299 Dest = mVarV / 2 * dT;
300
301
302
303 Dheader2 = ['Scaling to model...'];
304 Dhead2 = ['    D        Dn        L'];
305 Ddat2 = [D Dn L];
306
307 disp(' ')
310 disp(Ddat2)
311 fprintf('Estimation from velocities histogram:\n')
312 fprintf('Tested D = %.3g %s, compare to scaled Des value of %.3g %s\n', ...
313     Dest, [SPACE_UNITS '²/' TIME_UNITS], D, [SPACE_UNITS '²/' TIME_UNITS]);
314
315 % printf('D.psd target value was %.3g %s\n', ...
316 %     Dest, msdDpsd, [SPACE_UNITS '²/' TIME_UNITS]);
317
318 end
```

Text Text Text

Text Text Text

Text Text Text