傅里叶光束传播代码反向传播问题及正向实现修正咨询

我正在编写一段基于傅里叶传播的光束传播代码，但遇到了反向传播的问题——当设置传播距离dz=z_0时，光束轮廓反而收缩，看起来像是向束腰方向传播，而非正向传播。我尝试过交换fft和ifft、修改传播相位的符号，但交换变换函数没效果，改负号还出现了溢出和无效值的警告。想请教下为什么会出现反向传播，以及如何修正实现正向传播？

我的相关代码片段 #

import numpy as np
import scipy
from matplotlib import pyplot as plt
from numpy.fft import fft, ifft, fft2, ifft2, fftfreq, fftshift, ifftshift

exp=np.exp
pi=np.pi
sqrt=np.sqrt
tg=np.tan
arctg=np.arctan2
angle=np.angle
array=np.array

plt.rcParams["figure.autolayout"] = True

class beam:
    def __init__(self, z_0: float, compOnda: float):
        self.z_0 = z_0 #distância de Rayleigh
        self.compOnda = compOnda #comprimento de onda
        self.k = 2*pi / self.compOnda #número de onda
        self.W_0 = sqrt(self.compOnda * self.z_0 / pi)
        self.teta_0 = self.W_0 / self.z_0
        self.A_0 = sqrt(2 / (pi * self.W_0)) #amplitude em 0; definimos ela de forma a normalizar o pulso (gaussiano - preciso verificar se normaliza o laguerre gaussiano também)
        self.I_0 = self.A_0 * self.A_0.conjugate() #intensidade em 0

    def W(self,z: np.array):
        return self.W_0 * sqrt(1 + (z / self.z_0)**2)

    def R(self,z: np.array):
        z_=np.where(z==0.0,1e-20,z) #evitando divisão por zero
        R=z_ * (1 + (self.z_0 / z_)**2)
        R=np.where(z==0.0,np.inf,R)
        return R

    def zeta(self,z: np.array):
        return arctg(z, self.z_0)

class GaussianBeam(beam):
    def __init__(self, z_0: float, compOnda: float):
        super().__init__(z_0,compOnda)

    def Amplitude(self,z: array = 0.0, rho: array = 0.0, phi: array = 0.0) -> array:
        z=np.asarray(z,dtype=float)
        z_0=self.z_0
        A_0=self.A_0
        k=self.k
        W_0=self.W_0
        W=self.W(z)
        R=self.R(z)
        zeta=self.zeta(z)
        U=A_0*(W_0/W)*exp(-rho**2/W**2 + (-k*z-k*(rho**2)/(2*R)+zeta)*1.0j)
        return U

wavelength=800*(10**-9) #nm
z_0=pi*(10**(-3))/2 #m ->approx. 0.0015.
W_0=20.10⁻6
dz=0.006 #propagation distance (m)

def Propag(U,dz,gr,d):
    #frequency coordinates
    fx = fftfreq(gr, d=d)
    fy = fftfreq(gr, d=d)
    fx, fy = np.meshgrid(fx, fy)
    kx = 2*pi*fx
    ky = 2*pi*fy
    kz = sqrt(k**2 - kx**2 - ky**2 + 0j) #propagation constant in z
    #radial complex amplitude
    Fourier_U = fft2(U)
    #angular spectrum
    Propag_FU = Fourier_U * exp(1j * kz * dz)
    U_propag = ifft2(Propag_FU)
    return U_propag

fx=GaussianBeam(z_0,wavelength)
W_0=fx.W_0
gr=1000
k=fx.k
x = np.linspace(-4*W_0,4*W_0,gr) #grid size
y = np.linspace(-4*W_0,4*W_0,gr)
d=8*W_0/(gr-1)
X,Y = np.meshgrid(x,y)
rho=sqrt(X**2+Y**2)
phi=arctg(Y,X)

U=fx.Amplitude(z_0,rho,phi) #this is a bidimensional Gaussian
U_propag=Propag(U,dz,gr,d)

def plot(u):
    modulo = abs(u)
    fase = angle(u)
    # Gráfico
    plt.figure(figsize=(12, 5))
    plt.subplot(1, 2, 1)
    plt.pcolor(X, Y, modulo, shading='auto', cmap='viridis')
    plt.title('Módulo')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.colorbar()
    plt.subplot(1, 2, 2)
    plt.pcolor(X, Y, fase, shading='auto', cmap='twilight')
    plt.title('Fase')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.colorbar()
    plt.tight_layout()
    plt.show()

plot(U_propag)

问题原因分析 #

1. 初始场位置错误：你当前生成的初始场是z=z0处的高斯光束，而非束腰位置（z=0）的光束。当你尝试传播dz=z0时，相当于从z=z0往束腰方向（z=0）移动，自然会看到光束收缩，看起来像是反向传播。

2. 傅里叶频率轴与相位不匹配：numpy.fftfreq生成的频率轴默认是正频率在前、负频率在后，和光学角谱的常规中心对称顺序不一致，直接应用传播相位可能导致方向混淆。

3. 倏逝波导致数值溢出：当你修改相位符号为负时，kx²+ky² > k²的分量会让kz变成虚数，exp(-1j*kz*dz)会变成指数增长的形式，直接引发数值溢出。

修正方案 #

1. 修正初始场位置（最直接的解决方法）

生成束腰位置（z=0）的高斯光束，这样传播dz就是向z>0的正向传播：

# 替换原来的初始场生成代码
U=fx.Amplitude(0.0,rho,phi) # 生成z=0处的束腰高斯光束

此时传播dz=z0，光束宽度会变为W0*sqrt(2)，符合正向传播的扩散预期。

2. 调整傅里叶传播的频率轴与相位计算

如果需要保留z=z0的初始场，需要对齐频率轴顺序并过滤倏逝波：

def Propag(U,dz,gr,d):
    #frequency coordinates
    fx = fftfreq(gr, d=d)
    fy = fftfreq(gr, d=d)
    fx, fy = np.meshgrid(fx, fy)
    kx = 2*pi*fx
    ky = 2*pi*fy
    
    # 过滤倏逝波，只保留可传播的分量
    k_sq = k**2 - kx**2 - ky**2
    kz = np.sqrt(np.maximum(k_sq, 0) + 0j)
    
    # 对齐频率轴到中心对称顺序
    Fourier_U = fftshift(fft2(U))
    # 应用正向传播相位
    Propag_FU = Fourier_U * exp(1j * kz * dz)
    # 逆移位后做逆FFT
    U_propag = ifft2(ifftshift(Propag_FU))
    return U_propag

3. 修复语法错误

你的代码中W_0=20.10⁻6是语法错误，应该改为W_0=20*(10**-6)，否则数值会完全不符合预期。

验证修正效果 #

修正后，当设置dz=z0时，光束会呈现出正常的扩散效果，宽度变为束腰宽度的sqrt(2)倍，符合高斯光束正向传播的物理规律。

内容来源于stack exchange