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.

.. code:: ipython3

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

.. code:: ipython3

    # Parameters for the pupil function
    pupil_diameter = 0.019725 # m
    gap_size = 90e-6 # m
    num_rings = 3
    segment_flat_to_flat = (pupil_diameter - (2 * num_rings + 1) * gap_size) / (2 * num_rings + 1)
    focal_length = 1 # m
    
    # Parameters for the simulation
    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.

.. code:: ipython3

    # HCIPy grids and propagator
    pupil_grid = hcipy.make_pupil_grid(dims=num_pix, diameter=pupil_diameter)
    
    focal_grid = hcipy.make_focal_grid(sampling, num_airy,
                                       pupil_diameter=pupil_diameter,
                                       reference_wavelength=wavelength,
                                       focal_length=focal_length)
    focal_grid = focal_grid.shifted(focal_grid.delta / 2)
    
    prop = hcipy.FraunhoferPropagator(pupil_grid, focal_grid, focal_length)

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.

.. code:: ipython3

    aper, segments = hcipy.make_hexagonal_segmented_aperture(num_rings, 
                                                             segment_flat_to_flat, 
                                                             gap_size, 
                                                             starting_ring=1,
                                                             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, cmap='gray')




.. parsed-literal::

    <matplotlib.image.NonUniformImage at 0x11d868c32c8>




.. image:: output_6_1.png


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 still completely flat, we don’t actually have to apply
the SM, although of course we could.

.. code:: ipython3

    # 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.figure(figsize=(8, 8))
    plt.title('Wavefront intensity at HCIPy SM')
    hcipy.imshow_field(wf.intensity, cmap='gray')
    plt.colorbar()
    plt.show()



.. image:: output_8_0.png


Poppy SM: ``psm``
~~~~~~~~~~~~~~~~~

We’ll do the same for Poppy.

.. code:: ipython3

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

.. code:: ipython3

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




.. parsed-literal::

    <AxesSubplot:title={'center':'Transmissivity for Poppy SM'}, ylabel='[meters]'>




.. image:: output_11_1.png


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.

.. code:: ipython3

    # Apply SM to pupil plane wf
    wf_sm = hsm(wf)
    
    # Propagate from SM to image plane
    im_ref_hc = prop(wf_sm)

.. code:: ipython3

    # Display intensity and phase in image plane
    plt.figure(figsize=(8, 8))
    plt.suptitle('Image plane after HCIPy SM')
    
    # Get normalization factor for HCIPy reference image
    norm_hc = np.max(im_ref_hc.intensity)
    
    hcipy.imshow_psf(im_ref_hc, normalization='peak')




.. parsed-literal::

    <matplotlib.image.NonUniformImage at 0x11d8b23c648>




.. image:: output_14_1.png


Poppy reference image
~~~~~~~~~~~~~~~~~~~~~

For the Poppy propagation, we need to make an optical system of which we
then calculate the PSF. We match HCIPy’s image scale with Poppy.

.. code:: ipython3

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




.. parsed-literal::

    <poppy.poppy_core.Detector at 0x11d8b34b7c8>



.. code:: ipython3

    # Calculate the PSF
    psf = osys.calc_psf(wavelength)
    plt.figure(figsize=(8, 8))
    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)


.. parsed-literal::

    Poppy PSF shape: (160, 160)
    


.. image:: output_17_1.png


Both reference images side-by-side
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    plt.figure(figsize=(15,6))
    
    plt.subplot(1, 2, 1)
    hcipy.imshow_field(np.log10(im_ref_hc.intensity / norm_hc), vmin=-10, cmap='inferno')
    plt.title('HCIPy reference PSF')
    plt.colorbar()
    
    plt.subplot(1, 2, 2)
    plt.imshow(np.log10(im_ref_pop / norm_pop), origin='lower', vmin=-10, cmap='inferno')
    plt.title('Poppy reference PSF')
    plt.colorbar()




.. parsed-literal::

    <matplotlib.colorbar.Colorbar at 0x11d88fef548>




.. image:: output_19_1.png


.. code:: ipython3

    ref_dif = im_ref_pop / norm_pop - im_ref_hc.intensity.shaped / norm_hc
    
    lims = np.max(np.abs(ref_dif))
    
    plt.figure(figsize=(15, 6))
    plt.suptitle(f'Maximum relative error: {lims:0.2g} relative to the peak intensity')
    plt.subplot(1, 2, 1)
    plt.imshow(ref_dif, origin='lower', vmin=-lims, vmax=lims, cmap='RdBu')
    plt.title('Full image')
    plt.colorbar()
    plt.subplot(1, 2, 2)
    plt.imshow(ref_dif[60:100,60:100], origin='lower', vmin=-lims, vmax=lims, cmap='RdBu')
    plt.title('Zoomed in')
    plt.colorbar()




.. parsed-literal::

    <matplotlib.colorbar.Colorbar at 0x11d8916b808>




.. image:: output_20_1.png


Applying aberrations
--------------------

