pyoptics.raywave

Conversion between ray ensembles and sampled wave fields.

Geometrical optics describes a wave field by its rays: the eikonal S (the optical path, OPL) gives the phase k_0*S, the rays are the normals of its wavefronts, and the transport equation div(I grad S) = 0 says that the power in a tube of rays is conserved. This module goes back and forth between the two descriptions.

Rays to waves

Ray field summation (rays_to_field()). Rays launched on a structured grid (e.g. a grid of pupil coordinates) are traced to a plane; the grid of their intersections is cut into triangles, each a small tube of rays. A triangle carries the power of its launch triangle, so its irradiance in the plane is the launch irradiance times the ratio of the launch and target areas (the inverse Jacobian of the ray map, Born and Wolf ch. 3.1). Inside the triangle the phase is interpolated from the eikonal S and its gradient n*d_perp at the three rays (the mean of the linear interpolation and of the tangent planes, which is exact for quadratic wavefronts), the field of the triangle is

u = sqrt(I_launch*|A_launch/A_target|)*exp(i*k_0*S)*exp(-i*pi*m/2),

and the fields of all triangles that cover a sample are added. Where the ray map folds (behind a caustic), several sheets of triangles cover the same point and interfere; m, the Maslov (KMAH) index, counts the caustics that the tube has crossed: the number of negative eigenvalues of the 2 x 2 affine map from the launch triangle to the target triangle (1 at a line focus, 2 at a point focus: the Gouy phase -pi). This is the classical asymptotic field. It is exact far from caustics, where the field is locally a sum of plane waves with slowly varying amplitude, and it diverges on the caustic itself (A_target -> 0), where the rays touch and geometrical optics breaks down; in the shadow of a caustic it gives zero instead of the evanescent Airy tail (Kravtsov and Orlov 1993).

Gaussian beamlets (GaussianBeamlets). The standard remedy is to give every ray a finite width: the field in a launch plane is written as a sum of Gaussian beamlets centred on the rays, each with the local tilt and curvature of the wavefront (Červený, Popov and Pšenčík 1982; Arnaud 1985; the Gaussian beam decomposition of lens design programs). A frame of Gaussians exp(-r**2/w**2) on a grid of spacing Delta adds up to pi*w**2/Delta**2 (to better than 1e-8 for w >= 1.5*Delta), so the beamlet of a sample with complex amplitude A and the area Delta**2 has the amplitude A*Delta**2/(pi*w**2) and

u_j(x) = a_j*exp(i*(phi_j + k_j.xi + xi.Gamma_j.xi/2)),   xi = x - x_j,
Gamma_j = H_j + 2i/w_j**2,

with the transverse wave vector k_j and the Hessian H_j of the phase. In a homogeneous medium the angular spectrum of every beamlet is propagated in closed form with k_z expanded to second order about k_j (k_z ~ k_z0 + g.dk + dk.G.dk/2, any tilt):

u_j(x, z) = a_j/sqrt(det(I - z*Gamma_j*G))
            *exp(i*(phi_j + z*k_z0 + k_j.xi + eta.Gamma_z.eta/2)),
eta = xi + z*g = xi - z*k_j/k_z0,   Gamma_z = Gamma_j*(I - z*G*Gamma_j)^-1.

The centre moves along the ray, the width and curvature follow the complex q parameter, and the square root (taken on the continuous branch: the eigenvalues of Gamma*G have negative imaginary parts) gives the Gouy phase. The beamlets stay finite at foci and caustics, so their sum is the field everywhere, as long as the beamlets are narrow compared to the scale on which the wavefront curvature changes and wide compared to the wavelength. Beamlets can be launched from rays (GaussianBeamlets.from_rays(), curvature from neighbouring rays) or from a sampled field (GaussianBeamlets.from_field(), also from a QuadraticPhaseField of pyoptics.analytic_phase, whose analytic phase supplies tilt and curvature).

GaussianBeamlets.trace() carries the beamlets through the surfaces of a raytrace.System (Gaussian beam decomposition through real optics). Every beamlet is followed by its base ray and eight parabasal rays (centre and direction shifted in x and y), traced exactly. Their central differences in the output plane give the 4 x 4 ray transfer matrix M = [[A, B], [C, D]] of the neighbourhood of the base ray in the canonical coordinates (x, k_perp), which is symplectic, and the beamlet follows the complex ABCD law (a Gaussian beamlet is the family of complex rays k_perp = Gamma*xi, Arnaud 1985)

