Creating your own optical elements

We will present the internal workings of optical elements, and show how you can create your own optical elements in HCIPy.

Standard optical elements

We’ll start by importing all the necessary packages.

from hcipy import *
import numpy as np
import matplotlib.pyplot as plt
import traceback

An OpticalElement is an object that can propagate a Wavefront from one plane to another. This includes apodizers, deformable mirrors, focal-plane masks, coronagraphs or even complete optical systems. An optical element should in principle be able to handle any type of Wavefront passed to its forward() and backward() functions. This puts a large responsibility on the implementation of the optical element.

Any proper implementation of an optical element should therefore implement a number of functions to help signal its intent. Here, we are going to implement a very simple OpticalElement: a neutral-density filter. We’ll show the full implementation first, and then go in to the details for each implemented function.

class NeutralDensityFilter(OpticalElement):
    def __init__(self, transmittance):
        self.transmittance = transmittance

    def forward(self, wavefront):
        wf = wavefront.copy()

        wf.electric_field *= np.sqrt(self.transmittance)

        return wf

    def backward(self, wavefront):
        wf = Wavefront.copy()

        wf.electric_field *= np.sqrt(self.transmittance)

        return wf

A neutral-density filter is a filter that reduces the intensity of the transmitted light independent of wavelength or spatial position. Therefore, the above implementation can handle any incoming wavefront. The initializer of the class just stores the passed transmittance. In the initializer, usually we prefer to do as much computation as possible, to aliviate the burden on the forward() and backward() functions.

The forward() function propagates the wavefront through the filter. This function first creates a copy of the wavefront to work on, then modifies its electric field and returns it. Any implementation of forward() should not attempt to operate on a wavefront in-place and should always return a new Wavefront object.

The backward() function might need some extra explanation. It does not implement the propagation in the opposite direction through the optical element, or even the inverse or pseudoinverse of the forward propagation, but rather the adjoint. In absense of non-linear effects, any optical element can be thought of as a linear operator on a complex Hilbert space. This Hilbert space is the space of all possible wavefronts, and an optical element acts on a wavefront to produce another wavefront in the same space. This means that we can also construct the Hermitan adjoint of this operator. The backward() function implements this Hermitian adjoint propagation. In some cases, ie. when energy is conserved in the optical element, the backward() function is the inverse of the forward() function, but in general this is not the case.

Agnostic optical elements

Motivation

Sometimes we want more complicated behaviour. For instance, we might want the transmission of the optic to change as function of wavelength, or spatial position in the plane. As a case study, we take here an ColourFilter class, which implements the transmission through a colour filter. A naive implementation would be something along the lines of the following.

class NaiveColourFilter(OpticalElement):
    def __init__(self, filter_transmission):
        self.filter_transmission = filter_transmission

    def forward(self, wavefront):
        wf = wavefront.copy()

        transmission = self.filter_transmission(wavefront.wavelength)
        wf.electric_field *= np.sqrt(transmission)

        return wf

    def backward(self, wavefront):
        wf = wavefront.copy()

        transmission = self.filter_transmission(wavefront.wavelength)
        wf.electric_field *= np.sqrt(transmission)

        return wf

We would initialize this class with a transmission which is a function of wavelength.

central_wavelength = 700e-9 # m
spectral_bandwidth = 100e-9 # m

def block_filter(wavelength):
    lower = central_wavelength - spectral_bandwidth / 2
    upper = central_wavelength + spectral_bandwidth / 2

    return (lower < wavelength) and (wavelength < upper)

colour_filter = NaiveColourFilter(block_filter)

pupil_grid = make_pupil_grid(128)
wf = Wavefront(pupil_grid.ones())
post_filter = colour_filter(wf)

However, this implementation doesn’t support a float as the filter transmission, as a float cannot be called with a wavelength.

try:
    colour_filter = NaiveColourFilter(0.7)
    post_ilter = colour_filter(wf)
except Exception as e:
    traceback.print_exc()

Traceback (most recent call last):
  File "<ipython-input-7-f11581191e03>", line 3, in <module>
    post_ilter = colour_filter(wf)
  File "c:usersemiel porgithubhcipyhcipyopticsoptical_element.py", line 28, in __call__
    return self.forward(wavefront)
  File "<ipython-input-5-d848e712e35d>", line 8, in forward
    transmission = self.filter_transmission(wavefront.wavelength)
TypeError: 'float' object is not callable

This narrow behaviour is fine for quick development, that is, in cases where we expect ourselves to be the only ones using the optical element and/or there are only a few use cases when we would want to use it. For more general optical elements, however, we would prefer the objects to be able to handle a wide variety of input arguments.

Take for example the LinearRetarder optical element. The apodization pass as the argument to its initializer/constructor can be 1. a scalar, 2. a Field, 3. a Field generator, 4. a function of wavelength, returning either a Field or scalar, 5. a function of both a Grid and a wavelength, returning a Field.

If they apodization will be evaluated at the Grid and wavelength of the incoming Wavefront. All of this is done seamlessly and is hidden from view for the user. Another example is that of the LinearRetarder, which should be able to handle (at least) a phase retardance which is either constant with wavelength (ie. a scalar) or a function of wavelength. Manually distinguising between the above cases can become quite cumbersome, so HCIPy offers a simpler solution.

Implementation details

Let’s look at the implementation of Apodizer as an example:

