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

# How to calculate the pair correlation function g(r)

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

• Pick a value of dr
• Loop over all values of r that you care about:
1. Consider each particle you have in turn. Count all particles that are a distance between r and r + dr away from the particle you're considering. You can think of this as all particles in a spherical shell surrounding the reference particle. The shell has a thickness dr.
2. Divide your total count by N, the number of reference particles you considered -- probably the total number of particles in your data.
3. Divide this number by 4 pi r^2 dr, the volume of the spherical shell (the surface area 4 pi r^2, multiplied by the small thickness dr). This accounts for the fact that as r gets larger, for trivial reasons you find more particles with the given separation.
4. Divide this by the particle number density. This ensures that g(r)=1 for data with no structure. In other words, if you just had an arbitrarily placed spherical shell of inner radius r and outer radius r+dr, you'd expect to find about rho * V particles inside, where rho is the number density and V is the volume of that shell.
• In 2D, follow the algorithm as above but divide by 2 pi r dr instead of step #3 above.

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.

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