Particle tracking using IDL -- John C. Crocker and
Eric R. Weeks
Home |
Download software |
Tutorial |
Extra software
msd.pro
This program calculates the mean square displacement of data in a
track array:
IDL> m = msd(t)
IDL> plot,m(0,*),m(3,*),/xlog,/ylog; this is <x^2>
IDL> oplot,m(0,*),m(4,*),linestyle=1; this is <y^2>
The output of msd has several columns, they are:
- m(0,*) = lag time delta-t
- m(1,*) = < delta x >
- m(2,*) = < delta y >
- m(3,*) = < (delta x)^2 >
- m(4,*) = < (delta y)^2 >
- m(5,*) = < (delta r)^2 > = m(3,*) + m(4,*)
- m(6,*) = N
The angle brackets indicate averages over all particles and
all starting times. The last column, N, is an estimate of the
independent number of data points in the average. This is
a little tricky to explain. Start with N' being the actual
number of dx's that are computed. For dt = 1 frame, N' = N.
For larger dt's, N = 2 N' / dt. The rationale is that for
large dt's, there's a lot of redundant data, so less truly
independent measurements (see below). Anyway, if you want to get back N',
just do something like:
IDL> m = msd(t)
IDL> m2 = m
IDL> m2(6,*) = round(m2(6,*)*m2(0,*)/2.0)
IDL> m2(6,0) = m(6,0)
There are several keywords you can use:
- DIM = 2; sets the dimension of the data, default is 2. If
you set it to DIM=3, then the columns of the array are:
- m(0,*) = lag time delta-t
- m(1:3,*) = < delta x >, < delta y >, < delta z >
- m(4:6,*) = < (delta x)^2 >, < (delta y)^2 >, <
(delta z)^2 >
- m(7,*) = < (delta r)^2 > = m(4,*) + m(5,*) + m(6,*)
- m(8,*) = N
- MICPERPIX = 1.0; conversion factor, to convert from
pixels to microns.
- TIMESTEP = 1.0; conversion factor, to convert from frame
number to seconds (or minutes, or whatever). In other words, it
rescales m(0,*).
- MAXTIME = 100.0; maximum lag time to calculate to. The default
is to not compute any lag times for which there
is insufficient data. The program judges this somewhat
conservatively, which means you can trust what it gives you.
However, you may have data that would cover longer times and
want to force the program to be less conservative, so you can use
this keyword. The MAXTIME should be specified in dimensional
time, in other words, if you're using TIMESTEP then this should
be the maximum time in seconds, and otherwise, in frames. (What
I mean by data at longer times: for example, you might have lots
of trajectories that are 100 s in duration, and a few
trajectories that last for 1000 s. The program will probably
stop after a lag time of about 100 s. You can set MAXTIME=1000.0 to
force it to calculate the mean square displacement up to lag
times of 1000 s, despite the fact that you have very little data
to determine the MSD at such long lag times. Use with caution!)
- OUTFILE = "msd.gdf"; you can use this to specify a file to
save the output to.
- MYDTS = [array of delta-t]; the program's default mode is to
generate a list of logrithmically spaced delta-t's. You can
provide specific ones instead, if you'd like.
- /QUIET; turns off printing of messages
- /KEEPDRIFT; does not subtract off <x >^2 from <x^2
> (likewise
for y and z).
Calculating the mean square displacement
A very brief tutorial. Define:
dx(t,tau) = x(t+tau) - x(t)
Then, define
where the mean is taken over all initial times t, and all
particles. This is the mean displacement.
Likewise,
< dx^2 > = mean(dx(t,tau)^2)
is the mean square displacement. You can perhaps see why we
define N the way we do, above. For example, imagine you have a
trajectory that is 1000 time steps long, and you are calculating
the mean square displacement for dt = 998. This is:
((x[998] - x[0])^2 + (x[999] - x[1])^2) / 2
The average is taken over the initial times (t=0,1) and thus
there are two dx's that we're averaging. However, these two dx's
are nearly the same, as x[0] and x[1] are approximately the same,
and likewise x[998] and x[999]. Thus, there's approximately one
independent measurement being used, and we define N = (2 * 2) /
999 = 0.004 as a very small number indicating that this is a
rather poor average. Currently the default is that MSD is
calculated as long as N > 1000; this conservatively stops the MSD
calculation a little early. To override this, use the MAXTIME
keyword as described above, or use the secret keyword which you
learn because you've read all the way to the bottom of this
webpage:
- minN = 1000; set the minimum allowed value of N
One final note: For a dilute sample of colloidal particles, you
should find:
< dx^2 > = < dy^2 > = 2 D dt
That is, the MSD grows linearly with the lag time, with a
proportionality constant being related to the diffusion constant
D.
This program was written by John, and improved by
Victor
Breedveld, who fixed a DOS-related bug and made some other
minor improvements. This web page was written by Eric.
Click here to return to microrheology
summary page
Contact us