Gamma' = (C + D Gamma)(A + B Gamma)^-1,   a' = a/sqrt(det(A + B Gamma)),

with the phase k_0*OPL of the base ray. det(A + B Gamma) never vanishes, but its argument grows by pi through every line focus: the square root is taken on the branch that follows the beamlet (its Gouy phase), accumulated in closed form along every straight segment of the base ray. Beamlets whose rays are vignetted are dropped, so an aperture inside a system is resolved to the beamlet size; for a sharp cut, sum the beamlets in the aperture plane and decompose the field again with from_field().

Waves to rays

Phase gradient (field_to_rays()). The time averaged energy flux of a scalar field, the Poynting vector of the paraxial vector field, is proportional to Im(conj(u)*grad u) = |u|**2 grad(phase): the local ray direction is the phase gradient, k_perp = grad(phase), completed with k_z = sqrt(n**2 k_0**2 - |k_perp|**2), each sample becomes a ray with the power |u|**2 dx dy and the optical path phase/k_0. The gradient is computed from the phase differences of neighbouring samples (no unwrapping needed for the directions). The result can be traced with raytrace.System.trace(). It is the right picture where the field is locally a single plane wave; inside a caustic, where several sheets overlap, the local phase gradient is an average direction weighted with the interference pattern and need not describe any ray.

Local spectra and Wigner functions (local_spectrum_rays(), wigner()). Phase space distributions assign several directions to one point. The Wigner function

W(x, k) = int u(x + s/2) conj(u(x - s/2)) exp(-i*k*s) ds

(here in 1D, for a section of a field) has the intensity and the spectrum as marginals, is concentrated on the Lagrangian manifold of the ray family, k = dS/dx, which folds at caustics, and oscillates (negative values) between the sheets (Bastiaans 1979; Alonso 2011). Its smoothed, positive version is the spectrogram (Husimi function): the spectrum of the field in a Gaussian window. local_spectrum_rays() samples it in 2D: every window centre and every spectral component above a threshold is a ray, with the power of that component. It resolves crossing ray sheets at the price of the uncertainty relation (position to the window size, direction to 2*pi/window), and the rays carry no phase (an incoherent ray set).

Helpers: straight_rays() moves rays to another plane of a homogeneous medium, ray_irradiance() bins ray powers into an irradiance map.

Conventions as in pyoptics.propagators and pyoptics.raytrace: fields vary as exp(i(k.r - omega*t)), the phase of the field is +k_0*OPL (vacuum wavenumber times the optical path, which includes the refractive index), the light travels in +z, ray directions are unit vectors. Fields are arrays of shape (Ny, Nx) on linear coordinates x, y. |u|**2 is the irradiance on the plane (power per area of the plane, not per area normal to the ray).

References

M. Born and E. Wolf, “Principles of Optics”, 7th ed., Cambridge University Press (1999), ch. 3.

Yu. A. Kravtsov and Yu. I. Orlov, “Caustics, Catastrophes and Wave Fields”, Springer (1993).

V. P. Maslov and M. V. Fedoriuk, “Semi-Classical Approximation in Quantum Mechanics”, Reidel (1981).

V. Červený, M. M. Popov and I. Pšenčík, “Computation of wave fields in inhomogeneous media - Gaussian beam approach”, Geophys. J. R. Astr. Soc. 70, 109-128 (1982).

J. A. Arnaud, “Representation of Gaussian beams by complex rays”, Appl. Opt. 24, 538-543 (1985).

A. W. Greynolds, “Propagation of generally astigmatic Gaussian beams along skew ray paths”, Proc. SPIE 560, 33-50 (1985).

M. J. Bastiaans, “Wigner distribution function and its application to first-order optics”, J. Opt. Soc. Am. 69, 1710-1716 (1979).

M. A. Alonso, “Wigner functions in optics: describing beams as ray bundles and pulses as particle ensembles”, Adv. Opt. Photon. 3, 272-365 (2011).

G. W. Forbes and M. A. Alonso, “Using rays better. I. Theory for smoothly varying media”, J. Opt. Soc. Am. A 18, 1132-1145 (2001).

Classes

FieldRays	Rays from a field: start points and unit directions, shape (N, 3), optical path (N,) (NaN if the rays carry no phase) and power (N,).
GaussianBeamlets	A field written as a sum of Gaussian beamlets, propagated in closed form.

Functions

