我正在努力在 C++ 中实现傅里叶光束传播方法。我真的更像是一名程序员而不是物理学家，但我认为我对我正在尝试做的事情有很好的理解。这是我想要实现的。

两个问题。在流程图中他说你应该进行衍射，然后再进行折射。但是，从方程式看来，您首先应用折射，进行 FFT，然后应用衍射。同样在衍射公式中，我注意到 kz 是 omega、kx 和另一个 kx 的函数。这是他应该有 ky 而不是另一个 kx 的错字吗？

我在下面包括我的代码。我正在尝试通过圆形光圈创建衍射图案。我不确定单元格网格的合适大小。它应该是波长的数量级吗？我只是得到空白输出，逆变换中的大多数都是 NaN。我认为问题在于我如何计算、和，请参阅代码中的注释。任何帮助将不胜感激。$k_z$$k_x$$k_y$

#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;
}

这是我正在衍射的狭缝。

这是FFT。（看角落）

我没有费心包括衍射图案的最后一张图像，因为它只是黑色的。