class Apodizer(AgnosticOpticalElement):
    def __init__(self, apodization):
        self.apodization = apodization

        AgnosticOpticalElement.__init__(self, grid_dependent=True, wavelength_dependent=True, max_in_cache=11)

    def make_instance(self, instance_data, input_grid, output_grid, wavelength):
        instance_data.apodization = self.evaluate_parameter(self.apodization, input_grid, output_grid, wavelength)

    @property
    def apodization(self):
        return self._apodization

    @apodization.setter
    def apodization(self, apodization):
        self._apodization = apodization

        self.clear_cache()

    def get_input_grid(self, output_grid, wavelength):
        return output_grid

    def get_output_grid(self, input_grid, wavelength):
        return input_grid

    @make_agnostic_forward
    def forward(self, instance_data, wavefront):
        wf = wavefront.copy()
        wf.electric_field *= instance_data.apodization

        return wf

    @make_agnostic_backward
    def backward(self, instance_data, wavefront):
        wf = wavefront.copy()
        wf.electric_field *= np.conj(instance_data.apodization)

        return wf

The Apodizer class is derived from AgnosticOpticalElement rather than OpticalElement. The AgnosticOpticalElement class provides the ability of an optical element to more easily work correctly for a wide range of inputs. It allows the optical element to be written mostly expecting a single input Grid and wavelength. The data required for each set of grid and wavelength is cached, so expensive calculations are not redone each time a Wavefront is propagated through.

Setup code for an AgnosticOpticalElement is done in two different functions. The first is the initializer __init__(). Anything done in this function should not depend on wavelength or input grid. The above implmenetation just stores the value of apodization whatever that may be, and calls the __init__() of the parent class. This initializer has a few parameters, indicating if the AgnosticOpticalELlement class should create new instances when a Wavefront differs in Grid or in wavelength. For example, if our optical element is achromatic, but can still depend on the input grid, we would set wavelength_dependent to False. The last argument is the maximum number of instances to keep in the internal cache. If the class stores a lot of data per instance, you might run out of memory if too many instances are kept in memory. This value can be reduced if necessary (or increased if the set up time per instance is high, but memory usage is small).

The second initializer is make_instance(). This function should calculate any data that depends on input grid, output grid and/or wavelength. Any calculated data should be stored in instance_data which is the object that is stored in the cache. In the implmementation above, we evalute the apodization at the supplied input grid, output grid and wavelength. This is done with the evaluate_parameter() function, which does all the heavy lifting with figuring out what type of parameter it actually is (ie. is it dependent on grid, or wavelength, or both, or none), and evaluates with the appropriate signature.

The apodization is implemented in the class using a property. The reason for this is that we want to clear the cache when we change the apodization. If we wouldn’t do this, we would still use the old apodization during propagation on a grid and wavelength that we have seen before, and were in the cache already.

The next two functions let the AgnosticOpticaElement know what the corresponding input grid is, given an output grid and wavelength, and vice versa. This allows backwards propagation through an optical element, without having done a forward propagation first. (The exact reason for this is a bit harder to explain. When doing a backward propagation, it needs to look up the instance data in the cache. At that point, only the output grid is available, while the cache would have instance_datas listed by input grid. In this case it is necessary to know the input grid corresponding to previously calculated instances.)

Finally, we have the forward() and backward() functions. These are implemened with the use of the make_agnostic_forward() and make_agnostic_backward() decorators, which look up the correct instance_data corresponding to the supplied Wavefront and call the forward() or backward() functions with these as well. Then, inside the propagation functions themselves, you have access to the evaluated parameters via instance_data. In the case of the Apodizer, we have access to the instance_data.apodization, and we propagate through the apodizer in the usual way.

Adding properties based on other properties

Suppose now that we would want to know the phase of the apodization in the apodizer. If we would have the apodization as a Field or a scalar, it would be easy:

apodization = Field(np.random.randn(pupil_grid.size) + 1j * np.random.randn(pupil_grid.size), pupil_grid)
phase = np.angle(apodization)

But in case the apodization is a function, things are not that trivial. The AgnosticOpticalElement also provides a method for simplifying this operation. We will inherit from Apodizer and create a SpecialApodizer with this new method:

class SpecialApodizer(Apodizer):
    def __init__(self, apodization):
        Apodizer.__init__(self, apodization)

    @property
    def phase(self):
        return self.construct_function(np.angle, self.apodization)

The construct_function() method of AgnosticOpticalElement returns a function with the encompassing signature. This means that if one of the arguments of the function depends on input grid, then the returned function will also be a function of input grid. If at least one of the arguments of the function is a function of wavelength, then the returned function will also be a function of wavelength. The non of the arguments were functions, then the supplied function is evaluated with the supplied arguments.

The net effect is that, for this function, the phase will behave in the same way as apodization did.

Conclusion

Implementing your own optical elements can be as easy as deriving from the OpticalElement class and overloading (= rewriting) the forward() and backward() functions.

When arguments to your optical element depend on the grid or the wavelength of the propagated wavefront, it sometimes is easier to derive from AgnosticOpticalElement instead. This class provides several methods and a caching strategy to simplify otherwise complicated optical elements. The best way to really get an insight in the internals of HCIPy is to look in the source code. A good place to start is are the apodizers, Apodizer, PhaseApodizer and SurfaceApodizer. These classes provide an easy and uncomplicated starting point. If the user is familiar with polarization and Jones matrices, the JonesMatrixOpticalElement and PhaseRetarder classes are also good examples of relatively simple agnostic optical elements.