pyoptics.abcd

ABCD (ray transfer) matrices for paraxial ray optics and Gaussian beams.

Rays

A paraxial ray in a meridional plane is described by its height y and its reduced angle u = n*theta, where theta is the angle to the optical axis and n the refractive index of the medium the ray is in. An optical system maps the ray vector (y, u) in its input plane onto the output plane by a real 2x2 matrix

[y_out]   [A  B] [y_in]
[u_out] = [C  D] [u_in].

With reduced angles, every lossless system has det M = AD - BC = 1, and a flat interface is the identity. In air, u is simply the angle. Lengths inside a medium of index n enter as reduced distances d/n.

Matrices are written for the light travelling from left to right (+z). Matrices of consecutive elements multiply from the right to the left: a system made of M_1 followed by M_2 has M = M_2 @ M_1 (see compose()). Mirrors are unfolded, i.e. light continues in +z behind them.

Sign conventions

• Radii of curvature R > 0 if the centre of curvature lies to the right of the surface (in +z). A biconvex lens thus has R1 > 0 and R2 < 0.

• A converging thin lens has f > 0.

• Mirrors: R > 0 for a concave (focusing) mirror, which acts like a thin lens with f = R/2.

• Ray angles are positive for rays going upwards (+y) while travelling in +z.

Gaussian beams

A fundamental Gaussian beam is described by its complex beam parameter q, with

1/q = 1/R - i*lambda_0/(pi*n*w**2),   q = z + i*z_R at distance z
                                       behind the waist,

where w is the 1/e**2 intensity radius, R the radius of curvature of the wavefront, lambda_0 the vacuum wavelength and z_R = pi*n*w_0**2/lambda_0 the Rayleigh length. This is the q of most textbooks (Siegman, Kogelnik & Li). The ABCD law holds for the reduced parameter q/n:

q_out/n_out = (A*q_in/n_in + B)/(C*q_in/n_in + D).

All functions of this module take and return the true q (not q/n) and the refractive index separately.

In the exp(i(kz - omega*t)) convention of pyoptics, the transverse field of the beam is proportional to exp(+i*k*r**2/(2*R))*exp(-r**2/w**2) (see gaussian_field()), which is exp(i*k*r**2/(2*conj(q))). The Gouy phase psi enters as exp(-i*psi), as in pyoptics.beams. Wavelengths are always vacuum wavelengths.

Contents

• element matrices: free_space(), thin_lens(), refraction(), mirror(), thick_lens(), grin(), compose()

• OpticalSystem: builds a system element by element and traces rays and Gaussian beams through it (also inside the elements), computes cardinal points and images and draws a schematic layout

• Gaussian beams: q_from_waist(), q_from_w_R(), beam_radius(), curvature_radius(), waist(), transform_q(), gouy_phase(), gaussian_field(), coupling_efficiency()

• resonators: stability_parameter(), is_stable(), eigenmode(), round_trip_gouy_phase(), g_parameters()

• imaging: cardinal_points(), image()

• wave optics: collins(), the Collins integral (Fresnel diffraction through any ABCD system)

References

H. Kogelnik and T. Li, “Laser beams and resonators”, Appl. Opt. 5, 1550-1567 (1966).

1. 5. Siegman, “Lasers”, University Science Books (1986).

S. A. Collins, “Lens-system diffraction integral written in terms of matrix optics”, J. Opt. Soc. Am. 60, 1168-1177 (1970).

B. E. A. Saleh and M. C. Teich, “Fundamentals of Photonics”, 2nd ed., Wiley (2007).

Classes

BeamTrace	Gaussian beam traced through an OpticalSystem: beam parameter q, radius w, wavefront radius R, accumulated Gouy phase (unwrapped) and refractive index n at the positions z.
CardinalPoints	Cardinal points of an optical system.
OpticalSystem	Paraxial optical system built element by element along the z axis.

Functions

