This applet demonstrates the use of Bruce Miller's jnt.FFT Fast Fourier Transform library, and thus by analogy, the FFT libraries from the GNU Scientific Libraries collection. The particular interest is how exactly to get 3-dimensional x-ray crystallographic data into and out of these functions.

There are two directions for the transform, forward and backward, which differ by only a sign -- the same sign difference between e-i2pi(hx+ky+lz) and e+i2pi(hx+ky+lz) in the standard crystallography equations. Just don't ask me which is which, because I'm not sure. Actually, I know that these libraries call the one with the minus sign the forward direction, but I'm not sure how this relates to crystallography or whether it even matters.

What I have figured out is that real space corresponds to the time domain, and reciprocal space corresponds to the frequency domain. This impacts the data storage format. I'm going to discuss the 1-D case first, where everything is stored in a simple array. Real space data (electron density) is a series of sample points along an axis of the unit cell, indexed as real numbers from 0 to 1. These are stored in order, as you would expect. For example:
    [r0.00, r0.01, r0.02, ... , r0.49, r0.50, r0.51, ..., r0.98, r0.99, r1.00]
where rx is the electron density (usually rho, here r) at position x (representing the fraction of the way across the unit cell).

The frequency domain is stored in "wrap-around" order, where the negative-indexed F's are stored after all the rest, and the most negative precedes the least negative. Thus:
    [F0, F+1, F+2, ... , F+(N-1), F+N, F-N, F-(N-1), ... , F-3, F-2, F-1]
where the F's are the crystallographic structure factors, which would usually have three indices instead of one: Fhkl. As you can see, it's as though the expected order (-N ... -1 0 +1 ... +N) has been "rotated" right by half the length of the array (round down). This brings up the point that my examples have had an odd number of entries. There is no change to the real space data for an even number of entries, but for reciprocal space, it apparently turns out that F+N = F-N, so the two center entries "collapse" to one entry in the array.

To get the sort of diffraction pattern you expect in crystallography, you have to rearrange the data so that zero is in the middle:
    [F-N, F-(N-1), ... , F-3, F-2, F-1, F0, F+1, F+2, ... , F+(N-1), F+N]
Again, this is just a "rotation" of the data.

One last complication. All Fourier transforms deal with complex numbers on both sides. If you want to pretend that the input is real-valued (i.e., that all the imaginary components are zero), as is sometimes the case for crystallography, some optimizations are possible, but in my opinion it's not worth the complications it introduces in the data format. However, we do have to have somewhere to store these complex numbers. Thus, each entry in the arrays above is really a pair of numbers, real component followed by imaginary component. Thus the arrays are allocated twice as long as the logical number of entries, because each entry uses two floating-point values.

Remember that these real and imaginary components are to structure factors and phases as Cartesian coordinates are to polar ones. Thus:
    F = sqrt(Re*Re + Im*Im)
    phi = atan2(Im, Re)
Or equivalently:
    Re = F cos(phi)
    Im = F sin(phi)

For multi-dimensional cases, all the values are stored in one array. If x is considered the "rows" index, y is "columns", and z is "slices" (for lack of a better term), then the array indexes of the real and imaginary components are as follows:
    Re(x, y, z) = data[ x*numColumns*numSlices + y*numSlices + z ]
    Im(x, y, z) = data[ x*numColumns*numSlices + y*numSlices + z + 1]

Because this code uses libraries that descend from a GNU product, this project is hereby released under the GNU Public License, as required.