Particle tracking using IDL -- John C. Crocker and Eric R. Weeks
Home | Download software | Tutorial | Extra software

How to calculate the pair correlation function g(r)

This explanation is for three-dimensional data. To calculate g(r), do the following:

One caveat: For experimental data, you often have edges to your sample that are artificial. For example, you take a picture of particles but at the edges of your picture, the system extends further outwards. Thus when calculating g(r) based on reference particles near the edge of your image, you have a problem. You'll have to modify step #3 above with the correct volume/area that actually lies within the image you're looking at. The routines I've written take care of that.

I wrote IDL routines to calculate g(r) in 2D and 3D. These do some special tricks to deal with experimental data, in other words, to cleverly deal with the edge effects.

2D program: For a particle near the edge of a rectangular image, when I'm counting particles a distance r away from it, for many values of r the circle of radius r extends outside of the image. I did some math to figure out how to determine how much angular extent of the circle lies within the image for these cases (in the subroutine checkquadrant). Thus the edges are correctly accounted for.

3D program: For a particle near the edge of a rectangular image box, when I'm counting particles a distance r away from it, we have the same problem described for 2D data. However, calculating the resulting solid angle of the sphere contained within the box is more than I could handle mathematically, for arbitrary box dimensions and arbitrary locations within the box. So I did a trick. My image boxes tend to be short and wide, that is, very narrow in Z and large in X and Y.

My program does the following. I'm calculating g(r) for r < rmax, with a default rmax = 10. I consider only reference particles that are more than rmax away from the horizontal edges of the box, so I never have to worry about overlaps in X and Y. The resulting formula, to worry about overlaps in Z alone, turns out to be quite reasonable. (That is, I calculate the surface area of a spherical hemisphere of radius R that is cut off at a finite height H < R. It turns out this formula is simple: 2 pi R H.)

So, an important warning: Don't choose rmax to be more than half of the width of your data, in X and Y. This restriction is only for the 3D program.

See here for additional comments on my "special tricks".

  1. What is the pair correlation function?
  2. How to calculate g(r) (you are here)
  3. IDL routines to calculate g(r)
  4. Extra g(r) routines -- unsupported

Contact us