Фурье-метод распространения: как развернуть гауссов пучок без численного взрыва

Гауссов пучок «сжимается» при Фурье-пропагации методом углового спектра. Разбираем знаки в exp(1j*kz*dz) и решение через комплексное сопряжение на практике

Фурье-метод распространения пучка — компактный способ продвигать комплексные поля в пространстве, однако соглашения о знаках и фазовые множители легко приводят к неверному направлению. Типичный симптом — гауссов профиль, который сужается там, где ожидается расширение. Ниже — краткий разбор проблемы и практическая поправка в рамках исходной численной постановки.

Постановка задачи

Цель — распространить гауссов пучок на расстояние dz методом углового спектра. Вычисление выполняется, но визуально пучок идёт «назад»: профиль сжимается около dz = z0, словно поле движется к своей перетяжке. Перестановка fft/ifft результата не изменила. Смена знака в комплексной фазе exp(±1j · kz · dz) приводила к предупреждениям о переполнении и некорректным значениям FFT.

import numpy as np
from numpy.fft import fft2, ifft2, fftfreq
from matplotlib import pyplot as plt
EXP = np.exp
PI = np.pi
SQRT = np.sqrt
ATAN2 = np.arctan2
PHASE = np.angle
plt.rcParams["figure.autolayout"] = True
class RayField:
    def __init__(self, z_r: float, lam: float):
        self.z_r = z_r
        self.lam = lam
        self.k0 = 2 * PI / self.lam
        self.w0 = SQRT(self.lam * self.z_r / PI)
        self.theta0 = self.w0 / self.z_r
        self.a0 = SQRT(2 / (PI * self.w0))
        self.i0 = self.a0 * self.a0.conjugate()
    def w(self, z: np.array):
        return self.w0 * SQRT(1 + (z / self.z_r) ** 2)
    def R(self, z: np.array):
        zsafe = np.where(z == 0.0, 1e-20, z)
        curv = zsafe * (1 + (self.z_r / zsafe) ** 2)
        curv = np.where(z == 0.0, np.inf, curv)
        return curv
    def gouy(self, z: np.array):
        return ATAN2(z, self.z_r)
class GaussianMode(RayField):
    def __init__(self, z_r: float, lam: float):
        super().__init__(z_r, lam)
    def field(self, z=0.0, rho=0.0, phi=0.0):
        z = np.asarray(z, dtype=float)
        k = self.k0
        w0 = self.w0
        a0 = self.a0
        w = self.w(z)
        R = self.R(z)
        g = self.gouy(z)
        U = a0 * (w0 / w) * EXP(-rho ** 2 / w ** 2 + (-k * z - k * (rho ** 2) / (2 * R) + g) * 1.0j)
        return U
lam = 800 * (10 ** -9)
z_r = PI * (10 ** -3) / 2
step = 0.006
def propagate_spectrum(U, dz, npts, spacing):
    fx = fftfreq(npts, d=spacing)
    fy = fftfreq(npts, d=spacing)
    FX, FY = np.meshgrid(fx, fy)
    kx = 2 * PI * FX
    ky = 2 * PI * FY
    kz = SQRT(k ** 2 - kx ** 2 - ky ** 2 + 0j)
    spec = fft2(U)
    spec_out = spec * EXP(1j * kz * dz)  # приводит к «обратному» распространению луча
    U_out = ifft2(spec_out)
    return U_out
model = GaussianMode(z_r, lam)
N = 1000
k = model.k0
W0 = model.w0
xs = np.linspace(-4 * W0, 4 * W0, N)
ys = np.linspace(-4 * W0, 4 * W0, N)
dx = 8 * W0 / (N - 1)
X, Y = np.meshgrid(xs, ys)
RHO = SQRT(X ** 2 + Y ** 2)
PHI = ATAN2(Y, X)
U0 = model.field(z_r, RHO, PHI)
U1 = propagate_spectrum(U0, step, N, dx)
def show(u):
    mag = abs(u)
    ang = PHASE(u)
    plt.figure(figsize=(12, 5))
    plt.subplot(1, 2, 1)
    plt.pcolor(X, Y, mag, shading='auto', cmap='viridis')
    plt.title('Magnitude')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.colorbar()
    plt.subplot(1, 2, 2)
    plt.pcolor(X, Y, ang, shading='auto', cmap='twilight')
    plt.title('Phase')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.colorbar()
    plt.tight_layout()
    plt.show()
show(U1)

Что происходит

В реализации методом углового спектра направление распространения задаётся комплексным фазовым множителем, которым умножается спектр. В этой схеме простая замена 1j на -1j в экспоненте приводила к переполнениям и некорректным значениям FFT. Такое поведение согласуется с фазовой сменой, превращающей часть членов в экспоненциально растущие компоненты, что численно «взрывается». На практике видно, что перестановка fft и ifft не меняла картину, тогда как прямое изменение знака в экспоненте вызывало предупреждения времени выполнения. Итог — поле выглядит так, будто движется к перетяжке, а не от неё.

Рабочее решение

Комплексное сопряжение фазового множителя распространения разворачивает направление без переполнений, которые возникают при грубой смене знака. Иными словами, ядро оставляем как exp(1j · kz · dz) и перед умножением спектра берём от него комплексное сопряжение. После этого расчёт идёт корректно, и пучок распространяется от своей перетяжки.

def propagate_spectrum(U, dz, npts, spacing):
    fx = fftfreq(npts, d=spacing)
    fy = fftfreq(npts, d=spacing)
    FX, FY = np.meshgrid(fx, fy)
    kx = 2 * PI * FX
    ky = 2 * PI * FY
    kz = SQRT(k ** 2 - kx ** 2 - ky ** 2 + 0j)
    spec = fft2(U)
    phase = EXP(1j * kz * dz).conj()  # безопасно разворачивает направление
    spec_out = spec * phase
    U_out = ifft2(spec_out)
    return U_out

Поправка локальная и не требует менять остальную цепочку. Использование сопряжения избавляет от численного «взрыва», наблюдавшегося при смене знака внутри экспоненты.

Почему это важно

В фурье-пропагации даже небольшие расхождения в фазовых соглашениях накапливаются и приводят к заметным артефактам. Неверное направление может проявляться как «прояснение» профиля на неправильной дистанции и уводить к ошибочным выводам о фокусировке. Попытки менять знаки «на глаз» тоже могут дестабилизировать расчёт. Комплексное сопряжение фазового множителя распространения — минимальный и надёжный способ коррекции в той же формулировке.

Выводы

Если кажется, что фурье-распространение идёт в обратную сторону, проверьте фазовый множитель, который умножает угловой спектр. В этой схеме умножение на exp(1j · kz · dz).conj() исправляет направление и избегает переполнений, возникающих при замене 1j на -1j. Это сохраняет устойчивость численного метода и согласует картину с физикой пучка, удаляющегося от перетяжки.