Intro to Poke: Polarization Ray Tracing

written by Jaren N. Ashcraft

The first physics submodule of poke is for polarization ray tracing (PRT). All of the physics are done in the poke.polarization module, and everything else is just ray data. PRT is an expression of the Fresnel Equations for thin-film polarization in three dimensions. This allows for the propagation of polarization-dependent performance through a ray trace of an optical system

The desireable data product is a Jones Pupil, which is the 3x3 PRT matrix rotated into a local coordinate system. Poke does this using the double pole coordinate system descibed in Chipman, Lam, and Young (2018) Chapter 11.4. This coordinate system is robust to polarization singularities that arise in using the s- and p- basis.

Initializing a Rayfront as Polarized

So you want a Jones pupil of an optical system, this section will describe how we set up optical system parameters for a Rayfront to be built. First, we start with the system properties:

• path to sequential ray trace file

• thin film information

n_film = 1.200 - 1j*7.260  # Al at 600nm from CODE V's cocatings directory
pth_to_lens = '/Users/Work/Desktop/poke/test_files/hubble_test.len'

That wasn’t too bad. Note that we only specify a thin film index, which means that the software assumes the substrate is made of entirely silver. Poke also supports thin film stacks, but we will cover that in another tutorial. Now we must specify the surface data. Poke handles surface data with dictionaries as a low-level “user interface”, and stores them in a list in the order that they appear in the raytrace.

# The Primary Mirror, 2, 4
s1 = {
    "surf" : 1,         # surface number in zemax
    "coating" : n_film, # refractive index of surface
    "mode" : "reflect"  # compute in reflection ("reflect") or transmission ("transmit")
}

# The Secondary Mirror
s2 = {
    "surf" : 2,
    "coating" : n_film,
    "mode" : "reflect"
}

# The Fold Mirror
s3 = {
    "surf" : 3,
    "coating" : n_film,
    "mode" : "reflect"
}

surflist = [s1,s2,s3]

Now that we have the surface information, we can initialize a Rayfront and pass this surface data to it. When ray tracing with CODE V files, it is presently necessary to slightly undersize the normalized pupil radius, so that we don’t get vignetting errors.

from poke.poke_core import Rayfront
from poke.poke_math import np

# rayfront parameters
number_of_rays = 20 # across the entrance pupil
wavelength = 0.6e-6
pupil_radius = 1.2 # semi-aperture of Hubble
max_field_of_view = 0.08 # degrees

rays = Rayfront(number_of_rays,wavelength,pupil_radius,max_field_of_view,circle=False)

norm fov =  [0. 0.]
base ray shape  (4, 400)

Now we have a standard Rayfront, which is just a bundle of un-traced rays now. To enable the physical optics capabilities, we need to call the Rayfront.as_polarized() method and pass it the surface list.

rays.as_polarized(surflist)

Then we can propagate it through the optical system with the rays.trace_rayset() method by supplying the path specified earlier

import zosapi
rays.trace_rayset(pth_to_lens)

res /Users/Work/Desktop/poke/test_files/hubble_test.len
CODE V warning: Warning: Buffer number 0 does not exist.  Nothing deleted.
CODE V warning: Warning: Solves may be affected by a change in the reference wavelength
global coordinate reference set to surface 1
maxrays =  400
CODE V warning: Warning: Buffer number 1 does not exist.  Nothing deleted.
1 Raysets traced through 3 surfaces

Now the rays have position and angle. This tells us a little bit about how Rayfronts are constructed. They have some attribute Rayfront._Data that holds on to the coordinate _. The following are accessible:

• xData: position in x axis

• yData: position in y axis

• zData: position in z axis

• lData: direction cosine in x axis

• mData: direction cosine in y axis

• nData: direction cosine in z axis

• l2Data: surface normal direction cosine in x axis

• m2Data: surface normal direction cosine in x axis

• n2Data: surface normal direction cosine in x axis

Each of these are numpy arrays which have shape [raybundle,surface,coordinate]. We can plot the position and direction cosines on the primary mirror:

import matplotlib.pyplot as plt
plt.figure(figsize=[9,3])
plt.subplot(121)
plt.title('Position on Primary')
plt.scatter(rays.xData[0,0],rays.yData[0,0],c=rays.zData[0,0], cmap='inferno')
plt.colorbar()
plt.subplot(122)
plt.title('Direction Cosine on Primary')
plt.scatter(rays.lData[0,0],rays.mData[0,0],c=rays.nData[0,0],cmap='RdBu')
plt.colorbar()
plt.show()

Turns out all we need is angle of incidence, direction cosines, and refractive index data to compute the polarized exit pupil. If your final axis isn’t aligned with the z-axis this is slightly more involved but for now let’s keep it simple:

rays.compute_jones_pupil(aloc=np.array([0., -1., 0.]))

And we can use poke’s built-in plotting routine to display it. Turns out Silver is pretty good in the infrared!

import poke.plotting as plot
plot.jones_pupil(rays, coordinates='polar')

Verifying the Physics

Now we have illustrated that we can return a Jones pupil, but is it accurate? To check this, we can compare the results to the ouput of POLGRID.seq. This is a CODE V macro that does the exact same polarization ray tracing the Poke does. It has a macro-only interface (as far as I know), so we will use the following CODE V seqence to run POLGRID.

