Ins and Outs of Polynomials

This document serves as a reference for how prysm is set up to work with polynomials, in the context of OPD or surface figure error. Much of what differentiates prysm’s API in this area has to do with the fact that it expects the grid to exist at the user level, but there are some deep and consequential implementation differences, too. Before we get into those, we will create a working grid and a mask for visualization:

import numpy as np

from prysm.coordinates import make_xy_grid, cart_to_polar
from prysm.geometry import circle

from matplotlib import pyplot as plt

x, y = make_xy_grid(256, diameter=2)
r, t = cart_to_polar(x, y)
mask = circle(1,r) == 0

This is a long document, so you may wish to search for your preferred polynomial flavor:

• Hopkins

• Zernike

• Jacobi

• Chebyshev

• Legendre

• Qs

Note that all polynomial types allow evaluation for arbitrary order.

Hopkins

The simplest polynomials are Hopkins’:

\[\text{OPD} = W_{abc} \left[\cos\left(a\cdot\theta\right) \cdot \rho^b \cdot H^c \right]\]

for some set of coefficients. The usage of this should not be surprising, for \(W_{131}\), coma one can write:

from prysm.polynomials import hopkins
cma = hopkins(1, 3, 1, r, t, 1)
cma[mask]=np.nan
plt.imshow(cma)

<matplotlib.image.AxesImage at 0x7f04103385d0>

Note that we defined our grid to have a radius of 1, but often you may hold two copies of r, one which is normalized by some reference radius for polynomial evaluation, and one which is not for pupil geometry evaluation. There is no further complexity in using Hopkins’ polynomials.

Zernike

prysm has a fairly granular implementation of Zernike polynomials, and expects its users to assemble the pieces to synthesize higher order functionality. The basic building block is the zernike_nm function, which takes azimuthal and radial orders n and m, as in \(Z_n^m\). For example, to compute the equivalent “primary coma” Zernike mode as the hopkins example, one would:

from prysm.polynomials import zernike_nm
cmaZ = zernike_nm(3,1, r,t, norm=True)
cmaZ[mask]=np.nan
plt.imshow(cmaZ)

<matplotlib.image.AxesImage at 0x7f040e226ed0>

Note that the terms can be orthonormalized (given unit RMS) or not, based on the norm kwarg. The order m can be negative to give access to the sinusoidal terms instead of cosinusoidal. If you wish to work with a particular ordering scheme, prysm supports Fringe, Noll, and ANSI out of the box, all of which start counting at 1.

from prysm.polynomials import noll_to_nm, fringe_to_nm, ansi_j_to_nm

n, m = fringe_to_nm(9)
sphZ = zernike_nm(n, m, r, t, norm=False)
sphZ[mask]=np.nan
plt.imshow(sphZ)

<matplotlib.image.AxesImage at 0x7f040e1b6990>

These functions are not iterator-aware and should be used with, say, a list comprehension.

If you wish to compute Zernikes much more quickly, the underlying implementation in prysm allows the work in computing lower order terms to be used to compute the higher order terms. The Zernike polynomials are fundamentally two “pieces” which get multiplied. The radial basis is where much of the work lives, and most programs that do not type out closed form solutions use Rodrigues’ technique to compute the radial basis:

\[R_n^m (\rho) = \sum_{k=0}^{\frac{n-m}{2}} \frac{(-1)^k (n-k)!}{k!(\frac{n+m}{2}-k)!(\frac{n-m}{2}-k)!}\rho^{n-2k} \tag{1}\]

prysm does not do this, and instead uses the fact that the radial polynomial is a Jacobi polynomial under a change-of-basis:

\[R_n^m (\rho) = P_\frac{n-m}{2}^\left(0,|m|\right)\left(2\rho^2 - 1\right) \tag{2}\]

And the jacobi polynomials can be computed using a recurrence relation:

