Segmented deformable mirrors

We will use segmented deformable mirrors and simulate the PSFs that result from segment pistons and tilts. We will compare this functionality against Poppy, another optical propagation package.

First we’ll import all packages.

%matplotlib inline
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None, 'figsize': (8, 6)}

import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (8, 6)
mpl.rcParams['figure.dpi'] = 150
mpl.rcParams['savefig.dpi'] = 150

import os
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
import hcipy
import poppy

# Parameters - these will depend on the aperture function you use
PUP_DIAMETER = 0.019725   # m
GAPSIZE = 90e-6  # m
FLATTOFLAT = PUP_DIAMETER/7.   # m

# these will not
num_pix = 1024
wavelength = 638e-9
num_airy = 20
sampling = 4
norm = False

Instantiate the segmented mirrors

HCIPy SM: hsm

We need to generate a pupil grid for the aperture, and a focal grid and propagator for the focal plane images after the DM.

# HCIPy grids and propagator
pupil_grid = hcipy.make_pupil_grid(dims=num_pix, diameter=PUP_DIAMETER)
#focal_grid = hcipy.make_focal_grid(pupil_grid, sampling, num_airy, wavelength=wavelength)
focal_grid = hcipy.make_uniform_grid([2*num_airy*sampling]*2, 2*num_airy * wavelength / PUP_DIAMETER)
prop = hcipy.FraunhoferPropagator(pupil_grid, focal_grid)

We generate a segmented aperture for the segmented mirror. For convenience, we’ll use the HiCAT pupil without spiders. We’ll use supersampling to better resolve the segment gaps.

aper, segments = hcipy.make_hicat_aperture(normalized=norm, with_spiders=False, return_segments=True)
aper = hcipy.evaluate_supersampled(aper, pupil_grid, 1)
segments = hcipy.evaluate_supersampled(segments, pupil_grid, 1)

plt.title('HCIPy aperture')
hcipy.imshow_field(aper)

<matplotlib.image.NonUniformImage at 0x7fa7f4528690>

Now we make the segmented mirror. In order to be able to apply the SM to a plane, that plane needs to be a Wavefront, which combines a Field - here the aperture - with a wavelength, here wavelength.

In this example here, since the SM doesn’t have any extra effects on the pupil since it’s completely flat still, we don’t actually have to apply the SM, although of course we could.

# Instantiate the segmented mirror
hsm = hcipy.SegmentedDeformableMirror(segments)

# Make a pupil plane wavefront from aperture
wf = hcipy.Wavefront(aper, wavelength)

# Apply SM if you want to
wf = hsm(wf)

plt.title('Wavefront intensity at HCIPy SM')
hcipy.imshow_field(wf.intensity)

<matplotlib.image.NonUniformImage at 0x7fa7b4e69390>

Poppy SM: psm

We’ll do the same for Poppy.

psm = poppy.dms.HexSegmentedDeformableMirror(name='Poppy SM',
                                             rings=3,
                                             flattoflat=FLATTOFLAT*u.m,
                                             gap=GAPSIZE*u.m,
                                             center=False)

# Display the transmission and phase of the poppy sm
plt.figure(figsize=(16, 8))
psm.display(what='both')

(<matplotlib.axes._subplots.AxesSubplot at 0x7fa7e08e0b50>,
 <matplotlib.axes._subplots.AxesSubplot at 0x7fa7e08e8410>)

Create reference images

HCIPy reference image

We need to apply the SM to the wavefront in the pupil plane and then propagate it to the image plane.

# Apply SM to pupil plane wf
wf_sm = hsm(wf)

# Propagate from SM to image plane
im_ref_hc = prop(wf_sm)

# Display intensity and phase in image plane
plt.figure(figsize=(18, 9))
plt.suptitle('Image plane after HCIPy SM')

# Get normalization factor for HCIPy reference image
norm_hc = np.max(im_ref_hc.intensity)

plt.subplot(1, 2, 1)
hcipy.imshow_field(np.log10(im_ref_hc.intensity/norm_hc))
plt.title('Intensity')
plt.colorbar()

plt.subplot(1, 2, 2)
hcipy.imshow_field(im_ref_hc.phase, cmap='RdBu')
plt.title('Phase')
plt.colorbar()

print('HCIPy PSF shape: {}'.format(im_ref_hc.intensity.shaped.shape))