beam_radius	Beam radius of a Gaussian beam with parameter q.
cardinal_points	Cardinal points of a system.
collins	Propagate a field through a paraxial ABCD system (Collins integral).
compose	System matrix of elements given in the order the light passes them.
coupling_efficiency	Power coupling efficiency of two coaxial fundamental Gaussian beams.
curvature_radius	Wavefront radius of curvature R = 1/Re(1/q), np.inf at a waist.
eigenmode	Self-reproducing beam parameter of a round trip matrix M.
free_space	Propagation over a distance d in a homogeneous medium of index n.
g_parameters	Resonator g parameters of a two mirror resonator.
gaussian_field	Sample a fundamental Gaussian beam with beam parameter q.
gouy_phase	Gouy phase psi accumulated by a fundamental Gaussian beam in a system.
grin	Section of length d of a graded index medium n(r) = n0*(1 - g**2*r**2/2).
image	Image of an object at distance s in front of the input plane.
is_stable	True if the round trip matrix M describes a stable resonator.
mirror	Spherical mirror (unfolded), R > 0 for a concave (focusing) mirror.
q_from_w_R	Complex beam parameter from beam radius and wavefront radius.
q_from_waist	Complex beam parameter at distance z behind a waist.
rayleigh_length	Rayleigh length z_R = Im(q).
refraction	Refraction at a spherical interface from index n1 to index n2.
round_trip_gouy_phase	Gouy phase of the eigenmode per round trip (rotationally symmetric resonator; per transverse axis it is half of that): arccos((A + D)/2) for B > 0 and 2*pi - arccos((A + D)/2) for B < 0, i.e. in [0, 2*pi).
stability_parameter	m = (A + D)/2 of a round trip matrix; a resonator is stable for |m| < 1.
thick_lens	Thick lens in a medium of index n (on both sides).
thin_lens	Thin lens of focal length f (f > 0 converging, in air).
transform_q	Transform a beam parameter with the ABCD law.
waist	Waist radius w0 and the waist position relative to the plane of q.

pyoptics.abcd.free_space(d, n=1.0)

Propagation over a distance d in a homogeneous medium of index n.

Parameters:

dnumber

Geometrical distance.

nnumber, optional

Refractive index of the medium.

Returns:

Marray

[[1, d/n], [0, 1]] (reduced distance d/n).

pyoptics.abcd.thin_lens(f)

Thin lens of focal length f (f > 0 converging, in air).

Parameters:

fnumber

Focal length. np.inf gives the identity.

Returns:

Marray

[[1, 0], [-1/f, 1]].

pyoptics.abcd.refraction(R, n1, n2)

Refraction at a spherical interface from index n1 to index n2.

Parameters:

Rnumber

Radius of curvature, R > 0 if the centre of curvature lies behind the interface (in +z). np.inf for a flat interface.

n1, n2numbers

Refractive indices before and behind the interface.

Returns:

Marray

[[1, 0], [-(n2 - n1)/R, 1]]. A flat interface is the identity in reduced angles (Snell’s law n1*theta1 = n2*theta2).

pyoptics.abcd.mirror(R, theta=0.0, plane='tangential', n=1.0)

Spherical mirror (unfolded), R > 0 for a concave (focusing) mirror.

Parameters:

Rnumber

Radius of curvature. np.inf for a plane mirror.

thetanumber, optional

Angle of incidence in radians. Under oblique incidence, the effective radius is R*cos(theta) in the tangential plane (the plane of incidence) and R/cos(theta) in the sagittal plane.

plane{“tangential”, “sagittal”}, optional

nnumber, optional

Refractive index of the medium in front of the mirror.

Returns:

Marray

[[1, 0], [-2*n/R_eff, 1]], a thin lens with focal length R_eff/2 in air.

pyoptics.abcd.grin(d, n0, g)

Section of length d of a graded index medium n(r) = n0*(1 - g**2*r**2/2).

