pyoptics.analytic_phase

Propagation with analytic handling of the wavefront phase.

Behind a lens, near a point source or far from a waist, a field has a strongly curved wavefront. Sampled as a complex array, its phase must change by less than pi between neighbouring samples: a beam of width L with the radius of curvature R needs about L**2/(lambda*|R|) samples per axis, often thousands, although its amplitude is smooth and needs only a few dozen. Wyrowski and coworkers therefore write the field as

u(x, y) = A(x, y)*exp(i*psi(x, y)),

with a smooth, sampled amplitude A and a phase psi that is known analytically, and carry psi through the operators of the propagation (Fourier transforms, multiplication with the transfer function) in closed form (Wang et al. 2019, 2020; Wyrowski and Kuhn 2011). Here psi is a separable quadratic polynomial

psi = alpha_x*x**2 + alpha_y*y**2 + beta_x*x + beta_y*y + gamma,

i.e. a spherical or astigmatic wavefront (alpha = k/(2R)) with a tilt (beta, the carrier) and a constant. QuadraticPhaseField holds A, the grid, the coefficients and the wavelength; the same class holds spectra, with k_x, k_y in place of x, y.

Semi-analytical Fourier transform (SFT, Wang et al. 2019). Per axis, with kappa = k - beta and x_kappa = kappa/(2*alpha), completing the square gives the exact identity

int A(x) exp(i*(alpha*x**2 + beta*x)) exp(-i*k*x) dx
    = exp(-i*kappa**2/(4*alpha)) (A * c)(x_kappa),    c(u) = exp(i*alpha*u**2).

The spectrum is again “amplitude times quadratic phase”, with the analytic phase -kappa**2/(4*alpha) and the amplitude (A * c) on the grid k_j = beta + 2*alpha*x_j. The convolution with the chirp c is done with FFTs (its spectrum is known in closed form, sqrt(i*pi/alpha) exp(-i*nu**2/(4*alpha))). This is exact, and the sampling of the chirps is the only condition: with N samples of spacing dx, the chirp exp(i*alpha*x**2) is sampled in x if rho = |alpha|*N*dx**2/pi < 1 and its counterpart exp(-i*nu**2/(4*alpha)) in the frequency domain if rho > 1. The transforms here pick per axis: for rho < 1 the quadratic phase is multiplied into A and an ordinary (chirp-z) FFT is used, for rho >= 1 the SFT (“mode” forces one of them). In the limit rho >> 1 the chirp convolution becomes the identity: the SFT turns into the pointwise (homeomorphic) transform below.

Propagation (propagate()). The transfer function exp(i*k_z*z) is split into its paraxial part, exp(i*k*z)*exp(-i*z*(k_x**2 + k_y**2)/(2*k)), which is added to the analytic quadratic phase of the spectrum, and the non-paraxial remainder, which multiplies the sampled amplitude (paraxial=False, exact angular spectrum propagation) or is dropped (paraxial=True, exact Fresnel propagation). One forward and one inverse SFT, each with N samples per axis, propagate a converging or diverging beam to any distance; the output grid follows the beam (it scales with (R + z)/R like the geometric footprint), and near a focus the inverse transform switches to the FFT and can be evaluated on a chosen window (x_out, y_out, chirp-z transform).

Homeomorphic Fourier transform (HFT, Wang et al. 2020). If the phase dominates (rho >> 1), the Fourier integral is given by its stationary point x_s, grad psi(x_s) = k, alone

U(k) = A(x_s) 2*pi/sqrt(|det H|) exp(i*pi*sigma(H)/4) exp(i*(psi(x_s) - k.x_s)),

with the Hessian H of psi and its signature sigma: each sample is mapped to one point of the spectrum (homeomorphically), in O(N) operations and without any FFT. homeomorphic_FT() evaluates it for any smooth phase given as an unwrapped sampled array, geometric_propagation() chains it with the transfer function and the inverse transform, which reduces exactly to geometric optics with the transport of the amplitude: the sample at x moves to x + z*grad(psi)/k_z, its amplitude scales with 1/sqrt(|det(I - z H G)|) (G: Hessian of k_z) and it collects the phase z*k**2/k_z plus the Gouy phase -pi/2 per crossed caustic. It fails near foci and caustics, where the SFT (or the FFT) is required.