field_to_rays	Rays along the local phase gradient (energy flux) of a field.
local_spectrum_rays	Rays from the local spectra of a field (sampled spectrogram / Husimi function).
ray_irradiance	Irradiance of a ray ensemble: the powers binned on a grid, per area.
rays_to_field	Sampled field from rays on a structured grid (ray field summation).
straight_rays	Move rays along straight lines to the plane z of a homogeneous medium.
wigner	Wigner distribution of a 1D field (e.g. a section of a 2D field).

class pyoptics.raywave.FieldRays(points, directions, opl, power)

Bases: tuple

Rays from a field: start points and unit directions, shape (N, 3), optical path (N,) (NaN if the rays carry no phase) and power (N,). Has the attributes of raytrace.Rays, so it can be passed to System.trace() directly.

directions

Alias for field number 1

opl

Alias for field number 2

points

Alias for field number 0

power

Alias for field number 3

pyoptics.raywave.straight_rays(points, directions, opl, z, n=1.0)

Move rays along straight lines to the plane z of a homogeneous medium.

Parameters:

points, directionsarrays, shape (…, 3)

Positions and unit directions.

oplarray, shape (…)

Optical path at the positions.

znumber

z of the target plane.

nnumber, optional

Refractive index of the medium.

Returns:

points, oplarrays

Intersections with the plane and the optical path there (the path is negative for rays that travel backwards to reach the plane).

pyoptics.raywave.ray_irradiance(points, power, x_edges, y_edges)

Irradiance of a ray ensemble: the powers binned on a grid, per area.

Parameters:

pointsarray, shape (…, 2) or (…, 3)

Ray positions in the plane (x, y; z is ignored).

powerarray or number

Power of each ray.

x_edges, y_edgesarrays

Bin edges.

Returns:

Earray, shape (len(y_edges) - 1, len(x_edges) - 1)

Power per area of the bins (the geometrical intensity).

pyoptics.raywave.rays_to_field(x, y, points, directions, opl, wavelength, launch=None, irradiance=1.0, n=1.0, gouy=True, max_pairs=4000000)

Sampled field from rays on a structured grid (ray field summation).

Every cell of the ray grid is split into two triangles (tubes of rays). A triangle carries the power of its launch triangle, its irradiance in the target plane is the launch irradiance times the area ratio, its phase is interpolated from the eikonal and the ray directions of its corners, and the fields of all triangles covering a sample are added (several sheets behind a caustic). See the module docstring.

Parameters:

x, yarrays

Linear, equidistant coordinates of the output grid.

pointsarray, shape (Nv, Nu, 2) or (Nv, Nu, 3)

Ray intersections with the target plane, on the structured launch grid (neighbouring rays are neighbours in the array). NaN for lost rays; their triangles are skipped.

directionsarray, shape (Nv, Nu, 3)

Unit ray directions at the target plane.

oplarray, shape (Nv, Nu)

Optical paths at the target plane.

wavelengthnumber

Vacuum wavelength.

launcharray, shape (Nv, Nu, 2), optional

Positions of the rays in the launch plane (e.g. entrance pupil coordinates in length units, or the start points). Default: the array indices (u, v), so that irradiance is the power per grid cell.

irradiancenumber or array (Nv, Nu), optional

Irradiance at the launch positions (power per launch area).

nnumber, optional

Refractive index at the target plane.

gouybool, optional

Include the phase -pi/2 per crossed caustic (Maslov index).

max_pairsint, optional

Triangle-sample pairs evaluated at once (memory).

Returns:

uarray, shape (len(y), len(x))

Complex field with |u|**2 the geometrical irradiance and phase k_0*OPL.

class pyoptics.raywave.GaussianBeamlets(centres, k_perp, phase, amplitude, Gamma, wavelength, n=1.0, z0=0.0)

Bases: object

A field written as a sum of Gaussian beamlets, propagated in closed form.

Beamlet j in the launch plane z0

u_j = a_j*exp(i*(phase_j + k_j.xi + xi.Gamma_j.xi/2)),  xi = (x, y) - centre_j,

with the transverse wave vector k_j and the complex symmetric matrix Gamma_j = H_j + 2i/w_j**2 (H: Hessian of the wavefront phase, w: 1/e amplitude radius). field() evaluates the sum in any plane z of the homogeneous medium (see the module docstring). Usually built with from_rays() or from_field().

Parameters:

centresarray, shape (N, 2)

Positions (x, y) in the launch plane.

k_perparray, shape (N, 2)

Transverse wave vectors (k_x, k_y), |k_perp| < n*k_0.

phasearray, shape (N,)

Phase at the centres.

amplitudearray, shape (N,)

Complex amplitudes a_j.