Rays oscillate with the period 2*pi/g (pitch). For g -> 0 this becomes free_space(d, n0).

Parameters:

dnumber

Length of the section.

n0number

Refractive index on the axis.

gnumber

Gradient constant (1/length). g = 0 is a homogeneous medium.

Returns:

Marray

[[cos(g*d), sin(g*d)/(n0*g)], [-n0*g*sin(g*d), cos(g*d)]].

pyoptics.abcd.compose(*matrices)

System matrix of elements given in the order the light passes them.

Parameters:

*matricesarrays

2x2 matrices, first element first.

Returns:

Marray

compose(M1, M2, M3) = M3 @ M2 @ M1; the identity for no arguments.

pyoptics.abcd.thick_lens(R1, R2, d, n_lens, n=1.0)

Thick lens in a medium of index n (on both sides).

The matrix maps the vertex of the first surface onto the vertex of the second surface.

Parameters:

R1, R2numbers

Radii of curvature of the first and second surface. For a biconvex lens R1 > 0 and R2 < 0; np.inf for a flat surface.

dnumber

Centre thickness.

n_lensnumber

Refractive index of the lens.

nnumber, optional

Refractive index of the surrounding medium.

Returns:

Marray

The 2x2 matrix of refraction, glass and refraction.

pyoptics.abcd.q_from_waist(w0, z, wavelength, n=1.0)

Complex beam parameter at distance z behind a waist.

Parameters:

w0number or array

Waist radius (1/e**2 intensity).

znumber or array

Distance from the waist; z < 0 lies in front of the waist.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

Returns:

qcomplex or array

q = z + i*z_R with the Rayleigh length z_R = pi*n*w0**2/wavelength.

pyoptics.abcd.q_from_w_R(w, R, wavelength, n=1.0)

Complex beam parameter from beam radius and wavefront radius.

Parameters:

wnumber or array

Beam radius (1/e**2 intensity).

Rnumber or array

Wavefront radius of curvature: np.inf for a flat wavefront, R > 0 for a diverging beam.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

Returns:

qcomplex or array

1/q = 1/R - i*wavelength/(pi*n*w**2).

pyoptics.abcd.beam_radius(q, wavelength, n=1.0)

Beam radius of a Gaussian beam with parameter q.

Parameters:

qcomplex or array

Complex beam parameter.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

Returns:

wnumber or array

1/e**2 intensity radius w = sqrt(-wavelength/(pi*n*Im(1/q))).

pyoptics.abcd.curvature_radius(q)

Wavefront radius of curvature R = 1/Re(1/q), np.inf at a waist.

R > 0 for a diverging beam (behind the waist).

pyoptics.abcd.rayleigh_length(q)

Rayleigh length z_R = Im(q).

pyoptics.abcd.waist(q, wavelength, n=1.0)

Waist radius w0 and the waist position relative to the plane of q.

Parameters:

qcomplex or array

Complex beam parameter.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

Returns:

w0number

Waist radius.

z_waistnumber

Position of the waist relative to the current plane, -Re(q): positive if the beam is converging towards a waist that lies ahead (in +z).

pyoptics.abcd.transform_q(M, q, n_in=1.0, n_out=1.0)

Transform a beam parameter with the ABCD law.

q_out/n_out = (A*q/n_in + B)/(C*q/n_in + D).

Parameters:

Marray

2x2 matrix of the system.

qcomplex or array

Beam parameter in the input plane.

n_in, n_outnumbers, optional

Refractive indices in the input and output planes.

Returns:

q_outcomplex or array

Beam parameter in the output plane.

pyoptics.abcd.gouy_phase(M, q, n_in=1.0)

Gouy phase psi accumulated by a fundamental Gaussian beam in a system.

