Calculating diffraction patterns using FFT 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?
 A: 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$.

Now, I'm not sure how to discretise the wavevector and how this discretization will correspond to the lmn indices of the FFT: 

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})$.

An other detail I'm unsure of is how to get a 2D image corresponding to the diffraction pattern.

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).
A: Not sure if this question is still open or if my answer will exactly fit. I have developed a calculator/simulator for fringe patterns associated with straight edge, slit and multiple slit experiments. You will find Fringe pattern calculators on page 20 of my paper.
https://singleedgecertainty.files.wordpress.com/2015/03/pdf-single-edge-certainty.pdf
This is an Excel based program and will calculate fringe patterns in nano-meters and simulate in graph form what the fringe patterns will look like. The program only represents the right side of the pattern starting from the center and going as far out as the program is set for. The shaded and highlighted cells can be manipulated to set any of the distances in nano-meters. The chart will correspond with a simulated version of an actual fringe pattern showing the actual distances and spacing's again in nano-meters.
Here is a sample of my two slit calculator.
https://www.dropbox.com/s/7znc6u8vzw7iju1/Alsept%20%202%20Slit%20Calculator.xlsx 
After opening click download, direct download or save then open Excel program. You will need to enable editing.
The background of the program shows all the calculations and is where you inter the variables such as slit width, separation, distance from slits to detector etc.
The chart/graph shows the pattern with fringe spacing's. You can zoom in and out on this to see the whole effect. 
Please let me know if you have any questions or would just like to discuss. thanks  
