# 5.2.3. Example 5: Optical elements polarimetry¶

Optical elements polarimetry consists on measuring the influence of an optical element upon a light wave polarization, i.e., its Jones or Mueller matrix of the optical element.

There are different ways of implementing a polarimeter for optical elements. In this example we will describe one of the common ones: the four-axis polarimeter, a a polarimeter composed of two rotating quarter-wave plates (linear retarders with 90º retardance), two rotating perfect polarizers and a photodetector:

Most optical element polarimeters, like this one, are usually divided into two parts: the polarization states generator (PSG) and the polarization states analyzer (PSA). The polarization states generator is the part of the polarimeter capable of generating many polarization states. In our case, it is composed by the first rotating polarizer and waveplate. If both of them are perfect elements, this PSG can generate all the possible polarization states. If the incident light wave is circularly polarized, all the generated states will have the same intensity. The second part is the polarization states analyzer (PSA). This part is the same as the light polarimeter described in the previous example.

[2]:

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

from py_pol import degrees
from py_pol.jones_vector import Jones_vector
from py_pol.jones_matrix import Jones_matrix, create_Jones_matrices
from py_pol.stokes import Stokes, create_Stokes
from py_pol.mueller import Mueller, create_Mueller


## 5.2.3.1. Jones formalism¶

During the previous example it was explained that, as the Jones formalism works using electric field and the detector measures intensity. This means that Jones optical element polarimetry can be performed using the two methods described for light polarimetry:

1. Make 7 measurements and relate them to the Jones matrix parameters.
2. Make $$N\geq7$$ measurements and use a fitting algorithm.

The first method is thoroughly described in the paper “Phase and amplitude modulation of elliptic polarization states by nonabsorbing anisotropic elements: application to liquid-crystal devices”, J. Nicolás, J. Campos and M. J. Yzuel, Vol. 19, No. 5/May 2002/J. Opt. Soc. Am. A. The second method is similar to the one described in the ligh polarimetry example and is left for the user as an interesting exercise.

## 5.2.3.2. Mueller-Stokes formalism¶

The method for determining the Mueller matrix of an optical element is similar to the one for Stokes light polarimetry. The k-th measurement of the light polarimeter experiment is:

$$S_{f}^{k}=\left[\begin{array}{cccc} A_{00}^{k} & A_{01}^{k} & A_{02}^{k} & A_{03}^{k}\\ \bullet & \bullet & \bullet & \bullet\\ \bullet & \bullet & \bullet & \bullet\\ \bullet & \bullet & \bullet & \bullet \end{array}\right]\left[\begin{array}{cccc} M_{00} & M_{01} & M_{02} & M_{03}\\ M_{10} & M_{11} & M_{12} & M_{13}\\ M_{20} & M_{21} & M_{22} & M_{23}\\ M_{30} & M_{31} & M_{32} & M_{33} \end{array}\right]\left[\begin{array}{c} S_{0}^{k}\\ S_{1}^{k}\\ S_{2}^{k}\\ S_{3}^{k} \end{array}\right]=\left[\begin{array}{c} I^{k}\\ \bullet\\ \bullet\\ \bullet \end{array}\right]$$.

Being $$A^{k}=M_{D2}(\theta_{D2}^{k})*M_{R2}(\theta_{R2}^{k})$$ the Mueller matrix of the PSA, $$S^k=M_{R1}(\theta_{R1}^{k})*M_{D1}(\theta_{D1}^{k})*S_i$$ the Stokes vector generated by the PSG and $$M_{ij}$$ the elements of the Mueller matrix of the sample. Due to the fact that the detector only measures the intensity, only some parts of the Mueller matrix of the PSA are important (the rest of the parameters are replaced by $$\bullet$$). Then, if we consider the 16 required measurements:

$$\left[\begin{array}{ccc} A_{00}^{0}S_{0}^{0} & \ldots & A_{03}^{0}S_{3}^{0}\\ \vdots & \ddots & \vdots\\ A_{00}^{16}S_{0}^{16} & \ldots & A_{03}^{16}S_{3}^{16} \end{array}\right]\left[\begin{array}{c} M_{00}\\ \vdots\\ M_{44} \end{array}\right]=\left[\begin{array}{ccc} b_{0,0} & \ldots & b_{0,16}\\ \vdots & \ddots & \vdots\\ b_{16,0} & \ldots & b_{16,16} \end{array}\right]\left[\begin{array}{c} M_{00}\\ \vdots\\ M_{44} \end{array}\right]=\left[\begin{array}{c} I^{0}\\ \vdots\\ I^{16} \end{array}\right]$$.

