# 2D Fourier transform of a few (4) disjoint discs on a plane

I'd really appreciate some advice.

Short Version

I'm trying to calculate the following $$\psi(X,Y,z=d)=\underset{aperature}{\iint}\psi(x,y,z=0)e^{-i\frac{k}{d}\left(xX+yY\right)}dxdy\\\psi(x,y,z=0)={\sum}\left[\delta\left(x-x_{i}\right)\delta\left(y-y_{i}\right)e^{ih_{i}k}\right]\otimes circ(\sqrt{x_{i}^{2}+y_{i}^{2}})$$

Long version

For an optics project, i'm trying to use Mathematica to perform a 2D Fourier Transform of 4 discs where each disc has a complex exponent representing its constant, relative phase (discs = functions who are 1*e^(i*phase) if inside the disc and 0 else). Disc centers are at different points in the plane, and the discs are disjoint.

I'm very new to Mathematica and not finding a better way around, have been trying to use NIntegrate to numerically asses and plot (using Plot3D) the FT of the discs, having expressed them in Cartesian coordinates.

Besides being a nightmare to write, I can't get the code to finish running since the integrals are very oscillatory and Mathematica keep sending me these notifications:

"Numerical integration converging too slowly; suspect one of the \ following: singularity, value of the integration is 0, highly \ oscillatory integrand, or WorkingPrecision too small."

