Wavefront Aberrations with Zernike Polynomials¶

What you will learn:

What Zernike polynomials are and how they describe optical aberrations
How to inject WFE (Wavefront Error) into a telescope model
How individual and combined aberrations alter the PSF
How to visualise the wavefront OPD (Optical Path Difference) map

Prerequisites: 02_getting_started.ipynb

Background¶

Real telescopes are never perfect. Thermal gradients, mirror surface errors, and mechanical deformations introduce wavefront aberrations — deviations of the actual wavefront from a perfect sphere.

PSFCraft uses Zernike polynomials (Noll indexing) to describe these errors.
Each polynomial corresponds to a named optical aberration:

Noll index	Name
1	Piston
2	Tip (x-tilt)
3	Tilt (y-tilt)
4	Defocus
5	Oblique astigmatism
6	Vertical astigmatism
7	Vertical coma
8	Horizontal coma
9	Vertical trefoil
10	Oblique trefoil
11	Primary spherical

Coefficients are given in metres (RMS wavefront error).
A typical specification might be a few tens of nanometres RMS.

In [1]:

Copied!





%matplotlib inline
import psfcraft
from psfcraft.utils import Optical_aberration_names
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import psfcraft
from psfcraft.utils import Optical_aberration_names
import numpy as np
import matplotlib.pyplot as plt

pysynphot is not installed. Please install it to use PSF generation with polychromatic sources.

1. Visualising Individual Zernike Modes¶

We isolate each aberration by setting only one coefficient at a time to 100 nm RMS, and display the resulting OPD map alongside the PSF.

In [2]:

Copied!





import poppy

# Noll indices 2–6  (skip piston=1, which has no observable PSF effect)
noll_indices = [2, 3, 4, 5, 6]   # Tip, Tilt, Defocus, Oblique Astig, Vertical Astig
amplitude_nm = 100                # 100 nm RMS each
wavelength = 1e-6                 # 1 µm

fig, axes = plt.subplots(2, len(noll_indices), figsize=(15, 6))
plt.subplots_adjust(hspace=0.3, wspace=0.25)

for col, noll in enumerate(noll_indices):
    # Build coefficient vector: only the selected Noll index is non-zero
    coefs = np.zeros(15)
    coefs[noll - 1] = amplitude_nm * 1e-9   # convert nm → m

    tel = psfcraft.NewtonianTelescope(version="0", wfe_coefficients=coefs)
    tel.pixelscale = 0.05
    psf = tel.calc_psf(monochromatic=wavelength, fov_pixels=32)

    # Row 0: OPD map
    # Build a fresh optical system to access the ZernikeWFE plane (planes[1]).
    # tel.optsys (stored by calc_psf) omits hidden planes; _get_optical_system includes them.
    osys = tel._get_optical_system(fov_pixels=32)
    wfe_plane = osys.planes[1]   # ZernikeWFE element
    wf_tmp = poppy.Wavefront(wavelength=wavelength, npix=256, diam=2 * tel.primary_radius)
    # get_phasor returns the complex pupil; extract OPD from its phase
    phasor = wfe_plane.get_phasor(wf_tmp)
    opd = np.angle(phasor) * wavelength / (2 * np.pi)   # radians → metres

    im = axes[0, col].imshow(
        opd * 1e9,          # display in nm
        cmap="RdBu_r",
        origin="lower",
        vmin=-amplitude_nm, vmax=amplitude_nm,
    )
    axes[0, col].set_title(Optical_aberration_names[noll - 1], fontsize=10)
    axes[0, col].set_xticks([])
    axes[0, col].set_yticks([])

    # Row 1: PSF
    ax = psfcraft.display_psf(psf, ax=axes[1, col], title="", colorbar=False,
                              return_ax=True, crosshairs=True)
    ax.lines[0].set_color("white")
    ax.lines[1].set_color("white")

axes[0, 0].set_ylabel("WFE  (OPD map, nm)")
axes[1, 0].set_ylabel("PSF")
plt.suptitle(f"Individual Zernike modes  —  {amplitude_nm} nm RMS each", fontsize=13)
plt.show()
import poppy

# Noll indices 2–6  (skip piston=1, which has no observable PSF effect)
noll_indices = [2, 3, 4, 5, 6]   # Tip, Tilt, Defocus, Oblique Astig, Vertical Astig
amplitude_nm = 100                # 100 nm RMS each
wavelength = 1e-6                 # 1 µm

fig, axes = plt.subplots(2, len(noll_indices), figsize=(15, 6))
plt.subplots_adjust(hspace=0.3, wspace=0.25)