If $$b$$ is invertible ($$\det(b)\neq0$$) then $$M_{ij}$$ can be easily calculated as $$M_{ij}=b^{-1}I'$$. This condition is satisfied if there are no repeated set of four angles.

Let us reproduce this numerically. Let us define a function that simulates the experiment and calculates the $$b$$ matrix. Also, we will include an error amplitude which we will use in the future.

[20]:

def Mueller_polarimeter_measurement(Msample, Mpol, S_ini, angles, Ierror=0):
"""This function simulates the light polarimeter measurements.

Parameters:
Msample (Mueller): Mueller matrix of the sample.
Mpol (list): List with the four Mueller objects of the polarimeter: Md1, Mr1, Mr2, Md2.
Sini (Stokes): Stokes object of the illumination.
angles (numpy.ndarray): 4xN array with the angles used to perform the measurement.
Ierror (float): Intensity error amplitude.

Returns:
Ik: Measured intensities vector.
b: b matrix."""
# Extract the information from the input parameters
N = angles.shape[1]
Md1, Mr1, Mr2, Md2 = Mpol
# Rotate the optical elements
Md1_rotated = Md1.rotate(angle=angles[0,:], keep=True)
Mr1_rotated = Mr1.rotate(angle=angles[1,:], keep=True)
Mr2_rotated = Mr2.rotate(angle=angles[2,:], keep=True)
Md2_rotated = Md2.rotate(angle=angles[3,:], keep=True)
# Multiply the objects as required
Sg = Mr1_rotated * Md1_rotated * S_ini
A = Md2_rotated * Mr2_rotated
S_final = A * Msample * Sg
# Measure the intensity
Ik = S_final.parameters.intensity() + np.random.randn(N) * Ierror
# Calculate the b matrix
Ak = A.parameters.matrix()   # Extract the 4x4xN matrix of the Mueller object
Ak = Ak[0,:,:]   # Extract all the first rows
#     Ak = np.transpose(Ak)   # Transpose so Ak rows corrrespond to different measurements
Sk = Sg.parameters.matrix()   # Extract the 4xN matrix of the Stokes object
#     Sk = np.transpose(Sk)   # Transpose so Sk rows corrrespond to different measurements
b = np.zeros((N,16))
for ind in range(N):
b[ind,:] = np.outer(Ak[:,ind], Sk[:,ind]).flatten()
# Return
return Ik, b


Now, we can perform the experiment and see if the calculated Mueller matrix is correct. We will use as sample the product of a random diattenuator, a random polarizer and a random diagonal diattenuator. We will use a circularly polarized light wave for simplicity.

[23]:

# Start by defining the polarimeter system
S_ini = Stokes('Light source')
S_ini.circular_light(intensity=1)
Mpol = create_Mueller(('Diattenuator 1', 'Retarder 1', 'Retarder 2', 'Diattenuator 2'))
Mpol[0].diattenuator_perfect(azimuth=0)
Mpol[1].quarter_waveplate(azimuth=0)
Mpol[2].quarter_waveplate(azimuth=0)
Mpol[3].diattenuator_perfect(azimuth=0)

# Create the sample parameters
p1 = np.random.rand() * 0.5 + 0.5
p2 = np.random.rand() * 0.5
delayD = np.random.rand() * 360*degrees

R = np.random.rand() * 180*degrees
alphaR = np.random.rand() * 90*degrees
delayR = np.random.rand() * 360*degrees

d = [np.random.rand(), np.random.rand(), np.random.rand()]

# Create the sample
Md, Mr, Mp = create_Mueller(N=3)
Mr.retarder_charac_angles(R=R, alpha=alphaR, delay=delayR)
Mp.depolarizer_diagonal(d=d)

Msample = Md*Mr*Mp
Msample.name = 'Sample'

# Now, define the rotation angles for the optical elements
angles = np.random.rand(4,16) * 180*degrees

# Make the measurements
Ik, b = Mueller_polarimeter_measurement(Msample=Msample, Mpol=Mpol, S_ini=S_ini, angles=angles)