psi = -arg(A + B*n_in/q), wrapped to (-pi, pi]. The field is multiplied by exp(-i*psi)/|A + B*n_in/q| (rotationally symmetric beam in a rotationally symmetric system). In one transverse dimension (or for each axis of an astigmatic system) the phase is psi/2; a Hermite-Gaussian HG_lm acquires (l + 1/2)*psi_x + (m + 1/2)*psi_y.

For free space behind a waist, psi = arctan(z/z_R), consistent with pyoptics.beams. To follow psi beyond +-pi, sample the system densely and unwrap, as OpticalSystem.trace_beam() does.

Parameters:

Marray

2x2 matrix of the system.

qcomplex or array

Beam parameter in the input plane.

n_innumber, optional

Refractive index in the input plane.

Returns:

psinumber or array

Gouy phase in radians.

pyoptics.abcd.gaussian_field(x, y, q, wavelength, n=1.0, q_y=None, x0=0.0, y0=0.0)

Sample a fundamental Gaussian beam with beam parameter q.

Parameters:

x, yarrays

qcomplex

Beam parameter (in x, and in y unless q_y is given).

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

q_ycomplex, optional

Beam parameter in y for simple astigmatic beams.

x0, y0numbers, optional

Centre of the beam.

Returns:

Earray

sqrt(2/(pi*w_x*w_y))*exp(i*k*(x**2/(2*R_x) + y**2/(2*R_y))) *exp(-x**2/w_x**2 - y**2/w_y**2), normalized to unit power. The Gouy phase is not included (see gouy_phase()).

pyoptics.abcd.coupling_efficiency(q1, q2, wavelength, n=1.0)

Power coupling efficiency of two coaxial fundamental Gaussian beams.

Both beams are rotationally symmetric and given in the same plane

eta = |<E1|E2>|**2 = 4/((w1/w2 + w2/w1)**2
                        + (pi*n*w1*w2/wavelength)**2*(1/R1 - 1/R2)**2).

This is the mode matching efficiency into a fibre or a resonator mode. In one transverse dimension the efficiency is sqrt(eta).

Parameters:

q1, q2complex or arrays

Beam parameters of the two beams.

wavelengthnumber

Vacuum wavelength.

nnumber, optional

Refractive index.

Returns:

etanumber or array

Coupling efficiency between 0 and 1.

pyoptics.abcd.stability_parameter(M)

m = (A + D)/2 of a round trip matrix; a resonator is stable for |m| < 1.

pyoptics.abcd.is_stable(M)

True if the round trip matrix M describes a stable resonator.

pyoptics.abcd.eigenmode(M, n=1.0)

Self-reproducing beam parameter of a round trip matrix M.

Solves q = (A*q + B)/(C*q + D) for a confined beam (Im q > 0):

n/q = (D - A)/(2*B) - i*sqrt(1 - m**2)/|B|,   m = (A + D)/2.

Combine with beam_radius() (with the wavelength) to get the mode size. Raises ValueError for an unstable resonator (|m| >= 1).

Parameters:

Marray

Round trip matrix starting and ending in the reference plane.

nnumber, optional

Refractive index in the reference plane.

pyoptics.abcd.round_trip_gouy_phase(M)

Gouy phase of the eigenmode per round trip (rotationally symmetric resonator; per transverse axis it is half of that): arccos((A + D)/2) for B > 0 and 2*pi - arccos((A + D)/2) for B < 0, i.e. in [0, 2*pi). It sets the spacing of the transverse mode frequencies, Delta nu = (Delta psi/(2*pi))*FSR per unit of mode order.

pyoptics.abcd.g_parameters(L, R1, R2)

Resonator g parameters of a two mirror resonator.

Parameters:

Lnumber

Resonator length.

R1, R2numbers

Radii of curvature of the mirrors, R > 0 for concave mirrors; np.inf for a plane mirror.

Returns:

g1, g2numbers

g_i = 1 - L/R_i. The resonator is stable for 0 < g1*g2 < 1.

class pyoptics.abcd.CardinalPoints(f, f_prime, F, F_prime, H, H_prime, N, N_prime, ffl, bfl)

