TransWikia.com

Help with Fourier beam propagation method

Computational Science Asked on July 17, 2021

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.

enter image description here

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.
Slit (Electric field input)

Here is the FFT. (Look at the corners)
enter image description here

I didn’t bother including the last image of the diffraction pattern because it is just black.

One Answer

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.

Answered by lemon on July 17, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP