Introduction to polarization
============================

We will introduce how polarization is simulated in HCIPy, both for full
and partially polarized light. We will show which polarization optical
elements are implemented and how these are used.

We’ll start by implementing the required packages and set up a circular
pupil.

.. code:: ipython3

    from hcipy import *
    import numpy as np
    import matplotlib.pyplot as plt
    
    pupil_grid = make_pupil_grid(256)
    
    aperture = make_magellan_aperture(True)(pupil_grid)

Types of polarized Wavefronts
-----------------------------

HCIPy supports three different types of wavefronts, which indicate
different types of polarization. These are:

1. Scalar wavefronts
2. Fully polarized wavefronts
3. Partially polarized wavefronts

Each of these different types of wavefronts are implemented in a
``Wavefront`` object. There are a number of ways to get the polarization
state of the ``Wavefront``. These include (components of) the Stokes
vector, the degree of linear/circular polarization, etc… The different
types of wavefronts are distinguished by the tensor shape of the
``electric_field`` that is in the ``Wavefront``.

When you are not using polarization, you are most likely mostly using
scalar wavefronts. These contain a scalar ``electric_field``.

.. code:: ipython3

    scalar_wavefront = Wavefront(aperture)

Fully polarized wavefronts use Jones vectors as their electric fields.
This means that the electric field will be a vector ``Field`` with two
components:

.. math:: \begin{bmatrix}E_x \\ E_y\end{bmatrix}

.

.. code:: ipython3

    jones_vector = Field([aperture, 1j * aperture], pupil_grid)
    
    fully_polarized_wavefront = Wavefront(jones_vector)

Partially polarized wavefronts use Jones matrices as their electric
fields. The electric field will be a 2x2 tensor ``Field``:

.. math:: \begin{bmatrix} J_{xx} & J_{xy} \\ J_{yx} & J_{yy}\end{bmatrix}.

These Jones matrices describe the change in polarization state and also
phase and amplitude, from the initial wavefront to the current
wavefront. The polarization state itself is set as the Stokes vector of
the initial wavefront. Note that, in this way, HCIPy doesn’t support
spatially varying degrees of polarization.

While you can create a full Jones matrix ``Field`` and pass that as the
electric field argument of the ``Wavefront``, like this:

.. code:: ipython3

    jones_matrix = Field([[aperture, aperture * 1j], [aperture * 1j, aperture]], pupil_grid)
    stokes_vector = [1, 0.3, -0.2, 0.1]
    
    partially_polarized_wavefront = Wavefront(jones_matrix, input_stokes_vector=stokes_vector)

in most cases, we are trying to create a ``Wavefront`` at the start of
our optical system. In this case we can pass a scalar electric field,
but add an input Stokes vector as well. This creates a Jones matrix
``Field`` which has for each pixel an identity matrix multiplied by the
scalar wavefront at that position. For example, creating a partially
polarized ``Wavefront`` at our telescope pupil, we would do:

.. code:: ipython3

    stokes_vector = [1, 0.3, -0.2, 0.1]
    
    partially_polarized_wavefront = Wavefront(aperture, input_stokes_vector=stokes_vector)

Polarization attributes of Wavefronts
-------------------------------------

We can look at the Stokes components of the wavefront, and show them:

.. code:: ipython3

    def imshow_stokes_vector(wavefront):
        I_max = wavefront.I.max()
        
        plt.subplot(2, 2, 1)
        plt.title('I')
        imshow_field(wavefront.I, cmap='inferno')
        plt.colorbar()
        
        plt.subplot(2, 2, 2)
        plt.title('Q')
        imshow_field(wavefront.Q, cmap='RdBu', vmin=-I_max, vmax=I_max)
        plt.colorbar()
        
        plt.subplot(2, 2, 3)
        plt.title('U')
        imshow_field(wavefront.U, cmap='RdBu', vmin=-I_max, vmax=I_max)
        plt.colorbar()
        
        plt.subplot(2, 2, 4)
        plt.title('V')
        imshow_field(wavefront.V, cmap='RdBu', vmin=-I_max, vmax=I_max)
        plt.colorbar()

.. code:: ipython3

    plt.suptitle('Scalar wavefront')
    imshow_stokes_vector(scalar_wavefront)
    plt.show()



.. image:: output_12_0.png


.. code:: ipython3

    plt.suptitle('Fully polarized wavefront')
    imshow_stokes_vector(fully_polarized_wavefront)
    plt.show()



.. image:: output_13_0.png


.. code:: ipython3

    plt.suptitle('Partially polarized wavefront')
    imshow_stokes_vector(partially_polarized_wavefront)
    plt.show()



.. image:: output_14_0.png


We can also show the polarization fraction:

.. code:: ipython3

    fig, axes = plt.subplots(nrows=1, ncols=3)
    
    axes[0].set_title('Scalar wavefront')
    imshow_field(scalar_wavefront.degree_of_polarization, vmin=0, vmax=1, cmap='inferno', mask=aperture, ax=axes[0])
    
    axes[1].set_title('Fully polarized wavefront')
    imshow_field(fully_polarized_wavefront.degree_of_polarization, vmin=0, vmax=1, cmap='inferno', mask=aperture, ax=axes[1])
    
    axes[2].set_title('Partially polarized wavefront')
    imshow_field(partially_polarized_wavefront.degree_of_polarization, vmin=0, vmax=1, cmap='inferno', mask=aperture, ax=axes[2])
    
    plt.colorbar(ax=axes.tolist(), orientation="horizontal")