Bases: tuple

Cardinal points of an optical system.

f, f_primefront and rear effective focal length (f = n_in/P,

f_prime = n_out/P with the power P = -C)

F, F_prime, H, H_prime, N, N_primez positions of the focal, principal and

nodal points, measured from the input plane (z = 0); the output plane is at z = length

ffl, bflfront and back focal lengths, i.e. distances of F in front of the

input plane and of F_prime behind the output plane

F

Alias for field number 2

F_prime

Alias for field number 3

H

Alias for field number 4

H_prime

Alias for field number 5

N

Alias for field number 6

N_prime

Alias for field number 7

bfl

Alias for field number 9

f

Alias for field number 0

f_prime

Alias for field number 1

ffl

Alias for field number 8

pyoptics.abcd.cardinal_points(M, length=0.0, n_in=1.0, n_out=1.0)

Cardinal points of a system.

Parameters:

Marray

2x2 matrix from the input plane (z = 0) to the output plane.

lengthnumber, optional

z position of the output plane.

n_in, n_outnumbers, optional

Refractive indices in front of and behind the system.

Returns:

CardinalPoints

Focal lengths and positions of the focal, principal and nodal points.

Raises:

ValueError

For an afocal system (C = 0).

pyoptics.abcd.image(M, s, n_in=1.0, n_out=1.0)

Image of an object at distance s in front of the input plane.

Parameters:

Marray

2x2 matrix of the system.

snumber

Object distance in front of the input plane (s > 0: real object).

n_in, n_outnumbers, optional

Refractive indices in front of and behind the system.

Returns:

s_primenumber

Distance of the image behind the output plane (negative: virtual image in front of it). np.inf if the image lies at infinity.

mnumber

Lateral magnification.

class pyoptics.abcd.BeamTrace(z, q, w, R, gouy, n)

Bases: tuple

Gaussian beam traced through an OpticalSystem: beam parameter q, radius w, wavefront radius R, accumulated Gouy phase (unwrapped) and refractive index n at the positions z.

R

Alias for field number 3

gouy

Alias for field number 4

n

Alias for field number 5

q

Alias for field number 1

w

Alias for field number 2

z

Alias for field number 0

class pyoptics.abcd.OpticalSystem(n=1.0)

Bases: object

Paraxial optical system built element by element along the z axis.

The builder methods return the system itself and can be chained:

>>> system = OpticalSystem().free_space(100).thin_lens(50).free_space(100)
>>> system.matrix   # 2f-2f imaging: [[-1, 0], [-0.02, -1]]

Rays and beams can be evaluated at any z, also inside thick elements. Positions z < 0 lie in the input medium in front of the system and z > length in the output medium behind it. Thin elements act at their z position; quantities evaluated exactly there refer to the plane just behind the element.

Parameters:

nnumber, optional

Refractive index of the input medium.

property length

Axial length of the system.

property n_out

Refractive index behind the last element.

free_space(d, n=None, label=None)

Append a homogeneous medium.

A change of the refractive index implies a flat interface.

Parameters:

dnumber

Length.

nnumber, optional

Refractive index. Defaults to the current medium.

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

thin_lens(f, label=None)

Append a thin lens.

Parameters:

fnumber

Focal length, f > 0 for a converging lens.

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

refraction(R, n2, label=None)

Append a spherical or flat interface.

Parameters:

Rnumber

Radius of curvature (see refraction()); np.inf for a flat interface.

n2number

Refractive index behind the interface.

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

mirror(R, theta=0.0, plane='tangential', label=None)

Append a spherical mirror (unfolded).

Parameters:

Rnumber

Radius of curvature, R > 0 for a concave mirror.

thetanumber, optional

Angle of incidence in radians.

plane{“tangential”, “sagittal”}, optional

Plane for oblique incidence, see mirror().

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