Tilted planes (tilted_plane()). Propagation onto a plane with a rotated normal is a coordinate transformation of the spectrum: the plane wave with wave vector k reaches the tilted plane with k^ = R^T k (Matsushima, Schimmel and Wyrowski 2003, see pyoptics.tilted). The interpolation of the spectrum onto a uniform k^ grid is the expensive part of the standard method, because the spectrum of a curved beam is a strongly oscillating chirp. Here the spectrum’s phase is analytic: only the smooth amplitude is interpolated, the phase psi~(k(k^)) + k.r0 is evaluated exactly at the new points, its best separable quadratic (incl. the linear carrier, cf. Zhang et al. 2016) is kept analytic and only the small remainder is sampled, followed by an inverse SFT in the plane’s own coordinates. So yes: the approach carries over to tilted planes, with the same small N as for parallel planes.

Functions and classes:

• QuadraticPhaseField: amplitude, grid and analytic phase of a field or of a spectrum; field(), phase(), intensity(), resample(), required_samples()

• semi_analytical_FT(), inverse_semi_analytical_FT(): the transforms

• propagate(): propagation between parallel planes (angular spectrum or Fresnel)

• tilted_plane(): propagation onto a tilted plane

• homeomorphic_FT(): pointwise Fourier transform for any smooth phase

• geometric_propagation(): pointwise propagation (HFT, transfer function, inverse HFT)

Conventions as in pyoptics.propagators: fields vary as exp(i(k.r - omega*t)), U(k) = int u(x) exp(-i*k.x) d^2x and u(x) = 1/(4*pi**2) int U(k) exp(+i*k.x) d^2k (continuous transforms, i.e. including dx*dy). Grids must be linear and equidistant; they need not be centered.

References

Z. Wang, S. Zhang, O. Baladron-Zorita, C. Hellmann and F. Wyrowski, “Application of the semi-analytical Fourier transform to electromagnetic modeling”, Opt. Express 27, 15335-15350 (2019).

Z. Wang, O. Baladron-Zorita, C. Hellmann and F. Wyrowski, “Theory and algorithm of the homeomorphic Fourier transform for optical simulations”, Opt. Express 28, 10552-10571 (2020).

Z. Wang, O. Baladron-Zorita, C. Hellmann and F. Wyrowski, “Generalized Debye integral”, Opt. Express 28, 24459-24470 (2020).

F. Wyrowski and M. Kuhn, “Introduction to field tracing”, J. Mod. Opt. 58, 449-466 (2011).

S. Zhang, D. Asoubar, C. Hellmann and F. Wyrowski, “Propagation of electromagnetic fields between non-parallel planes: a fully vectorial formulation and an efficient implementation”, Appl. Opt. 55, 529-538 (2016).

K. Matsushima, H. Schimmel and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves”, J. Opt. Soc. Am. A 20, 1755-1762 (2003).

Classes

QuadraticPhaseField

A field (or spectrum) u = A*exp(i*psi) with a sampled amplitude A and the analytic phase psi = alpha_x*x**2 + alpha_y*y**2 + beta_x*x + beta_y*y + gamma.

Functions

geometric_propagation	Pointwise propagation: HFT, transfer function exp(i*k_z*z), inverse HFT.
homeomorphic_FT	Pointwise (homeomorphic) Fourier transform of A*exp(i*phase).
inverse_semi_analytical_FT	Field u(x) = 1/(4*pi**2) int U(k) exp(+i*k.x) d^2k of a spectrum; the inverse of semi_analytical_FT() (see there for the parameters).
propagate	Propagate a QuadraticPhaseField by z between parallel planes.
semi_analytical_FT	Spectrum U(k) = int u(x) exp(-i*k.x) d^2x of a QuadraticPhaseField.
tilted_plane	Propagate a QuadraticPhaseField from z = 0 onto a tilted plane.

class pyoptics.analytic_phase.QuadraticPhaseField(amplitude, x, y, wavelength, alpha=(0.0, 0.0), beta=(0.0, 0.0), gamma=0.0, n=1.0, domain='space')

Bases: object

A field (or spectrum) u = A*exp(i*psi) with a sampled amplitude A and the analytic phase psi = alpha_x*x**2 + alpha_y*y**2 + beta_x*x + beta_y*y + gamma.

Parameters:

amplitudearray

Complex amplitude A of shape (Ny, Nx). It must be smooth (sampled according to its own bandwidth); the wavefront goes into psi.

x, yarrays

Linear, ascending, equidistant coordinates (need not be centered). For a spectrum: the wavenumbers k_x, k_y.

wavelengthnumber

Vacuum wavelength.

alpha, betapairs of numbers, optional

Quadratic and linear coefficients of psi in x and y.

gammanumber, optional

Constant phase.

nnumber, optional

Refractive index of the medium.

