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.