.. image:: https://mybinder.org/badge_logo.svg
   :target: https://mybinder.org/v2/gh/ehpor/hcipy/HEAD?labpath=doc/tutorial_notebooks/CoreConceptsWavefrontsOpticalElements/CoreConceptsWavefrontsOpticalElements.ipynb
   :alt: Open in Binder

Core Concepts: Wavefronts and Optical Elements
==============================================

We’ll explore the ``Wavefront`` object and how optical elements
transform wavefronts through amplitude and phase modifications. We’ll
propagate wavefronts in different regimes and build a simple optical
system.

.. code:: ipython3

    import hcipy
    import numpy as np
    import matplotlib.pyplot as plt

A ``Wavefront`` stores the complex electric field on a grid, combining
both amplitude and phase information. Let’s start by setting up a
computational grid and defining a circular pupil aperture — these are
the basic building blocks we’ll use throughout.

.. code:: ipython3

    # Create grid and aperture
    pupil_diameter = 1
    pupil_grid = hcipy.make_pupil_grid(128, pupil_diameter * 1.2)
    aperture = hcipy.make_circular_aperture(1.0)(pupil_grid)

Now we’ll build a planar wavefront from this aperture. Looking at the
amplitude, phase, and intensity, we can see a perfect flat wavefront:
the amplitude is uniform inside the pupil and zero outside, the phase is
constant everywhere, and the intensity follows the same circular
footprint as the amplitude.

.. code:: ipython3

    # Create wavefront and inspect properties
    wavelength = 500e-9
    wf = hcipy.Wavefront(aperture, wavelength)
    
    print(f"Wavelength: {wf.wavelength*1e9:.0f} nm")
    print(f"Total power: {wf.total_power:.4f}")
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    hcipy.imshow_field(wf.amplitude, ax=axes[0])
    axes[0].set_title('Amplitude')
    hcipy.imshow_field(wf.phase, cmap='RdBu', ax=axes[1])
    axes[1].set_title('Phase')
    hcipy.imshow_field(wf.intensity, ax=axes[2])
    axes[2].set_title('Intensity')
    plt.show()


.. parsed-literal::

    Wavelength: 500 nm
    Total power: 0.7840



.. image:: output_5_1.png


Real optical systems introduce phase aberrations that distort the
wavefront. We’ll use Zernike polynomials to add two common aberrations:
defocus and astigmatism. Defocus produces a parabolic phase map (bright
center fading to dark edges), while astigmatism shows a distinctive
saddle-shaped pattern with two orthogonal lobes.

.. code:: ipython3

    # Defocus
    defocus_phase = hcipy.zernike(2, 0)(pupil_grid)
    wf_defocus = hcipy.Wavefront(aperture * np.exp(1j * defocus_phase), wavelength)
    
    # Astigmatism
    astig_phase = hcipy.zernike(2, 2)(pupil_grid)
    wf_astig = hcipy.Wavefront(aperture * np.exp(1j * astig_phase), wavelength)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    hcipy.imshow_field(wf_defocus.phase, cmap='RdBu', ax=axes[0])
    axes[0].set_title('Defocus Phase')
    hcipy.imshow_field(wf_astig.phase, cmap='RdBu', ax=axes[1])
    axes[1].set_title('Astigmatism Phase')
    plt.show()



.. image:: output_7_0.png


Next we’ll propagate the planar wavefront to the focal plane using
``FraunhoferPropagator``, which computes the far-field diffraction
pattern — the point spread function (PSF). The resulting Airy pattern
has a bright central core surrounded by concentric rings, with the first
dark ring at a radius of about 1.22 λ/D from the center.