domainstr, optional

“space” (default) or “spectrum”.

classmethod from_curvature(amplitude, x, y, wavelength, R=inf, carrier=(0.0, 0.0), n=1.0)

Field A*exp(i*k*(x**2/(2*R_x) + y**2/(2*R_y)))*exp(i*(c_x*x + c_y*y)).

R > 0 is a diverging, R < 0 a converging wave (focus at z = -R), R = inf a plane wavefront; R may be a pair (R_x, R_y). carrier is the tilt (c_x, c_y) as transverse wavenumbers.

classmethod from_field(field, x, y, wavelength, R=inf, carrier=(0.0, 0.0), n=1.0)

Split a sampled field into the analytic phase given by R and carrier (see from_curvature()) and the remaining amplitude.

property k

Wavenumber in the medium.

property radius_of_curvature

(R_x, R_y) = k/(2*alpha) of a field (inf for alpha = 0).

phase(x=None, y=None)

Analytic phase psi on the grid (or on the linear coordinates x, y).

field()

The complete field A*exp(i*psi) on the grid. Aliased where psi changes by more than pi between samples (see required_samples()).

intensity()

|A|**2: the intensity does not depend on psi.

required_samples()

Samples (N_x, N_y) that the complete field A*exp(i*psi) needs on the same window: the window length times the largest local frequency of psi, over pi, plus the samples of A.

chirp_ratio()

rho = |alpha|*N*d**2/pi per axis: the transforms use the SFT for rho >= 1 and the FFT (phase multiplied into A) for rho < 1.

resample(x, y)

The same field on other linear coordinates: the amplitude is interpolated (cubic splines, zero outside the grid), the phase is evaluated analytically. Accurate as long as A is smooth (a RuntimeWarning is issued if its phase changes by more than 0.25 rad between samples (weighted with the amplitude), e.g. behind an FFT step that multiplied the quadratic phase into A: pass the output grid to the transform instead).

pyoptics.analytic_phase.semi_analytical_FT(field, kx_out=None, ky_out=None, mode='auto', pad=2)

Spectrum U(k) = int u(x) exp(-i*k.x) d^2x of a QuadraticPhaseField.

Parameters:

fieldQuadraticPhaseField

Field in the space domain.

kx_out, ky_outarrays, optional

Output wavenumbers (linear, equidistant). Axes transformed with the FFT are evaluated there directly (chirp-z transform), axes transformed with the SFT are interpolated (zero outside the natural window of each). Default: the natural grid of each transform (N samples).

modestr, optional

“auto” (default: SFT where the chirp ratio rho >= 1, see the module docstring), “sft” or “fft”.

padnumber, optional

Zero padding factor of the chirp convolution in the SFT.

Returns:

spectrumQuadraticPhaseField

Amplitude and analytic phase of U on the wavenumbers x = k_x, y = k_y (domain “spectrum”).

pyoptics.analytic_phase.inverse_semi_analytical_FT(spectrum, x_out=None, y_out=None, mode='auto', pad=2)

Field u(x) = 1/(4*pi**2) int U(k) exp(+i*k.x) d^2k of a spectrum; the inverse of semi_analytical_FT() (see there for the parameters).

pyoptics.analytic_phase.propagate(field, z, paraxial=False, x_out=None, y_out=None, mode='auto', pad=2)

Propagate a QuadraticPhaseField by z between parallel planes.

The forward transform, the transfer function and the inverse transform all keep the quadratic phase analytic (see the module docstring), so that N samples per axis suffice at any distance if the amplitude is smooth. The output grid follows the beam: away from a focus it is the input grid scaled by about (R + z)/R; near a focus (chirp ratio < 1) the inverse transform is an FFT whose window can be chosen with x_out, y_out.

Parameters:

fieldQuadraticPhaseField

Field in the plane z = 0.

znumber

Propagation distance (negative: backwards).

paraxialboolean, optional

If True, propagate with the Fresnel transfer function exp(i*k*z)*exp(-i*z*(k_x**2 + k_y**2)/(2*k)), which is handled entirely analytically. If False (default), the non-paraxial remainder of exp(i*k_z*z) multiplies the sampled spectrum: exact angular spectrum (Rayleigh-Sommerfeld) propagation, as long as the remainder is sampled (a RuntimeWarning is issued if it changes by more than 0.25 rad between samples, weighted with the amplitude). Evanescent waves decay.

x_out, y_outarrays, optional

Output coordinates (linear, equidistant), see semi_analytical_FT().

mode, padoptional

See semi_analytical_FT().

Returns:

