I'm trying to write a piece of code that calculates a diffraction pattern similar to an X-ray experiment using a FFT.
From my knowledge, the diffraction pattern for point particles can be calculated from the following formula:
$$
S(\mathbf{k}) = \left| \int_V d\mathbf{r}\rho(\mathbf{r})e^{-i\mathbf{k}\cdot\mathbf{r}} \right|^2
$$
which is basically the structure factor and $\rho(\mathbf{r})$ is the density function.
My shot on this was to first discretise the above formula using a box of side $r_{max}$:
$$
S(\mathbf{k}) = \frac{r_{max}}{M} \sum_{lmn}^M \rho_{lmn} e^{-i\frac{r_{max}}{M}\left( k_x l + k_y m +k_z n \right)}
$$
Now, I'm not sure how to discretise the wavevector and how this discretization will correspond to the lmn indices of the FFT:
$$
A_{lmn} = \sum_{ijk}^M a_{ijk} e^{-i2\pi\frac{1}{M}\left( l i + mj +nk \right)}
$$
An other detail I'm unsure of is how to get a 2D image corresponding to the diffraction pattern. What I though of, was that $\mathbf{k} = \mathbf{k}_S – \mathbf{k}_I$, where S is for the scattered light and I for the incident and I can set $\mathbf{k}_I = (\frac{2\pi}{\lambda},0,0)$. In approximation, the scattered and incident waves have the same magnitude $\frac{2\pi}{\lambda}$ so for a pixel of the screen (which has a distance d from the sample) with coordinates $x,y$, I have:
$$
\mathbf{k}_S = \frac{2\pi}{\lambda} \frac{(d, x, y)}{\sqrt{d^2 + x^2 + y^2}}
$$
and consequently I can calculate a scattering wave-vector for each pixel.
As I mentioned I'm not so sure how to discretise my scattering wave-vector and how this corresponds to the FFT wave-vectors, and on the last step of finding the correct wave-vectors for each pixel. Could you perhaps help me out?
Best Answer
If you really have point particles, then $\rho(\mathbf{r})=\sum_i \rho_i \delta(\mathbf{r}-\mathbf{r}_i)$, and hence $S(\mathbf{k}) = \left| \sum_i\rho_i e^{-i\mathbf{k}\cdot\mathbf{r_i}} \right|^2$. Even if the $\mathbf{r}_i$ are not on a uniform grid, there exists non-uniform FFT techniques to evaluate such sums for a list of $\mathbf{k}_j$ using normal FFTs with appropriate oversampling + low-pass filters and interpolation kernels.
As you seem to assume that the $\mathbf{r}_i$ are on a uniform grid, you only need the interpolation kernels (and appropriate oversampling). The simplest proceeding is probably to compute the higher order derivatives by a FFT on a uniform grid, and use these derivatives to construct local polynomial interpolants which allow the evaluation of the spectrum at arbitrary points $\mathbf{k}_j$.
If asked like this, you just need to compare the coefficients of $l$, $m$ and $n$ in the exponent: The $k_x$, $k_y$ and $k_z$ must be such that $e^{-i\frac{r_{max}}{M}\left( k_x l + k_y m +k_z n \right)} = e^{-i2\pi\frac{1}{M}\left( l i + mj +nk \right)}$. This is equivalent to $r_{max}\left( k_x l + k_y m +k_z n \right) = 2\pi\left( l i + mj + nk + M\mathbb{Z} \right)$. If this should be true for arbitrary $l$, $m$ and $n$ ($\in \mathbb{Z}$) then $k_x = \frac{2\pi}{r_{max}}(i + M\mathbb{Z})$, $k_y = \frac{2\pi}{r_{max}}(j + M\mathbb{Z})$ and $k_z = \frac{2\pi}{r_{max}}(k + M\mathbb{Z})$.
I'm not sure, but I believe you need to select a diffraction theory like Rayleigh-Sommerfeld or Fresnel-Kirchhoff first. This will give you a decomposition of the diffracted field into a sum of plane waves. (Alternatively you can also solve the diffraction problem rigorously by Maxwell's equation to get this sum of plane waves.) You then must propagate these plane waves to the observation point (take care to treat evanescent plane waves properly...), and superpose them again (for example by another FFT).