\[a \cdot P_n^{(\alpha,\beta)} = b \cdot x \cdot P_{n-1}^{(\alpha,\beta)} - c \cdot P_{n-2}^{(\alpha,\beta)} \tag{3}\]

In other words, for a given \(m\), you can compute \(R\) for \(n=3\) from \(R\) for \(n=2\) and \(n=1\), and so on until you reach the highest value of N. Because the sum in the Rodrigues formulation is increasingly large as \(n,m\) grow, it has worse than linear time complexity. Because the recurrrence in Eq. (3) does not change as \(n,m\) grow it does have linear time complexity.

The use of this recurrence relation is hidden from the user in the zernike_nm function, and the recurrence relation is for a so-called auxiliary polynomial (\(R\)), so the Zernike polynomials themselves are not useful for this recurrence. You can make use of it by calling the zernike_nm_sequence function, a naming that will become familiar by the end of this reference guide. Consider the first 16 Fringe Zernikes:

from prysm.polynomials import zernike_nm_sequence

nms = [fringe_to_nm(i) for i in range(1,36)]

# zernike_nm_sequence returns a generator
%timeit polynomials = list(zernike_nm_sequence(nms, r, t)) # implicit norm=True

44.2 ms ± 118 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Compare the timing to not using the sequence flavored version:

%%timeit
for n, m in nms:
    zernike_nm(n, m, r, t)

144 ms ± 250 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

The sequence function returns a generator to leave the user in control of their memory usage. If you wished to compute 1,000 Zernike polynomials, this would avoid holding them all in memory at once while still improving performance. These is no benefit other than performance and plausibly reduced memory usage to the _sequence version of the function. A side benefit to the recurrence relation is that it is numerically stable to higher order than Rodrigues’ expression, so you can compute higher order Zernike polynomials without numerical errors. This is an especially useful property for using lower-precision data types like float32, since they will suffer from numerical imprecision earlier.

Jacobi

Of course, because the Zernike polynomials are related to them you also have access to the Jacobi polynomials:

from prysm.polynomials import jacobi, jacobi_sequence

x_ = x[0,:] # not required to be 1D, just for example
plt.plot(x_, jacobi(3,0,0,x_))

[<matplotlib.lines.Line2D at 0x7f040e0f9750>]

This shape may be familiar as the Zernike flavor of coma across one axis.

Chebyshev

Both types of Chevyshev polynomials are supported. They are both just special cases of Jacobi polynomials:

\[\text{cheby1} \equiv P_n^\left(-0.5,-0.5\right)(x) \quad / \quad P_n^\left(-0.5,-0.5\right)(1)\]

\[\text{cheby2} \equiv P_n^\left(0.5,0.5\right)(x) \quad / \quad P_n^\left(0.5,0.5\right)(1)\]

from prysm.polynomials import cheby1, cheby2, cheby1_sequence

plt.plot(x_, cheby1(3,x_), x_, cheby2(3,x_))
plt.legend(['first kind', 'second kind'])

<matplotlib.legend.Legend at 0x7f041156f850>

The most typical use of these polynomials in optics are as an orthogonal basis over some rectangular aperture. The calculation is separable in x and y, so it can be reduced from scaling by \(2(N\cdot M)\) to just \(N+M\). prysm will compute the mode for one column of x and one row of y, then broadcast to 2D to assemble the ‘image’. This introduces three new functions:

from prysm.polynomials import cheby1_2d_sequence, mode_1d_to_2d, sum_of_xy_modes # or cheby2_...

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
/tmp/ipykernel_367/3759324147.py in <module>
----> 1 from prysm.polynomials import cheby1_2d_sequence, mode_1d_to_2d, sum_of_xy_modes # or cheby2_...

ImportError: cannot import name 'cheby1_2d_sequence' from 'prysm.polynomials' (/home/docs/checkouts/readthedocs.org/user_builds/prysm/envs/v0.20/lib/python3.7/site-packages/prysm/polynomials/__init__.py)

