pyoptics.wpm

Wave propagation method (WPM) of Brenner and Singer.

The WPM propagates a scalar field through a medium with a spatially varying refractive index n(x, y, z) in thin slices delta_z. In each slice, the field is decomposed into plane waves, and every plane wave is propagated at every location (x, y) with the local k_z(x, y) = sqrt((n(x, y) k_0)^2 - k_x^2 - k_y^2):

U(x, y, z + delta_z) = sum_{k_x, k_y} A(k_x, k_y) exp(i (k_x x + k_y y))
                                      exp(i k_z(x, y) delta_z).

K.-H. Brenner and W. Singer, Appl. Opt. 32, 4984 (1993).

Two implementations of the same sum are available:

• “direct”: the double sum for every pixel, O(N_x^2 N_y^2) per step, parallelized with OpenMP.

• “levels”: for each distinct index value n_q in the slice, the whole spectrum is propagated with the homogeneous transfer function of n_q by FFT; each pixel takes the result of its own index value. This is the same sum, rearranged, and costs O(Q N_x N_y log(N_x N_y)) for Q distinct index values. Media made of a few materials (lenses, waveguides, gratings) are fast with it. This variant is known as the “fast WPM”.

The default method=”auto” chooses the cheaper one.

Convention: exp(i(k.r - omega*t)), i.e. propagation by delta_z multiplies a plane wave by exp(+i*k_z*delta_z). Absorbing media have n + i*kappa with kappa > 0. Evanescent waves decay (Im(k_z) >= 0).

Sampling: the plane wave decomposition is periodic in the window, and only plane waves up to abs(k_x), abs(k_y) <= pi/dx are represented. To cover all propagating waves in the densest medium, use dx <= wl/(2 n_max), see max_sampling_distance(). Like all FFT propagators, the WPM needs a guard band (or an absorbing index n + i*kappa at the window edges) so that light leaving the window does not re-enter on the other side.

Limitations: the WPM is a scalar method and neglects reflections; at index steps, the transmitted amplitude is not corrected by Fresnel coefficients.

References

K.-H. Brenner and W. Singer, “Light propagation through microlenses: a new simulation method”, Appl. Opt. 32, 4984-4988 (1993).

S. Schmidt, T. Tiess, S. Schröter, R. Hambach, M. Jäger, H. Bartelt, A. Tünnermann and H. Gross, “Wave-optical modeling beyond the thin-element-approximation”, Opt. Express 24, 30188 (2016).

Functions

averaging_z_sampling	Average func(x, y, z) between z1 and z2 by computing the integral
max_sampling_distance	Return the largest sampling distance dx for which the plane wave decomposition contains all propagating waves in a medium with the refractive index n_max: dx = wl/(2 n_max).
wpm_propagate	Propagate a field through one slice with the WPM.
wpm_propagate_z	Propagate a field through the slices between consecutive z_vals and return the field at all z_vals.

pyoptics.wpm.averaging_z_sampling(func, x, y, z1, z2, N)

Average func(x, y, z) between z1 and z2 by computing the integral

1/(z2 - z1) int_{z1}^{z2} func(x, y, z) dz

with an N point Gauss-Legendre quadrature.

Parameters:

funccallable

func(x, y, z), returns an array.

x, yarrays

Passed on to func.

z1, z2number

Integration limits.

Nint

Number of quadrature nodes. N = 1 is the midpoint rule.

Returns:

avgarray

func averaged over z, same shape as func(x, y, z).

pyoptics.wpm.max_sampling_distance(wl, n_max)

Return the largest sampling distance dx for which the plane wave decomposition contains all propagating waves in a medium with the refractive index n_max: dx = wl/(2 n_max).

Parameters:

wlnumber

Vacuum wavelength.

n_maxnumber

Largest (real part of the) refractive index in the problem.

Returns:

dxnumber

pyoptics.wpm.wpm_propagate(n_func, field_z, x, y, z, delta_z, wl, n_sampling='mid', method='auto')

Propagate a field through one slice with the WPM.

Parameters:

n_funccallable

Refractive index n = n_func(X, Y, z) with gridded X, Y (as from np.meshgrid(x, y)). May be complex (n + i*kappa, kappa > 0 absorbs). n_func() is called once per slice (more often for n_sampling=int) with the whole grid, so that it can be vectorized.

field_zarray

Sampled field at z, shape (len(y), len(x)). A one-dimensional field of shape (len(x),) describes a problem invariant in y (pass y=None); n_func then receives one-dimensional X and Y = 0.

x, yarrays

Equidistant coordinate arrays. y is ignored for one-dimensional fields.

znumber

Current z position.

delta_znumber

Distance to propagate.

wlnumber

Vacuum wavelength.

n_sampling“mid”, “start” or int, optional

Where in the slice [z, z + delta_z] the index is sampled: in the middle (default), at the start, or averaged with an n_sampling point Gauss-Legendre quadrature.

method“auto”, “direct” or “levels”, optional

Implementation of the WPM sum, see the module docstring. Both give the same result; “auto” picks the faster one.

Returns:

field_propagatedarray

Field at z + delta_z, same shape as field_z.

Notes

The plane wave decomposition uses the centered DFT of pyoptics.fft, i.e. phases refer to the grid point x[N//2], y[M//2]. Any equidistant grid works, centered or not.

pyoptics.wpm.wpm_propagate_z(n_func, field, x, y, z_vals, wl, n_sampling='mid', method='auto')

Propagate a field through the slices between consecutive z_vals and return the field at all z_vals.

Parameters:

n_funccallable

n_func(X, Y, z), see wpm_propagate().

fieldarray

Field at z_vals[0], shape (len(y), len(x)), or (len(x),) for a problem invariant in y (pass y=None).

x, yarrays

Coordinate arrays.

z_valsarray

Increasing (or decreasing) z positions. Need not be equidistant.

wlnumber

Vacuum wavelength.

n_sampling, methodoptional

See wpm_propagate().

Returns:

fieldsarray

Fields at z_vals, shape (len(z_vals),) + field.shape. fields[0] is the input field.