HCIPy PSF shape: (160, 160)

Poppy reference image

For the Poppy propagation, we need to make an optical system of which we then calculate the PSF.

I will try to match the image resolution and size of the HCIPy image. I first adjust the pixelscale and fov_arcsec such that their ratio works and then I add a tweak factor fac to scale it to the HCIPy image. I also set oversample to something that matches the HCIPy sampling (it’s close enough). I keep reusing these numbers and the tweak factor later on in the notebook.

# Make an optical system with the Poppy SM and a detector
psm.flatten()
osys = poppy.OpticalSystem()
osys.add_pupil(psm)
fac = 5300
pxscle = 10 * np.degrees(wavelength / PUP_DIAMETER) * 3600.0 / sampling * 0.97
fovarc = pxscle * 160 / 10
osys.add_detector(pixelscale=pxscle, fov_arcsec=fovarc, oversample=10)

<poppy.poppy_core.Detector at 0x7fa7d033f710>

# Calculate the PSF
psf = osys.calc_psf(wavelength)
plt.figure(figsize=(10, 10))
poppy.display_psf(psf, vmin=1e-9, vmax=0.1)

# Get the PSF as an array
im_ref_pop = psf[0].data
print('Poppy PSF shape: {}'.format(im_ref_pop.shape))

# Get normalization from Poppy reference image
norm_pop = np.max(im_ref_pop)

---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

~/anaconda3/envs/hcipy-docs-v0.3.0/lib/python3.7/site-packages/numpy/core/function_base.py in linspace(start, stop, num, endpoint, retstep, dtype, axis)
    116     try:
--> 117         num = operator.index(num)
    118     except TypeError:


TypeError: 'numpy.float64' object cannot be interpreted as an integer


During handling of the above exception, another exception occurred:


TypeError                                 Traceback (most recent call last)

<ipython-input-13-25ea51ea132a> in <module>
      2 psf = osys.calc_psf(wavelength)
      3 plt.figure(figsize=(10, 10))
----> 4 poppy.display_psf(psf, vmin=1e-9, vmax=0.1)
      5
      6 # Get the PSF as an array


~/anaconda3/envs/hcipy-docs-v0.3.0/lib/python3.7/site-packages/poppy/utils.py in display_psf(HDUlist_or_filename, ext, vmin, vmax, scale, cmap, title, imagecrop, adjust_for_oversampling, normalize, crosshairs, markcentroid, colorbar, colorbar_orientation, pixelscale, ax, return_ax, interpolation, cube_slice)
    256             )
    257         if scale.lower() == 'log':
--> 258             ticks = np.logspace(np.log10(vmin), np.log10(vmax), np.log10(vmax / vmin) + 1)
    259             if colorbar_orientation == 'horizontal' and vmax == 1e-1 and vmin == 1e-8:
    260                 ticks = [1e-8, 1e-6, 1e-4, 1e-2, 1e-1]  # looks better


<__array_function__ internals> in logspace(*args, **kwargs)