# orders 1, 2, 3 in x and 4, 5, 6 in y
ns = [1, 2, 3]
ms = [4, 5, 6]
modesx, modesy = cheby1_2d_sequence(ns, ms, x, y)
plt.plot(x_, modesx[0]) # modes are 1D
plt.title('a single mode, 1D')
plt.figure()
# and can be expanded to 2D
plt.imshow(mode_1d_to_2d(modesx[0], x, y))
plt.title('a single mode, 2D')

Wx = [1]*len(modesx)
Wy = [1]*len(modesy)
im = sum_of_xy_modes(modesx, modesy, x, y, Wx, Wy)
plt.imshow(im)
plt.title('a sum of modes')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_367/2418987938.py in <module>
      2 ns = [1, 2, 3]
      3 ms = [4, 5, 6]
----> 4 modesx, modesy = cheby1_2d_sequence(ns, ms, x, y)
      5 plt.plot(x_, modesx[0]) # modes are 1D
      6 plt.title('a single mode, 1D')

NameError: name 'cheby1_2d_sequence' is not defined

As a final note, there is no reason you can’t just use the cheby1/cheby2 functions with 2D arrays, it is only slower:

im = cheby2(3, y)
plt.imshow(im)

<matplotlib.image.AxesImage at 0x7f040decde90>

Legendre

These polynomials are just a special case of Jacobi polynomials:

\[\text{legendre} \equiv P_n^\left(0,0\right)(x)\]

Usage follows from the Chebyshev exactly, except the functions are prefixed by legendre.

from prysm.polynomials import legendre, legendre_sequence, legendre_2d_sequence

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
/tmp/ipykernel_367/2660797986.py in <module>
----> 1 from prysm.polynomials import legendre, legendre_sequence, legendre_2d_sequence

ImportError: cannot import name 'legendre_2d_sequence' from 'prysm.polynomials' (/home/docs/checkouts/readthedocs.org/user_builds/prysm/envs/v0.20/lib/python3.7/site-packages/prysm/polynomials/__init__.py)

Qs

Qs are Greg Forbes’ Q polynomials, \(Q\text{bfs}\), \(Q\text{con}\), and \(Q_n^m\). Qbfs and Qcon polynomials are radial only, and replace the ‘standard’ asphere equation. The implementation of all three of these also uses a recurrence relation, although it is more complicated and outside the scope of this reference guide. Each includes the leading prefix from the papers:

• \(\rho^2(1-\rho^2)\) for \(Q\text{bfs}\),

• \(\rho^4\) for \(Q\text{con}\),

• the same as \(Q\text{bfs}\) for \(Q_n^m\) when \(m=0\) or \(\rho^m \cos\left(m\theta\right)\) for \(m\neq 0\)

The \(Q_n^m\) implementation departs from the papers in order to have a more Zernike-esque flavor. Instead of having \(a,b\) coefficients and \(a\) map to \(\cos\) and \(b\) to \(\sin\), this implementation uses the sign of \(m\), with \(\cos\) for \(m>0\) and \(\sin\) for \(m<0\).

There are six essential functions:

from prysm.polynomials import (
    Qbfs, Qbfs_sequence,
    Qcon, Qcon_sequence,
    Q2d, Q2d_sequence,
)

p = Qbfs(2,r)
p[mask]=np.nan
plt.imshow(p)

<matplotlib.image.AxesImage at 0x7f040de9d250>

p = Qcon(2,r)
p[mask]=np.nan
plt.imshow(p)

<matplotlib.image.AxesImage at 0x7f040ddc44d0>

p = Q2d(2, 2, r, t) # cosine term
p[mask]=np.nan
plt.imshow(p)

<matplotlib.image.AxesImage at 0x7f040dd95a50>

p2 = Q2d(2, -2, r, t) # sine term
p2[mask]=np.nan
plt.imshow(p2)

<matplotlib.image.AxesImage at 0x7f040dcc6490>