I know there are no singularities since when I calculate each integral separately from the rest with smaller phase (so it won't oscillate allot) I do manage to get a correct 3D plot. However, when all of the integrals are calculated together, the calculation takes forever (I quit the kernel after about half an hour of running).

It'd really help if someone could advise me on how to write this code or suggest a simpler way to do this (there must be one). As I mentioned, I'm very new to Mathematica and so thorough explanations would help allot and be very much appreciated.

I'm attaching my code so far for reference.

My code:

\*discs radius*\
ρ = 1;
\*constants that determine phase*\
k = (2 π)/(0.5 10^-6);
d = 2;
Subscript[h, 1] = 0.001;
Subscript[h, 2] = 0.002;
Subscript[h, 3] = 0.003;
Subscript[x, 1] = 0;
Subscript[y, 1] = 0.03;
Subscript[x, 2] = 0.035;
Subscript[y, 2] = 0.01;
Subscript[x, 3] = 0.025;
Subscript[y, 3] = 0.02;

Plot3D[(Abs[ -((I k )/d) π ρ^2 ((
2 BesselJ[1, k/d  ρ Sqrt[u^2 + v^2]])/(
k/d  ρ Sqrt[u^2 + v^2])) - (2 I)/
v  ( E^(I k/d (Subscript[h, 1] - v Subscript[y, 1])) (
NIntegrate[
E^(-I  u x k/
d) (E^(-I v k/
d  (ρ^2 - ( x - Subscript[x, 1])^2)^0.5) - 1), {x,
Subscript[x, 1], Subscript[x, 1] + ρ},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (E^(-I v k/
d  (ρ^2 - ( x - Subscript[x, 1])^2)^0.5) - 1), {x,
Subscript[x, 1] - ρ, Subscript[x, 1]},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (1 - E^(
I v k/d (ρ^2 - ( x - Subscript[x, 1])^2)^0.5)), {x,
Subscript[x, 1] -ρ, Subscript[x, 1]},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (1 - E^(
I v k/d (ρ^2 - ( x - Subscript[x, 1])^2)^0.5)), {x,
Subscript[x, 1], Subscript[x, 1] + ρ},
MaxRecursion -> 5000]) +
E^(I k/d (Subscript[h, 2] - v Subscript[y, 2])) (
NIntegrate[
E^(-I  u x k/
d) (E^(-I v k/
d  (ρ^2 - ( x - Subscript[x, 2])^2)^0.5) - 1), {x,
Subscript[x, 2], Subscript[x, 2] + ρ},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (E^(-I v k/
d  (ρ^2 - ( x - Subscript[x, 2])^2)^0.5) - 1), {x,
Subscript[x, 2] - ρ, Subscript[x, 2]},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (1 - E^(
I v k/d (ρ^2 - ( x - Subscript[x, 2])^2)^0.5)), {x,
Subscript[x, 2] -ρ, Subscript[x, 2]},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (1 - E^(
I v k/d (ρ^2 - ( x - Subscript[x, 2])^2)^0.5)), {x,
Subscript[x, 2], Subscript[x, 2] + ρ},
MaxRecursion -> 5000]) +
E^(I k/d (Subscript[h, 3] - v Subscript[y, 3])) (
NIntegrate[
E^(-I  u x k/
d) (E^(-I v k/
d  (ρ^2 - ( x - Subscript[x, 3])^2)^0.5) - 1), {x,
Subscript[x, 3], Subscript[x, 3] + ρ},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (E^(-I v k/
d  (ρ^2 - ( x - Subscript[x, 3])^2)^0.5) - 1), {x,
Subscript[x, 3] - ρ, Subscript[x, 3]},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (1 - E^(
I v k/d (ρ^2 - ( x - Subscript[x, 3])^2)^0.5)), {x,
Subscript[x, 3] -ρ, Subscript[x, 3]},
MaxRecursion -> 5000] +
NIntegrate[
E^(-I  u x k/
d) (1 - E^(
I v k/d (ρ^2 - ( x - Subscript[x, 3])^2)^0.5)), {x,
Subscript[x, 3], Subscript[x, 3] + ρ},
MaxRecursion -> 5000]))])^2, {u, -0.001,
0.001}, {v, -0.001, 0.001}]

• A couple points: first, using NIntegrate to calculate Fourier transforms is not a good idea, as it's very inefficient (a naive Riemann integration-based 2D Fourier transform is an $O(n^4)$ operation, where $n$ is the number of subdivisions per axis). It would be far more efficient to use 2D fast Fourier transforms, which are implemented using Fourier. For that, you will have to carefully look at the definition of Fourier to determine what your axes are showing you, but it will execute in milliseconds, rather than hours. Oct 13 '14 at 20:15
• Second, unless I'm missing something, the frequency domain function is available in closed form, since the Fourier transform of a disk is a first-order Bessel function $J_1$, since it's the Hankel transform of a square pulse. Translations simply amount to multiplication by plane waves, so you have a simple closed-form solution which is a finite sum of product of Bessel functions and plane waves. See the Wikipedia article on Fourier transforms to get the exact expression. As such, a numerical computation is unnecessary. Oct 13 '14 at 20:19

Here's a qualitative way to do the computation using FFTs. First, make some data (in this, the disks all have phase 1, but that can be easily fixed):

w1 = 600;
w2 = 800;
dat = Sum[
RotateRight[
DiskMatrix[
RandomInteger[{1, 150}], {w1,
w2}], {RandomInteger[{-1000, 1000}],
RandomInteger[{-1000, 1000}]}], {k, 6}];


Here's a plot of the data:

Image[0.5 dat, Magnification -> 1]


Then you can plot the frequency domain picture like this:

<< Developer
Image[RotateRight[
Abs[Fourier[dat]], {w1/2, w2/2}]\[TensorProduct]ToPackedArray[{1.0,
0.3, 0.1}], Magnification -> 1]


• thank you very much for the detailed answer! As I said, I kinda new to this and was not yet familiar with the in and outs of working with Mathematica. I'll go over your answer and will follow up with my results. Oct 14 '14 at 7:48
• I tried using the FT property for translations that you mentioned with this: k = (2 [Pi])/(0.5 10^-6); [Rho] = 1; d = 2; ContourPlot[(Re[ I [Pi] [Rho]^2 (( 2 BesselJ[1, [Rho] Sqrt[u^2 + v^2]])/( [Rho] Sqrt[ u^2 + v^2])) (1 + E^(-I (1 u + 20 v)) + E^(-I (20 u + 1 v)) + E^(-I (50 u + 10 v)))])^2, {u, -20, 20}, {v, -20, 20}] I didn't get a correct plot and don't really know why. Also, thank you for the detailed code you posted. I managed to use it to get the aperture I need but can't find a way to add phases to each DiskMatrix. Oct 14 '14 at 9:06
• @AsafMiron: To add phases to the disks, instead of using DiskMatrix[...] you can use Exp[I phi]DiskMatrix[...] where phi` is the phase you want. As for the analytical form, I don't remember all the Fourier transform details, but your professor can probably help you sort out the paper-based math. Oct 14 '14 at 13:03
• Thanks! One last thing, how do I choose a vector in some direction in the image? (I want to FT along lines in the image). Tried this to no avail: Select[mask, mask[{#}, {-(1/2) # + 450}]]] Oct 14 '14 at 13:11
• @AsafMiron: What does "FT along lines" mean? A 2D FT is a map $f:\mathbb{C}^2\rightarrow\mathbb{C}^2$, so I'm not sure what you mean by Fourier transform along a line. Do you mean that you want to compute the Fourier integral for one particular choice of wavevector? Oct 14 '14 at 13:18