.. code:: ipython3

    # 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)))
    
    # Poppy and HCIPy have a different way of indexing segments
    # Figure out which index to poke on which mirror
    poppy_index_to_hcipy_index = []
    for n in range(1, num_rings + 1):
        base = list(range(3 * (n - 1) * n + 1, 3 * n * (n + 1) + 1))
    
        poppy_index_to_hcipy_index.extend(base[2 * n::-1])
        poppy_index_to_hcipy_index.extend(base[:2 * n:-1])
    
    poppy_index_to_hcipy_index = {j: i for i, j in enumerate(poppy_index_to_hcipy_index) if j is not None}
    
    hcipy_index_to_poppy_index = {j: i for i, j in poppy_index_to_hcipy_index.items()}


.. parsed-literal::

    Aberration: 4.0 rad
    Aberration: 4.061634147705169e-07 m
    

.. code:: ipython3

    # Flatten both SMs just to be sure
    hsm.flatten()
    psm.flatten()
    
    # Poking segment 35 and 25
    for i in [35, 25]:
        hsm.set_segment_actuators(i, aber_to_opd(aber_rad, wavelength) / 2, 0, 0)
        psm.set_actuator(hcipy_index_to_poppy_index[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('OPD for HCIPy SM')
    hcipy.imshow_field(hsm.surface * 2, mask=aper, cmap='RdBu_r', vmin=-5e-7, vmax=5e-7)
    plt.colorbar()
    plt.show()
    
    # Poppy
    plt.figure(figsize=(8,8))
    psm.display(what='opd')
    plt.show()



.. image:: output_23_0.png



.. image:: output_23_1.png


Show focal plane images
~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    ### 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
    # 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=(15, 6))
    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), cmap='inferno', vmin=-9)
    plt.title('HCIPy pistoned pair')
    plt.colorbar()
    
    plt.subplot(1, 2, 2)
    plt.imshow(np.log10(im_pistoned_pop / norm_pop), origin='lower', cmap='inferno', vmin=-9)
    plt.title('Poppy pistoned pair')
    plt.colorbar()




.. parsed-literal::

    <matplotlib.colorbar.Colorbar at 0x11d87f5e988>




.. image:: output_25_1.png


A mix of piston, tip and tilt (PTT)
-----------------------------------

.. code:: ipython3

    aber_rad_tt = 200e-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(i, opd_piston / 2, 0, 0)
        psm.set_actuator(hcipy_index_to_poppy_index[i], opd_piston * u.m, 0, 0)
        
    for i in [3, 35, 30, 8]:
        hsm.set_segment_actuators(i, -0.5 * opd_piston / 2, 0, 0)
        psm.set_actuator(hcipy_index_to_poppy_index[i], -0.5 * opd_piston * u.m, 0, 0)
        
    for i in [14, 18, 1, 32, 12]:
        hsm.set_segment_actuators(i, 0.3 * opd_piston / 2, 0, 0)
        psm.set_actuator(hcipy_index_to_poppy_index[i], 0.3 * opd_piston * u.m, 0, 0)
        
    ## TIP and TILT
    for i in [2, 5, 11, 15, 22]:
        hsm.set_segment_actuators(i, 0, aber_rad_tt / 2, 0.3 * aber_rad_tt / 2)
        psm.set_actuator(hcipy_index_to_poppy_index[i], 0, aber_rad_tt, 0.3 * aber_rad_tt)
        
    for i in [4, 6, 26]:
        hsm.set_segment_actuators(i, 0, -aber_rad_tt / 2, 0)
        psm.set_actuator(hcipy_index_to_poppy_index[i], 0, -aber_rad_tt, 0)
        
    for i in [34, 31, 7]:
        hsm.set_segment_actuators(i, 0, 0, 1.3 * aber_rad_tt / 2)
        psm.set_actuator(hcipy_index_to_poppy_index[i], 0, 0, 1.3 * aber_rad_tt)

.. code:: ipython3

    # 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()
    plt.show()
    
    # Poppy
    plt.figure(figsize=(8,8))
    psm.display(what='opd')
    plt.show()



.. image:: output_28_0.png



.. image:: output_28_1.png


.. code:: ipython3

    ### Propagate to image plane
    ## HCIPy
    # Propagate from pupil plane through SM to image plane
    im_pistoned_hc = prop(hsm(wf)).intensity
    
    ## Poppy
    # Calculate the PSF
    psf = osys.calc_psf(wavelength)
    
    # Get the PSF as an array
    im_pistoned_pop = psf[0].data

.. code:: ipython3

    ### 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 / norm_hc), cmap='inferno', vmin=-9)
    plt.title('HCIPy random arangement')
    plt.colorbar()
    
    plt.subplot(1, 2, 2)
    plt.imshow(np.log10(im_pistoned_pop / norm_pop), origin='lower', cmap='inferno', vmin=-9)
    plt.title('Poppy tipped arangement')
    plt.colorbar()
    plt.show()



.. image:: output_30_0.png