# Help with Fourier beam propagation method

I am working on implementing the Fourier beam propagation method in C++. I am really more of a programmer than a physicist but I think I have a good understanding of what I am trying to do. Here is what I am trying to implement. Two questions. In the flow chart he says you should perform the diffraction and after that do the refraction. From the equation however it looks like you apply the refraction first, do the FFT and then apply the diffraction. Also in the formula for diffraction I notice that the has kz as a function of omega, kx, and another kx. Is that a typo where he should have ky instead of another kx?

I am including my code bellow. I am trying to create the diffraction pattern through a circular aperture. I am not really sure what is an appropriate size for the cell grids. Should it be on the order of a wavelength? I am just getting blank output and most in the inverse transform are NAN. I think the issue lies with how I calculate kz kx and ky. See comments in code. Any help would be appreciated in getting this working.

#include <complex>
#include <cmath>
#include <iostream>
#include <stddef.h>
#include "fourier.hpp"
#include <omp.h>
#include <fstream>
#include <fftw3.h>
#include <cstring>

const int N = 1000;

void save(const char* fn, std::complex<double>* g)
{
std::ofstream out(fn);
out << "P5 " << N << " " << N << " 255 ";

double max = 0.0;
for(int i = 0; i < N * N; i++)
{
max = fmax(std::abs(g[i]), max);
}

for(int i = 0; i < N * N; i++)
{
unsigned char c = 255.0 * std::abs(g[i]) / max;
out.write((const char*)&c, 1);
}

out.close();
}

int main()
{
typedef std::complex<double> complex;
const complex I(0.0, 1.0);

complex* e_prime = new std::complex<double>[N * N];
complex* e = new std::complex<double>[N * N];

memset(e, 0, sizeof(std::complex<double>) * N * N);

const double lambda = 500.0e-9;
const double dz = 1.0e-8; //forward step size
const double n = 1.0; //index of retraction (just a constant 1)
const double k0 = 2.0 * M_PI / lambda;
const double cell_size = lambda / 4.0; //I want the grid to have a 1/4 wavelength resolution

//create a circle in the electric field grid
for(int y = 0; y < N; y++)
{
for(int x = 0; x < N; x++)
{
int dx = x - N / 2;
int dy = y - N / 2;

if(dx * dx + dy * dy < 1000)
e[y * N + x] = 1.0;
else
e[y * N + x] = 0.0;

}
}

save("slit.ppm", e);

//apply the refraction step
for(int y = 0; y < N; y++)
{
for(int x = 0; x < N; x++)
{
e[y * N + x] *= std::exp(-I * n * dz * k0);
}
}

//Perform the fourier transform of the refracted electric field
fftw_plan p = fftw_plan_dft_2d(N, N, (fftw_complex*)e, (fftw_complex*)e_prime, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_NO_SIMD);
fftw_execute(p);
fftw_destroy_plan(p);

//save the fourier transform for reference
save("ft.ppm", e_prime);

//apply the diffraction step
for(int y = 0; y < N; y++)
{
for(int x = 0; x < N; x++)
{
/*
* This is where I am a bit confused.  So far I have not
* found where in the math to account for the size of the
* grid cells.  I figure since in the spatial domain you would
* multiply by the cell_size in the frequency domain I suppose you
* would divide.
*
* I think this may be the issue because I am getting kx and ky > k0
*/
double kx = x / cell_size;
double ky = y / cell_size;
std::cout << std::sqrt(complex(k0 * k0 - kx * kx - ky * ky)) << std::endl;
//I think you would calculate kz = sqrt(k0^2 - kx^2 - ky^2)
e_prime[y * N + x] *= std::exp(-I * dz * std::sqrt(complex(k0 * k0 - kx * kx - ky * ky)));
}
}

//do the inverse fft
p = fftw_plan_dft_2d(N, N, (fftw_complex*)e_prime, (fftw_complex*)e, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_NO_SIMD);
fftw_execute(p);
fftw_destroy_plan(p);

//this should be what the diffraction pattern looks like.
save("diffracted.ppm", e);

delete[] e;
delete[] e_prime;
return 0;
}


Here is the slit I am diffracting through. Here is the FFT. (Look at the corners) I didn't bother including the last image of the diffraction pattern because it is just black.

• "Is that a typo where he should have ky instead of another kx?"; yes, that's a typo. Whether to diffract or refract first, since it's your program you can try both orders, and compare the results -- step by step -- and note which yields a better result for various test cases. It will be a computational physics experiment! – Peter Diehr Mar 13 '16 at 13:00
• That's some good advice. I didn't realize I can reverse the order. Unfortunately I am getting black output right now so I will try that when I get it working. – chasep255 Mar 13 '16 at 13:01
• Would Computational Science be a better home for this question? – Qmechanic Mar 13 '16 at 13:21
• It might be. I did not know about that forum. Thanks. – chasep255 Mar 13 '16 at 13:22
• My last fourier transform solver is years ago, but I remember in FFTW you need to normalize the output: that is, you need to multiply the inverse FFT by a factor of $N$ -- not sure whether this is done in your code. – davidhigh Apr 19 '17 at 10:17

Your cell_size should be the length of the box, not the resolution of each cell, so use:

 const double cell_size = lambda / 4.0 * N;


instead. And I also think you should be doing:

    double kx = (x-N/2)*2*M_PI/cell_size;
double ky = (y-N/2)*2*M_PI/cell_size;


A couple of side notes:

When using FFTs you should pick your grid size ($N$) to be a power of 2 for optimal efficiency. Also, you should pick units of nanometres so that you don't have those significant exponents knocking around and potentially messing up precision.

Edit In case you copy and pasted this, I originally made a typo of (x-N/2) instead of (y-N/2). Now corrected.

• I made those changes and I am still getting all black output. double kx = (x - N / 2) * 2 * M_PI / cell_size; double ky = (y - N / 2) * 2 * M_PI / cell_size; double kz = sqrt(k0 * k0 - kx * kx - ky * ky); e_prime[y * N + x] *= std::exp(-I * dz * kz); – chasep255 Mar 13 '16 at 13:41
• @chasep255 If you comment out the diffraction step, does it still vanish? – lemon Mar 13 '16 at 13:45
• If I comment it out I just get a picture of the slit on the output which is what is expected. – chasep255 Mar 13 '16 at 13:49
• @chasep255 By the way, the fact that you start with a discontinuous $E$-field may cause some problems (at least with efficiency). Try a Gaussian instead: e[y * N + x]=std::exp(-(dx*dx+dy*dy)/60.); – lemon Mar 13 '16 at 16:54
• @chasep255 I've run your code (with the above changes) and using a Gaussian as the initial dist. I also decreased dz by a factor of 10 and iterated 1,000 times. The result has radial symmetry - indeed, the gaussian seems to be a stable solution. It does scale in height thought, but that seems to be the only problem... – lemon Mar 13 '16 at 17:23