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.