thick_lens(R1, R2, d, n_lens, label=None)

Append a thick lens surrounded by the current medium.

The lens is added as two refracting surfaces and the glass in between.

Parameters:

R1, R2numbers

Radii of curvature of the two surfaces, see thick_lens().

dnumber

Centre thickness.

n_lensnumber

Refractive index of the lens.

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

grin(d, n0, g, label=None)

Append a graded index section n(r) = n0*(1 - g**2*r**2/2).

Parameters:

dnumber

Length of the section.

n0number

Refractive index on the axis.

gnumber

Gradient constant (1/length).

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

element(M, length=0.0, n=None, label=None)

Append any element given by its matrix.

Inside an element of nonzero length, rays and beams are only available at its ends.

Parameters:

Marray

Full 2x2 matrix of the element.

lengthnumber, optional

Axial length of the element.

nnumber, optional

Refractive index behind the element. Defaults to the current medium.

labelstr, optional

Label for draw().

Returns:

OpticalSystem

The system itself (for chaining).

append(other)

Append all elements of another OpticalSystem.

property matrix

System matrix from z = 0 to z = length.

n_at(z)

Refractive index at the position z.

matrix_at(z)

Matrix from z = 0 to the position z (just behind thin elements at z).

trace(rays, z)

Trace rays through the system.

Parameters:

raysarray

Ray(s) (y, theta) in the plane z = 0 (or (2, N) for N rays), theta being the true angle in the input medium.

zarray

Positions at which to evaluate the rays.

Returns:

y, thetaarrays

Height and true angle, shape (len(z), N) (or (len(z),) for a single ray).

trace_beam(q, wavelength, z)

Trace a Gaussian beam through the system.

Parameters:

qcomplex

Beam parameter in the plane z = 0.

wavelengthnumber

Vacuum wavelength.

zarray

Positions at which to evaluate the beam.

Returns:

BeamTrace

q, w, R, n and the Gouy phase at the positions z (relative to z = 0, unwrapped along z, so z should be sampled densely enough for the phase to change by less than pi between samples).

cardinal_points()

Cardinal points, see cardinal_points().

image(s)

Image of an object at distance s in front of the system, see image().

draw(ax=None, height=1.0, color='0.3')

Draw a schematic layout of the system into a matplotlib axis.

Lenses are drawn as double arrows (pointing outwards for converging, inwards for diverging lenses), mirrors as thick lines, interfaces as arcs and media with n != 1 as shaded areas.

Parameters:

axmatplotlib axis, optional

Axis to draw into. Defaults to the current axis.

heightnumber, optional

Half height of the symbols in data units.

colorcolor, optional

Color of the symbols.

Returns:

axmatplotlib axis

pyoptics.abcd.collins(field, x, y, M, wavelength, M_y=None)

Propagate a field through a paraxial ABCD system (Collins integral).

The generalized Fresnel integral of S. A. Collins (JOSA 60, 1168 (1970)),

E_out(x) = (i*lambda*B)**(-1/2) * integral E_in(x')
           * exp(i*pi/(lambda*B)*(A*x'**2 - 2*x'*x + D*x**2)) dx'

per transverse axis (with the vacuum wavelength and reduced B), is evaluated by splitting M into a thin lens, free space B and another thin lens; the free space part uses the Fresnel transfer function. The output is sampled on the input grid, so the grid has to hold the output field and resolve both lens chirps (wrap-around is not prevented). For B = 0 (imaging planes) the field is scaled by interpolation. The constant phase exp(i*k*L) of the axial optical path is omitted.

Parameters:

fieldarray

x, yarrays

Centered coordinates (linear or 2d).

Marray

ABCD matrix (in x, and in y unless M_y is given) with det M = 1.

wavelengthnumber

Vacuum wavelength.

M_yarray, optional

Matrix in y for systems with cylindrical elements or tilted mirrors.

Returns:

Earray

Field in the output plane on the grid x, y.