Gammaarray, shape (N, 2, 2)

Complex curvature matrices (positive definite imaginary part).

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index of the medium.

z0number, optional

z of the launch plane.

property k

Wavenumber in the medium.

classmethod from_field(field, x, y, wavelength, n=1.0, step=4, overlap=1.5, phase=None, threshold=0.0001, z0=0.0)

Beamlets on a sub grid of a sampled field.

Parameters:

fieldarray (Ny, Nx) or analytic_phase.QuadraticPhaseField

The field. For a QuadraticPhaseField, tilt and curvature of the analytic phase are taken exactly and only the smooth amplitude is differentiated (x, y and wavelength are then taken from it and the arguments ignored).

x, yarrays

Linear coordinates.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

stepint, optional

Beamlet spacing in samples.

overlapnumber, optional

Beamlet radius w over the spacing (>= 1.5 for a flat frame).

phasearray, optional

Unwrapped phase of the field; default: unwrapped by integrating the wrapped phase differences.

thresholdnumber, optional

Drop beamlets with |u|**2 below threshold times the maximum.

z0number, optional

z of the field.

classmethod from_rays(points, directions, opl, wavelength, launch=None, irradiance=1.0, n=1.0, overlap=1.5)

One beamlet per ray of a structured ray grid in a plane.

The wavefront curvature at every ray comes from the directions of its neighbours (H = dk_perp/dx, finite differences on the grid), the beamlet radius from the local ray spacing, w = overlap*sqrt(area), the amplitude from the irradiance, carried from the launch plane like in rays_to_field().

Parameters:

pointsarray, shape (Nv, Nu, 3)

Ray intersections with the launch plane of the beamlets (all with the same z).

directionsarray, shape (Nv, Nu, 3)

Unit directions.

oplarray, shape (Nv, Nu)

Optical paths.

wavelengthnumber

Vacuum wavelength.

launch, irradianceoptional

As in rays_to_field(): positions and irradiance of the rays where they were launched (default: the irradiance in the beamlet plane is the power per grid cell over the cell area, i.e. the given irradiance is the power per cell).

nnumber, optional

Refractive index of the medium.

overlapnumber, optional

Beamlet radius over the ray spacing.

propagated(z)

Centres, widths and curvatures in the plane z (see field()).

Returns:

centresarray (N, 2)

Where the beamlet axes cross the plane.

Gamma_zarray (N, 2, 2)

Complex curvature matrices there.

factorarray (N,)

Complex amplitude factor a_j/sqrt(det(I - (z - z0) Gamma G)) times the phase exp(i*(phase_j + (z - z0)*k**2/k_z)) on the axis.

field(x, y, z=None, max_pairs=20000000)

The sum of the beamlets at the points (x, y) of the plane z.

Parameters:

x, yarrays

Coordinates, broadcast against each other (e.g. a meshgrid, or a line).

znumber, optional

Plane (default: the launch plane z0).

max_pairsint, optional

Beamlet-point pairs evaluated at once (memory).

Returns:

ucomplex array of the broadcast shape of x and y

trace(system, z=None, clip=True, delta=0.001, details=False)

Carry the beamlets through the surfaces of a raytrace.System.

Every beamlet is followed by its base ray and eight parabasal rays (its centre shifted by +-delta*w in x and y, its transverse wave vector by +-delta*2/w in k_x and k_y), all traced exactly through the system. In the output plane, central differences of the parabasal rays give the 4 x 4 ray transfer matrix M = [[A, B], [C, D]] of the neighbourhood of the base ray, in the canonical coordinates (x, y, k_x, k_y) of the launch plane and of the output plane, and the beamlet follows the complex ABCD law (Arnaud 1985)

Gamma' = (C + D Gamma)(A + B Gamma)^-1,
a' = a/sqrt(det(A + B Gamma)),  phase' = k_0*OPL of the base ray.

The square root is taken on the branch that follows the beamlet through the system (its Gouy phase): the argument of det(A + B Gamma) is accumulated in planes normal to the base ray, in closed form along every straight segment and by the (small) change of the projection at every surface. A beamlet is dropped if one of its rays is lost (vignetting with clip=True, total internal reflection, missed surfaces). Mirrors fold the light out of +z; the output plane must lie in a medium where the light travels in +z.

Parameters:

systemraytrace.System

The system, with the same wavelength; the beamlets lie in its object space (index system.n[0]), in front of the first surface.

znumber, optional

Output plane (global z) behind the last surface; default: the image plane of the system.

