Гауссов пучок и метод углового спектра: как выбрать шаг по z

Разбор ошибки в методе углового спектра на БПФ: слишком малый шаг по z накапливает численную погрешность. Показаны пример кода, метрика и оптимальный выбор dz.

Распространение углового спектра на основе БПФ — надежный способ эволюции поперечного оптического поля на расстоянии, но он может подвести, когда сетка, края и величина шага вступают в конфликт. Обычно добавляют поглощающие границы и разбивают путь на мелкие шаги, чтобы поле не упиралось в края области. Если при больших дистанциях профиль начинает расходиться с аналитическим решением, подозрение часто падает на фильтр. Но в рассмотренном ниже случае фильтр ни при чем: нестабильность возникает из‑за слишком малого шага по z, который заставляет выполнять слишком много пар БПФ/обратное БПФ и накапливает численную ошибку.

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

Цель — стартовать с гауссова поля, дискретизированного при z = z0, распространить его на расстояние d и сравнить численный результат в плоскости z = z0 + d с аналитическим гауссовым пучком в той же плоскости. Для удержания энергии вдали от границ используется пошаговое распространение с поглощающей маской. На относительно коротких дистанциях совпадение хорошее. На длинных — когда пучок приближается к краю — численное решение заметно отклоняется от аналитического, даже при гладкой маске.

Пример кода, воспроизводящий расхождение

Фрагмент ниже демонстрирует основную логику: определение гауссова пучка, одношаговый распространитель в пространственно‑частотной области, многократный запуск с поглощением у края и тест, сравнивающий результат в целевой плоскости с аналитикой. Имена переменных на английском, логика соответствует описанию.

import numpy as np
from numpy.fft import fft2, ifft2, fftfreq
from numpy import exp, pi, sqrt

arctg = np.arctan2

class BeamCore:
  def __init__(self, zR: float, lambda0: float):
    self.zR = zR
    self.lambda0 = lambda0
    self.k = 2 * pi / self.lambda0
    self.w0 = sqrt(self.lambda0 * self.zR / pi)
    self.theta0 = self.w0 / self.zR
    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.zR)**2)

  def R(self, z: np.array):
    z_safe = np.where(z == 0.0, 1e-20, z)
    rval = z_safe * (1 + (self.zR / z_safe)**2)
    rval = np.where(z == 0.0, np.inf, rval)
    return rval

  def gouy(self, z: np.array):
    return arctg(z, self.zR)

class GaussianBeam(BeamCore):
  def __init__(self, zR: float, lambda0: float):
    super().__init__(zR, lambda0)

  def field(self, z: np.array = 0.0, r: np.array = 0.0, phi: np.array = 0.0) -> np.array:
    z = np.asarray(z, dtype=float)
    w = self.w(z)
    R = self.R(z)
    zeta = self.gouy(z)
    U = self.a0 * (self.w0 / w) * exp(-r**2 / w**2 + (-self.k * z - self.k * (r**2) / (2 * R) + zeta) * 1.0j)
    return U

def step_propagate(U, dz, kz):
  UF = fft2(U)
  if dz > 0.0:
    UFprop = UF * exp(1j * kz * dz).conj()
  elif dz < 0.0:
    UFprop = UF * exp(1j * kz * (-dz))
  Uout = ifft2(UFprop)
  return Uout

def propagate_sliced(U, zdist, x1d, y1d, k0, edge_mask):
  nx, ny = U.shape
  dx = (x1d[-1] - x1d[0]) / (nx - 1)
  dy = (y1d[-1] - y1d[0]) / (ny - 1)
  fx = fftfreq(nx, d=dx)
  fy = fftfreq(ny, d=dy)
  FX, FY = np.meshgrid(fx, fy)
  kx = 2 * pi * FX
  ky = 2 * pi * FY
  kz = sqrt(k0**2 - kx**2 - ky**2 + 0j)

  Ustep = U
  while zdist != 0.0:
    U_faded = Ustep * edge_mask
    dz = 0.00005  # очень маленький шаг
    dz = min(dz, zdist)
    Ustep = step_propagate(U_faded, dz, kz)
    zdist -= dz
  return Ustep

# Параметры
lambda0 = 800 * (10**-9)
zR = pi * (10**(-3)) / 2
beam = GaussianBeam(zR, lambda0)
w0 = beam.w0
k0 = beam.k

