# !/usr/bin/env python3
# -*- coding: utf-8 -*-
# ------------------------------------
# Authors: Luis Miguel Sanchez Brea and Jesus del Hoyo
# Date: 2019/01/09 (version 1.0)
# License: GPL
# -------------------------------------
"""
Mueller objects describe optical polarization elements in the Mueller-Stokes formalism.
**Class fields:**
* **M**: 4x4xN array containing all the Mueller matrices.
* **name**: Name of the object for print purposes.
* **shape**: Shape desired for the outputs.
* **size**: Number of stored Mueller matrices.
* **ndim**: Number of dimensions for representation purposes.
* **no_rotation**: If True, rotation method do not act upon the object. Useful for objects that shouldn't be rotated as mirrors.
* **type**: Type of the object ('Mueller'). This is used for determining the object class as using isinstance may throw unexpected results in .ipynb files.
* **parameters**: parameters of the Mueller matrices.
* **checks**: checks of the Mueller matrices.
* **analysis**: analysis of the Mueller matrices.
**Parent methods**
* **clear**: Removes data and name form Jones vector.
* **copy**: Creates a copy of the Jones_vector object.
* **stretch**: Stretches a Jones vector of size 1.
* **shape_like**: Takes the shape of another object to use as its own.
* **reshape**: Changes the shape of the object.
* **flatten**: Transforms N-D objects into 1-D objects (0-D if only 1 element).
* **flip**: Flips the object along some dimensions.
* **get_list**: Creates a list with single elements.
* **from_list**: Creates the object from a list of single elements.
* **concatenate**: Canocatenates several objects into a single one.
* **draw**: Draws the components of the object.
* **clear**: Clears the information of the object.
**Generation methods**
* **from_components**: Creates a Mueller matrix directly from the 16 $M_{ij}$ elements.
* **from_matrix**: Creates a Mueller object directly from a 4x4xN matrix.
* **from_normalized**: Creates a Mueller matrix directly from a normalized 4x4 matrix ($M_{norm} = M/M_{00}$).
* **from_Jones**: Creates a Mueller Matrix equivalent to a Jones matrix.
* **from_blocks**: Creates a Mueller matrix from the blocks of its decomposition.
* **from_covariance**: Creates a Mueller matrix from the equivalent covariant matrix.
* **from_inverse**: Creates a Mueller matrix from the inverse matrix.
* **from_list**: Creates a Jones_matrix object directly from a list of 4x4 numpy arrays.
* **vacuum**: Creates the matrix for vacuum.
* **mirror**: Creates the matrix for a mirror. NOTE: This matrix mus not be rotated.
* **filter_amplifier**: Creates the matrix for a neutral filter or amplifier element.
* **depolarizer_perfect**: Creates a perfect depolarizer.
* **depolarizer_diagonal**: Creates a depolarizer with elements just in the diagonal.
* **depolarizer_states**: Creates a general depolarizer from the diattenuation, polarizance and eigenstate vectors.
* **diattenuator_perfect**: Creates a perfect linear polarizer.
* **diattenuator_linear**: Creates a real diattenuator with perpendicular axes.
* **diattenuator_charac_angles**: Creates the most general homogeneous diattenuator with orthogonal eigenstates from the characteristic angles of the main eigenstate.
* **diattenuator_azimuth_ellipticity**: Creates the most general homogenous diattenuator from the characteristic angles of the main eigenstate.
* **diattenuator_vector**: Creates the most general homogenous diattenuator from the diattenuation vector.
* **quarter_waveplate**: Creates a perfect retarder with 90º retardance.
* **half_waveplate**: Creates a perfect retarder with 180º retardance.
* **retarder_linear**: Creates a linear retarder.
* **retarder_charac_angles**: Creates the most general homogeneous retarder from the characteristic angles of the fast eigenstate.
* **retarder_azimuth_ellipticity**: Creates the most general homogeneous retarder from the characteristic angles of the fast eigenstate.
* **retarder_from_vector**: Creates the most general homogeneous retarder from the retardance vector.
* **diattenuator_retarder_linear**: Creates an homogeneous linear diattenuator retarder with the same axes for diattenuation and retardance.
* **diattenuator_retarder_azimuth_ellipticity**: Creates the most general homogeneous diattenuator retarder with the same axes for diattenuation and retardance from the azimuth and ellipticity angle.
* **diattenuator_retarder_charac_angles**: Creates the most general homogeneous diattenuator retarder with the same axes for diattenuation and retardance from the characteristic angles.
* **general_eigenstates**: Generates the most general pure optical element from its eigenstates.
**Manipulation methods**
* **rotate**: Rotates the Mueller matrix.
* **sum**: Calculates the summatory of the Jones matrices in the object.
* **prod**: Calculates the product of the Jones matrices in the object.
* **remove_global_phase**: Removes the phase introduced by the optical element.
* **add_global_phase**: Increases the phase introduced by the optical element.
* **set_global_phase**: Sets the phase introduced by the optical element.
* **reciprocal**: Flips the optical element so the light transverses it in the opposite direction.
* **transpose**: Transposes the Mueller matrix of the element.
* **inverse**: Calculates the inverse matrix of the Mueller matrix.
* **covariant_matrix**: This method calculates the covariant matrix of the Mueller matrix of the object.
**Parameters subclass methods**
* **matrix**: Gets a numpy array with all the matrices.
* **components**: Extracts the four components of the Mueller matrix.
* **global_phase**: Extracts the global phase introduced by the object.
* **blocks**: Method that divides a mueller matrix in their blocks: mean transmission ($M_{00}$), diattenuation and polarizance vectors and small matrix m.
* **diattenuation_vector**: Extracts the 3xN array of diattenuation vectors.
* **polarizance_vector**: Extracts the 3xN array of polarizance vectors.
* **small_matrix**: Extracts the 3x3xN array of small matrix m.
* **retardance_vector**: Extracts the 3xN array of retardance vectors (if exists).
* **mean_transmission**: Calculates the mean transmission coefficient.
* **transmissions**: Calculates the maximum and minimum transmissions.
* **inhomogeneity**: Calculates the inhomogeneity parameter.
* **diattenuation**: Calculates the diattenuation of a Mueller matrix.
* **diattenuation_linear**: Calculates the linear diattenuation of a Mueller matrix.
* **diattenuation_circular**: Calculates the circular diattenuation of a Mueller matrix.
* **polarizance**: Calculates the polarizance of a Mueller matrix.
* **polarizance_linear**: Calculates the linear polarizance of a Mueller matrix.
* **polarizance_circular**: Calculates the delay of the matrix.
* **degree_polarizance**: Calculates the degree of polarizance.
* **spheric_purity**: Calculates the spheric purity grade.
* **retardance**: Calculates the retardance (also refered as delay) of the Mueller matrix of a pure retarder.
* **polarimetric_purity**: Calculates the degree of polarimetric purity of a Mueller matrix.
* **depolarization_index**: Calculates the depolarization_index of a Mueller matrix.
* **polarimetric_purity_indices**: Calculates the polarimetric purity indices of a Mueller matrix.
* **eig**: Calculates the eigenvalues and eigenstates (eigenvectors) of the Mueller matrices.
* **eigenvalues**: Calculates the eigenvalues and of the Mueller matrices.
* **eigenvectors**: Calculates the eigenvectors of the Mueller matrices.
* **eigenstates**: Calculates the eigenstates (Stokes vectors of the eigenvectors) of the Mueller matrices.
* **det**: Calculates the determinant and of the Mueller matrices.
* **trace**: Calculates the trace of the Mueller matrices.
* **norm**: Calculates the norm of the Mueller matrices.
* **get_all**: Returns a dictionary with all the parameters of the object.
**Checks subclass methods**
* **is_physical**: Conditions of physical realizability.
* **is_non_depolarizing / is_pure**: Checks if matrix is non-depolarizing.
* **is_homogeneous**: Checks if the matrix is homogeneous (eigenstates are orthogonal). It is implemented in two different ways.
* **is_retarder**: Checks if the matrix M corresponds to a pure retarder.
* **is_diattenuator / is_polarizer**: Checks if the matrix M corresponds to a pure homogeneous diattenuator.
* **is_depolarizer**: Checks if the object corresponds to a depolarizer.
* **is_singular**: Checks if the matrix is singular (at least one of its eigenvalues is 0).
* **is_symmetric**: Checks if the Mueller matrices are symmetric.
* **is_eigenstate**: Checks if a given light state is an eigenstate of the objct.
* **get_all**: Returns a dictionary with all the checks of the object.
**Analysis subclass methods**
* **diattenuator**: Calculates all the parameters from the Mueller Matrix of a diattenuator.
* **polarizer**: Calculates all the parameters from the Mueller Matrix of a diattenuator using the polarizance vector. If the polarizer is homogeneous, this is equivalent to the previous method.
* **retarder**: Calculates all the parameters from the Mueller Matrix of a retarder.
* **depolarizer**: Calculates some of the parameters from the Mueller matrix of a diattenuator.
* **filter_physical_conditions**: Method that filters experimental errors by forcing the Mueller matrix M to fulfill the conditions necessary for a matrix to be physicall.
* **filter_purify_number**: Purifies a Mueller matrix by choosing the number of eigenvalues of the covariant matrix that will be made 0.
* **filter_purify_threshold**: Purifies a Mueller matrix by making 0 the eigenvalues of the covariant matrix lower than a certain threshold.
* **decompose_pure**: Polar decomposition of a pure Mueller matrix in a retarder and a diattenuator.
* **decompose_polar**: Polar decomposition of a general Mueller matrix in a depolarizer, retarder and a diattenuator.
"""
from .utils import *
from .py_pol import Py_pol
from .stokes import Stokes, create_Stokes
from .jones_matrix import Jones_matrix
from .jones_vector import Jones_vector
from . import degrees, eps, num_decimals, number_types
from sympy.functions.special.tensor_functions import Eijk
from numpy.linalg import inv
from numpy import arctan2, array, cos, exp, matrix, pi, sin, sqrt
import numpy as np
from functools import wraps
import warnings
from copy import deepcopy
warnings.filterwarnings('ignore')
tol_default = eps
counter_max = 20
N_print_list = 5
print_list_spaces = 3
empty_matrix = np.zeros((4, 4, 1), dtype=float)
change_names = True
tol_default = eps
tol_default_filter = 1e-10
unknown_phase = False
default_phase = 0
zero_D = np.zeros(3)
zero_m = np.zeros((3, 3))
# Create a list with the base of matrices and its kronecker product
S = [
np.eye(2),
array([[1, 0], [0, -1]]),
array([[0, 1], [1, 0]]),
array([[0, -1j], [1j, 0]])
]
S_kron = []
for i in range(4):
for j in range(4):
S_kron.append(np.kron(S[i], np.conj(S[j])))
#############################################################################
# Methods
#############################################################################
[docs]
def create_Mueller(name='M', N=1, out_object=True):
"""Method that creates several Mueller_matrix objects at the same time from a list of names or a number.
Parameters:
name (string or list): Name of the object for print purposes. Default: 'M'.
N (int): Number of created elements. Default: 1.
out_object (bool): If N=1 and out_object is True the output is a Jones_matrix instead of a list. Default: True.
Attributes:
self.parameters (class): Class containing the measurable parameters of the Mueller matrices.
self.checks (class): Class containing the methods that check something about the Mueller matrices.
Returns:
J (list or Mueller_matrix): Result.
"""
J = []
if isinstance(name, list) or isinstance(name, tuple):
for n in name:
J.append(Mueller(n))
else:
for _ in range(N):
J.append(Mueller(name))
if len(J) == 1 and out_object:
J = J[0]
return J
[docs]
def set_printoptions(N_list=None, list_spaces=None):
"""Method that modifies the global print options parameters.
Parameters:
N_print_list (int): Number of matrices that will be printed as a list if the shape of the object is 1D. Default: None
print_list_spaces (int): Number of spaces between matrices if they are printed as a list. Default: None
"""
global N_print_list, print_list_spaces
if list_spaces is not None:
print_list_spaces = list_spaces
if N_list is not None:
N_print_list = N_list
################################################################################
# Main class
################################################################################
[docs]
class Mueller(Py_pol):
"""Class for Mueller matrices
Parameters:
name (str): name of Mueller matrix, for string representation
Attributes:
M (np.ndarray): 4x4xN array of floats containing all the Mueller matrices.
m00 (np.ndarray): 1xN array containing the $M_{00}$ element of all Mueller matrices.
D (np.ndarray): 3xN array containing the diattenuation vectors of the Mueller matrices.
P (np.ndarray): 3xN array containing the polarizance vectors of the Mueller matrices.
m (np.ndarray): 3xN array containing the rest of the Mueller matrices.
global_phase (np.ndarray): 1xN array storing the global phase introduced by the optical objects.
name (string): Name of the object for print purposes.
shape (tuple or list): Shape desired for the outputs.
size (int): Number of stored Mueller matrices.
ndim (int): Number of dimensions for representation purposes.
no_rotation (bool): If True, rotation method do not act upon the object. Useful for objects that shouldn't be rotated as mirrors.
type (string): Type of the object ('Jones_matrix'). This is used for determining the object class as using isinstance may throw unexpected results in .ipynb files.
Attributes:
self.parameters (class): parameters of the Mueller matrices.
self.checks (class): checks of the Mueller matrices.
self.analysis (class): analysis of the Mueller matrices.
"""
__array_priority__ = 30000
############################################################################
# Operations
############################################################################
def __init__(self, name='M'):
super().__init__(name=name, _class="Mueller")
self.no_rotation = False
self._global_phase = default_phase
self.parameters = Parameters_Mueller(self)
self.analysis = Analysis_Mueller(self)
self.checks = Check_Mueller(self)
def __add__(self, other):
"""Adds two Mueller matrices.
Parameters:
other (Mueller or Jones_vector): 2nd matrix to add.
Returns:
M2 (Mueller): Result.
"""
try:
if other.type in ('Jones_matrix', 'Mueller'):
# Transform other to Mueller if necessary
M2 = Mueller()
if other.type == 'Jones_matrix':
M1 = Mueller(other.name)
M1.from_Jones(other)
else:
M1 = other
# Calculate the new global phase
if self.global_phase is None or M1.global_phase is None:
phase = None
elif np.all(self.global_phase == M1.global_phase):
phase = self.global_phase
else:
phase = None
# Calculate the new matrix
M2.from_matrix(self.M + M1.M, global_phase=phase)
M2.shape = take_shape((self, M1))
if change_names:
M2.name = self.name + " + " + M1.name
return M2
else:
raise ValueError('other is {} instead of Jones_matrix.'.format(
other.type))
except:
raise ValueError('other is not a py_pol object but {}'.format(
type(other)))
def __sub__(self, other):
"""Substracts two Mueller matrices.
Parameters:
other (Mueller or Jones_vector): 2nd matrix to substract.
Returns:
M3 (Mueller): Result.
"""
M3 = self + ((-1) * other)
if change_names:
M3.name = self.name + " - " + other.name
return M3
def __rmul__(self, other):
"""Multiplies a Mueller matrix by a number. If the number is complex or real negative, the absolute value is used and the global phase is updated acordingly.
Parameters:
other (float, complex or numpy.ndarray): number to multiply.
Returns:
(Mueller): Result.
"""
M3 = Mueller()
# Save the Number of elements, and then flatten
if isinstance(other, number_types):
N = 1
other2 = other
elif isinstance(other, np.ndarray):
N = other.size
other2 = other.flatten()
else:
raise ValueError('Other is not a number or a numpy array')
# Calculate components
components = self.parameters.components(shape=False)
# Check that the multiplication can be performed
if N == self.size or self.size == 1 or N == 1:
# Calculate the absolute value and complex phase of the number
mod, phase = (np.abs(other2), np.angle(other2))
for ind in range(16):
components[ind] = components[ind] * mod
# Create the object
M3.from_components(components, global_phase=self.global_phase)
M3.add_global_phase(phase)
if isinstance(other, np.ndarray):
M3.shape = take_shape((self, other.shape))
elif isinstance(other, number_types):
M3.shape = self.shape
else:
raise ValueError(
'The number of elements in other ({}) and {} ({}) is not the same'
.format(N, self.name, self.size))
if isinstance(other, number_types) and change_names:
M3.name = str(other) + " * " + self.name
else:
M3.name = str(other)
return M3
def __mul__(self, other):
"""Multiplies the Mueller matrix by a number, an array of numbers, a Stokes or Jones vector, or another Mueller or Jones matrix.
Parameters:
other (float, numpy.ndarray, Stokes, Jones_vector, Mueller or Jones_matrix): 2nd object to multiply.
Returns:
(Mueller): Result.
"""
# Multiplication by numbers or arrays is already implemented in __mul__
if isinstance(other, number_types) or isinstance(other, np.ndarray):
return other * self
# Multiply by py_pol objects
else:
other_Jones = None
# try:
# Transform Jones vectors if required
if other.type == 'Jones_vector':
new_other = Stokes()
new_other.from_Jones(other)
other_Jones = other
elif other.type == 'Jones_matrix':
new_other = Mueller()
new_other.from_Jones(other)
other_Jones = other
else:
new_other = other
try:
if other.type == 'Mueller':
other_Jones = Jones_matrix().from_Mueller(other)
else:
other_Jones = Jones_vector().from_Stokes(other)
except:
other_Jones = Jones_matrix().vacuum()
try:
self_Jones = Jones_matrix().from_Mueller(self)
except:
self_Jones = Jones_matrix().vacuum()
# Prepare variables
new_self, new_other = expand_objects([self, new_other], copy=True)
if new_other.type == 'Stokes':
M3 = Stokes()
elif new_other.type == 'Mueller':
M3 = Mueller()
else:
raise ValueError(
'other is not a correct py_pol object, but {}.'.format(
other.type))
# Multiply
Mf = matmul_pypol(new_self.M, new_other.M)
Jf = self_Jones * other_Jones
phase = Jf.parameters.global_phase()
M3.from_matrix(Mf, global_phase=phase)
#M3.from_matrix(Mf)
if change_names:
M3.name = self.name + " * " + other.name
# except:
# raise ValueError(
# 'other is not number, numpy.ndarray or py_pol object, but {}.'.format(type(other)))
M3.shape = take_shape((self, other))
if isinstance(other, number_types) and change_names:
M3.name = str(other) + " * " + self.name
return M3
def __truediv__(self, other):
"""Divides a Mueller matrix by a number. If the number is complex or real negative, the absolute value is used and the global phase is updated acordingly.
Parameters:
other (float or numpy.ndarray): Divisor.
Returns:
(Mueller): Result.
"""
M3 = self * (other**(-1))
if isinstance(other, number_types) and change_names:
M3.name = self.name + " / " + str(other)
return M3
def __repr__(self):
"""prints information about class"""
# Extract the components
components = self.parameters.components()
# If the object is empty, say it
if self.size == 0 or np.all(self.M == 0):
return '{} is empty\n'.format(self.name)
# If the object is 0D or 1D, print it like a list or inline
elif self.size == 1 or self.shape is None or self.ndim < 2:
# Short enough objects can be printed as a list of matrices
if self.size <= N_print_list:
list = self.get_list(out_number=False)
str = "{} = \n".format(self.name)
str = str + PrintMatrices(list, print_list_spaces)
# Print the rest using a single line for each component
else:
str = "{} M00 = {}\n".format(self.name, components[0])
for ind1 in range(4):
for ind2 in range(4):
if ind1 + ind2 > 0:
str = str + " " * \
len(self.name) + " M{}{} = {}\n".format(ind1,
ind2, components[ind1 * 4 + ind2])
# Print higher dimensionality as pure arrays
else:
str = ""
for ind1 in range(4):
for ind2 in range(4):
str = str + \
"{} M{}{} =\n {}\n".format(
self.name, ind1, ind2, components[ind1 * 4 + ind2])
return str
def __getitem__(self, index):
"""
Implements object extraction from indices.
"""
if change_names:
E = Mueller(self.name + '_picked')
else:
E = Mueller(self.name)
# If the indices are 1D, act upon the matrix directly
if isinstance(index, (int, slice)) and self.ndim > 1:
E.from_matrix(self.M[:, :, index])
elif isinstance(index,
np.ndarray) and index.ndim == 1 and self.ndim > 1:
E.from_matrix(self.M[:, :, index])
# If not, act upon the components
else:
components = self.parameters.components(out_number=False)
for ind in range(16):
components[ind] = components[ind][index]
E.from_components(components)
return E
def __setitem__(self, index, data):
"""
Implements object inclusion from indices.
"""
# Check that data is a correct pypol object
if data.type == 'Jones_matrix':
data2 = Mueller(data.name)
data2.from_Jones(data)
elif data.type == 'Mueller':
data2 = data
else:
raise ValueError(
'data is type {} instead of Jones_vector or Stokes.'.format(
data.type))
# Expand phase if required
if self.global_phase is None:
self.global_phase = np.zeros(self.size) * np.nan
if data2.global_phase is None:
data2.global_phase = np.nan
# If the indices are 1D, act upon the matrix directly
if isinstance(index, int) and self.ndim > 1:
self.M[:, :, index] = np.squeeze(data2.M)
# Add global phase
self.global_phase[index] = data2.global_phase
elif isinstance(index, slice) and self.ndim > 1:
if data2.size == 1:
if index.step is None:
step = 1
else:
step = index.step
N = int((index.stop - index.start) / step)
data3 = data2.stretch(length=N, keep=True)
else:
data3 = data2
self.M[:, :, index] = np.squeeze(data3.M)
# Add global phase
self.global_phase[index] = data3.global_phase
elif isinstance(index,
np.ndarray) and index.ndim == 1 and self.ndim > 1:
self.M[:, :, index] = data2.M
# If not, act upon the components
else:
# Extract phase and components
components = self.parameters.components(out_number=False)
phase = self.parameters.global_phase(out_number=False)
components_new = data2.parameters.components(out_number=False)
for ind in range(16):
components[ind][index] = components_new[ind]
phase_new = data2.parameters.global_phase(out_number=False)
phase[index] = phase_new
# Set the new values
self.from_components(components, global_phase=phase)
def __eq__(self, other):
"""
Implements equality operation.
"""
try:
# Calculate the difference object
if other.type == 'Jones_matrix':
M2 = Mueller()
M2.from_Jones(other)
elif other.type == 'Mueller':
M2 = other
else:
return False
except:
return False
# Stretch if required
if self.size == 1 and M2.size > 1:
M1 = self.stretch(length=M2.size, keep=True)
else:
M1 = self
if M1.size > 1 and M2.size == 1:
M2 = M2.stretch(length=M1.size, keep=True)
shape = take_shape((self, other))
# Compare matrices
norm = np.linalg.norm(M1.M - M2.M, axis=(0, 1))
cond1 = norm < tol_default
if shape is not None:
cond1 = np.reshape(cond1, shape)
# Compare phases
if (M1.global_phase is None) or (M1.global_phase is None):
cond2 = np.nan
else:
phase1 = M1.parameters.global_phase(shape=shape)
phase2 = M2.parameters.global_phase(shape=shape)
cond2 = phase1 == phase2
# Merge conditions
cond = cond1 * cond2
return cond
################
## PROPERTIES
################
@property
def global_phase(self):
return self._global_phase
@global_phase.setter
def global_phase(self, value):
if value is None or isinstance(value, number_types):
self._global_phase = value
elif value.size == 1:
self._global_phase = np.ones(self.size) * value
elif value.size != self.size:
raise ValueError("Global phase has diffrent size ({}) than object ({})".format(value.size, self.size))
else:
self._global_phase = value
#####################
## MANIPULATION
#####################
[docs]
def add_global_phase(self,
phase=0,
unknown_as_zero=unknown_phase,
keep=False):
"""Method that adds a phase to the Mueller object.
Parameters:
phase (float or np.ndarray): Phase to be added to the Mueller object. Default: 0.
unknown_as_zero (bool): If True, takes unknown phase as zero. Default: False.
keep (bool): If True, self is not updated. Default: False.
Returns:
(Mueller): Recalculated Mueller object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Prepare variables
phase, new_obj, new_shape = prepare_variables([phase],
expand=[True],
obj=new_obj,
give_shape=True)
phase = phase.flatten()
# Add the phase
if self.global_phase is None:
if unknown_as_zero:
self.global_phase = phase
else:
self.global_phase = self.global_phase + phase
# End
new_obj.shape, _ = select_shape(new_obj, new_shape)
return new_obj
[docs]
def set_global_phase(self,
phase=0,
keep=False,
shape_like=None,
shape=None):
"""Method that sets the phase to the Mueller object.
Parameters:
phase (float or np.ndarray): Phase to be added to the Mueller object. Default: 0.
keep (bool): If True, self is not updated. Default: False.
shape_like (float or numpy.ndarray): Use the shape of this array. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Recalculated Mueller object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# If None or Nan, skip some steps
if phase is None or np.all(np.isnan(phase)):
self.phase = None
else:
# Prepare variables
phase, new_obj, new_shape = prepare_variables([phase],
expand=[True],
obj=new_obj,
give_shape=True)
phase = phase.flatten()
new_obj.shape, _ = select_shape(new_obj,
shape_var=new_shape,
shape_fun=shape,
shape_like=shape_like)
self.global_phase = phase
# End
return new_obj
[docs]
def remove_global_phase(self, keep=False):
"""Method that removes the phase to the Mueller object.
Parameters:
keep (bool): If True, self is not updated. Default: False.
Returns:
(Mueller): Recalculated Mueller object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Set the phase
self.global_phase = default_phase
# End
return new_obj
# @_actualize_
[docs]
def rotate(self, angle=0, keep=False, change_name=change_names):
"""Rotates a Mueller matrix a certain angle
M_rotated = R(-angle) * self.M * R(angle)
Parameters:
angle (float): angle of rotation in radians. Default: 0
keep (bool): If True, the original element is not updated. Default: False.
change_name (bool): If True and angle is of size 1, changes the object name adding @ XX deg, being XX the total rotation angle. Default: True.
Returns:
(Mueller): Mueller object rotated
"""
if self.no_rotation:
print('Warning: Tried to rotate {}, which must not be rotated.'.
format(self.name))
return self
else:
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Prepare variables
angle, new_obj, new_shape = prepare_variables([angle],
expand=[True],
obj=new_obj,
give_shape=True)
# Calculate the rotation objects
Jneg, Jpos = create_Mueller(('-', '+'))
Jneg.from_matrix(rotation_matrix_Mueller(-angle))
Jpos.from_matrix(rotation_matrix_Mueller(angle))
# Rotate
other = Jneg * (new_obj * Jpos)
new_obj.from_matrix(other.M)
# Update name
if change_name and angle.size == 1:
if np.abs(angle) > tol_default:
new_obj.name = new_obj.name + \
" @ {:1.2f} deg".format(angle[0] / degrees)
# Return
new_obj.shape, _ = select_shape(obj=new_obj,
shape_var=new_shape)
return new_obj
[docs]
def covariance_matrix(self,
keep=False,
shape_like=None,
shape=None,
change_name=change_names):
"""Calculates the covariance matrix of a Mueller matrix.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016) pp 171.
Notes:
The base of matrices S is used in an uncommon order.
In order to obtain the same result as in the book, the formula must be:
.. math:: H=0.25\sum(m[i,j]\,kron([S(i),S^{*}(j))].
Parameters:
keep (bool): If True, the original element is not updated. Default: False.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
change_name (bool): If True, changes the object name adding Recip. of at the beggining of the name. Default: True.
Returns:
(Mueller): Modified object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Calculate components
components = self.parameters.components(shape=False, out_number=False)
# Calculate the covariance matrices
H = np.zeros((4, 4, new_obj.size), dtype=complex)
for ind in range(16):
H = H + np.multiply.outer(S_kron[ind], components[ind])
new_obj.from_matrix(H / 4)
new_obj.shape, _ = select_shape(shape_var=self.shape,
shape_like=shape_like,
shape_fun=shape)
# Fix the name if required
if change_name:
new_obj.name = 'Covariant of ' + new_obj.name
return new_obj
[docs]
def inverse(self,
keep=False,
shape_like=None,
shape=None,
change_name=change_names):
"""Calculates the inverse matrix of the Mueller matrix.
Parameters:
keep (bool): If True, the original element is not updated. Default: False.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
change_name (bool): If True, changes the object name adding Recip. of at the beggining of the name. Default: True.
Returns:
(Mueller): Modified object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Caluclate inverse
new_obj.from_matrix(inv_pypol(self.M),
shape_like=shape_like,
shape=shape)
new_obj.shape, _ = select_shape(shape_var=self.shape,
shape_like=shape_like,
shape_fun=shape)
# Fix the name if required
if change_name:
new_obj.name = 'Inverse of ' + new_obj.name
return new_obj
# @_actualize_
[docs]
def reciprocal(self,
keep=False,
shape_like=None,
shape=None,
change_name=change_names):
"""Calculates the reciprocal of the optical element, so the light tranverses it in the opposite direction. In Mueller formalism, it is calculated as:
.. math:: M^{r}=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & -1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]M^{T}\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & -1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 111.
Parameters:
keep (bool): If True, the original element is not updated. Default: False.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
change_name (bool): If True, changes the object name adding Recip. of at the beggining of the name. Default: True.
Returns:
(Mueller): Result.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Calculate components
components = new_obj.transpose(keep=True).parameters.components(
shape=False)
components[2] = -components[2]
components[6] = -components[6]
components[8] = -components[8]
components[9] = -components[9]
components[11] = -components[11]
components[14] = -components[14]
# Create the object
new_obj.from_components(components)
new_obj.shape, _ = select_shape(shape_var=self.shape,
shape_like=shape_like,
shape_fun=shape)
# Fix the name if required
if change_name:
new_obj.name = 'Reciprocal of ' + new_obj.name
return new_obj
[docs]
def transpose(self,
keep=False,
shape_like=None,
shape=None,
change_name=change_names):
"""Calculates the transposed matrices of the Mueller matrices.
Parameters:
keep (bool): If True, the original element is not updated. Default: False.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
change_name (bool): If True, changes the object name adding Recip. of at the beggining of the name. Default: True.
Returns:
(Jones_matrix): Modified object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Caluclate inverse
new_obj.from_matrix(np.transpose(new_obj.M, axes=(1, 0, 2)))
new_obj.shape, _ = select_shape(shape_var=self.shape,
shape_like=shape_like,
shape_fun=shape)
# Fix the name if required
if change_name:
new_obj.name = 'Transpose of ' + new_obj.name
return new_obj
[docs]
def sum(self, axis=None, keep=False, change_name=change_names):
"""Calculates the sum of Mueller matrices stored in the object.
Parameters:
axis (int, list or tuple): Axes along which the summatory is performed. If None, all matrices are summed. Default: None
keep (bool): If True, the original element is not updated. Default: False.
change_name (bool): If True, changes the object name adding Recip. of at the beggining of the name. Default: True.
Returns:
(Mueller): Modified object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Simple case
if axis is None or new_obj.ndim <= 1:
M = np.sum(new_obj.M, axis=2)
# Complicated case
else:
# Calculate maximum axis
if isinstance(axis, int):
axis = axis + 2
m = axis
else:
axis = np.array(axis) + 2
m = np.max(axis)
# Check that the axes are correct
if m >= new_obj.ndim + 2:
raise ValueError(
'Axis {} greater than the number of dimensions of {}, which is {}'
.format(m, new_obj.name, new_obj.ndim))
# Reshape M to fit the current shape
shape = [4, 4] + new_obj.shape
M = np.reshape(new_obj.M, shape)
# check if the axis is int or not
if isinstance(axis, int):
M = np.sum(M, axis=axis)
else:
M = np.sum(M, axis=tuple(axis))
# Create the object and return it
new_obj.from_matrix(M)
if change_names:
new_obj.name = 'Sum of ' + new_obj.name
return new_obj
[docs]
def prod(self, axis=None, keep=False, change_name=change_names):
"""Calculates the product of Mueller matrices stored in the object.
Parameters:
axis (int, list or tuple): Axes along which the product is performed. If None, all matrices are multiplied. Default: None
keep (bool): If True, the original element is not updated. Default: False.
change_name (bool): If True, changes the object name adding Recip. of at the beggining of the name. Default: True.
Returns:
(Mueller): Modified object.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.copy()
else:
new_obj = self
# Simple case
if axis is not None:
N_axis = np.array(axis).size
if axis is None or new_obj.ndim <= 1 or new_obj.ndim == N_axis:
M = new_obj.M[:, :, 0]
for ind in range(1, new_obj.size):
M = M @ new_obj.M[:, :, ind]
# Complicated case
else:
# Calculate maximum axis
if isinstance(axis, int):
m = axis + 2
else:
axis = np.array(axis)
m = np.max(axis) + 2
# Check that the axes are correct
if m >= new_obj.ndim + 2:
raise ValueError(
'Axis {} greater than the number of dimensions of {}, which is {}'
.format(m, new_obj.name, new_obj.ndim))
# Calculate shapes, sizes and indices
if isinstance(axis, int):
shape_removed = new_obj.shape[axis]
else:
shape_removed = np.array(new_obj.shape)[axis]
N_removed = np.prod(shape_removed)
ind_removed = combine_indices(
np.unravel_index(np.array(range(N_removed)), shape_removed))
shape_matrix = np.delete(new_obj.shape, axis)
N_matrix = np.prod(shape_matrix)
ind_matrix = combine_indices(
np.unravel_index(np.array(range(N_matrix)), shape_matrix))
shape_final = [4, 4] + list(shape_matrix)
axes_aux = np.array(range(2, new_obj.ndim + 2))
shape_orig = [4, 4] + list(new_obj.shape)
# Make the for loop of the matrix to be calculated
M_orig = np.reshape(new_obj.M, shape_orig)
M = np.zeros(shape_final)
for indM in range(N_matrix):
# Make the multiplication loop
indices = merge_indices(ind_matrix[indM], ind_removed[0], axis)
aux = multitake(M_orig, indices, axes_aux)
for indR in range(1, N_removed):
indices = merge_indices(ind_matrix[indM],
ind_removed[indR], axis)
aux = aux @ multitake(M_orig, indices, axes_aux)
# Store the result
for i1 in range(4):
for i2 in range(4):
ind_aux = tuple([i1, i2] + list(ind_matrix[indM]))
M[ind_aux] = aux[i1, i2]
# Create the object and return it
new_obj.from_matrix(M)
if change_names:
new_obj.name = 'Prod of ' + new_obj.name
return new_obj
####################################################################
# Creation
####################################################################
[docs]
def from_components(self,
components,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates the Mueller matrix object form the arrays of its 16 components components.
Parameters:
components (tuple or list): A 4 element tuple containing the 6 components of the Mueller matrices (M00, M01, ..., M32, M33).
global_phase (float or numpy.ndarray): Adds a global phase to the object. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller_matrix): Created object.
"""
# Prepare variables
components = list(components)
components.append(global_phase)
expand = 16 * [True] + [False]
(components), new_shape = prepare_variables(vars=components,
expand=expand,
length=length,
give_shape=True)
# Create the matrix
self.M = np.array(
[[components[0], components[1], components[2], components[3]],
[components[4], components[5], components[6], components[7]],
[components[8], components[9], components[10], components[11]],
[components[12], components[13], components[14], components[15]]])
# Rest of operations
self.no_rotation = False
# self.size = components[0].size
self.shape, _ = select_shape(self,
shape_var=new_shape,
shape_fun=shape,
shape_like=shape_like)
self.set_global_phase(phase=global_phase)
self.no_rotation = False
return self
[docs]
def from_matrix(self,
M,
global_phase=default_phase,
length=1,
shape_like=None,
shape=None):
"""Create a Mueller object from an external array.
Parameters:
M (numpy.ndarray): New matrix. At least two dimensions must be of size 4.
global_phase (numpy.ndarray): Adds a global phase to the Stokes object. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Check if the matrix is of the correct Size
M = np.array(M)
s = M.size
# 1D and 2D
if M.ndim == 1 or M.ndim == 2:
if M.size % 16 == 0:
M = np.reshape(M, (4, 4, int(M.size / 16)))
else:
raise ValueError(
'M must have a number of elements multiple of 16.')
if M.size == 16:
sh = None
else:
sh = [int(M.size / 16)]
# 3D or more
elif M.ndim > 2:
sh = np.array(M.shape)
N = np.sum(sh == 4)
if N > 1:
# Find the matrix indices and the final shape
ind1 = np.argmin(~(sh == 4))
sh = np.delete(sh, ind1)
ind2 = np.argmin(~(sh == 4))
sh = np.delete(sh, ind2)
ind2 = ind2 + 1
# Calculate the components and construct the matrix from them
M = np.array([[
multitake(M, [0, 0], [ind1, ind2]).flatten(),
multitake(M, [0, 1], [ind1, ind2]).flatten(),
multitake(M, [0, 2], [ind1, ind2]).flatten(),
multitake(M, [0, 3], [ind1, ind2]).flatten()
],
[
multitake(M, [1, 0], [ind1, ind2]).flatten(),
multitake(M, [1, 1], [ind1, ind2]).flatten(),
multitake(M, [1, 2], [ind1, ind2]).flatten(),
multitake(M, [1, 3], [ind1, ind2]).flatten()
],
[
multitake(M, [2, 0], [ind1, ind2]).flatten(),
multitake(M, [2, 1], [ind1, ind2]).flatten(),
multitake(M, [2, 2], [ind1, ind2]).flatten(),
multitake(M, [2, 3], [ind1, ind2]).flatten()
],
[
multitake(M, [3, 0], [ind1, ind2]).flatten(),
multitake(M, [3, 1], [ind1, ind2]).flatten(),
multitake(M, [3, 2], [ind1, ind2]).flatten(),
multitake(M, [3, 3], [ind1, ind2]).flatten()
]])
else:
raise ValueError(
'M must have four elements in at least two dimensions. Instead, it has shape = {}'
.format(M.shape))
else:
raise ValueError('M can not be empty')
# Increase length if required
if M.size == 16 and length > 1:
M = np.multiply.outer(np.squeeze(M), np.ones(length))
# End operations
# self.size = int(M.size / 16)
self.M = M
self.no_rotation = False
self.shape, _ = select_shape(self,
shape_var=sh,
shape_fun=shape,
shape_like=shape_like)
self.set_global_phase(phase=global_phase)
return self
[docs]
def from_normalized(self,
m,
M00=None,
global_phase=default_phase,
shape_like=None,
shape=None):
"""Creates a Mueller object directly from the normalized matrix $m = M/M_{00}$, and $M_{00}$.
Parameters:
m (4x4 numpy.matrix): Mueller matrix
M00 (float): [0, 1] Mean transmission coefficient. Default: maximum possible (None).
global_phase (numpy.ndarray): Adds a global phase to the Stokes object. Default: 0.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Easy case, M00 is a number
if type(M00) in number_types:
M = m * M00
self.from_matrix(M,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
# Complicated case, M00 is an array
else:
# Check that the variables can be used
if m.size > 16 and m.size / 16 != M00.size:
raise ValueError(
'M00 of {} elements is incompatible with m of {} elements'.
format(M00.size, m.size / 16))
# Create the object as if the matrix were not normalized
self.from_matrix(M, shape=False)
components = self.parameters.components(shape=False)
# Denormalize
for ind in range(16):
components[ind] = components[ind] * M00
# Create the correct object
self.from_components(components,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
return self
[docs]
def from_blocks(self,
Dv=zero_D,
Pv=zero_D,
m=zero_m,
M00=1,
global_phase=default_phase,
length=1,
shape_like=None,
shape=None):
"""Method that creates a Mueller object from the block components of its matrix.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
Dv (numpy.ndarray): Diattenuation vector 3xN array.
Pv (numpy.ndarray): Polarizance vector 3xN array.
m (numpy.ndarray): Small matrix m 3x3xN array.
M00 (numpy.ndarray): Parameter of average intensity array of size N. Default: 1
global_phase (numpy.ndarray): Adds a global phase to the Stokes object. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Prepare the variables
M00, Dv, Pv, m, new_shape = prepare_variables_blocks(
M00,
Dv,
Pv,
m,
length=length,
give_shape=True,
multiply_by_M00=True)
# Build the object
components = [
M00, Dv[0, ], Dv[1], Dv[2], Pv[0], m[0, 0], m[0, 1], m[0, 2],
Pv[1], m[1, 0], m[1, 1], m[1, 2], Pv[2], m[2, 0], m[2, 1], m[2, 2]
]
self.from_components(components, global_phase=global_phase)
# Update object parameters
self.shape, _ = select_shape(self,
shape_var=new_shape,
shape_fun=shape,
shape_like=shape_like)
# self.M00, self.D, self.P, self.m = (M00, D, P, m)
return self
# @_actualize_
[docs]
def from_Jones(self, J, length=1, shape_like=None, shape=None):
"""Takes a Jones Matrix and converts into Mueller Matrix
.. math:: M(J)=\left[\begin{array}{cccc}
1 & 0 & 0 & 1\\
1 & 0 & 0 & -1\\
0 & 1 & 1 & 0\\
0 & i & -i & 0
\end{array}\right]\left(J\otimes J^{*}\right)\left[\begin{array}{cccc}
1 & 0 & 0 & 1\\
1 & 0 & 0 & -1\\
0 & 1 & 1 & 0\\
0 & i & -i & 0
\end{array}\right]^{-1}
References:
"Polarized light and the Mueller Matrix approach", 2nd Ed., J. J. Gil and E. Ossikowski, pp 107.
Parameters:
J (Jones_matrix): Jones matrix object.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (float or numpy.ndarray): Use the shape of this array. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Extract components from the Jones object and derivatives
J00, J01, J10, J11 = J.parameters.components(shape=False)
J00m = np.abs(J00)**2
J01m = np.abs(J01)**2
J10m = np.abs(J10)**2
J11m = np.abs(J11)**2
J00c = np.conj(J00)
J01c = np.conj(J01)
J10c = np.conj(J10)
J11c = np.conj(J11)
# Create the components of the Mueller object
components = [0] * 16
components[0] = (J00m + J01m + J10m + J11m) / 2
components[1] = (J00m - J01m + J10m - J11m) / 2
components[2] = (J00 * J01c + J10 * J11c).real
components[3] = (J00 * J01c + J10 * J11c).imag
components[4] = (J00m + J01m - J10m - J11m) / 2
components[5] = (J00m - J01m - J10m + J11m) / 2
components[6] = (J00 * J01c - J10 * J11c).real
components[7] = (J00 * J01c + J10c * J11).imag
components[8] = (J00 * J10c + J01 * J11c).real
components[9] = (J00 * J10c - J01 * J11c).real
components[10] = (J00 * J11c + J01 * J10c).real
components[11] = (J10 * J01c + J11c * J00).imag
components[12] = (J10 * J00c + J01c * J11).imag
components[13] = (J10 * J00c + J11c * J01).imag
components[14] = (J11 * J00c + J01c * J10).imag
components[15] = (J00 * J11c - J01 * J10c).real
phase = J.parameters.global_phase()
# Create the object
self.from_components(components,
global_phase=phase,
length=length,
shape=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def from_covariance(self,
H,
global_phase=default_phase,
length=1,
shape_like=None,
shape=None):
"""Calculates the Mueller matrix from the covariance matrix:
$M_{ij}=Trace\left[\left(\sigma_{i}\otimes\sigma_{j}^{*}\right)H\right]$
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
H (Mueller or numpy.ndarray): Covariance matrix.
global_phase (numpy.ndarray): Adds a global phase to the Stokes object. Default: 0.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Check if H is a py_pol object or an array
if isinstance(H, np.ndarray):
obj = Mueller(self.name)
obj.from_matrix(H)
else:
obj = H
old_shape = obj.shape
obj.shape = None
# Initialize the matrix
M = np.zeros((4, 4, obj.size), dtype=complex)
# Loop in elements
components = []
elem = Mueller()
for indI in range(4):
for indJ in range(4):
# Calculate the Sij matrix
Sij = np.multiply.outer(np.kron(S[indI], np.conj(S[indJ])),
np.ones(obj.size))
elem.from_matrix(Sij)
# Multiply it by H
elem = elem * obj
# Save the trace as the new component
components += [
np.array(elem.parameters.trace(shape=False), dtype=float)
]
# Create the new object
self.from_components(components, global_phase=global_phase)
# Reshape if necessary
self.shape, _ = select_shape(obj=self,
shape_var=old_shape,
shape_fun=shape,
shape_like=shape_like)
return self
# @_actualize_
# def from_inverse(self, M):
# """Calculates the Mueller matrix from the inverse matrix.
#
# Parameters:
# M (numpy.matrix 4x4 or Mueller object): Inverse matrix.
#
# Returns:
# (numpy.matrix): 4x4 matrix.
# """
# try:
# if M.type == 'Mueller':
# self.from_matrix(M.M.I)
# else:
# self.from_matrix(M.I)
# except:
# self.from_matrix(M.I)
# return self.M
[docs]
def vacuum(self, global_phase=0, length=1, shape_like=None, shape=None):
"""Creates the matrix for vacuum i.e., an optically neutral element.
Parameters:
global_phase (float or numpy.ndarray): Adds a global phase to the Jones matrix. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
self.filter_amplifier(D=1,
length=length,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
return self
[docs]
def filter_amplifier(self,
D=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates the Mueller object of neutral filters or amplifiers.
Parameters:
D (float or numpy.ndarray): Attenuation (gain if > 1). Default: 1.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Calculate
components = [D] + [0] * 4 + [D] + \
[0] * 4 + [D] + [0] * 4 + [D]
self.from_components(components=components,
global_phase=global_phase,
length=length,
shape=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def mirror(self,
ref=1,
ref_field=None,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Mueller matrix of a mirror.
Parameters:
ref (float or numpy.ndarray): Intensity reflectivity of the mirror. Default: 1.
ref_field (float or numpy.ndarray): Electric field reflectivity coefficient. If not None, it overrides REF. Default: None.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use the intensity reflectivity
if ref_field is not None:
ref = np.abs(ref_field)**2
# Prepare variables
(ref,
global_phase), new_shape = prepare_variables(vars=[ref, global_phase],
expand=[True, False],
length=length,
give_shape=True)
# Calculate
components = [ref] + [0] * 4 + [ref] + \
[0] * 4 + [-ref] + [0] * 4 + [-ref]
self.from_components(components=components,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.no_rotation = True
return self
# @_actualize_
[docs]
def depolarizer_perfect(self,
M00=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates a perfect depolarizer:
.. math:: M_{p}=\left[\begin{array}{cccc}
M_{00} & 0 & 0 & 0\\
0 & d_{1} & 0 & 0\\
0 & 0 & d_{2} & 0\\
0 & 0 & 0 & d_{3}
\end{array}\right]
Parameters:
M00 (float, default 1): Parameter of average intensity. Default: 1.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Calculate
components = [M00] + [0] * 15
self.from_components(components=components,
length=length,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def depolarizer_diagonal(self,
d,
M00=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates a diagonal depolarizer:
.. math:: M_{p}=\left[\begin{array}{cccc}
M_{00} & 0 & 0 & 0\\
0 & d_{1} & 0 & 0\\
0 & 0 & d_{2} & 0\\
0 & 0 & 0 & d_{3}
\end{array}\right]
Parameters:
d (list, float or numpy.ndarray): Absorption coefficients. If list, it must contain three float or numpy arrays, one for each diagonal value. If float or numpy.ndarray, $d_1 = d_2 = d_3$.
M00 (float, default 1): Parameter of average intensity. Default: 1.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Calculate
if isinstance(d, list):
components = [M00] + [0] * 4 + [d[0]] + [0] * 4 + [
d[1]
] + [0] * 4 + [d[2]]
else:
components = [M00] + [0] * 4 + [d] + [0] * 4 + [d] + [0] * 4 + [d]
self.from_components(components=components,
length=length,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def depolarizer_states(self,
d,
S,
Dv=zero_D,
Pv=zero_D,
M00=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates a general depolarizer form its three eigenstates and their depolarization factors (eigenvalues), plus its polarization or polarizance vector.
Parameters:
d (list, float or numpy.ndarray): Depolarization factors (eigenvalues of m). If list, it must contain three float or numpy arrays, one for each diagonal value. If float or numpy.ndarray, $d_1 = d_2 = d_3$.
S (list or Stokes): Principal states. If list, it must contain three Stokes objects. If Stokes, at least one dimension must have dimension 3.
Dv (numpy.ndarray): Diattenuation vector. If None, the polarizance vector is used instead. Default: None.
Pv (numpy.ndarray): Polarizance vector. Used only if Dv is None. If None, the depolarizer will have zero diattenuation and polarizance. Default: None.
M00 (float, default 1): Parameter of average intensity. Default: 1.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Prepare variables. First the extinction
if isinstance(d, list):
d_list = d
else:
sh = np.array(d.shape)
ind = np.argmin(~(sh == 3))
d_list = [
np.take(d, 0, axis=ind),
np.take(d, 1, axis=ind),
np.take(d, 2, axis=ind)
]
(d1, d2, d3), new_shape1 = prepare_variables(vars=d,
expand=[True] * 3,
length=length,
give_shape=True)
d1 = np.sqrt(d1.flatten())
d2 = np.sqrt(d2.flatten())
d3 = np.sqrt(d3.flatten())
# Now, the states
if isinstance(S, list):
S1, S2, S3 = S
else:
S1, S2, S3 = create_Stokes(N=3)
sh = np.array(S.shape)
ind = np.argmin(~(sh == 3))
M = S.parameters.matrix()
S1.from_matrix(np.take(M, 0, axis=ind + 1))
S2.from_matrix(np.take(M, 1, axis=ind + 1))
S3.from_matrix(np.take(M, 2, axis=ind + 1))
S1.normalize()
S2.normalize()
S3.normalize()
comp = [
S1.M[1, :], S1.M[2, :], S1.M[3, :], S2.M[1, :], S2.M[2, :],
S2.M[3, :], S3.M[1, :], S3.M[2, :], S3.M[3, :]
]
comp, new_shape2 = prepare_variables(vars=comp,
expand=[True] * 9,
length=length,
give_shape=True)
for ind in range(9):
comp[ind] = comp[ind].flatten()
new_shape = take_shape([new_shape1, new_shape2])
# Create the small matrix of the depolarizer
v1 = np.array([d1 * comp[0], d1 * comp[1], d1 * comp[2]])
v2 = np.array([d2 * comp[3], d2 * comp[4], d2 * comp[5]])
v3 = np.array([d3 * comp[6], d3 * comp[7], d3 * comp[8]])
m1 = kron_axis(v1, v1, axis=0)
m2 = kron_axis(v2, v2, axis=0)
m3 = kron_axis(v3, v3, axis=0)
m = m1 + m2 + m3
# Create the object
self.from_blocks(M00=M00, Dv=Dv, Pv=Pv, m=m)
self.shape, _ = select_shape(obj=self,
shape_var=new_shape,
shape_fun=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def diattenuator_perfect(self,
azimuth=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Mueller 4x4 matrix for a perfect diattenuator (polarizer).
Parameters:
azimuth (float or numpy.ndarray): Rotation angle of the high transmission polarizer axis. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
self.diattenuator_linear(p1=1,
p2=0,
azimuth=azimuth,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def diattenuator_linear(self,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
azimuth=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Mueller matrices of pure linear homogeneous diattenuators.
.. math:: M\left(\theta=0\right)=\frac{1}{2}\left[\begin{array}{cccc}
p_{1}^{2}+p_{2}^{2} & p_{1}^{2}-p_{2}^{2} & 0 & 0\\
p_{1}^{2}-p_{2}^{2} & p_{1}^{2}+p_{2}^{2} & 0 & 0\\
0 & 0 & 2p_{1}p_{2} & 0\\
0 & 0 & 0 & 2p_{1}p_{2}
\end{array}\right]
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), (4.79) - p. 143.
Handbook of Optics vol 2. 22.16 (Table 1).
Parameters:
p1 (float or numpy.ndarray): Field transmission of the transmission axis. Default: 1.
p2 (float or numpy.ndarray): Field transmission of the attenuation axis. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
azimuth (float or numpy.ndarray): Rotation angle of the high transmission polarizer axis. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use field transmission coefficients
if Tmax is None:
Tmax = p1**2
if Tmin is None:
Tmin = p2**2
# Prepare variables
(p1, p2,
azimuth), new_shape = prepare_variables(vars=[p1, p2, azimuth],
expand=[False, False, False],
length=length,
give_shape=True)
# Calculate intensity transmission coefficients
a = (Tmax + Tmin) / 2
b = (Tmax - Tmin) / 2
c = np.sqrt(Tmax * Tmin)
# Calculate the matrix
components = [a] + [b] + [0] * 2 + [b] + \
[a] + [0] * 4 + [c] + [0] * 4 + [c]
self.from_components(components=components,
global_phase=global_phase,
length=length)
self.rotate(azimuth)
self.shape, _ = select_shape(obj=self,
shape_var=new_shape,
shape_fun=shape,
shape_like=shape_like)
return self
# @_actualize_
[docs]
def diattenuator_charac_angles(self,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
alpha=0,
delay=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Method that calculates the most general homogenous diattenuator from diattenuator parameters with the intermediate step of calculating the diattenuation vector.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 142.
Parameters:
p1 (float or numpy.ndarray): Electric field transmission coefficient of the transmission eigenstate Default: 1
p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
alpha (float or numpy.ndarray): [0, pi/2]: tan(alpha) is the ratio between field amplitudes of X and Y components. Default: 0.
delay (float or numpy.ndarray): [0, 2*pi]: phase difference between X and Y field components. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use field transmission coefficients
if Tmax is None:
Tmax = p1**2
if Tmin is None:
Tmin = p2**2
# Prepare variables
(Tmax, Tmin, alpha, delay), new_shape = prepare_variables(
vars=[Tmax, Tmin, alpha, delay],
expand=[True, False, False, False],
length=length,
give_shape=True)
# Restrict parameter values to the correct interval
alpha = put_in_limits(alpha, "alpha")
delay = put_in_limits(delay, "delay")
# Calculate diattenuation vector and M00
cte = (Tmax - Tmin) / (Tmax + Tmin)
cond = np.isnan(cte)
if np.any(cond):
cte[cond] = 0
M00 = (Tmax + Tmin) / 2
Dv = [
cte * np.cos(2 * alpha), cte * np.sin(2 * alpha) * np.cos(delay),
cte * np.sin(2 * alpha) * np.sin(delay)
]
# Create the element
self.diattenuator_vector(Dv=Dv,
M00=M00,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
# @_actualize_
[docs]
def diattenuator_azimuth_ellipticity(self,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
azimuth=0,
ellipticity=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Method that calculates the most general homogenous diattenuator from
diattenuator parameters with the intermediate step of calculating the
diattenuation vector.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 142.
Parameters:
p1 (float or numpy.ndarray): Electric field transmission coefficient of the transmission eigenstate. Default: 1.
p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
azimuth (float): [0, pi]: Azimuth. Defaut: 0
ellipticity (float): [-pi/4, pi/4]: Ellipticity angle. Default: 0
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use field transmission coefficients
if Tmax is None:
Tmax = p1**2
if Tmin is None:
Tmin = p2**2
# Prepare variables
(Tmax, Tmin, azimuth, ellipticity), new_shape = prepare_variables(
vars=[Tmax, Tmin, azimuth, ellipticity],
expand=[True, False, False, False],
length=length,
give_shape=True)
# Restrict parameter values to the correct interval
azimuth = put_in_limits(azimuth, "azimuth")
ellipticity = put_in_limits(ellipticity, "ellipticity")
# Calculate the diattenuation vector and M00
cte = (Tmax - Tmin) / (Tmax + Tmin)
cond = np.isnan(cte)
if np.any(cond):
cte[cond] = 0
M00 = (Tmax + Tmin) / 2
Dv = [
cte * np.cos(2 * azimuth) * np.cos(2 * ellipticity),
cte * np.sin(2 * azimuth) * np.cos(2 * ellipticity),
cte * np.sin(2 * ellipticity)
]
# Create the element
self.diattenuator_vector(Dv=Dv,
M00=M00,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
# @_actualize_
[docs]
def diattenuator_vector(self,
Dv,
M00=None,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Method that calculates the most general homogenous diattenuator from the
Diattenuation or Polarizance vector.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 142.
Parameters:
Dv (3xN numpy.ndarray): Diattenuation or Polarizance vectors.
M00 (float or numpy.ndarray): Parameter of average intensity. If None, the maximum possible value is used. Default: None.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Prepare the variables
_, Dv, _, _, new_shape = prepare_variables_blocks(Dv=Dv,
give_shape=True)
# Calculate the diattenuation and make it physically realizable
d = np.linalg.norm(Dv, axis=0)
cond = d > 1
if np.any(cond):
Dv[0, cond] = Dv[0, cond] / d[cond]
Dv[1, cond] = Dv[1, cond] / d[cond]
Dv[2, cond] = Dv[2, cond] / d[cond]
d[cond] = 1
# Calculate maximum achievable m00
if M00 is None:
M00 = 1 / (1 + d)
# Calculate the m small matrix
skd = np.sqrt(1 - d**2)
m1 = np.multiply.outer(np.diag([1, 1, 1]), skd)
cte = (1 - skd) / d**2
cond = np.isnan(cte)
if np.any(cond):
cte[cond] = 0
m2 = kron_axis(Dv, Dv, axis=0) * np.multiply.outer(
np.ones((3, 3)), cte)
m = m1 + m2
# This equation fails if d=0
cond = d < tol_default
if np.any(cond):
m = np.array([np.eye(3)])
# Now we have all the necessary blocks
self.from_blocks(Dv,
Dv,
m,
M00,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
# @_actualize_
[docs]
def retarder_linear(self,
R,
azimuth=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Mueller matrix of homogeneous linear retarders.
.. math:: M\left(\theta=0\right)=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & \cos(\Delta) & \sin(\Delta)\\
0 & 0 & -\sin(\Delta) & \cos(\Delta)
\end{array}\right]
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), (4.31) - p. 132
Handbook of Optics vol 2. 22.16 (Table 1).
Parameters:
R (float or numpy.ndarray): [0, pi] Retardance introduced to the slow eigenstate respect to the fast eigenstate.
azimuth (float or numpy.ndarray): rotation angle of the high transmission polarizer axis. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Make sure that retardance is correct
R = put_in_limits(R, 'Retardance')
# Calculate the matrix components
components = [1] + [0] * 4 + [1] + [0] * 4 + \
[np.cos(R)] + [np.sin(R)] + [0] * 2 + [-np.sin(R)] + [np.cos(R)]
# Create the object
self.from_components(components,
global_phase=global_phase,
length=length,
shape=shape,
shape_like=shape_like)
self.rotate(azimuth)
return self
# @_actualize_
[docs]
def quarter_waveplate(self,
azimuth=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Mueller matrices of ideal quarter-wave retarder.
.. math:: M\left(\theta=0\right)=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & -1 & 0
\end{array}\right]
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), (4.32) - p. 132
Handbook of Optics vol 2. 22.16 (Table 1).
Parameters:
azimuth (float or numpy.ndarray): rotation angle of the high transmission polarizer axis. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Calculate the matrix components
components = [1] + [0] * 4 + [1] + [0] * 5 + \
[1] + [0] * 2 + [-1] + [0]
# Create the object
self.from_components(components,
global_phase=global_phase,
length=length,
shape=shape,
shape_like=shape_like)
self.rotate(azimuth)
return self
# @_actualize_
[docs]
def half_waveplate(self,
azimuth=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Mueller matrices of ideal half-wave retarders.
.. math:: M\left(\theta=0\right)=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & -1 & 0\\
0 & 0 & 0 & -1
\end{array}\right]
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), (4.32) - p. 132
Handbook of Optics vol 2. 22.16 (Table 1).
Parameters:
azimuth (float or numpy.ndarray): rotation angle of the high transmission polarizer axis. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Calculate the matrix components
components = [1] + [0] * 4 + [1] + [0] * 4 + [-1] + [0] * 4 + [-1]
# Create the object
self.from_components(components,
global_phase=global_phase,
length=length,
shape=shape,
shape_like=shape_like)
self.rotate(azimuth)
return self
# @_actualize_
[docs]
def retarder_charac_angles(self,
R,
alpha,
delay,
M00=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Method that calculates the most general homogeneous retarder from the characteristic angles of the fast eigenstate. The method calculates first the retardance vector, and uses it to calculate the Mueler matrix.
References:
"Polarized light and the Mueller Matrix approach", J. J. Gil, pp 125.
Parameters:
R (float or numpy.ndarray): [0, pi] Retardance introduced to the slow eigenstate respect to the fast eigenstate.
alpha (float or numpy.ndarray): [0, pi]: tan(alpha) is the ratio between amplitudes of the electric field of the fast eigenstate.
delay (float or numpy.ndarray): [0, 2*pi]: phase difference between both components of the electric field of the fast eigenstate.
M00 (float or numpy.ndarray): Parameter of average intensity. Default: 1
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Prepare variables
(R, alpha,
delay), new_shape = prepare_variables(vars=[R, alpha, delay],
expand=[False, True, True],
give_shape=True)
# Restrict parameter values to the correct interval
alpha = put_in_limits(alpha, "alpha")
delay = put_in_limits(delay, "delay")
# Calculate the normalized retardance vector
Rv = np.array([
np.cos(2 * alpha),
np.sin(2 * alpha) * np.cos(delay),
np.sin(2 * alpha) * np.sin(delay)
])
# Create the object
self.retarder_vector(Rv=Rv,
R=R,
kind='normalized',
M00=M00,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
# @_actualize_
[docs]
def retarder_azimuth_ellipticity(self,
R,
azimuth,
ellipticity,
M00=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Method that calculates the most general homogeneous retarder from azimuth and ellipticity of the fast eigenstate. The method calculates first the retardance vector, and uses it to calculate the Mueler matrix.
References:
"Polarized light and the Mueller Matrix approach", J. J. Gil, pp 125.
Parameters:
R (float or numpy.ndarray): [0, pi] Retardance introduced to the slow eigenstate respect to the fast eigenstate.
azimuth (float or numpy.ndarray): [0, pi]: Azimuth.
ellipticity (float or numpy.ndarray): [-pi/4, pi/4]: Ellipticity angle.
M00 (float or numpy.ndarray): Parameter of average intensity. Default: 1.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Prepare variables
(R, azimuth, ellipticity), new_shape = prepare_variables(
vars=[R, azimuth, ellipticity],
expand=[False, True, True],
give_shape=True)
# Restrict parameter values to the correct interval
azimuth[np.isnan(azimuth)] = 0 # Avoid nan values
azimuth = put_in_limits(azimuth, "azimuth")
ellipticity = put_in_limits(ellipticity, "ellipticity")
# Calculate the normalized retardance vector
Rv = np.array([
np.cos(2 * azimuth) * np.cos(2 * ellipticity),
np.sin(2 * azimuth) * np.cos(2 * ellipticity),
np.sin(2 * ellipticity)
])
# Create the object
self.retarder_vector(Rv=Rv,
R=R,
kind='normalized',
M00=M00,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
# @_actualize_
[docs]
def retarder_vector(self,
Rv,
R=90 * degrees,
kind='normalized',
M00=1,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Method that calculates the most general homogeneous retarder from the retardance vector.
References:
"Polarized light and the Mueller Matrix approach", J. J. Gil, pp 125.
Parameters:
Rv (3xN numpy.ndarray): Retardance vector.
R (float or numpy.ndarray): [0, pi] Retardance introduced to the slow eigenstate respect to the fast eigenstate. Default: 90 degrees.
kind (string): Identifies the type of retardance vector. There are three possibilities: NORMALIZED (also called Pauli vector), STRAIGHT or COMPLETE. Default: 'normalized'.
M00 (float or numpy.ndarray): Mean transmission coefficient. If different than 1, the object won't be a pure retarder. Default: 1.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Prepare the variables
_, Rv, _, _, new_shape = prepare_variables_blocks(Dv=Rv,
give_shape=True)
# Choose the right retardance
if kind in ('complete', 'COMPLETE', 'Complete'):
# Complete retardance vectors
R = np.linalg.norm(Rv, axis=0)
aux = np.multiply.outer(np.ones(3), R)
cond = R != 0
if np.any(cond):
Rv[:, cond] = Rv[:, cond] / aux[:, cond]
elif kind in ('straight', 'STRAIGHT', 'Straight'):
# Straight retardance vectors
R = np.linalg.norm(Rv, axis=0) * np.pi
aux = np.multiply.outer(np.ones(3), R / np.pi)
cond = R != 0
if np.any(cond):
Rv[:, cond] = Rv[:, cond] / aux[:, cond]
else:
# Normalized (Pauli) retardance vector
R = put_in_limits(R, 'Retardance')
norm = np.linalg.norm(Rv, axis=0)
aux = np.multiply.outer(np.ones(3), norm)
cond = (norm != 0) + (norm != 1)
if np.any(cond):
Rv[:, cond] = Rv[:, cond] / aux[:, cond]
# Reshape R
(R, M00), new_shape2 = prepare_variables(vars=[R, M00],
expand=[True, False],
length=Rv.shape[1],
give_shape=True)
new_shape = take_shape((new_shape, new_shape2))
# Calculate small m matrix
sR, cR = (np.sin(R), np.cos(R))
m = np.multiply.outer(np.eye(3), cR)
for indI in range(3):
for indJ in range(3):
m[indI, indJ, :] += Rv[indI, :] * Rv[indJ, :] * (1 - cR)
for indK in range(3):
m[indI, indJ, :] = m[indI, indJ, :] + \
Eijk(indI, indJ, indK) * Rv[indK, :] * sR
# Create the object
self.from_blocks(Dv=np.zeros(3),
Pv=np.zeros(3),
m=m,
M00=M00,
global_phase=global_phase,
length=length,
shape=shape,
shape_like=shape_like)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
# @_actualize_
[docs]
def diattenuator_retarder_linear(self,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
R=90 * degrees,
azimuth=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates the Mueller matrices of linear diattenuator retarders with the same eigenstates for diattenuation and retardance.
.. math:: M\left(\theta=0\right)=\frac{1}{2}\left[\begin{array}{cccc}
p_{1}^{2}+p_{2}^{2} & p_{1}^{2}-p_{2}^{2} & 0 & 0\\
p_{1}^{2}-p_{2}^{2} & p_{1}^{2}+p_{2}^{2} & 0 & 0\\
0 & 0 & 2p_{1}p_{2}\cos(\varDelta) & 2p_{1}p_{2}\sin(\varDelta)\\
0 & 0 & -2p_{1}p_{2}\sin(\varDelta) & 2p_{1}p_{2}\cos(\varDelta)
\end{array}\right]
References:
Handbook of Optics vol 2. 22.16 (Table 1).
Parameters:
p1 (float or numpy.ndarray): Field transmission of the transmission axis. Default: 1.
p2 (float or numpy.ndarray): Field transmission of the attenuation axis. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
R (float or numpy.ndarray): [0, pi] Retardance introduced to the slow eigenstate respect to the fast eigenstate. Default: 90 degrees.
azimuth (float or numpy.ndarray): rotation angle of the high transmission polarizer axis. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use field transmission coefficients
if Tmax is None:
Tmax = p1**2
if Tmin is None:
Tmin = p2**2
# Prepare variables
(p1, p2,
azimuth), new_shape = prepare_variables(vars=[p1, p2, azimuth],
expand=[False, False, False],
length=length,
give_shape=True)
# Calculate intensity transmission coefficients
a = (Tmax + Tmin) / 2
b = (Tmax - Tmin) / 2
c = np.sqrt(Tmax * Tmin)
# Calculate the matrix
components = [a] + [b] + [0] * 2 + [b] + [a] + [0] * 4 + \
[c * np.cos(R)] + [c * np.sin(R)] + [0] * 2 + \
[-c * np.sin(R)] + [c * np.cos(R)]
self.from_components(components=components,
global_phase=global_phase,
shape=shape,
shape_like=shape_like)
self.rotate(azimuth)
self.shape, _ = select_shape(obj=self, shape_var=new_shape)
return self
[docs]
def diattenuator_retarder_azimuth_ellipticity(self,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
R=0,
azimuth=0,
ellipticity=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates the most general homogenous diattenuator retarder from the azimuth and ellipticity of the fast eigenstate.
Parameters:
p1 (float or numpy.ndarray): Field transmission of the fast axis. Default: 1.
p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
R (float or numpy.ndarray): Retardance. Default: 0.
azimuth (float or numpy.ndarray): rotation angle of the high transmission polarizer axis. Default: 0.
ellipticity (float): [-pi/4, pi/]: Ellipticity angle. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use field transmission coefficients
if Tmax is not None:
p1 = np.sqrt(Tmax)
if Tmin is not None:
p2 = np.sqrt(Tmin)
# Prepare variables
(p1, p2, R,
azimuth, ellipticity), new_shape = prepare_variables(vars=[p1, p2, R, azimuth, ellipticity],
expand=[False, False, False, False, False],
length=length,
give_shape=True)
# Restrict retardance between 0 and 360 degrees
R = put_in_limits(R, "delay", out_number=False)
# Create the two objects
name = self.name
E1 = Mueller().diattenuator_azimuth_ellipticity(p1=p1,
p2=p2,
azimuth=azimuth,
ellipticity=ellipticity,
shape=shape,
shape_like=shape_like,
length=length)
# Orthogonal state for R > 180º
cond = R > np.pi
if cond.size > 0:
azimuth[cond] = azimuth[cond] + np.pi/2
ellipticity[cond] = -ellipticity[cond]
E2 = Mueller().retarder_azimuth_ellipticity(R=R,
azimuth=azimuth,
ellipticity=ellipticity,
shape=shape,
shape_like=shape_like,
length=length)
self = E1 * E2
self.name = name
return self
[docs]
def diattenuator_retarder_charac_angles(self,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
R=0,
alpha=0,
delay=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates the most general homogenous diattenuator retarder from the characteristic angles of the fast eigenstate.
Parameters:
p1 (float or numpy.ndarray): Field transmission of the fast axis. Default: 1.
p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
R (float or numpy.ndarray): Retardance. Default: 0.
alpha (float or numpy.ndarray): [0, pi/2]: tan(alpha) is the ratio between field amplitudes of X and Y components. Default: 0.
delay (float or numpy.ndarray): [0, 2*pi]: phase difference between X and Y field components. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Mueller): Created object.
"""
# Use field transmission coefficients
if Tmax is not None:
p1 = np.sqrt(Tmax)
if Tmin is not None:
p2 = np.sqrt(Tmin)
# Prepare variables
(p1, p2, R,
alpha, delay), new_shape = prepare_variables(vars=[p1, p2, R, alpha, delay],
expand=[False, False, False, False, False],
length=length,
give_shape=True)
# Restrict retardance between 0 and 360 degrees
R = put_in_limits(R, "delay", out_number=False)
# Create the two objects
name = self.name
E1 = Mueller().diattenuator_charac_angles(p1=p1,
p2=p2,
alpha=alpha,
delay=delay,
shape=shape,
shape_like=shape_like,
length=length)
# Orthogonal state for R > 180º
cond = R > np.pi
if cond.size > 0:
alpha[cond] = np.pi/2 - alpha[cond]
delay[cond] = delay[cond] + np.pi
E2 = Mueller().retarder_charac_angles(R=R,
alpha=alpha,
delay=delay,
shape=shape,
shape_like=shape_like,
length=length)
self = E1 * E2
self.name = name
return self
#
# def diattenuator_retarder_azimuth_ellipticity(self,
# p1=1,
# p2=1,
# Tmax=None,
# Tmin=None,
# R=0,
# azimuth=0,
# ellipticity=0,
# global_phase=0,
# length=1,
# shape_like=None,
# shape=None):
# """Creates the most general homogenous diattenuator retarder from the azimuth and ellipticity of the fast eigenstate.
#
# Parameters:
# p1 (float or numpy.ndarray): Field transmission of the fast axis. Default: 1.
# p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
# Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
# Tmax (float or numpy.ndarray): Minimum transmission. If not None, overrides p1. Default: None.
# R (float or numpy.ndarray): Retardance. Default: 0.
# azimuth (float or numpy.ndarray): rotation angle of the high transmission polarizer axis. Default: 0.
# ellipticity (float): [-pi/4, pi/]: Ellipticity angle. Default: 0.
# global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
# length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
# shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
# shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
#
# Returns:
# (Mueller): Created object.
# """
# # Use field transmission coefficients
# if Tmax is not None:
# p1 = np.sqrt(Tmax)
# if Tmin is not None:
# p2 = np.sqrt(Tmin)
# # Create the two objects
# E1 = Mueller()
# E1.diattenuator_azimuth_ellipticity(p1=p1,
# p2=p2,
# azimuth=azimuth,
# ellipticity=ellipticity,
# shape=shape,
# shape_like=shape_like,
# length=length)
# E2 = Mueller()
# E2.retarder_azimuth_ellipticity(R=R,
# azimuth=azimuth,
# ellipticity=ellipticity,
# shape=shape,
# shape_like=shape_like,
# length=length)
# # Multiply and extract
# new_obj = E1 * E2
# self.from_matrix(new_obj.M)
# self.shape, _ = new_obj.shape, new_obj.ndim
# # return self
#
# def diattenuator_retarder_charac_angles(self,
# p1=1,
# p2=1,
# Tmax=None,
# Tmin=None,
# R=0,
# alpha=0,
# delay=0,
# global_phase=0,
# length=1,
# shape_like=None,
# shape=None):
# """Creates the most general homogenous diattenuator retarder from the characteristic angles of the fast eigenstate.
#
# Parameters:
# p1 (float or numpy.ndarray): Field transmission of the fast axis. Default: 1.
# p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
# Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
# Tmax (float or numpy.ndarray): Minimum transmission. If not None, overrides p1. Default: None.
# R (float or numpy.ndarray): Retardance. Default: 0.
# alpha (float or numpy.ndarray): [0, pi/2]: tan(alpha) is the ratio between field amplitudes of X and Y components. Default: 0.
# delay (float or numpy.ndarray): [0, 2*pi]: phase difference between X and Y field components. Default: 0.
# length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
# shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
# shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
#
# Returns:
# (Mueller): Created object.
# """
# # Use field transmission coefficients
# if Tmax is not None:
# p1 = np.sqrt(Tmax)
# if Tmin is not None:
# p2 = np.sqrt(Tmin)
# # Create the two objects
# E1, E2 = create_Mueller(N=2)
# E1.diattenuator_charac_angles(p1=p1,
# p2=p2,
# alpha=alpha,
# delay=delay,
# shape=shape,
# shape_like=shape_like,
# length=length)
# E2.retarder_charac_angles(R=R,
# alpha=alpha,
# delay=delay,
# shape=shape,
# shape_like=shape_like,
# length=length)
# # Multiply and extract
# new_obj = E1 * E2
# self.from_matrix(new_obj.M)
# self.shape, _ = new_obj.shape, new_obj.ndim
# return self
[docs]
def general_eigenstates(self,
E1,
E2=None,
p1=1,
p2=0,
Tmax=None,
Tmin=None,
R=0,
global_phase=0,
length=1,
shape_like=None,
shape=None):
"""Creates the most general pure optical element from its eigenstates.
Parameters:
E1 (Jones_vector): First eigenstate.
E2 (Jones_vector): Second eigenstate. If None, E2 is taken as the perpendicular state to E1, so the optical object is homogenous. Default: None
p1 (float or numpy.ndarray): Field transmission of the fast axis. Default: 1.
p2 (float or numpy.ndarray): Electric field transmission coefficient of the extinction eigenstate. Default: 0.
Tmax (float or numpy.ndarray): Maximum transmission. If not None, overrides p1. Default: None.
Tmin (float or numpy.ndarray): Minimum transmission. If not None, overrides p2. Default: None.
R (float or numpy.ndarray): Retardance. Default: 0.
global_phase (float or numpy.ndarray): Global phase introduced by the optical element. Default: 0.
length (int): If final object is of size 1, it is stretched to match this size. Default: 1.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
Returns:
(Jones_matrix): Created object.
"""
# Use field transmission coefficients
if Tmax is not None:
p1 = np.sqrt(Tmax)
if Tmin is not None:
p2 = np.sqrt(Tmin)
# Main calculation
if E2 is None:
# Simple case: homogenous case
az, el = E1.parameters.azimuth_ellipticity()
self.diattenuator_retarder_azimuth_ellipticity(
p1=p1,
p2=p2,
R=R,
azimuth=az,
ellipticity=el,
shape=shape,
shape_like=shape_like,
length=length,
global_phase=global_phase)
else:
# Complicated case: inhomogeneous case. It must be done from Jones
E = Jones_vector(self.name)
E.general_eigenstates(E1=E1,
E2=E2,
p1=p1,
p2=p2,
R=R,
global_phase=global_phase,
length=length,
shape_like=shape_like,
shape=shape)
self.from_Jones(E)
return self
########################################################################
# Parameters
########################################################################
[docs]
class Parameters_Mueller(object):
"""Class for Mueller Matrix Parameters.
Parameters:
self.parent (Mueller_matrix): Parent object.
"""
def __init__(self, parent):
self.parent = parent
def __repr__(self):
"""print all parameters
TODO: print all as mueller_matrix"""
dict = self.get_all(verbose=True, draw=True)
return ''
[docs]
def get_all(self, verbose=False, draw=False):
"""Creates a dictionary with all the parameters of Mueller matrix.
Parameters:
verbose (bool): If True, print all parameters. Default: False.
draw (bool): If True, draw all plots/images of the parameters. Default: False.
"""
dict_params = {}
dict_params['Mean transmission'] = self.mean_transmission(
verbose=verbose, draw=draw)
dict_params['Inhomogeneity'] = self.inhomogeneity(verbose=verbose,
draw=draw)
dict_params['Diattenuation'] = self.diattenuation(verbose=verbose,
draw=draw)
dict_params['Linear diattenuation'] = self.diattenuation_linear(
verbose=verbose, draw=draw)
dict_params[
'Circular diattenuation'] = self.diattenuation_circular(
verbose=verbose, draw=draw)
dict_params['Polarizance'] = self.polarizance(verbose=verbose,
draw=draw)
dict_params['Linear polarizance'] = self.polarizance_linear(
verbose=verbose, draw=draw)
dict_params['Circular polarizance'] = self.polarizance_circular(
verbose=verbose, draw=draw)
dict_params['Spheric purity'] = self.spheric_purity(
verbose=verbose, draw=draw)
dict_params['Retardance'] = self.retardance(verbose=verbose,
draw=draw)
dict_params['Polarimetric purity'] = self.polarimetric_purity(
verbose=verbose, draw=draw)
dict_params['Depolarization degree'] = self.depolarization_index(
verbose=verbose, draw=draw)
dict_params[
'Polarimetric purity indices'] = self.polarimetric_purity_indices(
verbose=verbose, draw=draw)
dict_params['Transmissions'] = self.transmissions(kind='all',
verbose=verbose,
draw=draw)
dict_params['Retardance'] = self.retardance(verbose=verbose,
draw=draw)
dict_params['Eigenstates'], dict_params[
'Eigenvectors'] = self.eig(verbose=verbose, draw=draw)
dict_params['Determinant'] = self.det(verbose=verbose, draw=draw)
dict_params['Trace'] = self.trace(verbose=verbose, draw=draw)
dict_params['Norm'] = self.norm(verbose=verbose, draw=draw)
return dict_params
[docs]
def matrix(self, shape=None, shape_like=None):
"""Returns the numpy array of Mueller matrices.
Parameters:
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
Returns:
(numpy.array) 4x4xN numpy array.
"""
shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
if shape is not None and len(shape) > 1:
shape = tuple([4, 4] + list(shape))
M = np.reshape(self.parent.M, shape)
else:
M = self.parent.M
return M
[docs]
def components(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Extracts the matrix components of the Mueller matrix.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
M00 (numpy.ndarray): array of the 0, 0 element of the matrix.
M01 (numpy.ndarray): array of the 0, 1 element of the matrix.
...
M32 (numpy.ndarray): array of the 3, 2 element of the matrix.
M33 (numpy.ndarray): array of the 3, 3 element of the matrix.
"""
# Calculate the components
components = []
for ind1 in range(4):
for ind2 in range(4):
comp = self.parent.M[ind1, ind2, :]
if out_number and comp.size == 1:
comp = comp[0]
components.append(comp)
# Reshape if required
components = reshape(components,
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The matrix components of {} are:'.format(
self.parent.name)
title = ('M00', 'M01', 'M02', 'M03', 'M10', 'M11', 'M12', 'M13',
'M20', 'M21', 'M22', 'M23', 'M30', 'M31', 'M32', 'M33')
PrintParam(param=components,
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
# Return
return components
[docs]
def mean_transmission(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Gives the mean transmission coefficients.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): if True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
M00 = self.parent.M[0, 0, :]
# If the result is a number and the user asks for it, return a float
if out_number and M00.size == 1:
M00 = M00[0]
# Reshape if neccessary
M00 = reshape([M00],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The mean transmission of {} is:'.format(
self.parent.name)
PrintParam(param=M00,
shape=self.parent.shape,
title='Mean transmission',
heading=heading,
verbose=verbose,
draw=draw)
return M00
[docs]
def diattenuation_vector(self,
normalize=True,
as_Stokes=False,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Gives the diattenuation vector. The first dimension will always have size 3.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
normalize (bool): If True, normalizes the diattenuation vector to M00. Default: True.
as_Stokes (bool): If True, the vectors are inside a normlized Stokes object. Default: False.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): if True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or Stokes): Result.
"""
# Extract components
# TODO: No idea why deepcopy is required here
comp = deepcopy(
self.components(shape=shape,
shape_like=shape_like,
out_number=False))
# Normalize (avoiding dividing by 0)
if normalize:
cond = comp[0] > tol_default
if np.any(cond):
for ind in range(1, 4):
comp[ind][cond] = comp[ind][cond] / comp[0][cond]
# Transform in Stokes if needed
if as_Stokes:
S0 = np.sqrt(comp[1]**2 + comp[2]**2 + comp[3]**2)
D = Stokes().from_components([S0, comp[1], comp[2], comp[3]])
D.normalize()
else:
D = np.array([comp[1], comp[2], comp[3]])
# Print the result if required
if verbose or draw:
if as_Stokes:
D.parameters.components(verbose=verbose, draw=draw)
else:
heading = 'The diattenuation vector of {} is:'.format(
self.parent.name)
title = ('D[0]', 'D[1]', 'D[2]')
PrintParam(param=[D[0, :], D[1, :], D[2, :]],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
return D
[docs]
def polarizance_vector(self,
normalize=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Gives the polarizance vector. The first dimension will always have size 3.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
normalize (bool): If True, normalizes the diattenuation vector to M00. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray): Result.
"""
# Extract components
# TODO: No idea why deepcopy is required here
comp = deepcopy(
self.components(shape=shape,
shape_like=shape_like,
out_number=False))
# Normalize (avoiding dividing by 0)
if normalize:
cond = comp[0] > tol_default
if np.any(cond):
for ind in (4, 8, 12):
comp[ind][cond] = comp[ind][cond] / comp[0][cond]
P = np.array([comp[4], comp[8], comp[12]])
# Print the result if required
if verbose or draw:
heading = 'The polarizance vector of {} is:'.format(
self.parent.name)
title = ('P[0]', 'P[1]', 'P[2]')
PrintParam(param=[P[0, :], P[1, :], P[2, :]],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
return P
[docs]
def small_matrix(self,
normalize=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Gives the small matrix m. The first two dimensions will always have size 3.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
normalize (bool): If True, normalizes the diattenuation vector to M00. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray): Result.
"""
# Extract components
# TODO: No idea why deepcopy is required here
comp = deepcopy(
self.components(shape=shape,
shape_like=shape_like,
out_number=False))
# Normalize (avoiding dividing by 0)
if normalize:
cond = comp[0] > tol_default
for ind in (5, 6, 7, 9, 10, 11, 13, 14, 15):
comp[ind][cond] = comp[ind][cond] / comp[0][cond]
m = np.array([[comp[5], comp[6], comp[7]],
[comp[9], comp[10], comp[11]],
[comp[13], comp[14], comp[15]]])
# Print the result if required
if verbose or draw:
heading = 'The small matrix of {} is:'.format(self.parent.name)
title = ('m[0,0]', 'm[0,1]', 'm[0,2]', 'm[1,0]', 'm[1,1]',
'm[1,2]', 'm[2,0]', 'm[2,1]', 'm[2,2]')
PrintParam(param=[
m[0, 0, :], m[0, 1, :], m[0, 2, :], m[1, 0, :], m[1, 1, :],
m[1, 2, :], m[2, 0, :], m[2, 1, :], m[2, 2, :]
],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
return m
[docs]
def blocks(self,
normalize=True,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Method that gives the Mueller matrix block components: $M_{00}$ (mean transmission), $D$ (diattenuation vector), $P$ (polarizance vector) and $m$ (small matrix).
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
normalize (bool): If True, normalizes the diattenuation vector to M00. Default: True.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
m00 (float or numpy.ndarray): Average intensity.
D (numpy.ndarray): Diattenuation vectors 3xN.
P (numpy.ndarray): Diattenuation vector 3xN.
m (numpy.ndarray): Small m matrix 3x3xN.
"""
# Extract the variables from the parent object
M00 = self.mean_transmission(shape=shape,
shape_like=shape_like,
out_number=out_number)
D = self.diattenuation_vector(shape=shape, shape_like=shape_like)
P = self.polarizance_vector(shape=shape, shape_like=shape_like)
m = self.small_matrix(shape=shape, shape_like=shape_like)
# Print the result if required
if verbose or draw:
print('The block components of {} are:'.format(self.parent.name))
heading = ' - M_00 of {} is:'.format(self.parent.name)
title = ('M00')
PrintParam(param=[M00],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
heading = ' - Diattenuation vector of {} is:'.format(
self.parent.name)
title = ('D[0]', 'D[1]', 'D[2]')
PrintParam(param=[D[0, :], D[1, :], D[2, :]],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
heading = ' - Polarizance vector of {} is:'.format(
self.parent.name)
title = ('P[0]', 'P[1]', 'P[2]')
PrintParam(param=[P[0, :], P[1, :], P[2, :]],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
heading = ' - Small matrix of {} is:'.format(self.parent.name)
title = ('m[0,0]', 'm[0,1]', 'm[0,2]', 'm[1,0]', 'm[1,1]',
'm[1,2]', 'm[2,0]', 'm[2,1]', 'm[2,2]')
PrintParam(param=[
m[0, 0, :], m[0, 1, :], m[0, 2, :], m[1, 0, :], m[1, 1, :],
m[1, 2, :], m[2, 0, :], m[2, 1, :], m[2, 2, :]
],
shape=self.parent.shape,
title=title,
heading=heading,
verbose=verbose,
draw=draw)
return M00, D, P, m
[docs]
def global_phase(self,
give_nan=True,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the phase of J00 (which is the reference for global phase in py_pol model).
Parameters:
give_nan(bool): If False, NaN values are transformed into 0. Default: True.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
gp = self.parent.global_phase
# If the phase is None (unknown), give nans or zeros
if gp is None:
if self.parent.shape is None:
gp = np.array([0])
else:
gp = np.zeros(self.parent.size)
if give_nan:
gp = gp * np.nan
# If the result is a number and the user asks for it, return a float
if out_number and gp.size == 1:
gp = gp[0]
# Calculate Ez and reshape if required
gp = reshape([gp],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The global phase of {} is (deg):'.format(
self.parent.name)
PrintParam(param=(gp / degrees),
shape=self.parent.shape,
title=("Global phase (deg)"),
heading=heading,
verbose=verbose,
draw=draw)
return gp
[docs]
def inhomogeneity(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the inhomogeneity parameter. This parameter is 0 for homogeneous optical elements and 1 for totally inhomogeneous (degenerate) elements.
Note: The equation of the reference shows at least an incorrect result in the diattenuator retarders.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 119.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): if True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the components
c = self.components(shape=False, out_number=False)
det = self.det(shape=False, out_number=False)
tr = self.trace(shape=False, out_number=False)
# Calculate the parameter
cte = (tr + c[1] + c[4] + 1j * (c[11] + c[14])) / \
np.sqrt(2 * (c[0] + c[5] + c[1] + c[4]))
num = 4 * c[0] - np.abs(cte)**2 - np.abs(cte**2 - 4 * det**0.25)
den = 4 * c[0] - np.abs(cte)**2 + np.abs(cte**2 - 4 * det**0.25)
cond = np.abs(num) < tol_default
if np.any(cond):
num[cond] = 0 # Safety for -0 numeric error
eta = np.sqrt(num / den)
# Safety: 0/0 must be changed to 0
cond = (num < tol_default) and (den < tol_default)
if np.any(cond):
eta[cond] = 0
# If the result is a number and the user asks for it, return a float
if out_number and eta.size == 1:
eta = eta[0]
# Reshape if neccessary
eta = reshape([eta],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The inhomogeneity parameter of {} is:'.format(
self.parent.name)
PrintParam(param=eta,
shape=self.parent.shape,
title='Inhomogeneity parameter',
heading=heading,
verbose=verbose,
draw=draw)
return eta
[docs]
def diattenuation(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the diattenuation of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", pp. 200, CRC Press (2016)
Parameters:
out_number (bool): if True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the diattenuation
D = self.diattenuation_vector(shape=shape, shape_like=shape_like)
D = np.linalg.norm(D, axis=0)
# If the result is a number and the user asks for it, return a float
if out_number and D.size == 1:
D = D[0]
# Print the result if required
if verbose or draw:
heading = 'The diattenuation of {} is:'.format(self.parent.name)
PrintParam(param=D,
shape=self.parent.shape,
title='Diattenuation',
heading=heading,
verbose=verbose,
draw=draw)
return D
[docs]
def diattenuation_linear(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the linear diattenuation of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): if True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the diattenuation
D = self.diattenuation_vector(shape=shape, shape_like=shape_like)
D = np.linalg.norm(D[0:2, :], axis=0)
# If the result is a number and the user asks for it, return a float
if out_number and D.size == 1:
D = D[0]
# Print the result if required
if verbose or draw:
heading = 'The linear diattenuation of {} is:'.format(
self.parent.name)
PrintParam(param=D,
shape=self.parent.shape,
title='Linear diattenuation',
heading=heading,
verbose=verbose,
draw=draw)
return D
[docs]
def diattenuation_circular(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the circular diattenuation of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): if True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the diattenuation
D = self.diattenuation_vector(shape=shape, shape_like=shape_like)
D = np.abs(D[2, :])
# If the result is a number and the user asks for it, return a float
if out_number and D.size == 1:
D = D[0]
# Print the result if required
if verbose or draw:
heading = 'The circular diattenuation of {} is:'.format(
self.parent.name)
PrintParam(param=D,
shape=self.parent.shape,
title='Circular diattenuation',
heading=heading,
verbose=verbose,
draw=draw)
return D
[docs]
def polarizance(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the polarizance of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the polarizance
P = self.polarizance_vector(shape=shape, shape_like=shape_like)
P = np.linalg.norm(P, axis=0)
# If the result is a number and the user asks for it, return a float
if out_number and P.size == 1:
P = P[0]
# Reshape if neccessary
P = reshape([P],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The polarizance of {} is:'.format(self.parent.name)
PrintParam(param=P,
shape=self.parent.shape,
title='Polarizance',
heading=heading,
verbose=verbose,
draw=draw)
return P
[docs]
def polarizance_linear(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the linear polarizance of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the polarizance
P = self.polarizance_vector(shape=shape, shape_like=shape_like)
P = np.linalg.norm(P[0:2, :], axis=0)
# If the result is a number and the user asks for it, return a float
if out_number and P.size == 1:
P = P[0]
# Print the result if required
if verbose or draw:
heading = 'The linear polarizance of {} is:'.format(
self.parent.name)
PrintParam(param=P,
shape=self.parent.shape,
title='Linear polarizance',
heading=heading,
verbose=verbose,
draw=draw)
return P
[docs]
def polarizance_circular(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the linear polarizance of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): if True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the polarizance
P = self.polarizance_vector(shape=shape, shape_like=shape_like)
P = np.abs(P[2, :])
# If the result is a number and the user asks for it, return a float
if out_number and P.size == 1:
P = P[0]
# Print the result if required
if verbose or draw:
heading = 'The circular polarizance of {} is:'.format(
self.parent.name)
PrintParam(param=P,
shape=self.parent.shape,
title='Circular polarizance',
heading=heading,
verbose=verbose,
draw=draw)
return P
[docs]
def degree_polarizance(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the degree of polarizance of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the degree
P = self.polarizance(out_number=out_number,
shape=shape,
shape_like=shape_like)
D = self.diattenuation(out_number=out_number,
shape=shape,
shape_like=shape_like)
Dp = sqrt((P**2 + D**2) / 2)
# Print the result if required
if verbose or draw:
heading = 'The degree of polarizance of {} is:'.format(
self.parent.name)
PrintParam(param=Dp,
shape=self.parent.shape,
title='Degree of polarizance',
heading=heading,
verbose=verbose,
draw=draw)
return Dp
[docs]
def spheric_purity(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the degree of spheric purity of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016) pp 204.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the degree
m = self.small_matrix(shape=shape, shape_like=shape_like)
SP = np.linalg.norm(m, axis=(0, 1)) / sqrt(3)
# If the result is a number and the user asks for it, return a float
if out_number and SP.size == 1:
SP = SP[0]
# Print the result if required
if verbose or draw:
heading = 'The degree of spherical purity of {} is:'.format(
self.parent.name)
PrintParam(param=SP,
shape=self.parent.shape,
title='Spherical purity',
heading=heading,
verbose=verbose,
draw=draw)
return SP
# Similar to purity grades
[docs]
def retardance(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the retardance of the Mueller matrix of a pure retarder.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016) pp 129.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# TODO: (Jesus) We have to find a way to know if we are in the 180 -> 360 degrees region. Clue: In that case, opper triangular part is negative (except for ellipticity negative and phi lower than 90 deg).
# Calculate the retardance
m = self.small_matrix(shape=shape, shape_like=shape_like)
cosR = (np.trace(m) - 1) / 2
R = np.arccos(cosR)
# If the result is a number and the user asks for it, return a float
if out_number and R.size == 1:
R = R[0]
# Print the result if required
if verbose or draw:
heading = 'The retardance of {} is (deg):'.format(self.parent.name)
PrintParam(param=R / degrees,
shape=self.parent.shape,
title='Retardance (deg)',
heading=heading,
verbose=verbose,
draw=draw)
return R
# Polarization or despolarization
[docs]
def polarimetric_purity(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the degree of polarimetric purity of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the degree
Pp = self.degree_polarizance(out_number=out_number,
shape=shape,
shape_like=shape_like)
Ps = self.spheric_purity(out_number=out_number,
shape=shape,
shape_like=shape_like)
PP = sqrt(2. / 3. * Pp**2 + Ps**2)
# Print the result if required
if verbose or draw:
heading = 'The degree of polarimetric purity of {} is:'.format(
self.parent.name)
PrintParam(param=PP,
shape=self.parent.shape,
title='Polarimetric purity',
heading=heading,
verbose=verbose,
draw=draw)
return PP
[docs]
def depolarization_index(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the depolarization index of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
PP = self.polarimetric_purity(out_number=out_number,
shape=shape,
shape_like=shape_like)
DI = sqrt(1. - PP**2)
# Print the result if required
if verbose or draw:
heading = 'The depolarization index of {} is:'.format(
self.parent.name)
PrintParam(param=DI,
shape=self.parent.shape,
title='Depolarization index',
heading=heading,
verbose=verbose,
draw=draw)
return DI
# def depolarization_factors(self):
# """Calculates the Euclidean distance and depolarization factor
#
# References:
# Handbook of Optics vol 2. 22.49 (46 and 47)
#
# Parameters:
# M (4x4 numpy.matrix): Mueller matrix
#
# Returns:
# (float): Euclidean distance of the normalized Mueller matrix from an ideal depolarizer
# (float): Dep(M) depolarization of the matrix
# """
# # TODO: (Jesus) Check if Mnorm must be used instead of M
#
# M = self.M
# quadratic_sum = (array(M)**2).sum()
# euclidean_distance = sqrt(quadratic_sum - M[0, 0]**2) / M[0, 0]
# depolarization = 1 - euclidean_distance / sqrt(3)
#
# return euclidean_distance, depolarization
# Polarimetric purity
[docs]
def polarimetric_purity_indices(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates the polarimetric purity indices of the Mueller matrices.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 208.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
P1, P2, P3 (numpy.ndarray or float): Results.
"""
# Calculate the covariance matrix
H = self.parent.covariance_matrix(keep=True)
# Calculate parameters of covariance matrix
th = np.abs(H.parameters.trace(shape=False, out_number=out_number))
v1, v2, v3, v4 = H.parameters.eigenvalues(shape=False,
out_number=False)
# Order the eigenvalues
vals = np.array([np.abs(v1), np.abs(v2), np.abs(v3), np.abs(v4)])
vals = np.sort(vals, axis=0)
# Calculate indices
P1 = (vals[3, :] - vals[2, :]) / th
P2 = (vals[3, :] + vals[2, :] - 2 * vals[1, :]) / th
P3 = (vals[3, :] + vals[2, :] + vals[1, :] - 3 * vals[0, :]) / th
# Size 1 objects need some care
if P1.size == 1 and P1.ndim == 1:
P1, P2, P3 = (P1[0], P2[0], P3[0])
# Reshape if neccessary
P1, P2, P3 = reshape([P1, P2, P3],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
# Eigenvalues
heading = 'The polarimetric purity indices of {} are:'.format(
self.parent.name)
PrintParam(param=(P1, P2, P3),
shape=self.parent.shape,
title=('P1', 'P2', 'P3'),
heading=heading,
verbose=verbose,
draw=draw)
return P1, P2, P3
[docs]
def transmissions(self,
kind='Intensity',
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculate the maximum and minimum transmitance of an optical element.
References:
Handbook of Optics vol 2. 22.32 (eq.38)
Parameters:
kind (str): There are three options, FIELD, INTENSITY or ALL.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
T_max (numpy.ndarray): Maximum intensity transmission.
T_min (numpy.ndarray): Minimum intensity transmission.
p1 (numpy.ndarray): Maximum field transmission.
p2 (numpy.ndarray): Minimum field transmission.
"""
# Take the parameters
M00 = self.mean_transmission(shape=shape,
shape_like=shape_like,
out_number=out_number)
D = self.diattenuation(shape=shape,
shape_like=shape_like,
out_number=out_number)
# Calculate the transmissions
T_max = M00 * (1 + D)
T_min = M00 * (1 - D)
if kind.upper() in ('FIELD', 'ALL'):
p1 = np.sqrt(T_max)
p2 = np.sqrt(T_min)
# Print the result if required
if verbose or draw:
# Intensity
if kind.upper() in ('INTENSITY', 'ALL'):
heading = 'The intensity transmissions of {} are:'.format(
self.parent.name)
PrintParam(param=(T_max, T_min),
shape=self.parent.shape,
title=('Maximum (int.)', 'Minimum (int.)'),
heading=heading,
verbose=verbose,
draw=draw)
# Field
if kind.upper() in ('FIELD', 'ALL'):
heading = 'The field transmissions of {} are:'.format(
self.parent.name)
PrintParam(param=(p1, p2),
shape=self.parent.shape,
title=('Maximum (int.)', 'Minimum (int.)'),
heading=heading,
verbose=verbose,
draw=draw)
# Return
ret = []
if kind.upper() in ('INTENSITY', 'ALL'):
ret += [T_max, T_min]
if kind.upper() in ('FIELD', 'ALL'):
ret += [p1, p2]
return ret
[docs]
def retardance_vector(self,
kind='norm',
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculate the retardance of a vector.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 128-130.
Args:
kind (string): Identifies the type of retardance vector. There are three possibilities: NORM (also called Pauli vector), STRAIGHT or COMPLETE. Default: 'norm'.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray): 3xN array of the result.
"""
# Extract the necessary parameters
R = self.retardance(shape=False, out_number=False)
m = self.small_matrix(shape=False)
Rv = np.zeros((3, self.parent.size))
# Transform depending on the type of vector
if kind in ('complete', 'COMPLETE', 'Complete'):
cte = R
elif kind in ('straight', 'STRAIGHT', 'Straight'):
cte = R / np.pi
else:
cte = np.ones_like(R)
# General case: retardance different than 0 or 180º
cond_G = (R >= tol_default) * (R <= 180 * degrees - tol_default)
if np.any(cond_G):
sR = 2 * np.sin(R[cond_G]) / cte[cond_G]
Rv[0, cond_G] = (m[1, 2, cond_G] - m[2, 1, cond_G]) / sR
Rv[1, cond_G] = (m[2, 0, cond_G] - m[0, 2, cond_G]) / sR
Rv[2, cond_G] = (m[0, 1, cond_G] - m[1, 0, cond_G]) / sR
# Particular case: R = 0
cond_Z = R < tol_default
if np.any(cond_Z):
# If R = 0, we don't have a retarder after all, so any vector is valid
Rv[0, cond_Z] = cte[cond_Z]
# Particular case: R = 180º
cond_P = R > 180 * degrees - tol_default
if np.any(cond_P):
# This case is more tricky, we have to calculate Rv from its eigenstates (supposing that it is a pure retarder)
S1, _ = self.eigenstates(shape=False)
alpha, delay = S1.parameters.charac_angles(shape=False,
out_number=False)
Rv[0, cond_P] = cte[cond_P] * np.cos(2 * alpha[cond_P])
Rv[1, cond_P] = cte[cond_P] * \
np.sin(2 * alpha[cond_P]) * np.cos(2 * delay[cond_P])
Rv[2, cond_P] = cte[cond_P] * \
np.sin(2 * alpha[cond_P]) * np.sin(2 * delay[cond_P])
# Print the result if required
if verbose or draw:
# Extract the components
r0, r1, r2 = (Rv[0, :], Rv[1, :], Rv[2, :])
r0, r1, r2 = reshape([r0, r1, r2],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print
heading = 'The retardance vector ({}) components components of {} are:'.format(
kind, self.parent.name)
PrintParam(param=[r0, r1, r2],
shape=self.parent.shape,
title=('r0', 'r1', 'r2'),
heading=heading,
verbose=verbose,
draw=draw)
# Reshape the vector if neccessary
shape = take_shape([self.parent, shape_like, shape])
if shape is not None and len(shape) > 1:
shape = np.insert(shape, 0, 3)
Rv = np.reshape(Rv, shape)
return Rv
[docs]
def eig(self,
values_as_matrix=False,
vectors_as_matrix=False,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the eigenvalues and eigenvectors of the Mueller matrices.
Parameters:
values_as_matrix (bool): If True, the eigenvalues output is a numpy.ndarray instead of a list of arrays. Default: False.
vectors_as_matrix (bool): If True, the eigenvectors output is a numpy.ndarray instead of a list of arrays. Default: False.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
eigenvalues (list or numpy.ndarray): List with the four eigenvalues (if values_as_matrix is False) or 4xN array (if values_as_matrix is True).
eigenvectors (list or numpy.ndarray): List with the four eigenvectors (if vectors_as_matrix is False) or 4x4xN array (if vectors_as_matrix is True).
"""
# Calculate the eig
M = np.moveaxis(self.parent.M, -1, 0)
val, vect = np.linalg.eig(M)
# Order the eigenvalues in the py_pol way
if ~values_as_matrix or verbose or draw:
v1, v2, v3, v4 = (val[:, 0], val[:, 1], val[:, 2], val[:, 3])
# Size 1 objects need some care
if v1.size == 1 and v1.ndim > 1:
v1, v2, v3, v4 = (v1[0], v2[0], v3[0], v4[0])
if out_number and v1.size == 1:
v1, v2, v3, v4 = (v1[0], v2[0], v3[0], v4[0])
# Reshape if neccessary
eigenvalues = reshape([v1, v2, v3, v4],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Reshape eigenvalues if output as matrix
if values_as_matrix:
val = np.moveaxis(val, -1, 0)
new_shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
if new_shape is not None and len(new_shape) > 1:
new_shape = [4] + list(new_shape)
val = np.reshape(val, new_shape)
# Reshape eigenvectors in pypol way
if ~vectors_as_matrix or verbose or draw:
e1 = np.array(
[vect[:, 0, 0], vect[:, 1, 0], vect[:, 2, 0], vect[:, 3, 0]])
e2 = np.array(
[vect[:, 0, 1], vect[:, 1, 1], vect[:, 2, 1], vect[:, 3, 1]])
e3 = np.array(
[vect[:, 0, 2], vect[:, 1, 2], vect[:, 2, 2], vect[:, 3, 2]])
e4 = np.array(
[vect[:, 0, 3], vect[:, 1, 3], vect[:, 2, 3], vect[:, 3, 3]])
eigenvectors = []
for elem in (e1, e2, e3, e4):
for ind in range(4):
eigenvectors.append(elem[ind, :])
eigenvectors = reshape(eigenvectors,
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
new_shape = [4] + list(eigenvectors[0].shape)
if len(new_shape) > 2:
e1 = np.reshape(e1, new_shape)
e2 = np.reshape(e2, new_shape)
e3 = np.reshape(e3, new_shape)
e4 = np.reshape(e4, new_shape)
# Reshape eigenvectors if output as matrix
if vectors_as_matrix:
vect = np.moveaxis(vect, -2, 0)
vect = np.moveaxis(vect, -1, 1)
new_shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
if new_shape is not None and len(new_shape) > 1:
new_shape = [4, 4] + list(new_shape)
vect = np.reshape(vect, new_shape)
# Print the result if required
if verbose or draw:
# Eigenvalues
heading = 'The eigenvalues of {} are:'.format(self.parent.name)
PrintParam(param=eigenvalues,
shape=self.parent.shape,
title=('v1', 'v2', 'v3', 'v4'),
heading=heading,
verbose=verbose,
draw=draw)
# Eigenvectors
heading = 'The eigenvectors of {} are:'.format(self.parent.name)
PrintParam(param=eigenvectors,
shape=self.parent.shape,
title=('e1 I', 'e1 Q', 'e1 U', 'e1 V', 'e2 I', 'e2 Q',
'e2 U', 'e2 V', 'e3 I', 'e3 Q', 'e3 U', 'e3 V',
'e4 I', 'e4 Q', 'e4 U', 'e4 V'),
heading=heading,
verbose=verbose,
draw=draw)
# Return
if values_as_matrix:
eigenvalues = val
if vectors_as_matrix:
eigenvectors = vect
else:
eigenvectors = (e1, e2, e3, e4)
return eigenvalues, eigenvectors
[docs]
def eigenvectors(self,
vectors_as_matrix=False,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the eigenvectors of the Mueller matrices.
Parameters:
vectors_as_matrix (bool): If True, the eigenvectors output is a numpy.ndarray instead of a list of arrays. Default: False.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray): 4x4xN eigenvector matrix (if vectors_as_matrix is True).
e1, e2, e3, e4 (numpy.ndarray): 4xN eigenvector matrices (if vectors_as_matrix is False).
"""
# Calculate
M = np.moveaxis(self.parent.M, -1, 0)
_, vect = np.linalg.eig(M)
# Reshape eigenvectors in pypol way
if ~vectors_as_matrix or verbose or draw:
e1 = np.array(
[vect[:, 0, 0], vect[:, 1, 0], vect[:, 2, 0], vect[:, 3, 0]])
e2 = np.array(
[vect[:, 0, 1], vect[:, 1, 1], vect[:, 2, 1], vect[:, 3, 1]])
e3 = np.array(
[vect[:, 0, 2], vect[:, 1, 2], vect[:, 2, 2], vect[:, 3, 2]])
e4 = np.array(
[vect[:, 0, 3], vect[:, 1, 3], vect[:, 2, 3], vect[:, 3, 3]])
eigenvectors = []
for elem in (e1, e2, e3, e4):
for ind in range(4):
eigenvectors.append(elem[ind, :])
eigenvectors = reshape(eigenvectors,
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# new_shape = [4] + list(eigenvectors[0].shape)
# if len(new_shape) > 2:
# e1 = np.reshape(e1, new_shape)
# e2 = np.reshape(e2, new_shape)
# e3 = np.reshape(e3, new_shape)
# e4 = np.reshape(e4, new_shape)
# Reshape eigenvectors if output as matrix
if vectors_as_matrix:
vect = np.moveaxis(vect, -2, 0)
vect = np.moveaxis(vect, -1, 1)
new_shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
if new_shape is not None and len(new_shape) > 1:
new_shape = [4, 4] + list(new_shape)
vect = np.reshape(vect, new_shape)
# Print the result if required
if verbose or draw:
# Eigenvectors
# print_vect = []
heading = 'The eigenvectors of {} are:'.format(self.parent.name)
# for elem in eigenvectors:
# for ind in range(4):
# print_vect.append(elem[ind, :])
# print_vect = reshape(
# print_vect,
# shape_like=shape_like,
# shape_fun=shape,
# obj=self.parent)
PrintParam(param=eigenvectors,
shape=self.parent.shape,
title=('e1 I', 'e1 Q', 'e1 U', 'e1 V', 'e2 I', 'e2 Q',
'e2 U', 'e2 V', 'e3 I', 'e3 Q', 'e3 U', 'e3 V',
'e4 I', 'e4 Q', 'e4 U', 'e4 V'),
heading=heading,
verbose=verbose,
draw=draw)
# Return
if vectors_as_matrix:
eigenvectors = vect
return eigenvectors
[docs]
def eigenstates(self,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the eigenstates of the optical object. It must be done in Jones formalism, so it is only valid for pure Mueller matrices.
Parameters:
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
S1, S2, S3, S4 (Stokes): eigenstates.
"""
# Calculate the eigenstates in Jones
J = Jones_matrix(self.parent.name)
J.from_Mueller(self.parent, shape_like=self.parent)
E1, E2 = J.parameters.eigenstates(shape=shape, shape_like=shape_like)
# Transfrom to Stokes
S1, S2 = create_Stokes(('1st eigenstate', '2nd eigenstate'))
S1.from_Jones(E1)
S2.from_Jones(E2)
# Print the result if required
if verbose or draw:
# Extract the info
heading = 'The eigenstates of {} are:'.format(self.parent.name)
c1 = S1.parameters.components()
c2 = S2.parameters.components()
components = list(c1) + list(c2)
# Print
PrintParam(param=components,
shape=self.parent.shape,
title=('S1 I', 'S1 Q', 'S1 U', 'S1 V', 'S2 I', 'S2 Q',
'S2 U', 'S2 V'),
heading=heading,
verbose=verbose,
draw=draw)
return S1, S2
[docs]
def eigenvalues(self,
values_as_matrix=False,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the eigenvalues and eigenstates of the Jones object.
Parameters:
values_as_matrix (bool): If True, the eigenvalues output is a numpy.ndarray instead of a list of arrays. Default: False.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
v (numpy.ndarray): 4xN eigenvalues matrix (if values_as_matrix is True).
v1, v2, v3, v4 (numpy.ndarray or float): Individual eigenvalues (if values_as_matrix is False).
"""
# Calculate the eigenvalues
M = np.moveaxis(self.parent.M, -1, 0)
val = np.linalg.eigvals(M)
# Order the eigenvalues in the py_pol way
if ~values_as_matrix or verbose or draw:
v1, v2, v3, v4 = (val[:, 0], val[:, 1], val[:, 2], val[:, 3])
# Size 1 objects need some care
if v1.size == 1 and v1.ndim > 1:
v1, v2, v3, v4 = (v1[0], v2[0], v3[0], v4[0])
if out_number and v1.size == 1:
v1, v2, v3, v4 = (v1[0], v2[0], v3[0], v4[0])
# Reshape if neccessary
eigenvalues = reshape([v1, v2, v3, v4],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Reshape eigenvalues if output as matrix
if values_as_matrix:
val = np.moveaxis(val, -1, 0)
new_shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
if new_shape is not None and len(new_shape) > 1:
new_shape = [4] + list(new_shape)
val = np.reshape(val, new_shape)
else:
v1, v2, v3, v4 = reshape([v1, v2, v3, v4],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
# Eigenvalues
heading = 'The eigenvalues of {} are:'.format(self.parent.name)
PrintParam(param=(v1, v2, v3, v4),
shape=self.parent.shape,
title=('v1', 'v2', 'v3', 'v4'),
heading=heading,
verbose=verbose,
draw=draw)
# Return
if values_as_matrix:
eigenvalues = val
else:
eigenvalues = (v1, v2, v3, v4)
return eigenvalues
[docs]
def det(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the determinants of the Mueller matrices.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray, float or complex): Result.
"""
# Calculate the eigenstates
M = np.moveaxis(self.parent.M, -1, 0)
det = np.linalg.det(M)
# If the result is a number and the user asks for it, return a float
if out_number and det.size == 1:
det = det[0]
# Reshape if neccessary
det = reshape([det],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The determinant of {} is:'.format(self.parent.name)
PrintParam(param=det,
shape=self.parent.shape,
title='Determinant',
heading=heading,
verbose=verbose,
draw=draw)
return det
[docs]
def trace(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the trace of the Mueller matrices.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or float): Result.
"""
# Calculate the eigenstates
trace = np.trace(self.parent.M)
# If the result is a number and the user asks for it, return a float
if out_number and trace.size == 1:
trace = trace[0]
# Reshape if neccessary
trace = reshape([trace],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The trace of {} is:'.format(self.parent.name)
PrintParam(param=trace,
shape=self.parent.shape,
title='Trace',
heading=heading,
verbose=verbose,
draw=draw)
return trace
[docs]
def norm(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Calculates the Frobenius norm of the Mueller matrices.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray, float or complex): Result.
"""
# Calculate the eigenstates
norm = np.linalg.norm(self.parent.M, axis=(0, 1))
# If the result is a number and the user asks for it, return a float
if out_number and norm.size == 1:
norm = norm[0]
# Reshape if neccessary
norm = reshape([norm],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = 'The norm of {} is:'.format(self.parent.name)
PrintParam(param=norm,
shape=self.parent.shape,
title='Norm',
heading=heading,
verbose=verbose,
draw=draw)
return norm
[docs]
class Analysis_Mueller(object):
"""Class for Analysis of Mueller Analysis
Parameters:
mueller_matrix (Mueller_matrix): Mueller Matrix
Attributes:
self.marent (Mueller): parent object.
"""
def __init__(self, parent):
self.parent = parent
[docs]
def depolarizer(self,
transmissions='all',
angles="all",
depolarization='all',
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates some of the parameters from the Mueller matrix of a diattenuator.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
transmissions (string): Determines the type of transmission output, FIELD, INTENSITY or ALL. Default: All.
angles (string): Determines the type of angles output, CHARAC (characteristic angles), AZIMUTH (azimuth and ellipticity) or ALL. Default: All.
depolarization (string): Determines the type of depolarization information: INDEX, FACTORS or ALL. Default: All.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
trans_D, transP (list). Transmissions calculated from the diattenuation and polarizance vectors (in that order). Depending on the input parameter transmissions, it can contain the field transmissions, the intensity transmissions or both.
angles_D, angles_P (list). Angles of the transmission eigenstate calculated from the diattenuation and polarizance vectors respectively. Depending on the input parameter angles, it can contain the characteristic angles, the azimuth and ellipticity, or all of them.
depol (numpy.ndarray or float): Depolarization index.
S (list): List with the three principal states.
"""
# Extract the information
D = self.parent.parameters.diattenuation(shape=shape,
shape_like=shape_like,
out_number=out_number)
Dv = self.parent.parameters.diattenuation_vector(shape=shape,
shape_like=shape_like)
P = self.parent.parameters.polarizance(shape=shape,
shape_like=shape_like,
out_number=out_number)
Pv = self.parent.parameters.polarizance_vector(shape=shape,
shape_like=shape_like)
M00 = self.parent.parameters.mean_transmission(shape=shape,
shape_like=shape_like,
out_number=out_number)
m = self.parent.parameters.small_matrix(shape=False)
depol = self.parent.parameters.depolarization_index(
shape=shape, shape_like=shape_like, out_number=out_number)
# Calculate transmissions
Tmax = M00 * (1 + D)
Tmin = M00 * (1 - D)
p1, p2 = (np.sqrt(Tmax), np.sqrt(Tmin))
trans_D, title_trans = ([], [])
if transmissions in ('ALL', 'All', 'all', 'INTENSITY', 'Intensity',
'intensity'):
trans_D += [Tmax, Tmin]
title_trans += ['Max. transmission', 'Min. transmission']
if transmissions in ('ALL', 'All', 'all', 'FIELD', 'Field', 'field'):
trans_D += [p1, p2]
title_trans += ['p1', 'p2']
Tmax = M00 * (1 + P)
Tmin = M00 * (1 - P)
p1, p2 = (np.sqrt(Tmax), np.sqrt(Tmin))
trans_P = []
if transmissions in ('ALL', 'All', 'all', 'INTENSITY', 'Intensity',
'intensity'):
trans_P += [Tmax, Tmin]
if transmissions in ('ALL', 'All', 'all', 'FIELD', 'Field', 'field'):
trans_P += [p1, p2]
# Calculate angles
azimuth, ellipticity = extract_azimuth_elipt(Dv, use_nan=False)
alpha, delay = azimuth_elipt_2_charac_angles(azimuth, ellipticity)
ang_D, title_ang = ([], [])
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac', 'charac'):
ang_D += [alpha, delay]
title_ang += ['Alpha', 'Delay']
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth', 'azimuth'):
ang_D += [azimuth, ellipticity]
title_ang += ['Azimuth', 'Ellipticity angle']
azimuth, ellipticity = extract_azimuth_elipt(Pv, use_nan=False)
alpha, delay = azimuth_elipt_2_charac_angles(azimuth, ellipticity)
ang_P = []
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac', 'charac'):
ang_P += [alpha, delay]
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth', 'azimuth'):
ang_P += [azimuth, ellipticity]
# Calculate depolarization factors and eigenstates
m = np.moveaxis(m, -1, 0)
val, vect = np.linalg.eig(m)
# val, vect = order_eig(val, vect, 'reverse')
d1, d2, d3 = (val[:, 0], val[:, 1], val[:, 2])
S1, S2, S3 = create_Stokes(
('First principal state', 'Second principal state',
'Third principal state'))
S1.from_components([1, vect[:, 0, 0], vect[:, 1, 0], vect[:, 2, 0]])
S2.from_components([1, vect[:, 0, 1], vect[:, 1, 1], vect[:, 2, 1]])
S3.from_components([1, vect[:, 0, 2], vect[:, 1, 2], vect[:, 2, 2]])
# Reshaoe them
d1, d2, d3 = reshape([d1, d2, d3],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
S1.shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
S2.shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
S3.shape, _ = select_shape(obj=self.parent,
shape_fun=shape,
shape_like=shape_like)
# Save them
principal_states = []
depolar = []
if depolarization.upper() in ('ALL', 'INDEX'):
depolar += [depol]
if depolarization.upper() in ('ALL', 'FACTORS'):
depolar = [d1, d2, d3]
principal_states += [S1, S2, S3]
# Print the result if required
if verbose or draw:
# Transform angles to degrees for representation
angles_rep_D = []
for a in ang_D:
angles_rep_D.append(a / degrees)
angles_rep_P = []
for a in ang_P:
angles_rep_P.append(a / degrees)
print('\nAnalysis of {} as depolarizer:'.format(self.parent.name))
# Depolarization index
if depolarization.upper() in ('ALL', 'INDEX'):
heading = '- Depolarization index of {} is:'.format(
self.parent.name)
PrintParam(param=[depol],
shape=self.parent.shape,
title=('Depolarization index'),
heading=heading,
verbose=verbose,
draw=draw)
# Depolarization factors and eigenstates
if depolarization.upper() in ('ALL', 'FACTORS'):
heading = '- First depolarization factor of {} is:'.format(
self.parent.name)
PrintParam(param=[d1],
shape=self.parent.shape,
title=('1st depolarization factor'),
heading=heading,
verbose=verbose,
draw=draw)
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac',
'charac'):
_ = S1.parameters.charac_angles(verbose=verbose, draw=draw)
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth',
'azimuth'):
_ = S1.parameters.azimuth_ellipticity(verbose=verbose,
draw=draw)
heading = '- Second depolarization factor of {} is:'.format(
self.parent.name)
PrintParam(param=[d2],
shape=self.parent.shape,
title=('2nd depolarization factor'),
heading=heading,
verbose=verbose,
draw=draw)
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac',
'charac'):
_ = S2.parameters.charac_angles(verbose=verbose, draw=draw)
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth',
'azimuth'):
_ = S2.parameters.azimuth_ellipticity(verbose=verbose,
draw=draw)
heading = '- Third depolarization factor of {} is:'.format(
self.parent.name)
PrintParam(param=[d3],
shape=self.parent.shape,
title=('3rd depolarization factor'),
heading=heading,
verbose=verbose,
draw=draw)
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac',
'charac'):
_ = S3.parameters.charac_angles(verbose=verbose, draw=draw)
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth',
'azimuth'):
_ = S3.parameters.azimuth_ellipticity(verbose=verbose,
draw=draw)
# Depolarizers usually have only diattenuation or polarizance.
if np.any(D > tol_default):
heading = '- Transmissions of {} from diattenuation are:'.format(
self.parent.name)
PrintParam(param=trans_D,
shape=self.parent.shape,
title=title_trans,
heading=heading,
verbose=verbose,
draw=draw)
heading = '- Angles of {} from polarizance are:'.format(
self.parent.name)
PrintParam(param=angles_rep_P,
shape=self.parent.shape,
title=title_ang,
heading=heading,
verbose=verbose,
draw=draw)
else:
print('- {} has no diattenuation.'.format(self.parent.name))
if np.any(P > tol_default):
heading = '- Transmissions of {} from polarizance are:'.format(
self.parent.name)
PrintParam(param=trans_P,
shape=self.parent.shape,
title=title_trans,
heading=heading,
verbose=verbose,
draw=draw)
heading = '- Angles of {} from diattenuation are:'.format(
self.parent.name)
PrintParam(param=angles_rep_D,
shape=self.parent.shape,
title=title_ang,
heading=heading,
verbose=verbose,
draw=draw)
else:
print('- {} has no polarizance.'.format(self.parent.name))
# Output
return trans_D, trans_P, ang_D, ang_P, depolar, principal_states
[docs]
def diattenuator(self,
transmissions='all',
angles="all",
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates all the parameters from the Mueller Matrix of a pure homogeneous diattenuator (using the diattenuance vector). If the object is not a pure homogenous diattenuator, some parameters may be wrong.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
transmissions (string): Determines the type of transmission output, FIELD, INTENSITY or ALL. Default: All.
angles (string): Determines the type of angles output, CHARAC (characteristic angles), AZIMUTH (azimuth and ellipticity) or ALL. Default: All.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
transmissions (list). Transmissions of the diattenuator. Depending on the input parameter transmissions, it can contain the field transmissions, the intensity transmissions or both.
angles (list). Angles of the transmission eigenstate. Depending on the input parameter angles, it can contain the characteristic angles, the azimuth and ellipticity, or all of them.
"""
# Extract the diattenuation and diattenuation vector
D = self.parent.parameters.diattenuation(shape=shape,
shape_like=shape_like,
out_number=out_number)
Dv = self.parent.parameters.diattenuation_vector(shape=shape,
shape_like=shape_like)
M00 = self.parent.parameters.mean_transmission(shape=shape,
shape_like=shape_like,
out_number=out_number)
# Calculate transmissions
Tmax = M00 * (1 + D)
Tmin = M00 * (1 - D)
p1, p2 = (np.sqrt(Tmax), np.sqrt(Tmin))
trans, title_trans = ([], [])
if transmissions in ('ALL', 'All', 'all', 'INTENSITY', 'Intensity',
'intensity'):
trans += [Tmax, Tmin]
title_trans += ['Max. transmission', 'Min. transmission']
if transmissions in ('ALL', 'All', 'all', 'FIELD', 'Field', 'field'):
trans += [p1, p2]
title_trans += ['p1', 'p2']
# Calculate angles
azimuth, ellipticity = extract_azimuth_elipt(Dv, use_nan=False)
alpha, delay = azimuth_elipt_2_charac_angles(azimuth, ellipticity)
ang, title_ang = ([], [])
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac', 'charac'):
ang += [alpha, delay]
title_ang += ['Alpha', 'Delay']
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth', 'azimuth'):
ang += [azimuth, ellipticity]
title_ang += ['Azimuth', 'Ellipticity angle']
# Print the result if required
if verbose or draw:
# Transform angles to degrees for representation
angles_rep = []
for a in ang:
angles_rep.append(a / degrees)
print('\nAnalysis of {} as diattenuator:\n'.format(
self.parent.name))
heading = '- Transmissions of {} are:'.format(self.parent.name)
PrintParam(param=trans,
shape=self.parent.shape,
title=title_trans,
heading=heading,
verbose=verbose,
draw=draw)
heading = '- Angles of {} are:'.format(self.parent.name)
PrintParam(param=angles_rep,
shape=self.parent.shape,
title=title_ang,
heading=heading,
verbose=verbose,
draw=draw)
# Output
return trans, ang
[docs]
def polarizer(self,
transmissions='all',
angles="all",
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates all the parameters from the Mueller Matrix of a pure homogeneous diattenuator (using the polarizance vector). If the object is not a pure homogenous diattenuator, some parameters may be wrong.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
transmissions (string): Determines the type of transmission output, FIELD, INTENSITY or ALL. Default: All.
angles (string): Determines the type of angles output, CHARAC (characteristic angles), AZIMUTH (azimuth and ellipticity) or ALL. Default: All.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
transmissions (list). Transmissions of the diattenuator. Depending on the input parameter transmissions, it can contain the field transmissions, the intensity transmissions or both.
angles (list). Angles of the transmission eigenstate. Depending on the input parameter angles, it can contain the characteristic angles, the azimuth and ellipticity, or all of them.
"""
# Extract the diattenuation and diattenuation vector
P = self.parent.parameters.polarizance(shape=shape,
shape_like=shape_like,
out_number=out_number)
Pv = self.parent.parameters.polarizance_vector(shape=shape,
shape_like=shape_like)
M00 = self.parent.parameters.mean_transmission(shape=shape,
shape_like=shape_like,
out_number=out_number)
# Calculate transmissions
Tmax = M00 * (1 + P)
Tmin = M00 * (1 - P)
p1, p2 = (np.sqrt(Tmax), np.sqrt(Tmin))
trans, title_trans = ([], [])
if transmissions in ('ALL', 'All', 'all', 'INTENSITY', 'Intensity',
'intensity'):
trans += [Tmax, Tmin]
title_trans += ['Max. transmission', 'Min. transmission']
if transmissions in ('ALL', 'All', 'all', 'FIELD', 'Field', 'field'):
trans += [p1, p2]
title_trans += ['p1', 'p2']
# Calculate angles
azimuth, ellipticity = extract_azimuth_elipt(Pv)
alpha, delay = azimuth_elipt_2_charac_angles(azimuth, ellipticity)
ang, title_ang = ([], [])
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac', 'charac'):
ang += [alpha, delay]
title_ang += ['Alpha', 'Delay']
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth', 'azimuth'):
ang += [azimuth, ellipticity]
title_ang += ['Azimuth', 'Ellipticity angle']
# Print the result if required
if verbose or draw:
# Transform angles to degrees for representation
angles_rep = []
for a in ang:
angles_rep.append(a / degrees)
print('\nAnalysis of {} as polarizer:\n'.format(self.parent.name))
heading = '- Transmissions of {} are:'.format(self.parent.name)
PrintParam(param=trans,
shape=self.parent.shape,
title=title_trans,
heading=heading,
verbose=verbose,
draw=draw)
heading = '- Angles of {} are:'.format(self.parent.name)
PrintParam(param=angles_rep,
shape=self.parent.shape,
title=title_ang,
heading=heading,
verbose=verbose,
draw=draw)
# Output
return trans, ang
[docs]
def retarder(self,
angles="all",
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Calculates all the parameters from the Mueller Matrix of a pure homogeneous retarder. If the object is not a pure homogenous retarder, some parameters may be wrong.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016)
Parameters:
angles (string): Determines the type of angles output, CHARAC (characteristic angles), AZIMUTH (azimuth and ellipticity) or ALL. Default: All.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
R (np.ndarray). Retardance.
angles (list). Angles of the transmission eigenstate. Depending on the input parameter angles, it can contain the characteristic angles, the azimuth and ellipticity, or all of them.
"""
# Extract the diattenuation and diattenuation vector
R = self.parent.parameters.retardance(shape=shape,
shape_like=shape_like,
out_number=out_number)
Rv = self.parent.parameters.retardance_vector(shape=shape,
shape_like=shape_like)
# Calculate angles
azimuth, ellipticity = extract_azimuth_elipt(Rv)
alpha, delay = azimuth_elipt_2_charac_angles(azimuth, ellipticity)
ang, title_ang = ([], [])
if angles in ('ALL', 'All', 'all', 'CHARAC', 'Charac', 'charac'):
ang += [alpha, delay]
title_ang += ['Alpha', 'Delay']
if angles in ('ALL', 'All', 'all', 'AZIMUTH', 'Azimuth', 'azimuth'):
ang += [azimuth, ellipticity]
title_ang += ['Azimuth', 'Ellipticity angle']
# Print the result if required
if verbose or draw:
# Transform angles to degrees for representation
angles_rep = []
for a in ang:
angles_rep.append(a / degrees)
print('\nAnalysis of {} as retarder:\n'.format(self.parent.name))
heading = '- Retardance of {} is:'.format(self.parent.name)
PrintParam(param=(R / degrees),
shape=self.parent.shape,
title=('Retardance'),
heading=heading,
verbose=verbose,
draw=draw)
heading = '- Angles of {} are:'.format(self.parent.name)
PrintParam(param=angles_rep,
shape=self.parent.shape,
title=title_ang,
heading=heading,
verbose=verbose,
draw=draw)
# Output
return R, ang
# # Matrix filtering
[docs]
def filter_purify_number(self, Neig=3, keep=False):
"""Method that filters experimental errors by making zero a certain number of eigenvalues of the covariance matrix.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016) pp 226.
Parameters:
Neig (int): Number of eigenvalues (1-3) of the covariant matrix to be made 0. Default: 3.
keep (bool): If True, the original object won't be altered. Default: False.
Returns:
(Mueller): Result.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.parent.copy()
else:
new_obj = self.parent
old_shape = new_obj.shape
# Calculate covariance matrix eigenvalues
new_obj.covariance_matrix()
val, vect = new_obj.parameters.eig(values_as_matrix=True,
vectors_as_matrix=True,
shape=False)
# Find the order of eigenvalues
order = np.argsort(val, axis=0)
# Make 0 the desired eigenvalues
val[order < Neig] = 0
# Recompose the matrix
H, diag, Ht = create_Mueller(N=3)
diag.from_components([val[0, :]] + 4 * [0] + [val[1, :]] + 4 * [0] +
[val[2, :]] + 4 * [0] + [val[3, :]])
H.from_matrix(vect)
Ht.from_matrix(np.conj(np.transpose(vect, axes=(1, 0, 2))))
result = H * diag * Ht
new_obj.from_covariance(result, shape=old_shape)
return new_obj
# # Matrix filtering
[docs]
def filter_purify_threshold(self, threshold=0.01, keep=False):
"""Method that filters experimental errors by making zero a certain number of eigenvalues of the covariance matrix.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016) pp 226.
Parameters:
threshold (float): If eigenvalues are lower than threshold, they will be make 0. Default: 0.01.
keep (bool): If True, the original object won't be altered. Default: False.
Returns:
(Mueller): Result.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.parent.copy()
else:
new_obj = self.parent
old_shape = new_obj.shape
# Calculate covariance matrix eigenvalues
new_obj.covariance_matrix()
val, vect = new_obj.parameters.eig(values_as_matrix=True,
vectors_as_matrix=True,
shape=False)
# Make 0 the desired eigenvalues
val[val < threshold] = 0
# Recompose the matrix
H, diag, Ht = create_Mueller(N=3)
diag.from_components([val[0, :]] + 4 * [0] + [val[1, :]] + 4 * [0] +
[val[2, :]] + 4 * [0] + [val[3, :]])
H.from_matrix(vect)
Ht.from_matrix(np.conj(np.transpose(vect, axes=(1, 0, 2))))
result = H * diag * Ht
new_obj.from_covariance(result, shape=old_shape)
return new_obj
[docs]
def filter_physical_conditions(self,
tol=tol_default_filter,
ignore_cond=None,
keep=False,
_counter=0):
"""Method that filters experimental errors by forcing the Mueller matrix to fulfill the conditions necessary for a matrix to be a real optical component.
Parameters:
tol (float): Tolerance in equalities.
ignore_cond (list): Conditions to ignore. If False or None, no condition is ignored. Default: None.
keep (bool): If True, the object is updated to the filtered result. If false, a new fresh copy is created. Default: True.
_counter (int): Auxiliar variable that shouldn't be used when calling this function from outside.
Returns:
Mf (4x4 matrix): Filtered Mueller matrix.
"""
# Act differently if we want to keep self intact
if keep:
new_obj = self.parent.copy()
else:
new_obj = self.parent
old_shape = new_obj.shape
# Start by calculating if the object is physically realizable
cond, partial = new_obj.checks.is_physical(give_all=True,
tol=tol,
ignore_cond=ignore_cond,
shape=False,
out_number=False)
# Act only if there is some violations or max numebr of iter is reached
if np.any(~cond) and _counter <= counter_max:
# Start by calculating the required information
comp = new_obj.parameters.components(shape=False, out_number=False)
D = new_obj.parameters.diattenuation(shape=False, out_number=False)
P = new_obj.parameters.polarizance(shape=False, out_number=False)
# Now solve the conditions one by one. However, in order to avoid problems, once we have acted in one matrix of the object, we won't modify it again until next iterations, except if the conditions are complementary
used = np.zeros(new_obj.size, dtype=bool)
# Condition 1: Elements must be real
cond = ~partial[0]
if np.any(cond):
for ind, elem in enumerate(comp):
comp[ind] = np.array(elem, dtype=float)
used = used + ~partial[0]
# Condition 2a: M00 must be greater than 0
cond = ~partial[1]
if np.any(cond):
comp[0][cond] = 0
used = used + cond
# Condition 2b: M00 must be lower than 1
cond = ~partial[2]
if np.any(cond):
M00 = comp[0][cond]
for ind, elem in enumerate(comp):
comp[ind][cond] = comp[ind][cond] / M00
used = used + cond
# Condition 3: All elements must be lower than M00
cond = ~partial[3]
if np.any(cond):
M00 = comp[0]
for ind, elem in enumerate(comp):
cond2 = np.abs(comp[ind]) > M00
comp[ind][cond * cond2] = np.sign(
comp[ind][cond * cond2]) * M00[cond2]
used = used + cond
# Condition 4a: Diattenuation can't be greater than 1
cond = ~partial[4] * ~used
if np.any(cond):
comp[1][cond] = comp[1][cond] / D[cond]
comp[2][cond] = comp[2][cond] / D[cond]
comp[3][cond] = comp[3][cond] / D[cond]
used = used + cond
# Condition 4b: Polarizance can't be greater than 1
cond = ~partial[5] * ~used
if np.any(cond):
comp[4][cond] = comp[4][cond] / P[cond]
comp[8][cond] = comp[8][cond] / P[cond]
comp[12][cond] = comp[12][cond] / P[cond]
used = used + cond
# Condition 5a: Total transmission can't be greater than 1
cond = ~partial[6] * ~used
if np.any(cond):
aux = 1 / comp[0][cond] - 1
comp[1][cond] = comp[1][cond] * aux
comp[2][cond] = comp[2][cond] * aux
comp[3][cond] = comp[3][cond] * aux
used = used + cond
# Condition 5b: Total reciproc transmission can't be greater than 1
cond = ~partial[7] * ~used
if np.any(cond):
aux = 1 / comp[0][cond] - 1
comp[4][cond] = comp[4][cond] * aux
comp[8][cond] = comp[8][cond] * aux
comp[12][cond] = comp[12][cond] * aux
used = used + cond
# Condition 6: Tr(M*M^T) <= 4*M00. Leave away M00, P and D components, as they should be fixed now
cond = ~partial[8] * ~used
if np.any(cond):
# Leave away M00, P and D components
cte1 = 2 * (comp[1][cond] * comp[4][cond] + comp[2][cond] *
comp[8][cond] + comp[3][cond] * comp[12][cond])
cte2 = comp[5][cond]**2 + comp[10][cond]**2 + comp[15][
cond]**2 + 2 * (comp[6][cond] * comp[9][cond] +
comp[7][cond] * comp[13][cond] +
comp[11][cond] * comp[14][cond])
aux = np.sqrt((3 * comp[0][cond]**2 - cte1) / cte2)
for ind, elem in enumerate(comp):
if not ind in (0, 1, 2, 3, 4, 8, 12):
comp[ind][cond] = comp[ind][cond] * aux
used = used + cond
# Cond 7a: Condition in m matrix
cond = ~partial[9] * ~used
if np.any(cond):
M00 = comp[0][cond]
# Calculate the row weights and signs
k1 = comp[5][cond] + comp[6][cond] + comp[7][cond]
k2 = comp[9][cond] + comp[10][cond] + comp[11][cond]
k3 = comp[13][cond] + comp[14][cond] + comp[15][cond]
s1 = np.sign(1 - k1 / (M00 * D[cond]))
s2 = np.sign(1 - k2 / (M00 * D[cond]))
s3 = np.sign(1 - k3 / (M00 * D[cond]))
kT = np.abs(k1) + np.abs(k2) + np.abs(k3)
# Fix
cond2 = (k1 != 0) * (comp[1][cond] != 0)
if np.any(cond2):
aux = (M00[cond2] * D[cond][cond2] / k1[cond2]) * (
1 - s1[cond2] * M00[cond2] * (1 - D[cond][cond2]) /
(np.abs(comp[1][cond][cond2] * k1[cond2] / kT[cond2])))
comp[5][cond][cond2] *= aux
comp[6][cond][cond2] *= aux
comp[7][cond][cond2] *= aux
cond2 = (k2 != 0) * (comp[2][cond] != 0)
if np.any(cond2):
aux = (M00[cond2] * D[cond][cond2] / k2[cond2]) * (
1 - s2[cond2] * M00[cond2] * (1 - D[cond][cond2]) /
(np.abs(comp[2][cond][cond2] * k2[cond2] / kT[cond2])))
comp[9][cond][cond2] *= aux
comp[10][cond][cond2] *= aux
comp[11][cond][cond2] *= aux
cond2 = (k3 != 0) * (comp[3][cond] != 0)
if np.any(cond2):
aux = (M00[cond2] * D[cond][cond2] / k3[cond2]) * (
1 - s3[cond2] * M00[cond2] * (1 - D[cond][cond2]) /
(np.abs(comp[3][cond][cond2] * k3[cond2] / kT[cond2])))
comp[13][cond][cond2] *= aux
comp[14][cond][cond2] *= aux
comp[15][cond][cond2] *= aux
used = used + cond
# Cond 7b: Condition in reciprocal m matrix
cond = ~partial[10] * ~used
if np.any(cond):
M00 = comp[0][cond]
# Calculate the row weights and signs
k1 = comp[5][cond] + comp[9][cond] + comp[13][cond]
k2 = comp[6][cond] + comp[10][cond] + comp[14][cond]
k3 = comp[7][cond] + comp[11][cond] + comp[15][cond]
s1 = np.sign(1 - k1 / (M00 * P[cond]))
s2 = np.sign(1 - k2 / (M00 * P[cond]))
s3 = np.sign(1 - k3 / (M00 * P[cond]))
kT = np.abs(k1) + np.abs(k2) + np.abs(k3)
# Fix
cond2 = (k1 != 0) * (comp[4][cond] != 0)
if np.any(cond2):
aux = (M00[cond2] * D[cond][cond2] / k1[cond2]) * (
1 - s1[cond2] * M00[cond2] * (1 - D[cond][cond2]) /
(np.abs(comp[4][cond][cond2] * k1[cond2] / kT[cond2])))
comp[5][cond][cond2] *= aux
comp[9][cond][cond2] *= aux
comp[13][cond][cond2] *= aux
cond2 = (k2 != 0) * (comp[2][cond] != 0)
if np.any(cond2):
aux = (M00[cond2] * D[cond][cond2] / k2[cond2]) * (
1 - s2[cond2] * M00[cond2] * (1 - D[cond][cond2]) /
(np.abs(comp[7][cond][cond2] * k2[cond2] / kT[cond2])))
comp[6][cond][cond2] *= aux
comp[10][cond][cond2] *= aux
comp[14][cond][cond2] *= aux
cond2 = (k3 != 0) * (comp[3][cond] != 0)
if np.any(cond2):
aux = (M00[cond2] * D[cond][cond2] / k3[cond2]) * (
1 - s3[cond2] * M00[cond2] * (1 - D[cond][cond2]) /
(np.abs(comp[10][cond][cond2] * k3[cond2] / kT[cond2]))
)
comp[7][cond][cond2] *= aux
comp[11][cond][cond2] *= aux
comp[15][cond][cond2] *= aux
used = used + cond
# Cond 8: Conditions in the eigenvalues of the covariance matrix
cond = (~partial[11] + ~partial[12] + ~partial[13]) * ~used
if np.any(cond):
# Calculate the covariance matrix and its eigenvalues
obj = new_obj.covariance_matrix(keep=True)
val, vect = obj.parameters.eig(values_as_matrix=True,
vectors_as_matrix=True,
shape=False)
# Condition 8a: Eigenvalues must be real
cond2 = ~partial[11] * ~used
if np.any(cond2):
val = np.real(val)
used = used + cond2
# Condition 8a: Eigenvalues must be greater than 0
cond2 = ~partial[12] * ~used
if np.any(cond2):
val[val < 0] = 0
# Condition 8a: Eigenvalues must be lower than 1
cond2 = ~partial[13] * ~used
if np.any(cond2):
val[val > 1] = 1
# Now, lets merge the fix in the last condition with the rest
H, diag, Ht = create_Mueller(N=3)
diag.from_components([val[0, :]] + 4 * [0] + [val[1, :]] +
4 * [0] + [val[2, :]] + 4 * [0] +
[val[3, :]])
H.from_matrix(vect)
Ht.from_matrix(np.conj(np.transpose(vect, axes=(1, 0, 2))))
result = H * diag * Ht
result.from_covariance(result)
comp2 = result.parameters.components(shape=False,
out_number=False)
for ind, elem in enumerate(comp2):
comp[ind][cond] = elem[cond]
# Finally, use the corrected values
new_obj.from_components(comp, shape=old_shape)
# Call again this function to solve all the problems iteratively through recurrence
new_obj.analysis.filter_physical_conditions(
tol=tol,
ignore_cond=ignore_cond,
keep=keep,
_counter=_counter + 1)
return new_obj
####################################################################
# Matrix decomposition
####################################################################
[docs]
def decompose_pure(
self,
decomposition='RP',
give_all=False,
tol=tol_default,
# filter=False,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False,
transmissions='all',
angles="all"):
"""Polar decomposition of a pure Mueller matrix in a retarder and a diattenuator.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 151.
S. Y. Lu, R. A. Chipman; "Interpretation of Mueller matrices based on polar decomposition"; J. Opt. Soc. Am. A/Vol. 13, No. 5 (1996)
Parameters:
decomposition (string): string with the order of the elements: retarder (R) or diattenuator/polarizer (D or P). Default: RP.
give_all (bool): If True, the dictionary of parameters will be given in the returned. Default: False.
tol (float): Tolerance in equalities.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (float): If true, the function prints out some information about the matrices.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
transmissions (string): Determines the type of transmission output, FIELD, INTENSITY or ALL. Default: All.
angles (string): Determines the type of angles output, CHARAC (characteristic angles), AZIMUTH (azimuth and ellipticity) or ALL. Default: All.
Returns:
Mr (Mueller): Mueller matrix of the retarder.
Md (Mueller): Mueller matrix of the diattenuator.
param (dictionary, optional): Dictionary with all the parameters of the decomposed elements.
"""
# In order to be efficient, start by printing soem things
Md = self.parent.copy()
Md.name = 'Diattenuator of ' + self.parent.name
if verbose or draw:
if decomposition[0] in ('Rr'):
decomposition = 'Mr * Md'
else:
decomposition = 'Md * Mr'
print("\n------------------------------------------------------")
print('Pure decomposition of {} as M = {}.'.format(
self.parent.name, decomposition))
# Filter the matrix if required TODO
# First step, extract the diattenuator parameters
if decomposition[0] in ('Rr'):
trans, ang = Md.analysis.diattenuator(transmissions='all',
angles='all',
out_number=out_number,
shape=shape,
shape_like=shape_like,
verbose=verbose,
draw=draw)
else:
trans, ang = Md.analysis.polarizer(transmissions='all',
angles='all',
out_number=out_number,
shape=shape,
shape_like=shape_like,
verbose=verbose,
draw=draw)
p1, p2 = trans[2:]
azD, elD = ang[2:]
# Now, we can create the diatenuator object
Md.diattenuator_azimuth_ellipticity(p1=p1,
p2=p2,
azimuth=azD,
ellipticity=elD)
# Now, the retarder. To do that, we have to invert the diattenuator
Mdi = Mueller()
Mdi.diattenuator_azimuth_ellipticity(p1=1 / p1,
p2=1 / p2,
azimuth=azD,
ellipticity=elD)
if decomposition[0] in ('Rr'):
Mr = self.parent * Mdi
else:
Mr = Mdi * self.parent
# In some cases, the diattenuator is not invertible. Solve those cases
cond = p2 < tol
if np.any(cond):
# If p1 was also 0, the original matrix is a total absorbr, so a retarder does not ahve sense
cond2 = p1 < tol
if np.any(cond2):
Mr[cond2] = self.parent[cond2]
# For the rest of cases, the retarder is not unique. Let's chose the one with the lowest retardance.
cond2 = ~cond2 * cond
if np.any(cond2):
# First, the retardance
Dv = self.parent.parameters.diattenuation_vector(
shape=shape, shape_like=shape_like)
Pv = self.parent.parameters.polarizance_vector(
shape=shape, shape_like=shape_like)
R = np.arccos(np.sum(Dv * Pv, axis=0))
# Now, the retardance vector
Rv = np.cross(Pv, Dv, axis=0)
norm = np.linalg.norm(Rv, axis=0)
Rv[0, :] = Rv[0, :] / norm
Rv[1, :] = Rv[1, :] / norm
Rv[2, :] = Rv[2, :] / norm
# Finally, construct the object
Mr2 = Mueller()
Mr2.retarder_vector(Rv=Rv,
R=R,
shape=shape,
shape_like=shape_like)
Mr[cond2] = Mr2[cond2]
Mr.name = 'Retarder of ' + self.parent.name
# Set the global phases
if self.parent.global_phase is not None:
Md.set_global_phase(self.parent.global_phase / 2,
shape=shape,
shape_like=shape_like)
Mr.set_global_phase(self.parent.global_phase / 2,
shape=shape,
shape_like=shape_like)
# Information dictionary and visualization
if give_all or verbose or draw:
# Diattenuator
trans, ang = Md.analysis.diattenuator(transmissions=transmissions,
angles=angles,
out_number=out_number,
verbose=verbose,
draw=draw)
parameters = {}
if transmissions.upper() == 'FIELD':
parameters['p1'], parameters['p2'] = trans
elif transmissions.upper() == 'INTENSITY':
parameters['Tmax'], parameters['Tmin'] = trans
else:
parameters['Tmax'], parameters['Tmin'], parameters[
'p1'], parameters['p2'] = trans
if angles.upper() == 'CHARAC':
parameters['alpha D'], parameters['delay D'] = ang
elif angles.upper() == 'ALL':
parameters['alpha D'], parameters['delay D'], parameters[
'azimuth D'], parameters['ellipticity D'] = ang
else:
parameters['azimuth D'], parameters['ellipticity D'] = ang
# Retarder
R, ang = Mr.analysis.retarder(angles=angles,
out_number=out_number,
shape=shape,
shape_like=shape_like,
verbose=verbose,
draw=draw)
parameters['R'] = R
if angles.upper() == 'CHARAC':
parameters['alpha R'], parameters['delay R'] = ang
elif angles.upper() == 'ALL':
parameters['alpha R'], parameters['delay R'], parameters[
'azimuth R'], parameters['ellipticity R'] = ang
else:
parameters['azimuth R'], parameters['ellipticity R'] = ang
# Error
if decomposition[0] in ('Rr'):
Mrec = Mr * Md
else:
Mrec = Md * Mr
error = self.parent.parameters.matrix(
shape=shape, shape_like=shape_like) - Mrec.parameters.matrix(
shape=shape, shape_like=shape_like)
error = np.linalg.norm(error, axis=(0, 1)) / 16
if ~out_number and self.parent.size > 1:
error = np.array([error])
parameters['Error'] = error
if verbose or draw:
heading = '{} decomposition mean square error:'.format(
self.parent.name)
PrintParam(param=[error],
shape=self.parent.shape,
title=['MSE'],
heading=heading,
verbose=verbose,
draw=draw)
print(
"------------------------------------------------------\n")
# Return
if give_all:
return Mr, Md, parameters
else:
return Mr, Md
[docs]
def decompose_polar(self,
decomposition='PRD',
give_all=False,
tol=tol_default,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False,
transmissions='all',
angles="all",
depolarization='all'):
"""Polar decomposition of a physically realizable Mueller matrix in a partial depolarizer, retarder and a diattenuator.
TODO: When the depolarizer is singular with 2 or 3 non-zero eigenvalues, the decomposed retarder is often erroneous (60-80% of the time).
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 257.
Parameters:
decomposition (string): string with the order of the elements: depolarizer (P), retarder (R) or diattenuator (D). There are six possible decompositions: PRD, RDP, PDR, RDP, DRP and DPR. Default: PRD.
give_all (bool): If true, the dictionary of parameters will be given in the returned. Default: False.
tol (float): Tolerance in equalities.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (float): If true, the function prints out some information about the matrices.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
transmissions (string): Determines the type of transmission output, FIELD, INTENSITY or ALL. Default: All.
angles (string): Determines the type of angles output, CHARAC (characteristic angles), AZIMUTH (azimuth and ellipticity) or ALL. Default: All.
depolarization (string): Determines the type of depolarization information: INDEX, FACTORS or ALL. Default: All.
Returns:
Mr (Mueller): Mueller matrix of the retarder.
Md (Mueller): Mueller matrix of the diattenuator/polarizer.
Mp (Mueller): Mueller matrix of the depolarizer.
param (dictionary, optional): Dictionary with all the parameters of the decomposed elements.
"""
# Some common parameters
M = self.parent.copy()
pure = M.checks.is_pure(shape=shape, shape_like=shape_like)
# # If M is completelly pure, there is no point in continuing in this path, go to the pure decomposition instead
# if np.all(pure):
# Mp = Mueller('Depolarizer of ' + self.parent.name)
# # Pure case (no depolarizer).
# new_dec = decomposition.replace('P', '')
# new_dec.replace('p', '')
# Mr, Md = self.decompose_pure(decomposition=new_dec, tol=tol)
# Mp.vacuum(length=M.size, shape=Mr.shape)
# else:
# Create objects and extract usefull info
Mp, Md = create_Mueller(N=2)
M00 = M.parameters.mean_transmission(shape=shape,
shape_like=shape_like)
singular = M.checks.is_singular(shape=shape, shape_like=shape_like)
det = self.parent.parameters.det(shape=False)
# Chhoose decomposition
if decomposition.upper() == 'PRD':
# Start by calculating the diattenuator
Dv = M.parameters.diattenuation_vector(shape=shape,
shape_like=shape_like)
Md.diattenuator_vector(Dv=Dv, M00=M00)
# Now, calculate the inverse of the diattenuator and extract it from the total matrix
(p1, p2), (alpha,
delay) = Md.analysis.diattenuator(transmissions='field',
angles='charac')
Mdi = Mueller()
Mdi.diattenuator_charac_angles(p1=1 / p1,
p2=1 / p2,
alpha=alpha,
delay=delay)
Mdi.M[np.isnan(Mdi.M)] = 0 # Allow continuing if p2 = 0
Mdi.M[np.isinf(Mdi.M)] = 0 # Allow continuing if p2 = 0
Mf = M * Mdi
# Now, we have the polarizance vector of Md. Extract it
_, _, Pvf, mf = Mf.parameters.blocks()
Mf.from_blocks(m=mf)
# Calculate eigenvalues and eigenvectors of the square matrix
Mf2 = Mf * Mf.transpose(keep=True)
val, vect = Mf2.parameters.eig(values_as_matrix=True,
vectors_as_matrix=True,
shape=False)
# print(val, vect)
# Recompose the matrix of depolarizer using the square root of the eigenvalues
H, diag, Ht = create_Mueller(('R2', 'Diag', 'R2T'))
diag.from_components([np.sqrt(val[0, :])] + 4 * [0] +
[np.sqrt(val[1, :])] + 4 * [0] +
[np.sqrt(val[2, :])] + 4 * [0] +
[np.sign(det) * np.sqrt(val[3, :])])
H.from_matrix(vect)
Ht.from_matrix(np.conj(np.transpose(vect, axes=(1, 0, 2))))
result = H * diag * Ht
# print(H, diag, Ht)
# print('result1', result)
mp = result.parameters.small_matrix()
Mp.from_blocks(Pv=Pvf, m=mp)
# Finally, the retarder can be calculated by using the inverse matrix of Mp
diag.from_components([1 / np.sqrt(val[0, :])] + 4 * [0] +
[1 / np.sqrt(val[1, :])] + 4 * [0] +
[1 / np.sqrt(val[2, :])] + 4 * [0] +
[np.sign(det) / np.sqrt(val[3, :])])
result = H * diag * Ht
# print('result2', result)
Mr = result * Mf
# print('Mf', Mf, result, Mr)
# Now, let's see if we have the special cases of singular matrices.
# First case: Mp is singular and Md not
cond = singular * (p2 > tol)
if np.any(cond):
Mp_aux, Mr_aux = create_Mueller(N=2)
# We have to wnow how many eigenvalues of Mf are different from 0
N = np.sum(val > tol, axis=0)
# print('N = ', N)
if Mr.ndim > 1:
N = np.reshape(N, Mr.shape)
condN1 = cond * (N == 1)
condN2 = cond * (N == 2)
condN3 = cond * (N == 3)
condN23 = condN2 + condN3
# Only one eigenvalue, No retarder at all. This is very simple
if np.any(condN1):
Mp_aux.from_blocks(Pv=Pvf, shape=Mp.shape)
Mr_aux.vacuum(length=Mr.size, shape=Mr.size)
Mr[condN1] = Mr_aux[condN1]
Mp[condN1] = Mp_aux[condN1]
# We start to need the right eigenvectors of Mf^2
if np.any(condN23):
# Left eigenvectors. Make eigenvale 1 = 0 to remove it
val_aux = val
val_aux[val == 1] = 0
val_ord, vect_ord = order_eig(val_aux, vect, 'reverse')
v1 = vect_ord[1:, 0, :]
# Left eigenvectors. Make eigenvale 1 = 0 to remove it
Mf2_b = Mf.transpose(keep=True) * Mf
val_b, vect_b = Mf2_b.parameters.eig(
values_as_matrix=True,
vectors_as_matrix=True,
shape=False)
val_b[val_b == 1] = 0
val_b_ord, vect_b_ord = order_eig(val_b, vect_b, 'reverse')
w1 = vect_b_ord[1:, 0, :]
# Two eigenvalues.
if np.any(condN2):
# TODO: Works 30% of the time
# Calculate Mp
trace_m2 = Mf2.parameters.trace(
) - Mf2.parameters.mean_transmission()
result = Mf2 / np.sqrt(trace_m2)
Mp_aux.from_blocks(
Pv=Pvf,
m=result.parameters.small_matrix(normalize=False),
shape=Mp.shape)
Mp[condN2] = Mp_aux[condN2]
# Calculate retardance
trace_m = Mf.parameters.trace(
) - Mf.parameters.mean_transmission()
R = np.arccos(trace_m / np.sqrt(trace_m2))
# Calculate normalized retardance vector
Rv = np.cross(v1, w1, axis=0)
norm = np.linalg.norm(Rv, axis=0)
Rv[0, :] = Rv[0, :] / norm
Rv[1, :] = Rv[1, :] / norm
Rv[2, :] = Rv[2, :] / norm
# Calculate the retarder
Mr_aux.retarder_vector(R=R, Rv=Rv)
Mr[condN2] = Mr_aux[condN2]
# Three eigenvalues
if np.any(condN3):
# TODO: Works 15% of the time
# Calculate Mp
suma = np.sqrt(val_ord[0, :]) + np.sqrt(val_ord[1, :])
producto = np.sqrt(val_ord[0, :] * val_ord[1, :])
result = Mf2 + producto * result.vacuum(length=M.size)
result.inverse()
result = suma * result * Mf2
Mp_aux.from_blocks(
Pv=Pvf,
m=result.parameters.small_matrix(normalize=False),
shape=Mp.shape)
Mp[condN3] = Mp_aux[condN3]
# Calculate the second set of vectors
v2 = vect_ord[1:, 1, :]
w2 = vect_b_ord[1:, 1, :]
prod_v = np.cross(v1, v2, axis=0)
norm = np.linalg.norm(prod_v, axis=0)
prod_v[0, :] /= norm
prod_v[1, :] /= norm
prod_v[2, :] /= norm
prod_w = np.cross(w1, w2, axis=0)
norm = np.linalg.norm(prod_w, axis=0)
prod_w[0, :] /= norm
prod_w[1, :] /= norm
prod_w[2, :] /= norm
# Calculate Mr
mr = kron_axis(v1, w1, axis=0) + \
kron_axis(v2, w2, axis=0) + kron_axis(prod_v, prod_w, axis=0)
Mr_aux.from_blocks(m=mr)
Mr[condN3] = Mr_aux[condN3]
# Now, the diattenuator is singular
cond = p2 <= tol
if np.any(cond):
# In this case, both depolarizer and retarder are not unique. We choose a diagonal depolarizer (easy to understand) and the retarder with minimum retardance.
Mp_aux, Mr_aux = create_Mueller(N=2)
# Start by the depolarizer
P = M.parameters.polarizance(shape=shape,
shape_like=shape_like)
Mp_aux.depolarizer_diagonal(d=P)
Mp[cond] = Mp_aux[cond]
# Now, the retarder
Pv = M.parameters.polarizance_vector(shape=shape,
shape_like=shape_like)
Dv = M.parameters.diattenuation_vector(shape=shape,
shape_like=shape_like)
R = np.arccos(np.sum(Pv * Dv, axis=0) / P)
Rv = np.cross(Pv, Dv, axis=0)
norm = np.linalg.norm(Rv, axis=0)
Rv[0, :] /= norm
Rv[1, :] /= norm
Rv[2, :] /= norm
Mr_aux.retarder_vector(R=R, Rv=Rv)
Mr[cond] = Mr_aux[cond]
elif decomposition.upper() == 'PDR':
Mr, Md, Mp = M.analysis.decompose_polar(tol=tol,
shape_like=shape_like,
shape=shape)
Md = Mr * Md * Mr.transpose(keep=True)
elif decomposition.upper() == 'RPD':
Mr, Md, Mp = M.analysis.decompose_polar(tol=tol,
shape_like=shape_like,
shape=shape)
Mp = Mr.transpose(keep=True) * Mp * Mr
elif decomposition.upper() == 'DRP':
# Same as above for the transposed matrix
M2 = M.transpose(keep=True)
Mr, Md, Mp = M2.analysis.decompose_polar(tol=tol,
shape_like=shape_like,
shape=shape)
Mr.transpose()
Mp.transpose()
elif decomposition.upper() == 'RDP':
Mr, Md, Mp = M.analysis.decompose_polar(tol=tol,
shape_like=shape_like,
shape=shape,
decomposition='DRP')
Md = Mr.transpose(keep=True) * Md * Mr
elif decomposition.upper() == 'DPR':
Mr, Md, Mp = M.analysis.decompose_polar(tol=tol,
shape_like=shape_like,
shape=shape,
decomposition='DRP')
Mp = Mr * Mp * Mr.transpose(keep=True)
else:
raise ValueError(
'Decomposition {} is not valid, must be a combination of R, P and D without repetition'
.format(decomposition))
# Correct names
Mp.name = 'Depolarizer of ' + self.parent.name
Md.name = 'Diattenuator of ' + self.parent.name
Mr.name = 'Retarder of ' + self.parent.name
# Set glpbal phase
if M.global_phase is not None:
gp = M.parameters.global_phase(shape=shape, shape_like=shape_like)
Md.set_global_phase(gp / 3)
Mp.set_global_phase(gp / 3)
Mr.set_global_phase(gp / 3)
# Finally, print options and extract the info
# print the heading
if verbose or draw or give_all:
if verbose or draw:
decomposition = ''
for elem in decomposition:
if elem in 'Rr':
decomposition += 'Mr'
elif elem in 'Dd':
decomposition += 'Md'
else:
decomposition += 'Mp'
if elem != decomposition[-1]:
decomposition += ' * '
print(
"\n------------------------------------------------------")
print('Polar decomposition of {} as M = {}.'.format(
self.parent.name, decomposition))
# Parameters
parameters = {}
# Diattenuator
trans, ang = Md.analysis.diattenuator(transmissions=transmissions,
angles=angles,
out_number=out_number,
verbose=verbose,
draw=draw)
if transmissions.upper() == 'FIELD':
parameters['p1'], parameters['p2'] = trans
elif transmissions.upper() == 'INTENSITY':
parameters['Tmax'], parameters['Tmin'] = trans
else:
parameters['Tmax'], parameters['Tmin'], parameters[
'p1'], parameters['p2'] = trans
if angles.upper() == 'CHARAC':
parameters['alpha D'], parameters['delay D'] = ang
elif angles.upper() == 'ALL':
parameters['alpha D'], parameters['delay D'], parameters[
'azimuth D'], parameters['ellipticity D'] = ang
else:
parameters['azimuth D'], parameters['ellipticity D'] = ang
# Retarder
R, ang = Mr.analysis.retarder(angles=angles,
out_number=out_number,
shape=shape,
shape_like=shape_like,
verbose=verbose,
draw=draw)
parameters['R'] = R
if angles.upper() == 'CHARAC':
parameters['alpha R'], parameters['delay R'] = ang
elif angles.upper() == 'ALL':
parameters['alpha R'], parameters['delay R'], parameters[
'azimuth R'], parameters['ellipticity R'] = ang
else:
parameters['azimuth R'], parameters['ellipticity R'] = ang
# Depolarizer
trans_D, trans_P, angles_D, angles_P, depol, principal_states = Mp.analysis.depolarizer(
transmissions=transmissions,
angles=angles,
depolarization=depolarization,
out_number=out_number,
verbose=verbose,
draw=draw)
if transmissions.upper() == 'FIELD':
parameters['p1 Depol D'], parameters['p2 Depol D'] = trans_D
parameters['p1 Depol P'], parameters['p2 Depol P'] = trans_P
elif transmissions.upper() == 'INTENSITY':
parameters['Tmax Depol D'], parameters[
'Tmin Depol D'] = trans_D
parameters['Tmax Depol P'], parameters[
'Tmin Depol P'] = trans_P
else:
parameters['Tmax Depol D'], parameters[
'Tmin Depol D'], parameters['p1 Depol D'], parameters[
'p2 Depol D'] = trans_D
parameters['Tmax Depol P'], parameters[
'Tmin Depol P'], parameters['p1 Depol P'], parameters[
'p2 Depol P'] = trans_P
if angles.upper() == 'CHARAC':
parameters['alpha Depol D'], parameters[
'delay Depol D'] = angles_D
parameters['alpha Depol P'], parameters[
'delay Depol P'] = angles_P
elif angles.upper() == 'ALL':
parameters['alpha Depol D'], parameters[
'delay Depol D'], parameters[
'azimuth Depol D'], parameters[
'ellipticity Depol D'] = angles_D
parameters['alpha Depol P'], parameters[
'delay Depol P'], parameters[
'azimuth Depol P'], parameters[
'ellipticity Depol P'] = angles_P
else:
parameters['azimuth Depol D'], parameters[
'ellipticity Depol D'] = angles_D
parameters['azimuth Depol P'], parameters[
'ellipticity Depol P'] = angles_P
# Error
Mrec = Mueller()
Mrec.vacuum(length=M.size)
Mrec.shape = Md.shape
for elem in decomposition:
if elem.upper() == 'R':
Mrec = Mrec * Mr
elif elem.upper() == 'D':
Mrec = Mrec * Md
else:
Mrec = Mrec * Mp
error = M.parameters.matrix(
shape=shape, shape_like=shape_like) - Mrec.parameters.matrix(
shape=shape, shape_like=shape_like)
error = np.linalg.norm(error, axis=(0, 1)) / 16
if out_number and M.size == 1:
error = error[0]
parameters['Error'] = error
if verbose or draw:
heading = '{} decomposition mean square error:'.format(
self.parent.name)
PrintParam(param=[error],
shape=self.parent.shape,
title=['MSE'],
heading=heading,
verbose=verbose,
draw=draw)
print(
"------------------------------------------------------\n")
# Return
if give_all:
return Mr, Md, Mp, parameters
else:
return Mr, Md, Mp
# # Calculate the diattenuator/polarizer
# p1, p2, alphaP, delayP, fiP, chiP = self.diattenuator()
# Mp.diattenuator_charac_angles_from_vector(
# p1, p2, alphaP, delayP)
# D = Mp.parameters.diattenuation()
# # Sometimes, due to numeric calculation, D may be slightly higher than 1. Fix it.
# if D > 1 and D < 1 + eps:
# D = 1
# elif D > 1 + eps:
# raise ValueError(
# "Mueller matrix is not real, diattenuation is > 1.")
# # Check if the matrix M is singular or not.
# singM = parent.checks.is_singular(tol=tol)
# singMp = Mp.checks.is_singular(tol=tol)
# nz = 0
# if singMp:
# # We have to determine if only Md is singular or not
# P = parent.parameters.polarizance()
# cond3 = np.abs(1 - P) <= tol
# if cond3:
# # Print type of decomposition
# if verbose:
# print(
# "Both the depolarizer and the polarizer are singular."
# )
# # Homogeneous case
# Md.from_matrix(np.identity(4))
# Mr, Mp = decompose_pure(M, decomposition='PR', tol=tol)
# else:
# # Print type of decomposition
# if verbose:
# print("The polarizer is singular.")
# # Calculate the depolarizer polarizance vector
# Pdv = parent.P
# Mr.from_matrix(np.identity(4))
# cero = np.matrix(np.zeros(3))
# ceroM = np.zeros((3, 3))
# Md.from_blocks(cero, Pdv, ceroM)
# else:
# # Calculate the depolarizer polarizance vector
# # Dv, Pv, m, m00 = divide_in_blocks(M)
# # Pdv = (Pv - m * Dv.T) / (1 - D**2)
# # For calculating the small matrix m of the depolarizer we need an
# # auxiliary matrix mf
# Gaux = matrix(np.diag([1, -1, -1, -1]))
# if singM:
# Mpinv = Gaux * Mp * Gaux * D**(-2)
# else:
# Mpinv = Gaux * Mp * Gaux * (1 - D**2)**(-1)
# Mf = parent * Mpinv
# _, Pdv, mf, _ = divide_in_blocks(Mf.M)
# md2 = mf * mf.T
# qi2, mr2 = np.linalg.eigh(md2)
# # check_eig(qi2, mr2, md2)
# qi = np.sqrt(qi2)
# cero = np.matrix(np.zeros(3))
# # Calculation method depends on Md being singular or not
# if singM: # If M is singular and Mp is not => Md is singular
# # Calculate the number of eigenvalues that are zero and order them
# qi2, mr2 = order_eig(qi2, mr2)
# qi = np.sqrt(qi2)
# nz = sum(qi < tol)
# # Calculate other auxiliary matrices and vectors
# md1 = mf.T * mf
# qi12, mr1 = np.linalg.eigh(md1)
# qi12, mr1 = order_eig(qi12, mr1)
# v1, v2, w1, w2 = (mr2[:, 0], mr2[:, 1], mr1[:, 0],
# mr1[:, 1])
# if nz == 3:
# # Print type of decomposition
# if verbose:
# print(
# "Depolarized matrix singular case with three null eigenvalues."
# )
# # Trivial case
# md = np.zeros([3, 3])
# Md.from_blocks(cero, Pdv, md)
# Mr.from_matrix(np.eye(4))
# elif nz == 2:
# # Print type of decomposition
# if verbose:
# print(
# "Depolarized matrix singular case with two null eigenvalues."
# )
# # Depolarizer
# md = mf * mf.T / sqrt(np.trace(mf * mf.T))
# Md.from_blocks(cero, Pdv, md)
# # Retarder. Note that it is not unique.
# # TODO: inexacto
# cR = np.trace(mf) / sqrt(np.trace(mf * mf.T))
# R = np.arccos(cR)
# x1 = np.cross(v1.T, w1.T)
# Mr.retarder_from_vector(R,
# x1[0] / np.linalg.norm(x1))
# else:
# # Print type of decomposition
# if verbose:
# print(
# "Depolarized matrix singular case with one null eigenvalue."
# )
# # Depolarizer
# mat_aux = mf * mf.T + qi[0] * qi[1] * np.eye(3)
# md = (qi[0] + qi[1]) * (mat_aux).I * mf * mf.T
# Md.from_blocks(cero, Pdv, md)
# # Retarder. Note that it is not unique.
# # TODO: No funciona de momento
# # mr1 = np.matrix(mr1)
# # mr2 = np.matrix(mr2)
# # print([qi, qi2])
# # print(mf - mr2 * np.matrix(np.diag(qi)) * mr1.T)
# # print(mf.T * mf - mr1 * np.diag(qi2) * mr1.T)
# # print(mf * mf.T - mr2 * np.diag(qi2) * mr2.T)
# # print("")
# # print(mr1 * mr1.T)
# (y1, y2) = (np.cross(v1.T, v2.T),
# np.cross(w1.T, w2.T))
# mr = v1 * w1.T + v2 * w2.T + y1 * y2.T / (
# np.linalg.norm(y1) * np.linalg.norm(y2))
# # mr = np.matrix(mr)
# # print(mr * mr.T)
# Mr.from_blocks(cero, cero.T, mr)
# else:
# # Print type of decomposition
# if verbose:
# print("General case.")
# # General case
# s = np.sign(np.linalg.det(M))
# md = np.diag([qi[0], qi[1], s * qi[2]])
# md = mr2 * md * mr2.T
# Md.from_blocks(cero, Pdv, md)
# # Calculate the retarder
# mdinv = mr2 * np.diag(
# [1 / qi[0], 1 / qi[1], s / qi[2]]) * mr2.T
# mr = mdinv * mf
# Mr.from_blocks(cero, cero.T, mr)
# if decomposition in ('DPR', 'dpr'):
# Mpure = Mr * Mp
# Mr, Mp = Mpure.analysis.decompose_pure(
# decomposition='PR', tol=tol)
#
# elif decomposition in ('PRD', 'prd', 'RPD', 'rpd'):
# # This procedure is simple, we just have to traspose
# Mtras = Mueller(parent.name)
# Mtras.from_matrix(M.T)
# Md, Mr, Mp = Mtras.analysis.decompose_polar(tol=tol)
# Mr.from_matrix(Mr.M.T)
# Md.from_matrix(Md.M.T)
# # Reverse order of polarizer and retarder if required
# if decomposition in ('RPD', 'rpd'):
# Mpure = Mp * Mr
# Mr, Mp = Mpure.analysis.decompose_pure(
# decomposition='RP', tol=tol)
# elif decomposition in ('PDR', 'pdr'):
# # We can calculate this one slightly differently
# Mp, Mr, Md = parent.analysis.decompose_polar(
# decomposition='PRD', tol=tol)
# Md_matrix = Mr.M * Md.M * Mr.M.I
# Md.from_matrix(Md_matrix)
# elif decomposition in ('RDP', 'rdp'):
# # We can calculate this one as the traspose of the previous one
# Mtras = Mueller(parent.name)
# Mtras.from_matrix(M.T)
# Mp, Md, Mr = Mtras.analysis.decompose_polar(
# decomposition='PDR', tol=tol)
# Mr.from_matrix(Mr.M.T)
# Md.from_matrix(Md.M.T)
# else:
# raise ValueError("Decomposition not yet implemented.")
# # Order the output matrices
# Mout = [0, 0, 0]
# for ind in range(3):
# if decomposition[ind] == 'D':
# Mout[ind] = Md
# elif decomposition[ind] == 'P':
# Mout[ind] = Mp
# else:
# Mout[ind] = Mr
# # Calculate parameters
# if verbose or give_all:
# R, alphaR, delayR, fiR, chiR = Mr.analysis.retarder()
# Pd = Md.parameters.polarizance()
# Desp = Md.parameters.depolarization_index()
# # Calculate error
# if give_all or verbose:
# if decomposition == 'DRP':
# Mt = Md * Mr * Mp
# D = np.abs(Mt.M - M)
# MeanErr = np.linalg.norm(D)
# MaxErr = D.max()
# # Print results
# if verbose:
# if decomposition == 'DRP':
# print("Polar decomposition of the matrix M = Mdesp * Mr * Mp:")
# for ind in range(3):
# print("")
# if decomposition[ind] == 'D':
# print("The depolarizer Mueller matrix is:")
# print(Md)
# print("Parameters:")
# print((" - Polarizance = {}.".format(Pd)))
# print((" - Depolarization degree = {}.".format(Desp)))
# elif decomposition[ind] == 'P':
# print("The diatenuator/polarizer Mueller matrix is:")
# if singM and nz < 3:
# print(
# "Warning: Retarder matrix may be slightly inaccurate"
# )
# print(Mp)
# print("Parameters:")
# print((" - p1 = {}; p2 = {}.".format(p1, p2)))
# print((
# " - Angle = {} deg; Delay between components = {} deg."
# .format((alphaP / degrees), (delayP / degrees))))
# print(
# (" - Azimuth = {} deg; Ellipticity = {} deg.".format(
# (fiP / degrees), (chiP / degrees))))
# Tmax, Tmin = Mp.parameters.transmissions()
# print((
# " - Max. transmission = {} %; Min. transmission = {} %."
# .format((Tmax * 100), (Tmin * 100))))
# else:
# print("The retarder Mueller matrix is:")
# print(Mr)
# print("Parameters:")
# print((" - Delay = {} deg.".format((R / degrees))))
# print((
# " - Angle = {} deg; Delay between components = {} deg."
# .format((alphaR / degrees), (delayR / degrees))))
# print(
# (" - Azimuth = {} deg; Ellipticity = {} deg.".format(
# (fiR / degrees), (chiR / degrees))))
# print("")
# print(("The mean square error in the decomposition is: {}".format(
# MeanErr)))
# print(
# ("The maximum error in the decomposition is: {}".format(MaxErr)
# ))
# print("------------------------------------------------------")
# # Dictionary of parameters
# if give_all:
# param = dict(
# Delay=R,
# AngleR=alphaR,
# AxisDelayR=delayR,
# AzimuthR=chiR,
# EllipticityR=fiR,
# p1=p1,
# p2=p2,
# AngleP=alphaP,
# AxisDelayP=delayP,
# AzimuthP=fiP,
# EllipticityP=chiP,
# DespPolarizance=Pd,
# DespDegree=Desp,
# MeanError=MeanErr,
# MaxError=MaxErr)
# # Output
# if give_all:
# return Mout[0], Mout[1], Mout[2], param
# else:
# return Mout[0], Mout[1], Mout[2]
[docs]
class Check_Mueller(object):
"""Class for Check of Mueller Matrices
Parameters:
parent (Mueller_matrix): Mueller Matrix
Attributes:
self.parent (Mueller): parent object.
"""
def __init__(self, parent):
self.parent = parent
def __repr__(self):
"""Print all parameters."""
self.get_all(verbose=True, draw=True)
return ''
[docs]
def get_all(self, verbose=False, draw=False):
"""Creates a dictionary with all the checks of Mueller matrix.
Parameters:
verbose (bool): If True, print all parameters. Default: False.
draw (bool): If True, draw all plots/images of the parameters. Default: False.
"""
dict_params = {}
dict_params['Physical'] = self.is_physical(verbose=verbose,
draw=draw)
dict_params['Pure'] = self.is_pure(verbose=verbose, draw=draw)
dict_params['Homogenous'] = self.is_homogeneous(verbose=verbose,
draw=draw)
dict_params['Retarder'] = self.is_retarder(verbose=verbose,
draw=draw)
dict_params['Diattenuator'] = self.is_diattenuator(
verbose=verbose, draw=draw)
dict_params['Depolarizer'] = self.is_depolarizer(verbose=verbose,
draw=draw)
dict_params['Singular'] = self.is_singular(verbose=verbose,
draw=draw)
dict_params['Symmetric'] = self.is_symmetric(verbose=verbose,
draw=draw)
return dict_params
[docs]
def is_physical(self,
tol=tol_default,
give_all=False,
ignore_cond=None,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""A Mueller matrix must fulfill several conditions in order to be physically realizable:
1. $M_{ij}\in\mathbb{R}$ for i, j = 0, 1, 2 and 3.
2a. $M_{00} \geq 0$.
2b. $M_{00} \leq 1$ (except in active media).
3. $M_{00}\geq abs(M_{ij})$ for i, j = 0, 1, 2 and 3.
4a. $D \leq 1$
4a. $P \leq 1$.
5a. $M_{00} (1 + D) \leq 1$ (except in active media).
5b. $M_{00} (1 + P) \leq 1$ (except in active media).
6. $Tr(M*M^T)\leq 4(M_{00})^2$.
7a. $M_{00}^{2}\left(1-D\right)^{2}\geq\mathop{\sum_{i=1}^{3}M_{0i}^{2}\left(1-\sum_{j=1}^{3}\frac{M_{ij}}{M_{00}D}\right)}^{2}$.
7b. $M_{00}^{2}\left(1-P\right)^{2}\geq\mathop{\sum_{i=1}^{3}M_{i0}^{2}\left(1-\sum_{j=1}^{3}\frac{M_{ji}}{M_{00}P}\right)}^{2}$.
8a. $\lambda_{i}\in\mathbb{R}$
8b. $\lambda_{i} \geq 0$.
8c. $\lambda_{i} \leq 1$ (except in active media).
Being D the diattenuation, P the polarizance and $\lambda_{i}$ the eigenvalues of the covariance matrix.
References:
Handbook of Optics vol 2. 22.34 (There is an errata in the equation equivalent to condition 7).
"Polarized light and the Mueller Matrix approach", J. J. Gil, pp 187.
Parameters:
tol (float): Tolerance in conditions. Default: tol_default.
give_all (bool): If True, the method also gives a list with the individual conditions. Default: False.
ignore_cond (list): Conditions to ignore. If False or None, no condition is ignored. Default: None.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
cond (numpy.ndarray or bool): Result.
partial_conditions (list): List with the partial conditions.
"""
# Prepare to ignore conditions
if ignore_cond is None or ignore_cond is False:
use_cond = [True] * 8
else:
use_cond = [True] * 8
for ind in range(8):
if ind in ignore_cond:
use_cond[ind] = False
# Calculate the parameters
M00 = self.parent.parameters.mean_transmission(shape=shape,
shape_like=shape_like,
out_number=False)
D = self.parent.parameters.diattenuation(shape=shape,
shape_like=shape_like,
out_number=False)
P = self.parent.parameters.polarizance(shape=shape,
shape_like=shape_like,
out_number=False)
comp = self.parent.parameters.components(shape=shape,
shape_like=shape_like,
out_number=False)
H = self.parent.covariance_matrix(keep=True)
vals = H.parameters.eigenvalues(shape=shape,
shape_like=shape_like,
out_number=False)
# Calculate the conditions
if use_cond[0]:
cond_1 = np.isreal(comp[0])
for ind in range(1, 16):
cond_1 *= np.isreal(comp[ind])
else:
cond_1 = np.ones_like(D, dtype=bool)
if use_cond[1]:
cond_2a = M00 >= -tol
cond_2b = M00 <= 1 + tol
else:
cond_2a = np.ones_like(D, dtype=bool)
cond_2b = cond_2a
if use_cond[2]:
cond_3 = np.ones_like(D, dtype=bool)
for ind in range(1, 16):
cond_3 *= comp[0] >= np.abs(comp[ind])
else:
cond_3 = np.ones_like(D, dtype=bool)
if use_cond[3]:
cond_4a = D <= 1 + tol
cond_4b = P <= 1 + tol
else:
cond_4a = np.ones_like(D, dtype=bool)
cond_4b = cond_4a
if use_cond[4]:
cond_5a = M00 * (1 + D) <= 1 + tol
cond_5b = M00 * (1 + P) <= 1 + tol
else:
cond_5a = np.ones_like(D, dtype=bool)
cond_5b = cond_5a
if use_cond[5]:
sum = comp[5]**2 + comp[10]**2 + comp[15]**2 + 2 * (
comp[1] * comp[4] + comp[2] * comp[8] + comp[3] * comp[12] +
comp[6] * comp[9] + comp[7] * comp[13] + comp[11] * comp[14])
cond_6 = sum <= 3 * M00**2 + tol
else:
cond_6 = np.ones_like(D, dtype=bool)
if use_cond[6]:
aux1 = M00**2 + (1 - D)**2
aux2 = (comp[1]**2 * (1 - (comp[5] + comp[6] + comp[7]) /
(M00 * D))**2 + comp[2]**2 *
(1 - (comp[9] + comp[10] + comp[11]) / (M00 * D))**2 +
comp[3]**2 * (1 - (comp[13] + comp[14] + comp[15]) /
(M00 * D))**2)
cond_7a = aux1 >= aux2
aux1 = M00**2 + (1 - P)**2
aux2 = (comp[4]**2 * (1 - (comp[5] + comp[9] + comp[13]) /
(M00 * P))**2 + comp[8]**2 *
(1 - (comp[6] + comp[10] + comp[14]) / (M00 * P))**2 +
comp[12]**2 * (1 - (comp[7] + comp[11] + comp[15]) /
(M00 * P))**2)
cond_7b = aux1 >= aux2
else:
cond_7a = np.ones_like(D, dtype=bool)
cond_7b = cond_5a
if use_cond[7]:
cond_8a = np.ones_like(D, dtype=bool)
cond_8b = np.ones_like(D, dtype=bool)
cond_8c = np.ones_like(D, dtype=bool)
for elem in vals:
# print('elements', elem)
cond_8a *= np.abs(np.imag(elem)) < tol
# print('result', cond_8a)
cond_8b *= elem >= -tol
cond_8c *= elem <= 1 + tol
else:
cond_8a = np.ones_like(D, dtype=bool)
cond_8b = cond_8a
cond_8c = cond_8a
# Merge de conditions and multiply
cond_partial = [
cond_1, cond_2a, cond_2b, cond_3, cond_4a, cond_4b, cond_5a,
cond_5b, cond_6, cond_7a, cond_7b, cond_8a, cond_8b, cond_8c
]
cond = np.ones_like(D, dtype=bool)
for elem in cond_partial:
cond *= elem
# If the result is a number and the user asks for it, return a float
if out_number and elem.size == 1:
elem = elem[0]
if out_number and cond.size == 1:
cond = cond[0]
# Print the result if required
if verbose or draw:
if give_all:
cond_plot = [cond] + cond_partial
titles = [
'Physical', 'Real elements', 'M00 >= 0', 'M00 <= 1',
'abs(Mij) <= M00', 'D <= 1', 'P <= 1', 'Tmax <= 1',
'Recip Tmax <= 1', 'Tr(M*M^T) <= 4*M00^2', 'm cond (D)',
'm cond (P)', 'Real eigenvalues', 'Eigenvalues >= 0',
'Eigenvalues <= 1'
]
else:
cond_plot = [cond]
titles = ['Physical']
heading = '{} is physically realizable:'.format(self.parent.name)
PrintParam(param=cond_plot,
shape=self.parent.shape,
title=titles,
heading=heading,
verbose=verbose,
draw=draw)
# Return
if give_all:
return cond, cond_partial
else:
return cond
[docs]
def is_non_depolarizing(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if matrix is pure, i.e., is non-depolarizing (the degree of polarimetric purity must be 1).
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or bool): Result.
"""
# Calculate the parameters
PP = self.parent.parameters.polarimetric_purity(shape=shape,
shape_like=shape_like,
out_number=out_number)
cond = 1 - PP <= tol_default
# Print the result if required
if verbose or draw:
heading = '{} is pure (non-depolarizing):'.format(self.parent.name)
PrintParam(param=cond,
shape=self.parent.shape,
title='Pure',
heading=heading,
verbose=verbose,
draw=draw)
return cond
[docs]
def is_pure(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if matrix is pure, i.e., is non-depolarizing (the degree of polarimetric purity must be 1).
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or bool): Result.
"""
# Calculate the parameters
PP = self.parent.parameters.polarimetric_purity(shape=shape,
shape_like=shape_like,
out_number=out_number)
cond = 1 - PP <= tol_default
# Print the result if required
if verbose or draw:
heading = '{} is pure (non-depolarizing):'.format(self.parent.name)
PrintParam(param=cond,
shape=self.parent.shape,
title='Pure',
heading=heading,
verbose=verbose,
draw=draw)
return cond
[docs]
def is_homogeneous(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if the matrix is homogeneous, i.e., its two eigenstates are perpendicular. If true, the inhomogeneity parameter must be 0.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 119.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or bool): Result.
"""
# Calculate the parameters
par = self.parent.parameters.inhomogeneity(shape=shape,
shape_like=shape_like,
out_number=out_number)
cond = par <= tol_default
# Print the result if required
if verbose or draw:
heading = '{} is homogeneous:'.format(self.parent.name)
PrintParam(param=cond,
shape=self.parent.shape,
title='Homogeneous',
heading=heading,
verbose=verbose,
draw=draw)
return cond
[docs]
def is_retarder(self,
give_all=False,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if the matrix M corresponds to a pure retarder.There are three
conditions:
1. Is physical
2. Is pure
3. Diatteunation = 0.
4. Polarizance = 0.
5. M must be unitary ($M^{T}=M^{-1}$).
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 129.
Parameters:
give_all (bool): If True, the method also gives a list with the individual conditions. Default: False.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
cond (numpy.ndarray or bool): Result.
partial_conditions (list): List with the partial conditions (if give_all is True).
"""
# Reality and purity
cond1 = self.is_physical()
cond2 = self.is_pure()
# Diattenuation and polarizance
par = self.parent.parameters.diattenuation(shape=shape,
shape_like=shape_like,
out_number=out_number)
cond3 = par <= tol_default
par = self.parent.parameters.polarizance(shape=shape,
shape_like=shape_like,
out_number=out_number)
cond4 = par <= tol_default
# Unitary
obj_T = self.parent.transpose(keep=True)
obj_inv = self.parent.inverse(keep=True)
obj = obj_T - obj_inv
cond5 = norm_1a = obj_error(obj_T,
obj_inv,
shape=shape,
shape_like=shape_like,
out_number=out_number) < tol_default
# Reshape
cond1, cond2, cond3, cond4, cond5 = reshape([cond1, cond2, cond3, cond4, cond5],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
cond = cond1 * cond2 * cond3 * cond4 * cond5
# Print the result if required
if verbose or draw:
if give_all:
cond_plot = [cond, cond1, cond2, cond3, cond4, cond5]
titles = ['Retarder', "Physical", "Pure", 'D = 0', 'P = 0', 'Unitary']
else:
cond_plot = [cond]
titles = ['Retarder']
heading = '{} is a retarder:'.format(self.parent.name)
PrintParam(param=cond_plot,
shape=self.parent.shape,
title=titles,
heading=heading,
verbose=verbose,
draw=draw)
# Return
if give_all:
return cond, [cond1, cond2, cond3, cond4, cond5]
else:
return cond
[docs]
def is_diattenuator(self,
give_all=False,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if the matrix M corresponds to a pure homogeneous diattenuator. It must fullfill several conditions:
1. Is physical
2. Is pure
3. Diattenuation > 0.
4. $M = M^T$.
5. The eigenstates of M are the Stokes vectors (1, D) and (1, -D).
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 142.
Parameters:
give_all (bool): If True, the method also gives a list with the individual conditions. Default: False.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
cond (numpy.ndarray or bool): Result.
partial_conditions (list): List with the partial conditions (if give_all is True).
"""
# Reality and purity
cond1 = self.is_physical()
cond2 = self.is_pure()
# Calculate the first two conditions
par = self.parent.parameters.diattenuation(shape=False,
out_number=False)
cond3 = par > tol_default
Dv = self.parent.parameters.diattenuation_vector(shape=False)
obj_T = self.parent.transpose(keep=True)
cond4 = obj_error(self.parent,
obj_T,
shape=shape,
shape_like=shape_like,
out_number=out_number) < tol_default
# The third one is more complicated. Start by calculating the diattenuation vector
D = self.parent.parameters.diattenuation_vector(shape=False, as_Stokes=True)
# Check if it is eigenstate
condA = self.is_eigenstate(D)
D.M[1:,:] = -D.M[1:,:]
condB = self.is_eigenstate(D)
cond5 = condA * condB
# Reshape
cond1, cond2, cond3, cond4, cond5 = reshape([cond1, cond2, cond3, cond4, cond5],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Total condition
cond = cond1 * cond2 * cond3 * cond4 * cond5
# Print the result if required
if verbose or draw:
if give_all:
cond_plot = [cond, cond1, cond2, cond3, cond4, cond5]
titles = [
'Diattenuator', "Physical", "Pure", 'D > 0', 'Symetric', 'Correct eigenstates'
]
else:
cond_plot = [cond]
titles = ['Diattenuator']
heading = '{} is a diattenuator:'.format(self.parent.name)
PrintParam(param=cond_plot,
shape=self.parent.shape,
title=titles,
heading=heading,
verbose=verbose,
draw=draw)
# Return
if give_all:
return cond, [cond1, cond2, cond3, cond4, cond5]
else:
return cond
[docs]
def is_polarizer(self, *args, **kwargs):
"""Checks if the matrix M corresponds to a pure homogeneous polarizer (diattenuator). It must fullfill several conditions:
1. Is physical
2. Is pure
3. Polarizance > 0.
4. $M = M^T$.
5. The eigenstates of M are the Stokes vectors (1, P) and (1, -P).
Note: This method is the same as is_diattenuator.
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 142.
Parameters:
give_all (bool): If True, the method also gives a list with the individual conditions. Default: False.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
cond (numpy.ndarray or bool): Result.
partial_conditions (list): List with the partial conditions.
"""
return self.transpose(keep=True).is_diattenuator(*args, **kwargs)
[docs]
def is_depolarizer(self,
give_all=False,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if the matrix M corresponds to a depolarizer. It must fullfill several conditions:
1. Is physical
2. Depolarization index > 0.
3. $m = m^T$ (m being the small m matrix).
References:
S. Y. Lu, R. A. Chipman; "Interpretation of Mueller matrices based on polar decomposition"; J. Opt. Soc. Am. A/Vol. 13, No. 5 (1996)
Parameters:
give_all (bool): If True, the method also gives a list with the individual conditions. Default: False.
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
cond (numpy.ndarray or bool): Result.
partial_conditions (list): List with the partial conditions.
"""
# Reality
cond1 = self.is_physical()
# Calculate the parameters
inhom = self.parent.parameters.depolarization_index(
shape=shape, shape_like=shape_like, out_number=False)
m = self.parent.parameters.small_matrix(shape=False)
# First condition
cond2 = inhom > tol_default
# Second condition
mT = np.transpose(m, axes=(1, 0, 2))
cond3 = np.linalg.norm(m - mT, axis=(0, 1)) < tol_default
# Reshape
cond1, cond2, cond3 = reshape([cond1, cond2, cond3],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Total condition
cond = cond1 * cond2 * cond3
# Draw
if verbose or draw:
if give_all:
cond_plot = [cond, cond1, cond2, cond3]
titles = ['Depolarizer', "Physical", 'Depol. index > 0', 'm is symetric']
else:
cond_plot = [cond]
titles = ['Depolarizer']
heading = '{} is a depolarizer:'.format(self.parent.name)
PrintParam(param=cond_plot,
shape=self.parent.shape,
title=titles,
heading=heading,
verbose=verbose,
draw=draw)
# Return
if give_all:
return cond, [cond1, cond2, cond3]
else:
return cond
[docs]
def is_singular(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""Checks if the matrix is singular (det(M) = 0).
References:
J.J. Gil, R. Ossikovsky "Polarized light and the Mueller Matrix approach", CRC Press (2016), pp 282.
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or bool): Result.
"""
# Calculate the parameters
par = self.parent.parameters.det(shape=shape,
shape_like=shape_like,
out_number=out_number)
cond = par <= tol_default
# Print the result if required
if verbose or draw:
heading = '{} is singular:'.format(self.parent.name)
PrintParam(param=cond,
shape=self.parent.shape,
title='Singular',
heading=heading,
verbose=verbose,
draw=draw)
return cond
[docs]
def is_symmetric(self,
out_number=True,
shape_like=None,
shape=None,
verbose=False,
draw=False):
"""
Determines if the object matrix is symmetric (i.e. $$M = M^T$$).
Parameters:
out_number (bool): If True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): If True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or bool): Result.
"""
# Calculate the condition
obj_T = self.parent.transpose(keep=True)
cond = obj_error(self.parent,
obj_T,
shape=shape,
shape_like=shape_like,
out_number=out_number) < tol_default
# Reshape if neccessary
cond = reshape([cond],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = '{} is symmetric:'.format(self.parent.name)
PrintParam(param=cond,
shape=self.parent.shape,
title='Symmetric',
heading=heading,
verbose=verbose,
draw=draw)
return cond
[docs]
def is_eigenstate(self, S, out_number=True, shape_like=None, shape=None, verbose=False, draw=False):
"""
Determines if the vector S is an eigenstate of the object.
Parameters:
S (Stokes or Jones_vector): State to test.
out_number (bool): if True and the result is a 1x1 array, return a number instead. Default: True.
shape_like (numpy.ndarray or py_pol object): Use the shape of this object. Default: None.
shape (tuple or list): If no shape_like array is given, use this shape instead. Default: None.
verbose (bool): if True prints the parameter. Default: False.
draw (bool): If True and the object is a 1D or 2D, plot it. Default: False.
Returns:
(numpy.ndarray or bool): Result.
"""
# Multiply
Sout = self.parent * S
# Check if S and Sout are proportional
prop = Sout.M[0,:] / S.M[0,:]
cond1 = np.abs(Sout.M[1,:] / S.M[1,:] - prop) < tol_default
cond2 = np.abs(Sout.M[2,:] / S.M[2,:] - prop) < tol_default
cond3 = np.abs(Sout.M[3,:] / S.M[3,:] - prop) < tol_default
cond = cond1 * cond2 * cond3
# Reshape if neccessary
cond = reshape([cond],
shape_like=shape_like,
shape_fun=shape,
obj=self.parent)
# Print the result if required
if verbose or draw:
heading = '{} is an eigenstate of {}:'.format(S.name, self.parent.name)
PrintParam(param=cond,
shape=self.parent.shape,
title='Eigenstate',
heading=heading,
verbose=verbose,
draw=draw)
return cond