.. parsed-literal::

    c:\users\emiel por\github\hcipy\hcipy\optics\wavefront.py:253: RuntimeWarning: invalid value encountered in true_divide
      return np.sqrt(self.Q**2 + self.U**2 + self.V**2) / self.I
    



.. parsed-literal::

    <matplotlib.colorbar.Colorbar at 0x21e70614bc8>




.. image:: output_16_2.png


We can see that the scalar wavefront is considered an unpolarized
wavefront, the fully polarized wavefront is 100% polarized, and the
partially polarized wavefront permits a partial polarization fraction,
which is exactly as we expected.

Optical propagation with polarized wavefronts
---------------------------------------------

Optical propagations through optical elements with polarized wavefronts
are seamless. Take for example a ``FraunhoferPropagator`` propagation.

.. code:: ipython3

    focal_grid = make_focal_grid(3, 8)
    prop = FraunhoferPropagator(pupil_grid, focal_grid)
    
    scalar_image = prop(scalar_wavefront)
    fully_polarized_image = prop(fully_polarized_wavefront)
    partially_polarized_image = prop(partially_polarized_wavefront)

.. code:: ipython3

    plt.suptitle('Partially polarized wavefront')
    imshow_stokes_vector(partially_polarized_image)
    plt.show()



.. image:: output_20_0.png


As the ``FraunhoferPropagator`` does not act differently depending on
the polarization of the incoming light, we see a normal Airy core which
has the same polarization state as we put in. Note that the call to the
propagator does not change depending on the type of the ``Wavefront``
that we propagate through it. All calls look exactly the same. This
means that if we create a full optical system, consisting of many
optical elements in a row, that any ``Wavefront`` can be propgated
through, and we do not need to worry about having to rewrite the optical
system when we want to switch to (partially) polarized light.

HCIPy contains many polarization optical elements, such as half-wave
plates, quarter-wave plates, general phase retarders, polarizers,
polarizing beam splitters and more. Let’s simulate a high-contrast
imaging scenario of an unpolarized star and a partially polarized planet
next to it.

.. code:: ipython3

    stokes_vector_planet = [1, 0.1, 0, 0]
    contrast = 3e-4
    angular_separation = 5 # lambda / D
    num_photons_star = 1e9
    
    wf_star = Wavefront(aperture)
    wf_star.total_power = num_photons_star
    
    e_planet = aperture * np.exp(2j * np.pi * pupil_grid.x * angular_separation)
    wf_planet = Wavefront(e_planet, input_stokes_vector=stokes_vector_planet)
    wf_planet.total_power = num_photons_star * contrast
    
    # Calculate science image for the star and planet
    science_image = prop(wf_star).power
    science_image += prop(wf_planet).power
    
    # Simulate photon noise
    science_image = large_poisson(science_image)
    
    imshow_field(np.log10(science_image / science_image.max() + 1e-20), vmin=-5, vmax=0)
    plt.colorbar()
    plt.show()



.. image:: output_22_0.png


In this intensity image it is impossible to see the planet. It is
completely buried below the light from the star. We need to look at a
polarized light image. The simplest way, in fact to the point of
oversimplification, is to split the light with a linear polarizing beam
splitter and look at the difference between the two images. If the
planet is polarized, it will show up in the difference image.

Let’s create a linear polarizing beam splitter object, and propagate the
star and planet ``Wavefront``\ s through it. This creates two
wavefronts, one for each port of the beam splitter. The image for each
port is calculated.

.. code:: ipython3

    pbs = LinearPolarizingBeamSplitter(0)
    
    img_star_1, img_star_2 = pbs(prop(wf_star))
    img_planet_1, img_planet_2 = pbs(prop(wf_planet))
    
    science_image_1 = img_star_1.power + img_planet_1.power
    science_image_2 = img_star_2.power + img_planet_2.power
    
    science_image_1 = large_poisson(science_image_1)
    science_image_2 = large_poisson(science_image_2)

We can now show these two images and their difference.

.. code:: ipython3

    fig, axes = plt.subplots(figsize=(10, 6), nrows=1, ncols=2)
    
    axes[0].set_title('Image for port 1')
    imshow_field(np.log10(science_image_1 / science_image_1.max() + 1e-20), vmin=-5, vmax=0, ax=axes[0])
    
    axes[1].set_title('Image for port 2')
    imshow_field(np.log10(science_image_2 / science_image_2.max() + 1e-20), vmin=-5, vmax=0, ax=axes[1])
    
    plt.colorbar(ax=axes.tolist(), orientation="horizontal")
    plt.show()
    
    plt.title('Difference')
    imshow_field((science_image_1 - science_image_2) / science_image_1.max(), cmap='RdBu', vmin=-1e-4, vmax=1e-4)
    plt.annotate('planet', xy=(5, -1), xytext=(3, -7), arrowprops={'arrowstyle': '-|>'})
    plt.colorbar()
    plt.show()



.. image:: output_26_0.png



.. image:: output_26_1.png


In the difference image we can clearly see the planet light up. It has a
contrast of :math:`3\times10^{-5}`, which is partly due to the contrast
of the planet (:math:`3\times10^{-4}`) and the degree of polarization of
the planet (10%). While this is clearly a simplified example (it can
only detect the Q-component of the planet), it shows the simplicity and
power of HCIPy for propagating polarized light.

This introduction falls short of the many other optical elements that
act on polarization. Other tutorials showcase some more of these optical
elements. These include, but are not limited to, `a simple
polarimeter <../SimplePolarimeter/SimplePolarimeter.rst>`__ and `imaging
with a vector apodizing phase plate
coronagraph <../VectorApodizingPhasePlate/VectorApodizingPhasePlate.rst>`__.