Particle tracking using IDL -- John C. Crocker and Eric R. Weeks

# 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

< dx > = mean(dx(t,tau))

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.