field_zQuadraticPhaseField

Field in the plane z.

pyoptics.analytic_phase.tilted_plane(field, z0, R=None, angles=None, x_out=None, y_out=None, N_out=None, k_margin=0.05, mode='auto', pad=2, threshold=1e-06)

Propagate a QuadraticPhaseField from z = 0 onto a tilted plane.

Geometry as in pyoptics.tilted: the target plane passes through r0 = (0, 0, z0) (or r0 = z0 for a 3-vector), its local axes x^, y^ and normal are the columns of the rotation matrix R, and the field is returned in the local coordinates (x^, y^).

Steps (see the module docstring): SFT of the source; a uniform grid of target wavenumbers k^ covering the rotated spectrum; k = R k^ and the Jacobian k_z/k^_z; cubic interpolation of the smooth spectral amplitude only; the exact phase psi~(k) + k.r0 at the new points, of which the best separable quadratic (incl. the carrier) is kept analytic; inverse SFT in the plane’s coordinates. Only plane waves with k_z > 0 and k^_z > 0 are kept.

Parameters:

fieldQuadraticPhaseField

Field in the plane z = 0.

z0number or array

Distance of the target plane’s origin, or its position (x0, y0, z0).

Rarray, optional

Rotation matrix (e.g. geometry.R_y(theta)).

anglestuple, optional

(theta_x, theta_y) instead of R, see tilted.rotation_from_angles().

x_out, y_outarrays, optional

Output coordinates in the plane (linear, equidistant), see semi_analytical_FT(). Default: the natural grid of the inverse transform.

N_outint or pair, optional

Samples of the target spectrum (and field) per axis. Default: twice the samples of the source (the remainder of the phase that is not quadratic in k^ must be sampled; a RuntimeWarning is issued otherwise).

k_marginnumber, optional

Relative margin of the target spectrum around the rotated support.

mode, padoptional

See semi_analytical_FT().

thresholdnumber, optional

Relative intensity of the spectrum that defines its support.

Returns:

field_tQuadraticPhaseField

The field on the tilted plane, in its local coordinates.

pyoptics.analytic_phase.homeomorphic_FT(amplitude, phase, x, y)

Pointwise (homeomorphic) Fourier transform of A*exp(i*phase).

Every sample x_j is mapped to the wavenumber k_j = grad phase(x_j), where it is the stationary point of the Fourier integral, and

U(k_j) = A_j 2*pi/sqrt(|det H_j|) exp(i*pi*sigma_j/4) exp(i*(phase_j - k_j.x_j))

(H: Hessian of the phase, sigma: its signature). Valid if the phase dominates (many oscillations across the beam, no caustic: det H != 0 and A slowly varying on the scale of the stationary zone).

Parameters:

amplitudearray

Smooth amplitude A, shape (Ny, Nx).

phasearray

Smooth, unwrapped phase on the same grid (it need not be sampled below pi per sample: unwrapped values are smooth).

x, yarrays

Linear coordinates.

Returns:

kx, kyarrays

Wavenumbers (non-uniform, shape (Ny, Nx)) of the spectrum samples.

U_amplitudearray

Complex amplitude (incl. 2*pi/sqrt(|det H|) and the signature phase).

U_phasearray

Unwrapped phase phase - k.x of the spectrum.

pyoptics.analytic_phase.geometric_propagation(amplitude, phase, x, y, z, wavelength, n=1.0)

Pointwise propagation: HFT, transfer function exp(i*k_z*z), inverse HFT.

The chain of the two stationary phase evaluations is evaluated in closed form (no division by det H, so plane wavefronts are allowed): the sample at x moves along its ray to x + z*k_perp/k_z with k_perp = grad phase, and

u = A/sqrt(|det M|) exp(i*(phase + z*k**2/k_z)) exp(-i*pi*n_neg/2),
M = I - z*H*G,

with the Hessians H of the phase and G of k_z(k_x, k_y), and the number n_neg of negative eigenvalues of M (foci crossed: the Gouy phase). This is geometrical optics with amplitude transport; it diverges at foci and caustics (det M = 0) and ignores diffraction at the edges of A.

Parameters:

amplitude, phasearrays

Smooth amplitude and unwrapped phase, shape (Ny, Nx).

x, yarrays

Linear coordinates.

znumber

Distance.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

Returns:

x_z, y_zarrays

Positions (shape (Ny, Nx), non-uniform) of the samples in the plane z.

u_amplitudearray

Complex amplitude there (incl. the Gouy phase).

u_phasearray

Unwrapped phase there, phase + z*k**2/k_z.