~/anaconda3/envs/hcipy-docs-v0.3.0/lib/python3.7/site-packages/numpy/core/function_base.py in logspace(start, stop, num, endpoint, base, dtype, axis)
    270
    271     """
--> 272     y = linspace(start, stop, num=num, endpoint=endpoint, axis=axis)
    273     if dtype is None:
    274         return _nx.power(base, y)


<__array_function__ internals> in linspace(*args, **kwargs)


~/anaconda3/envs/hcipy-docs-v0.3.0/lib/python3.7/site-packages/numpy/core/function_base.py in linspace(start, stop, num, endpoint, retstep, dtype, axis)
    119         raise TypeError(
    120             "object of type {} cannot be safely interpreted as an integer."
--> 121                 .format(type(num)))
    122
    123     if num < 0:


TypeError: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.

Both reference images side-by-side

plt.figure(figsize=(20,10))

plt.subplot(1, 2, 1)
hcipy.imshow_field(np.log10(im_ref_hc.intensity / norm_hc))
plt.title('HCIPy reference PSF')
plt.colorbar()

plt.subplot(1, 2, 2)
plt.imshow(np.log10(im_ref_pop / norm_pop), origin='lower')
plt.title('Poppy reference PSF')
plt.colorbar()

ref_dif = im_ref_pop / norm_pop - im_ref_hc.intensity.shaped / norm_hc

plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
plt.imshow(ref_dif, origin='lower')
plt.title('Full image')
plt.colorbar()
plt.subplot(1, 2, 2)
plt.imshow(ref_dif[60:100,60:100], origin='lower')
plt.title('Zoomed in')
plt.colorbar()

Applying aberrations

# Define function from rad of phase to m OPD
def aber_to_opd(aber_rad, wavelength):
    aber_m = aber_rad * wavelength / (2 * np.pi)
    return aber_m

aber_rad = 4.0

print('Aberration: {} rad'.format(aber_rad))
print('Aberration: {} m'.format(aber_to_opd(aber_rad, wavelength)))

# Flatten both SMs just to be sure
hsm.flatten()
psm.flatten()

poppy_to_hcipy_index = {
    1: 2, 2: 1, 3: 0, 4: 5, 5: 4, 6: 3,
    7: 10, 8: 9, 9: 8, 10: 7, 11: 6, 12: 17, 13: 16, 14: 15, 15: 14, 16: 13, 17: 12, 18: 11,
    19: 24, 20: 23, 21: 22, 22: 21, 23: 20, 24: 19, 25: 18,
    26: 35, 27: 34, 28: 33, 29: 32, 30: 31, 31: 30, 32: 29, 33: 28, 34: 27, 35: 26, 36: 25}

for i in [35, 25]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], aber_to_opd(aber_rad, wavelength) / 2, 0, 0)
    psm.set_actuator(i, aber_to_opd(aber_rad, wavelength) * u.m, 0, 0)

# Display both segmented mirrors in OPD

# HCIPy
plt.figure(figsize=(8,8))
plt.title('Phase for HCIPy SM')
hcipy.imshow_field(hsm.surface * 2, mask=aper, cmap='RdBu_r', vmin=-5e-7, vmax=5e-7)
plt.colorbar()

# Poppy
plt.figure(figsize=(8,8))
psm.display(what='opd')

Show focal plane images

### HCIPy
# Apply SM to pupil plane wf
wf_fp_pistoned = hsm(wf)

# Propagate from SM to image plane
im_pistoned_hc = prop(wf_fp_pistoned)

### Poppy
# Make an optical system with the Poppy SM and a detector
osys = poppy.OpticalSystem()
osys.add_pupil(psm)
pxscle = 0.0031*fac      # I'm tweaking pixelscale and fov_arcsec to match the HCIPy image
fovarc = 0.05*fac
osys.add_detector(pixelscale=pxscle, fov_arcsec=fovarc, oversample=10)

# Calculate the PSF
psf = osys.calc_psf(wavelength)
plt.figure(figsize=(10, 10))

# Get the PSF as an array
im_pistoned_pop = psf[0].data

### Display intensity of both cases image plane
plt.figure(figsize=(18, 9))
plt.suptitle('Image plane after SM for $\phi$ = ' + str(aber_rad) + ' rad')

plt.subplot(1, 2, 1)
hcipy.imshow_field(np.log10(im_pistoned_hc.intensity / norm_hc))
plt.title('HCIPy pistoned pair')
plt.colorbar()

plt.subplot(1, 2, 2)
plt.imshow(np.log10(im_pistoned_pop / norm_pop), origin='lower')
plt.title('Poppy pistoned pair')
plt.colorbar()

Image degradation as function of rms piston errors

We will plot the

# Aberration range
aber_array = np.linspace(0, 2*np.pi, 50, True)
print('Aber in rad: \n{}'.format(aber_array))
print('Aber in m: \n{}'.format(aber_to_opd(aber_array, wavelength)))

# Apply pistons
hc_ims = []
pop_ims = []
for aber_rad in aber_array:
    # Flatten both SMs
    hsm.flatten()
    psm.flatten()

    # HCIPy
    for i in [34, 25]:
        opd = aber_to_opd(aber_rad, wavelength)
        hsm.set_segment_actuators(poppy_to_hcipy_index[i], opd / 2, 0, 0)   # hcipy uses SURFACE, not OPD
        psm.set_actuator(i, opd * u.m, 0, 0)

    # Propagate to image plane

    # HCIPy:
    # Propagate from pupil plane through SM to image plane
    im_pistoned_hc = prop(hsm(wf))

    # Poppy:
    # Make an optical system with the Poppy SM and a detector
    osys = poppy.OpticalSystem()
    osys.add_pupil(psm)
    pxscle = 0.0031 * fac  # I'm tweaking pixelscale and fov_arcsec to match the HCIPy image
    fovarc = 0.05 * fac
    osys.add_detector(pixelscale=pxscle, fov_arcsec=fovarc, oversample=10)

    # Calculate the PSF
    psf = osys.calc_psf(wavelength)

    # Get the PSF as an array
    im_pistoned_pop = psf[0].data

    hc_ims.append(im_pistoned_hc.intensity.shaped / im_pistoned_hc.intensity.max())
    pop_ims.append(im_pistoned_pop / im_pistoned_pop.max())

hc_ims = np.array(hc_ims)
pop_ims = np.array(pop_ims)

### Quantify with image sums
sum_hc = np.sum(hc_ims, axis=(1,2))
sum_hc /= sum_hc[0]
sum_pop = np.sum(pop_ims, axis=(1,2))
sum_pop /= sum_pop[0]

plt.suptitle('Image degradation of SMs')
plt.plot(aber_array, sum_hc, '-', label='HCIPy SM')
plt.plot(aber_array, sum_pop, '--', label='Poppy SM')
plt.xlabel('phase aberration (rad)')
plt.ylabel('image sum')
plt.legend()
plt.show()

A mix of piston, tip and tilt (PTT)

aber_rad_tt = 500e-6
aber_rad_p = 1.8

opd_piston = aber_to_opd(aber_rad_p, wavelength)

### Put aberrations on both SMs
# Flatten both SMs
hsm.flatten()
psm.flatten()

## PISTON
for i in [19, 28, 23, 16]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], opd_piston / 2, 0, 0)
    psm.set_actuator(i, opd_piston * u.m, 0, 0)

for i in [3, 35, 30, 8]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], 0.5 * opd_piston / 2, 0, 0)
    psm.set_actuator(i, 0.5 * opd_piston * u.m, 0, 0)

for i in [14, 18, 1, 32, 12]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], 0.3 * opd_piston / 2, 0, 0)
    psm.set_actuator(i, 0.3 * opd_piston * u.m, 0, 0)

## TIP and TILT
for i in [2, 5, 11, 15, 22]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], 0, aber_rad_tt / 2, 0.3 * aber_rad_tt / 2)
    psm.set_actuator(i, 0, aber_rad_tt, 0.3 * aber_rad_tt)

for i in [4, 6, 36]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], 0, aber_rad_tt / 2, 0)
    psm.set_actuator(i, 0, aber_rad_tt, 0)

for i in [34, 31, 7]:
    hsm.set_segment_actuators(poppy_to_hcipy_index[i], 0, 0, 1.3 * aber_rad_tt / 2)
    psm.set_actuator(i, 0, 0, 1.3 * aber_rad_tt)

# Display both segmented mirrors in OPD

# HCIPy
plt.figure(figsize=(8,8))
plt.title('OPD for HCIPy SM')
hcipy.imshow_field(hsm.surface * 2, mask=aper, cmap='RdBu_r', vmin=-5e-7, vmax=5e-7)
plt.colorbar()

# Poppy
plt.figure(figsize=(8,8))
psm.display(what='opd')

### Propagate to image plane
## HCIPy
# Propagate from pupil plane through SM to image plane
im_pistoned_hc = prop(hsm(wf))

## Poppy
# Make an optical system with the Poppy SM and a detector
osys = poppy.OpticalSystem()
osys.add_pupil(psm)
pxscle = 0.0031 * fac  # I'm tweaking pixelscale and fov_arcsec to match the HCIPy image
fovarc = 0.05 * fac
osys.add_detector(pixelscale=pxscle, fov_arcsec=fovarc, oversample=10)

# Calculate the PSF
psf = osys.calc_psf(wavelength)

# Get the PSF as an array
im_pistoned_pop = psf[0].data

### Display intensity of both cases image plane
plt.figure(figsize=(18, 9))
plt.suptitle('Image plane after SM forrandom arangement')

plt.subplot(1, 2, 1)
hcipy.imshow_field(np.log10(im_pistoned_hc.intensity/norm_hc))
plt.title('HCIPy random arangement')
plt.colorbar()

plt.subplot(1, 2, 2)
plt.imshow(np.log10(im_pistoned_pop/norm_pop), origin='lower')
plt.title('Poppy tipped arangement')
plt.colorbar()
plt.show()