clipbool, optional

Drop beamlets that hit a semi diameter.

deltanumber, optional

Size of the parabasal offsets relative to the beamlet.

detailsbool, optional

Also return the transfer matrices and Gouy phases.

Returns:

beamletsGaussianBeamlets

In the plane z, in the medium behind the last surface.

Marray (N, 4, 4)

Only with details=True: transfer matrices of the kept beamlets.

gouyarray (N,)

Only with details=True: -arg(det(A + B Gamma))/2, unwrapped (for a round beamlet through a focus: -pi/2 per line focus).

keeparray (N,) of bool

Only with details=True: which input beamlets were kept.

pyoptics.raywave.field_to_rays(field, x=None, y=None, wavelength=None, n=1.0, z=0.0, step=1, threshold=0.001, phase=None, structured=False)

Rays along the local phase gradient (energy flux) of a field.

Every sample (of a sub grid with the given step) with |u|**2 above threshold times the maximum (and so are its four neighbours) becomes a ray at (x, y, z) with the direction (k_x, k_y, k_z)/(n*k_0), k_perp = grad(phase), the power |u|**2 times the area per ray and the optical path phase/k_0.

Parameters:

fieldarray (Ny, Nx) or analytic_phase.QuadraticPhaseField

Field. For a QuadraticPhaseField the analytic phase is differentiated exactly (x, y, wavelength and n are taken from it).

x, yarrays

Linear coordinates.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

znumber, optional

z of the plane (for the start points).

stepint, optional

Take every step-th sample along x and y.

thresholdnumber, optional

Relative intensity below which no rays are made.

phasearray, optional

Unwrapped phase; default: unwrapped from the field (for the optical path; the directions never need unwrapping).

structuredbool, optional

If True, keep the grid: the arrays have the shape (Ny’, Nx’, …) of the (sub sampled) grid, with NaN for the dropped samples, so that the rays (also after System.trace() and a reshape) can be turned into a field again with rays_to_field() (launch: the start points, irradiance: |u|**2).

Returns:

FieldRays

Can be traced with raytrace.System.trace(). Samples whose gradient exceeds n*k_0 (evanescent) are dropped.

pyoptics.raywave.local_spectrum_rays(field, x, y, wavelength, n=1.0, z=0.0, window=32, step=None, threshold=0.01, pad=2, jitter=None)

Rays from the local spectra of a field (sampled spectrogram / Husimi function).

The field is cut into Gaussian windowed patches (1/e**2 intensity radius window/4 samples) centred every step samples. Every component of the (zero padded) patch spectrum with a power above threshold times the largest one becomes a ray at the patch centre with the direction of its wave vector and the power of that component, scaled such that all rays together carry the power of the field. Several ray directions per point: crossing sheets of a caustic are resolved. The rays carry no phase (opl NaN).

Parameters:

fieldarray (Ny, Nx)

Sampled field.

x, yarrays

Linear coordinates.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

znumber, optional

z of the plane.

windowint, optional

Patch size in samples (even).

stepint, optional

Distance of the patch centres in samples (default window//4).

thresholdnumber, optional

Relative power below which spectral components are dropped.

padint, optional

Zero padding factor of the patch FFT (finer direction sampling).

jitternumpy.random.Generator or int, optional

If given (a generator or a seed), every ray is moved to a random point of its phase space cell (uniformly within the step in x and y and the spectral bin in k_x and k_y). Without, the rays sit on a lattice of positions and directions, which shows as moire patterns when the rays are binned in other planes.

Returns:

FieldRays

The power is normalized such that its sum equals sum(|u|**2)*dx*dy.

pyoptics.raywave.wigner(u, x, upsample=2)

Wigner distribution of a 1D field (e.g. a section of a 2D field).

W(x, k) = int u(x + s/2) conj(u(x - s/2)) exp(-i*k*s) ds, evaluated on the sample positions x and the wavenumbers k (the field is interpolated by Fourier zero padding first, by the factor upsample, so that the wavenumbers cover the full band of the input sampling without aliasing).

Its marginals are the intensity, int W dk/(2*pi) = |u(x)|**2, and the power spectrum; it is real and can be negative (interference between ray sheets).

Parameters:

uarray (N,)

Field.

xarray (N,)

Linear coordinates.

upsampleint, optional

Interpolation factor (>= 2).

Returns:

karray (N*upsample,)

Wavenumbers, ascending, from -pi/dx to pi/dx.

Warray (len(k), N)

Wigner distribution, W[i, j] at (x[j], k[i]).