for col, noll in enumerate(noll_indices):
    # Build coefficient vector: only the selected Noll index is non-zero
    coefs = np.zeros(15)
    coefs[noll - 1] = amplitude_nm * 1e-9   # convert nm → m

    tel = psfcraft.NewtonianTelescope(version="0", wfe_coefficients=coefs)
    tel.pixelscale = 0.05
    psf = tel.calc_psf(monochromatic=wavelength, fov_pixels=32)

    # Row 0: OPD map
    # Build a fresh optical system to access the ZernikeWFE plane (planes[1]).
    # tel.optsys (stored by calc_psf) omits hidden planes; _get_optical_system includes them.
    osys = tel._get_optical_system(fov_pixels=32)
    wfe_plane = osys.planes[1]   # ZernikeWFE element
    wf_tmp = poppy.Wavefront(wavelength=wavelength, npix=256, diam=2 * tel.primary_radius)
    # get_phasor returns the complex pupil; extract OPD from its phase
    phasor = wfe_plane.get_phasor(wf_tmp)
    opd = np.angle(phasor) * wavelength / (2 * np.pi)   # radians → metres

    im = axes[0, col].imshow(
        opd * 1e9,          # display in nm
        cmap="RdBu_r",
        origin="lower",
        vmin=-amplitude_nm, vmax=amplitude_nm,
    )
    axes[0, col].set_title(Optical_aberration_names[noll - 1], fontsize=10)
    axes[0, col].set_xticks([])
    axes[0, col].set_yticks([])

    # Row 1: PSF
    ax = psfcraft.display_psf(psf, ax=axes[1, col], title="", colorbar=False,
                              return_ax=True, crosshairs=True)
    ax.lines[0].set_color("white")
    ax.lines[1].set_color("white")

axes[0, 0].set_ylabel("WFE  (OPD map, nm)")
axes[1, 0].set_ylabel("PSF")
plt.suptitle(f"Individual Zernike modes  —  {amplitude_nm} nm RMS each", fontsize=13)
plt.show()

No description has been provided for this image

Interpretation:

Tip / Tilt — shift the PSF centroid; they are equivalent to pointing errors.
Defocus — broadens the PSF into a ring structure.
Astigmatism — elongates the PSF along a preferred axis.

2. Combining Aberrations¶

In practice, many modes are present simultaneously.
Here we set a realistic WFE budget similar to space telescope requirements.

In [3]:

Copied!





# Realistic multi-mode WFE (Noll indexing, in nm)
# index:   1    2    3    4    5    6    7    8    9
wfe_nm = [ 0,   0,   0,  10,  50,  50, 100,  25,  25]
wfe_m  = [c * 1e-9 for c in wfe_nm]

tel = psfcraft.NewtonianTelescope(version="1_3", wfe_coefficients=wfe_m)
tel.pixelscale = 0.05

psf_aberrated = tel.calc_psf(monochromatic=1e-6, fov_pixels=64)

# Compare perfect vs aberrated
tel_perfect = psfcraft.NewtonianTelescope(version="1_3")
tel_perfect.pixelscale = 0.05
psf_perfect = tel_perfect.calc_psf(monochromatic=1e-6, fov_pixels=64)

# OPD map: build fresh optsys to get ZernikeWFE plane, extract via get_phasor
wavelength = 1e-6
osys = tel._get_optical_system(fov_pixels=64)
wfe_plane = osys.planes[1]   # ZernikeWFE element
wf_tmp = poppy.Wavefront(wavelength=wavelength, npix=256, diam=2 * tel.primary_radius)
phasor = wfe_plane.get_phasor(wf_tmp)
opd_map = np.angle(phasor) * wavelength / (2 * np.pi) * 1e9   # nm

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
psfcraft.display_psf(psf_perfect,   ax=axes[0], title="Perfect (no WFE)", colorbar=False)
psfcraft.display_psf(psf_aberrated, ax=axes[1], title="Aberrated",        colorbar=False)

im = axes[2].imshow(opd_map, cmap="RdBu_r", origin="lower")
axes[2].set_title("WFE map")
axes[2].set_xticks([])
axes[2].set_yticks([])
fig.colorbar(im, ax=axes[2], label="OPD (nm)")

plt.tight_layout()
plt.show()

# Print RMS WFE
rms_nm = np.sqrt(np.sum(np.array(wfe_nm[1:])**2))
print(f"Total RMS WFE (quadratic sum, modes 2–9): {rms_nm:.1f} nm")
# Realistic multi-mode WFE (Noll indexing, in nm)
# index:   1    2    3    4    5    6    7    8    9
wfe_nm = [ 0,   0,   0,  10,  50,  50, 100,  25,  25]
wfe_m  = [c * 1e-9 for c in wfe_nm]

tel = psfcraft.NewtonianTelescope(version="1_3", wfe_coefficients=wfe_m)
tel.pixelscale = 0.05

psf_aberrated = tel.calc_psf(monochromatic=1e-6, fov_pixels=64)

# Compare perfect vs aberrated
tel_perfect = psfcraft.NewtonianTelescope(version="1_3")
tel_perfect.pixelscale = 0.05
psf_perfect = tel_perfect.calc_psf(monochromatic=1e-6, fov_pixels=64)