# Поглощающая маска на большем поле
pad = 500
grid = 1000 + 2 * pad
X = np.linspace(-8 * w0, 8 * w0, grid)
Y = np.linspace(-8 * w0, 8 * w0, grid)
XX, YY = np.meshgrid(X, Y)
r_mask = sqrt(XX**2 + YY**2)

def gauss_edge(r, pad_size):
  sigma = 1
  r_shift = (r - 4 * w0) / w0
  return exp(-(r_shift**2) / (2 * sigma**2))

edge_mask = np.where(r_mask > 4 * w0, gauss_edge(r_mask, pad), 1)

# Рабочая сетка для поля и распространения
x = np.linspace(-8 * w0, 8 * w0, 2000)
y = np.linspace(-8 * w0, 8 * w0, 2000)
Xg, Yg = np.meshgrid(x, y)
r = sqrt(Xg**2 + Yg**2)
phi = arctg(Yg, Xg)

# Начальное поле при z = zR
U0 = beam.field(zR, r, phi)

# Распространяем на 5 * zR, чтобы попасть в z = 6 * zR
U_num = propagate_sliced(U0, 5 * zR, x, y, k0, edge_mask)

# Аналитическое поле при z = 6 * zR
U_ref = beam.field(6 * zR, r, phi)

# Простая метрика расстояния по пикселям
def pix_distance(a: np.array, b: np.array):
  v = abs((a - b).flatten())
  return np.dot(v, v) / len(v)

err = pix_distance(U_num, U_ref)
print(err)

С столь малым шагом отклонение на больших расстояниях велико. Поле стартует корректно, края приглушены, но итоговый профиль не совпадает с аналитическим гауссианом, как должен бы.

Что на самом деле идет не так

Проблема не в поглощающем слое по краю. Ломается всё из‑за слишком малого шага распространения. Каждый шаг делает БПФ и обратное БПФ. Когда шагов неоправданно много, повторные преобразования накапливают численную погрешность. Эта погрешность суммируется на всей траектории и начинает доминировать, когда пучок расширяется.

Решение: увеличить шаг распространения

Увеличение шага уменьшает число пар БПФ/обратное БПФ, и накопленная ошибка падает. В этой конфигурации использование dz = 0.001 дает заметно лучшее совпадение с аналитикой в целевой плоскости.

def propagate_sliced(U, zdist, x1d, y1d, k0, edge_mask):
  nx, ny = U.shape
  dx = (x1d[-1] - x1d[0]) / (nx - 1)
  dy = (y1d[-1] - y1d[0]) / (ny - 1)
  fx = fftfreq(nx, d=dx)
  fy = fftfreq(ny, d=dy)
  FX, FY = np.meshgrid(fx, fy)
  kx = 2 * pi * FX
  ky = 2 * pi * FY
  kz = sqrt(k0**2 - kx**2 - ky**2 + 0j)

  Ustep = U
  while zdist != 0.0:
    U_faded = Ustep * edge_mask
    dz = 0.001  # больший шаг здесь улучшает результат
    dz = min(dz, zdist)
    Ustep = step_propagate(U_faded, dz, kz)
    zdist -= dz
  return Ustep

Оценить улучшение можно той же функцией «пиксельного» расстояния. При dz = 0.001 расстояние около 2.35, тогда как при dz = 0.00005 оно подскакивает примерно до 64.49. Меняйте dz и смотрите, как ведет себя метрика на вашей дальности.

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

В методе углового спектра шаг — это не просто удобный параметр. Он определяет, сколько спектральных преобразований вы выполните. Если шаг слишком мал, метод тратит точность на избыточные обновления. Даже при гладкой поглощающей маске и достаточно большой сетке слишком мелкий шаг может определить весь бюджет ошибки, просто вынуждая делать слишком много циклов БПФ/обратное БПФ.

Вывод

Если пошаговый распространитель на БПФ на дальних дистанциях уходит от аналитического гауссиана, не спешите винить поглощение на краю. Сначала проверьте шаг по z. В приведенной настройке увеличение dz с 0.00005 до 0.001 резко сокращает рассогласование, что подтверждает простая метрика по пикселям. Начните с более крупного шага, сверяйтесь с аналитическим профилем и уже затем уточняйте параметры под вашу дистанцию и сетку.