.. code:: ipython3

    wf_in = hcipy.Wavefront(aperture, wavelength)
    
    # Create focal grid and Fraunhofer propagator
    focal_grid = hcipy.make_focal_grid(q=4, num_airy=20, spatial_resolution=wavelength / pupil_diameter)
    
    fraunhofer = hcipy.FraunhoferPropagator(pupil_grid, focal_grid)
    wf_focal = fraunhofer(wf_in)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    hcipy.imshow_field(wf_in.amplitude, ax=axes[0])
    axes[0].set_title('Pupil Plane Amplitude')
    
    # Log-scale PSF
    img = wf_focal.intensity / wf_focal.intensity.max()
    hcipy.imshow_field(np.log10(img + 1e-20), cmap='inferno', vmin=-5, ax=axes[1])
    axes[1].set_title('Focal Plane PSF (Airy Pattern)')
    plt.show()



.. image:: output_9_0.png


For shorter propagation distances, the Fraunhofer approximation breaks
down, so we’ll use ``FresnelPropagator`` instead. This accounts for the
quadratic phase factor that accumulates during propagation. With a
Fresnel number of about 200, we’re firmly in the near-field regime where
the intensity pattern still resembles the pupil shape with some edge
diffraction effects.

.. code:: ipython3

    # Propagate with FresnelPropagator
    distance = 10000.0  # meters
    fresnel = hcipy.FresnelPropagator(pupil_grid, distance)
    wf_near = fresnel(wf_in)
    
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    hcipy.imshow_field(wf_near.intensity, cmap='inferno', ax=axes[0])
    axes[0].set_title(f'Intensity at {distance}m')
    hcipy.imshow_field(wf_near.phase, cmap='RdBu', ax=axes[1])
    axes[1].set_title(f'Phase at {distance}m')
    plt.show()
    
    print(f"Fresnel number: {pupil_diameter**2 / (wavelength * distance):.2f}")



.. image:: output_11_0.png


.. parsed-literal::

    Fresnel number: 200.00


Let’s add some aberrations to see how wavefront errors degrade image
quality. We’ll create a phase screen — a sinusoidal ripple modulated by
a Gaussian envelope — and apply it to the wavefront.

.. code:: ipython3

    # Create phase screen
    phase_screen = 0.5 * np.sin(2 * np.pi * pupil_grid.x) * np.exp(-pupil_grid.as_('polar').r**2)
    
    # Apply to wavefront
    wf_aberrated = hcipy.Wavefront(wf_in.electric_field * np.exp(1j * phase_screen), wavelength)

Now we propagate both the ideal and aberrated wavefronts to the focal
plane and compare their PSFs. The Strehl ratio compares the peak
focal-plane intensity of the aberrated PSF to the ideal one — values
near 1 mean diffraction-limited performance, while our value will be
noticeably lower, reflecting the phase errors we introduced.

.. code:: ipython3

    wf_focal_ideal = fraunhofer(wf_in)
    wf_focal_aberrated = fraunhofer(wf_aberrated)
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    hcipy.imshow_field(phase_screen, cmap='RdBu', ax=axes[0,0])
    axes[0,0].set_title('Phase Screen (radians)')
    
    hcipy.imshow_field(np.log10(wf_focal_ideal.intensity / wf_focal_ideal.intensity.max() + 1e-20), cmap='inferno', vmin=-5, ax=axes[0, 1])
    axes[0,1].set_title('Ideal PSF')
    
    hcipy.imshow_field(wf_aberrated.phase, cmap='RdBu', ax=axes[1,0])
    axes[1,0].set_title('Aberrated Wavefront Phase')
    
    hcipy.imshow_field(np.log10(wf_focal_aberrated.intensity / wf_focal_aberrated.intensity.max() + 1e-20), cmap='inferno', vmin=-5, ax=axes[1, 1])
    axes[1,1].set_title('Aberrated PSF')
    
    plt.tight_layout()
    plt.show()
    
    # Strehl ratio
    strehl = wf_focal_aberrated.intensity.max() / wf_focal_ideal.intensity.max()
    print(f"Strehl ratio: {strehl:.4f}")



.. image:: output_15_0.png


.. parsed-literal::

    Strehl ratio: 0.9544