! Define input variables
num ^input_array(2,400) ^output_array(43,401) ^input_ray(4)
num ^success
^nrays == 400
^nrays_across == sqrt(^nrays)
^x_start == -1
^x_end == 1
^y_start == -1
^y_end == 1
! compute step size and construct input array
^step_size == absf(^y_start-^y_end)/(^nrays_across-1)
^x_cord == ^x_start
^y_cord == ^y_start
^next_row == 0
FOR ^iter 1 ^nrays
    ^input_array(1,^iter) == ^x_cord
    ^input_array(2,^iter) == ^y_cord
    ! update the x_cord
    ^x_cord == ^x_cord + ^step_size
    IF modf(^iter,^nrays_across) = 0
        ^next_row == 1
    END IF
    IF ^next_row = 1
        ! update y_cord
        ^y_cord == ^y_cord + ^step_size
        ! reset x_cord
        ^x_cord == ^x_start
        ! reset next_row
        ^next_row == 0
    END IF
END FOR
! Run POLGRID
^success == POLGRID(1, 1, 1, 0, ^nrays, ^input_array, ^output_array)
wri ^success
! Write output_array to text file
BUF DEL B1
^result == ARRAY_TO_BUFFER(^output_array,1,0)
BUF EXP B1 SEP "polgrid_output.txt"

Before executing this, we will need to set up our coated HST for comparison. We will use the REFL_AL_450nm_700nm.mul coating for comparison, which comes installed with CODE V in the “coating” directory. We then save the macro above as run_polgrid.seq to perform the polarization ray trace. You can find both the coated HST model and the polgrid macro in the experiments/physics_validation/ directory of poke.

# Compare to the polgrid output
from poke.poke_math import np
pth = '../../experiments/physics_validation/polgrid_output.txt'
polgrid_data = np.genfromtxt(pth, delimiter=' ')

# jones matrix indices, see CODE V docs for the remaining rows
jxx_r = polgrid_data[20,1:]
jxx_i = polgrid_data[21,1:]
jxy_r = polgrid_data[22,1:]
jxy_i = polgrid_data[23,1:]
jyx_r = polgrid_data[24,1:]
jyx_i = polgrid_data[25,1:]
jyy_r = polgrid_data[26,1:]
jyy_i = polgrid_data[27,1:]

Jmat = np.zeros([jxx_r.shape[0], 2, 2], dtype=np.complex128)
Jmat[..., 0, 0] = jxx_r + 1j*jxx_i
Jmat[..., 0, 1] = jxy_r + 1j*jxy_i
Jmat[..., 1, 0] = jyx_r + 1j*jyx_i
Jmat[..., 1, 1] = jyy_r + 1j*jyy_i

x, y = rays.xData[0, 0], rays.yData[0, 0]

fig, axs = plt.subplots(figsize=[12, 5], nrows=2, ncols=4)
plt.suptitle("Jones Pupil - CODE V")
for j in range(2):
    for k in range(2):
        ax = axs[j, k]
        ax.set_title("|J{j}{k}|".format(j=j, k=k))
        sca = ax.scatter(x, y, c=np.abs(Jmat[..., j, k]), cmap="inferno")
        fig.colorbar(sca, ax=ax)

        # turn off the ticks
        if j != 1:
            ax.xaxis.set_visible(False)
        if k != 0:
            ax.yaxis.set_visible(False)

# theres a phase offset
offset = 0

for j in range(2):
    for k in range(2):

        if j == 1:
            if k == 1:
                offset = 0
        ax = axs[j, k + 2]
        ax.set_title(r"$\angle$" + "J{j}{k}".format(j=j, k=k))
        sca = ax.scatter(x, y, c=np.angle(Jmat[..., j, k]) + offset, cmap="coolwarm")
        fig.colorbar(sca, ax=ax)

        # turn off the ticks
        if j != 1:
            ax.xaxis.set_visible(False)

        ax.yaxis.set_visible(False)
plt.show()

poke_jones = rays.jones_pupil[-1]
fig, axs = plt.subplots(figsize=[12, 5], nrows=2, ncols=4)
plt.suptitle("Jones Pupil - CODE V")
for j in range(2):
    for k in range(2):
        ax = axs[j, k]
        ax.set_title("|J{j}{k}|".format(j=j, k=k))
        sca = ax.scatter(x, y, c=np.abs(Jmat[..., j, k]) - np.abs(poke_jones[..., j, k]), cmap="inferno")
        fig.colorbar(sca, ax=ax)

        # turn off the ticks
        if j != 1:
            ax.xaxis.set_visible(False)
        if k != 0:
            ax.yaxis.set_visible(False)

# theres a phase offset
offset = 0

for j in range(2):
    for k in range(2):

        if k == 1:
            offset = 0
        else:
            offset = 0
        ax = axs[j, k + 2]
        ax.set_title(r"$\angle$" + "J{j}{k}".format(j=j, k=k))
        sca = ax.scatter(x, y, c=np.angle(Jmat[..., j, k]) + offset - np.angle(poke_jones[..., j, k]), cmap="coolwarm")
        fig.colorbar(sca, ax=ax)

        # turn off the ticks
        if j != 1:
            ax.xaxis.set_visible(False)

        ax.yaxis.set_visible(False)
plt.show()

The differences in simulation are dominated by a couple numerical factors.

• The on-diagonal amplitudes appear limited by some geometric transformation at the \(10^{-5}\%\) level. This discrepancy can arise from the differences in local coordinate system chosen to conduct the PRT matrix to Jones pupil conversion

• The on-diagonal phases are limited by numerical precision, given the order of magnitude of the difference and apparent randomness of the data

• The off-diagonal phases and amplitudes appear limited by the precise location of the zero point in amplitude. This is not due to the PRT algorithm, but more a difference in how POLGRID and Poke define the initial rayset. What we see here are rays landing in slightly different locations.

Further physics validation experiments are conducted in the poke/experiments/physics_validation/ directory, if you would like to view them or play around yourself!