# Calculate the original Mueller matrix
b_inv = np.linalg.inv(b)
Ik = np.reshape(Ik, (16,1))
M_calc = Mueller('Calculated')
M_calc.from_matrix(M=b_inv@Ik)

# Compare the result
print(Msample, M_calc)

Sample =
[+0.444 -0.074 +0.200 -0.181]
[+0.004 +0.054 +0.042 -0.083]
[+0.041 +0.004 +0.102 +0.115]
[-0.401 +0.083 -0.211 +0.212]
Calculated =
[+0.444 -0.074 +0.200 -0.181]
[+0.004 +0.054 +0.042 -0.083]
[+0.041 +0.004 +0.102 +0.115]
[-0.401 +0.083 -0.211 +0.212]



As explained in the previous example, experiments are subject to errors. This means that the calculated Mueller matrices will be wrong. However, the error in the calculated matrix can be reduced by oversampling, i.e., taking more than the minimum 16 measurements. In this case, $$b$$ will be a Nx16 matrix and $$I'$$ a Nx1 vector. Then, $$M_{ij}$$ can be calculated using the pseudoinverse of $$b$$:

$$b^+=(b^T*b)^{-1}*b^T$$.

Let us define a function that performs one polarimetry experiment.

[41]:

def Mueller_polariemter_experiment(N, Ierror):
"""This funcion simulates a polarimetry experiment with errors in the detection.

Parameters:
N (int): Number of measurements.
Ierror (float): Intensity error amplitude.

Returns:
dif (numpy.ndarray): Difference between the calculated and the original array."""
# Create the sample parameters
p1 = np.random.rand() * 0.5 + 0.5
p2 = np.random.rand() * 0.5
delayD = np.random.rand() * 360*degrees

R = np.random.rand() * 180*degrees
alphaR = np.random.rand() * 90*degrees
delayR = np.random.rand() * 360*degrees

d = [np.random.rand(), np.random.rand(), np.random.rand()]

# Create the sample
Md, Mr, Mp = create_Mueller(N=3)
Mr.retarder_charac_angles(R=R, alpha=alphaR, delay=delayR)
Mp.depolarizer_diagonal(d=d)

Msample = Md*Mr*Mp
Msample.name = 'Sample'

# Now, define the rotation angles for the optical elements
angles = np.random.rand(4,N) * 180*degrees

# Make the measurements
Ik, b = Mueller_polarimeter_measurement(Msample=Msample, Mpol=Mpol, S_ini=S_ini, angles=angles, Ierror=Ierror)

# Calculate the Stokes vector
b_T = np.transpose(b)
b_inv = np.linalg.inv(b_T @ b)
b_inv = b_inv @ b_T
Ik = np.reshape(Ik, (N,1))
M_calc = b_inv @ Ik

# Calculate the difference
dif = np.linalg.norm(Msample.M.flatten() - M_calc.flatten())

# Return
return dif


Now, lets do the experiment many times, changing the number of measurements and the error amplitude.

[42]:

# Define some variables
Nstat = 100      # Number of experiments done under the same condition.
I_errors = np.array([0, 0.005, 0.01, 0.02])    # Array of error factors
Ns = np.arange(20, 100, 10)    # Array of the number of measurements
tol = 1e-12
I_wave = 1    # We will maintain this parameter constant so we can compare the errors easily

# Create the result variables
error = np.zeros((Ns.size, I_errors.size))
legend = []

# Make the loops changing the variables
for indI, Ierror in enumerate(I_errors):   # Loop in the error amplitude
for indN, N in enumerate(Ns):   # Loop in the number of measurements
# Create the temporal variables
aux = np.zeros((Nstat))
# Repeat the same conditions several times
for indS in range(Nstat):
# Do the experiment
dif = Mueller_polariemter_experiment(N=N, Ierror=Ierror);
# Calculate the errors
aux[indS] = dif
# Calculate the mean square error
error[indN, indI] = np.linalg.norm(aux) / Nstat
# Create the legend for the plots
legend += ['Ierror = {} %'.format(Ierror*100)]

[43]:

# Plot the results
plt.figure(figsize=(8,8))
plt.plot(Ns, error[:,:])
plt.legend(legend)
plt.xlabel('Number of measurements')
plt.ylabel('MSE')

[43]:

Text(0, 0.5, 'MSE')


These plots show that the error in the determiantion of the Stokes vectors is proportional to the initial error amplitude, and can be decreased increasing the number of measurements.