# OPD map: build fresh optsys to get ZernikeWFE plane, extract via get_phasor
wavelength = 1e-6
osys = tel._get_optical_system(fov_pixels=64)
wfe_plane = osys.planes[1]   # ZernikeWFE element
wf_tmp = poppy.Wavefront(wavelength=wavelength, npix=256, diam=2 * tel.primary_radius)
phasor = wfe_plane.get_phasor(wf_tmp)
opd_map = np.angle(phasor) * wavelength / (2 * np.pi) * 1e9   # nm

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
psfcraft.display_psf(psf_perfect,   ax=axes[0], title="Perfect (no WFE)", colorbar=False)
psfcraft.display_psf(psf_aberrated, ax=axes[1], title="Aberrated",        colorbar=False)

im = axes[2].imshow(opd_map, cmap="RdBu_r", origin="lower")
axes[2].set_title("WFE map")
axes[2].set_xticks([])
axes[2].set_yticks([])
fig.colorbar(im, ax=axes[2], label="OPD (nm)")

plt.tight_layout()
plt.show()

# Print RMS WFE
rms_nm = np.sqrt(np.sum(np.array(wfe_nm[1:])**2))
print(f"Total RMS WFE (quadratic sum, modes 2–9): {rms_nm:.1f} nm")

Total RMS WFE (quadratic sum, modes 2–9): 127.9 nm

3. Sign Sensitivity of Aberrations¶

Some aberrations (e.g. astigmatism, coma) produce different PSF shapes for positive and negative coefficients. This matters for wavefront sensing applications where sign reconstruction is ambiguous.

In [4]:

Copied!





modes_to_test = {
    "Tip (2)":         2,
    "Astigmatism (5)": 5,
    "Coma (7)":        7,
}
amplitude_nm = 100
wavelength = 5e-7  # 500 nm — shorter λ makes aberrations more visible

fig, axes = plt.subplots(2, len(modes_to_test), figsize=(12, 6))
plt.subplots_adjust(hspace=0.35, wspace=0.2)

for col, (label, noll) in enumerate(modes_to_test.items()):
    for row, sign in enumerate([-1, +1]):
        coefs = np.zeros(15)
        coefs[noll - 1] = sign * amplitude_nm * 1e-9

        tel = psfcraft.NewtonianTelescope(version="0", wfe_coefficients=coefs)
        tel.pixelscale = 0.05
        psf = tel.calc_psf(monochromatic=wavelength, fov_pixels=32)

        sign_str = "−" if sign < 0 else "+"
        ax = psfcraft.display_psf(psf, ax=axes[row, col], colorbar=False,
                                   title=f"{label}  α {sign_str}", return_ax=True,
                                   crosshairs=True)
        ax.lines[0].set_color("white")
        ax.lines[1].set_color("white")

plt.suptitle(f"Sign sensitivity of Zernike aberrations  ({amplitude_nm} nm RMS)", fontsize=13)
plt.show()
modes_to_test = {
    "Tip (2)":         2,
    "Astigmatism (5)": 5,
    "Coma (7)":        7,
}
amplitude_nm = 100
wavelength = 5e-7  # 500 nm — shorter λ makes aberrations more visible

fig, axes = plt.subplots(2, len(modes_to_test), figsize=(12, 6))
plt.subplots_adjust(hspace=0.35, wspace=0.2)

for col, (label, noll) in enumerate(modes_to_test.items()):
    for row, sign in enumerate([-1, +1]):
        coefs = np.zeros(15)
        coefs[noll - 1] = sign * amplitude_nm * 1e-9

        tel = psfcraft.NewtonianTelescope(version="0", wfe_coefficients=coefs)
        tel.pixelscale = 0.05
        psf = tel.calc_psf(monochromatic=wavelength, fov_pixels=32)

        sign_str = "−" if sign < 0 else "+"
        ax = psfcraft.display_psf(psf, ax=axes[row, col], colorbar=False,
                                   title=f"{label}  α {sign_str}", return_ax=True,
                                   crosshairs=True)
        ax.lines[0].set_color("white")
        ax.lines[1].set_color("white")

plt.suptitle(f"Sign sensitivity of Zernike aberrations  ({amplitude_nm} nm RMS)", fontsize=13)
plt.show()

Key Takeaways¶

wfe_coefficients is a list of Zernike coefficients in metres, indexed from Noll 1.
Set a coefficient to amplitude_nm * 1e-9 to inject amplitude_nm nm RMS of that mode.
Tip and tilt shift the PSF; defocus broadens it; astigmatism and coma break symmetry.
The OPD map (planes[1].display(what='opd')) lets you verify the injected aberration visually.
Sign matters for non-symmetric modes (coma, trefoil, …) — relevant for wavefront sensing.

Next: 05_encircled_energy.ipynb — quantify PSF quality with Encircled Energy.