"""
Useful routines for managing polarization with the Jones calculus.
The routines are broken into four broad categories: creating Jones vectors,
creating Jones Matrices, interpreting Jones Vectors, and converting to the
Mueller / Stokes matrix calculus.
Creating Jones vectors for specific polarization states::
* field_linear(angle)
* field_left_circular()
* field_right_circular()
* field_horizontal()
* field_vertical()
* field_ellipsometry(tanpsi, Delta)
* field_elliptical(azimuth, elliptic_angle)
* field_components(Ex, Ey) [raw-component constructor]
Creating Jones Matrices for polarizing elements::
* op_linear_polarizer(angle)
* op_retarder(fast_axis_angle, retardance)
* op_attenuator(transmittance)
* op_mirror()
* op_rotation(angle)
* op_quarter_wave_plate(fast_axis_angle)
* op_half_wave_plate(fast_axis_angle)
* op_fresnel_reflection(index_of_refraction, angle)
* op_fresnel_transmission(index_of_refraction, angle)
Interpreting the polarization state::
* use_alternate_convention(boolean)
* interpret(jones_vector)
* intensity(jones_vector)
* phase(jones_vector)
* ellipse_orientation(jones_vector) [alias]
* ellipse_azimuth(jones_vector)
* ellipse_axes(jones_vector)
* ellipticity(jones_vector)
* ellipticity_angle(jones_vector)
* amplitude_ratio(jones_vector)
* amplitude_ratio_angle(jones_vector)
* normalize(jones_vector, preserve_global_phase=True)
* polarization_variable(jones_vector)
Converting to Mueller formalism::
* jones_op_to_mueller_op(jones_matrix)
* jones_to_stokes(jones_vector)
"""
import numpy as np
import pypolar.fresnel
__all__ = (
"use_alternate_convention",
"op_linear_polarizer",
"op_retarder",
"op_attenuator",
"op_mirror",
"op_rotation",
"op_quarter_wave_plate",
"op_half_wave_plate",
"op_fresnel_reflection",
"op_fresnel_transmission",
"field_linear",
"field_left_circular",
"field_right_circular",
"field_horizontal",
"field_vertical",
"field_ellipsometry",
"field_elliptical",
"field_components",
"interpret",
"intensity",
"phase",
"ellipse_orientation",
"ellipse_azimuth",
"ellipse_axes",
"ellipticity",
"ellipticity_angle",
"amplitude_ratio",
"amplitude_ratio_angle",
"normalize",
"polarization_variable",
"jones_op_to_mueller_op",
"jones_to_stokes",
)
alternate_sign_convention = False
[docs]
def use_alternate_convention(state):
"""
Change sign convention used for Jones calculus.
Read the documentation about the different conventions possible.
The default convention is to assume the wave function is
represented by exp(j*(omega * t - k * z)) and that the perspective
when viewing a sectional pattern is to look back along the
optical axis towards the source. This is the most commonly used
convention, but there are noteable exceptions (Wikipedia, Fowler, and
Hecht (sometimes).
Call this function once at the beginning and everything should
be just fine.
"""
global alternate_sign_convention
alternate_sign_convention = state
[docs]
def op_linear_polarizer(theta):
"""
Jones matrix operator for a rotated linear polarizer.
The polarizer has been rotated around a normal to its surface.
Args:
theta: rotation angle measured from the horizontal plane [radians]
"""
return np.array(
[[np.cos(theta) ** 2, np.sin(theta) * np.cos(theta)], [np.sin(theta) * np.cos(theta), np.sin(theta) ** 2]]
)
[docs]
def op_retarder(theta, delta):
"""
Jones matrix operator for an rotated optical retarder.
The retarder has been rotated around a normal to its surface.
Args:
theta: rotation angle between fast - axis and the horizontal plane [radians]
delta: phase delay introduced between fast and slow - axes [radians]
"""
if alternate_sign_convention:
theta *= -1
P = np.exp(+delta / 2 * 1j)
Q = np.exp(-delta / 2 * 1j)
D = np.sin(delta / 2) * 2j
C = np.cos(theta)
S = np.sin(theta)
retarder = np.array([[C * C * P + S * S * Q, C * S * D], [C * S * D, C * C * Q + S * S * P]])
if alternate_sign_convention:
return np.conjugate(retarder)
return retarder
[docs]
def op_attenuator(t):
"""
Jones matrix operator for an isotropic optical attenuator.
The transmittance t=I / I_0 is the fraction of light getting
through the attenuator or absorber.
Args:
t: fraction of intensity getting through attenuator [---]
"""
f = np.sqrt(t)
return np.array([[f, 0], [0, f]])
[docs]
def op_mirror():
"""Jones matrix operator for a perfect mirror."""
return np.array([[1, 0], [0, -1]])
[docs]
def op_rotation(theta):
"""
Jones matrix operator to rotate light around the optical axis.
Args:
theta : angle of rotation around optical axis [radians]
Returns:
2x2 matrix of the rotation operator [-]
"""
return np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]])
[docs]
def op_quarter_wave_plate(theta):
"""
Jones matrix operator for an rotated quarter - wave plate.
The QWP had been rotated around a normal to its surface.
Args:
theta : angle from fast - axis to horizontal plane [radians]
Returns:
2x2 matrix of the quarter - wave plate operator [-]
"""
return op_retarder(theta, np.pi / 2)
[docs]
def op_half_wave_plate(theta):
"""
Jones matrix operator for a rotated half - wave plate.
The half wave plate has been rotated around a normal to
the surface of the plate.
Args:
theta : angle from fast - axis to horizontal plane [radians]
Returns:
2x2 matrix of the half - wave plate operator [-]
"""
return op_retarder(theta, np.pi)
[docs]
def op_fresnel_reflection(m, theta, n_i=1):
"""
Jones matrix operator for Fresnel reflection at angle theta.
Args:
m : complex index of refraction [-]
theta : angle from normal to surface [radians]
n_i: real refractive index of incident medium [-]
Returns:
2x2 matrix of the Fresnel reflection operator [-]
"""
return np.array(
[
[pypolar.fresnel.r_par_amplitude(m, theta, n_i=n_i), 0],
[0, pypolar.fresnel.r_per_amplitude(m, theta, n_i=n_i)],
]
)
[docs]
def op_fresnel_transmission(m, theta, n_i=1):
"""
Jones matrix operator for Fresnel transmission at angle theta.
This operator acts on field amplitudes in the p/s basis and returns the
transmitted field amplitudes. It does not include irradiance
normalization factors.
Args:
m : complex index of refraction [-]
theta : angle from normal to surface [radians]
n_i: real refractive index of incident medium [-]
Returns:
2x2 Fresnel transmission operator [-]
"""
tper = pypolar.fresnel.t_per_amplitude(m, theta, n_i=n_i)
tpar = pypolar.fresnel.t_par_amplitude(m, theta, n_i=n_i)
return np.array([[tpar, 0], [0, tper]], dtype=complex)
[docs]
def field_linear(theta):
"""Jones vector for linear polarized light at angle theta from horizontal plane."""
return np.array([np.cos(theta), np.sin(theta)]).T
[docs]
def field_right_circular():
"""Jones Vector for right circular polarized light."""
J = 1 / np.sqrt(2) * np.array([1, 1j])
if alternate_sign_convention:
return np.conjugate(J)
return J
[docs]
def field_left_circular():
"""Jones Vector for left circular polarized light."""
J = 1 / np.sqrt(2) * np.array([1, -1j])
if alternate_sign_convention:
return np.conjugate(J)
return J
[docs]
def field_horizontal():
"""Jones Vector for horizontal polarized light."""
return np.array([1, 0])
[docs]
def field_vertical():
"""Jones Vector for vertical polarized light."""
return np.array([0, 1])
[docs]
def field_ellipsometry(tanpsi, Delta):
"""
Jones vector for using ellipometer parameters.
This creates a Jones vector for the specific set of ellipsometry
parameters tanpsi and Delta. See Fujiwara table 3.1 for example
Args:
tanpsi: abs(E_x / E_y) [-]
Delta: angle(E_x) - angle(E_y) [radians]
Returns:
Jones vector with specified characteristics
"""
psi = np.arctan(tanpsi)
J = np.array([np.sin(psi) * np.exp(1j * Delta), np.cos(psi)]).T
if alternate_sign_convention:
return np.conjugate(J)
return J
[docs]
def field_elliptical(azimuth, elliptic_angle, phi_x=0, E_0=1):
"""
Jones vector for elliptically polarized light.
Uses Azzam's equation 1.75
Args:
azimuth: tilt angle of ellipse from x - axis [radians]
elliptic_angle: arctan(minor - axis / major - axis) [radians]
phi_x: phase for E field in x - direction [radians]
E_0: amplitude of field
Returns:
Jones vector with specified characteristics
"""
ce = np.cos(elliptic_angle)
se = np.sin(elliptic_angle)
ca = np.cos(azimuth)
sa = np.sin(azimuth)
J = E_0 * np.array([ca * ce - sa * se * 1j, sa * ce + ca * se * 1j])
J *= np.exp(1j * (phi_x - np.angle(J[0])))
if alternate_sign_convention:
return np.conjugate(J)
return J
[docs]
def field_components(Ex, Ey):
"""
Build a Jones vector directly from x/y complex field components.
Args:
Ex: Complex electric-field component along x.
Ey: Complex electric-field component along y.
Returns:
Jones vector `[Ex, Ey]`, conjugated when the alternate sign convention
is active.
"""
J = np.array([Ex, Ey])
if alternate_sign_convention:
return np.conjugate(J)
return J
[docs]
def interpret(J):
"""
Interpret a Jones vector.
Args:
J: Jones vector with two complex entries.
Examples:
interpret([1, -1j]) -> "Right circular polarization"
interpret([0.5, 0.5]) -> "Linear polarization at 45.000000 degrees CCW from x-axis"
interpret(np.array([exp(-1j * pi), exp(-1j * pi / 3)]))
-> "Left elliptical polarization, rotated with respect to the axes"
Returns:
Human-readable interpretation string. Malformed or unphysical inputs are
reported with diagnostic strings.
"""
try:
arr = np.asarray(J, dtype=complex)
except (TypeError, ValueError):
message = "Malformed input: Jones vector must contain two numeric elements"
return message
if arr.size != 2:
message = "Malformed input: Jones vector must have exactly two elements"
return message
JJ = np.reshape(arr, 2)
if not np.all(np.isfinite(JJ.real)) or not np.all(np.isfinite(JJ.imag)):
message = "Malformed input: Jones vector contains NaN or infinite values"
return message
inten = intensity(JJ)
if inten <= 1e-12:
message = "Unphysical Jones vector: zero intensity (polarization state is undefined)"
return message
# Normalize interpretation to a single internal convention. Jones vectors
# created by this module flip handedness when `alternate_sign_convention`
# changes, so conjugate here in the default mode to keep naming consistent.
if alternate_sign_convention:
JI = JJ
else:
JI = np.conjugate(JJ)
mag1, mag2 = abs(JI[0]), abs(JI[1])
eps = 1e-9
delta = phase(JI)
phaze = np.degrees(delta)
azi = np.degrees(ellipse_azimuth(JI))
ell = ellipticity(JI)
ell_angle_rad = ellipticity_angle(JI)
ell_angle = np.degrees(ell_angle_rad)
def add_ellipticity_metrics(text):
"""Append ellipticity metrics with correct units."""
text += " ellipticity angle = %.1f°\n" % ell_angle
text += " ellipticity (b/a) = %.3f" % ell
return text
s = "Intensity is %.3f\n" % inten
s += "Phase is %.1f°\n" % phaze
if abs(ell_angle_rad) < eps:
ang = np.arctan2(mag2, mag1) * 180 / np.pi
return s + "Linear polarization at %f degrees CCW from x - axis" % ang
if abs(abs(ell_angle_rad) - np.pi / 4) < eps:
if ell_angle_rad < 0:
return s + "Right circular polarization"
return s + "Left circular polarization"
if ell_angle_rad < 0:
hand = "Right"
non_rotated = abs(delta + np.pi / 2) < eps
else:
hand = "Left"
non_rotated = abs(delta - np.pi / 2) < eps
if non_rotated:
s += "%s elliptical polarization, non - rotated\n" % hand
else:
s += "%s elliptical polarization\n" % hand
s += " rotated %.1f° respect to the axes\n" % azi
return add_ellipticity_metrics(s)
[docs]
def normalize(J, preserve_global_phase=True):
"""
Normalize a Jones vector by dividing by its Euclidean norm.
Args:
J: Jones vector.
preserve_global_phase: If `True` (default), preserve the global phase.
If `False`, remove global phase by making the first non-zero
component real and non-negative.
Returns:
Normalized Jones vector. A zero vector is returned unchanged.
"""
arr = np.asarray(J, dtype=complex)
if arr.size != 2:
raise ValueError("Jones vector must have exactly two elements")
vec = np.reshape(arr, 2)
norm = np.linalg.norm(vec)
if norm == 0:
return vec
out = vec / norm
if preserve_global_phase:
return out
eps = 1e-15
if abs(out[0]) > eps:
phi = np.angle(out[0])
else:
phi = np.angle(out[1])
return out * np.exp(-1j * phi)
[docs]
def intensity(J):
"""Return the intensity."""
inten = abs(J[..., 0]) ** 2 + abs(J[..., 1]) ** 2
return inten
[docs]
def phase(J):
"""Return the phase."""
gamma = np.angle(J[..., 1]) - np.angle(J[..., 0])
# Keep the phase difference on the principal branch [-pi, pi]
# so equivalent Jones vectors do not differ by 2*pi multiples.
return np.angle(np.exp(1j * gamma))
def _ellipticity_handedness(J):
"""Return handedness sign term proportional to the Stokes S3 component."""
return np.imag(np.conjugate(J[..., 0]) * J[..., 1])
[docs]
def ellipse_azimuth(J):
"""
Return the angle between the major semi - axis and the x - axis.
The polarization ellipse is rotated by this angle (called
the azimuth) relative to the laboratory frame.
"""
Ex0, Ey0 = np.abs(J)
delta = phase(J)
numer = 2 * Ex0 * Ey0 * np.cos(delta)
denom = Ex0**2 - Ey0**2
psi = 0.5 * np.arctan2(numer, denom)
return psi
[docs]
def ellipse_orientation(J):
"""
Return the orientation angle of the polarization ellipse.
This is a naming-parity alias for :func:`ellipse_azimuth`.
"""
return ellipse_azimuth(J)
[docs]
def ellipse_axes(J):
"""
Return the semi - major and semi - minor radii of the ellipse.
Twice these values will be the semi - major or semi - minor diameters.
"""
Ex0, Ey0 = np.abs(J)
alpha = ellipse_azimuth(J)
delta = phase(J)
C = np.cos(alpha)
S = np.sin(alpha)
asqr = (Ex0 * C) ** 2 + (Ey0 * S) ** 2 + 2 * Ex0 * Ey0 * C * S * np.cos(delta)
bsqr = (Ex0 * S) ** 2 + (Ey0 * C) ** 2 - 2 * Ex0 * Ey0 * C * S * np.cos(delta)
a = np.sqrt(abs(asqr))
b = np.sqrt(abs(bsqr))
if a < b:
return b, a
return a, b
[docs]
def ellipticity(J):
"""
Return the ellipticity of the polarization ellipse.
This is the ratio of semi - minor to semi - major radii.
The ellipticity is a measure of the fatness of the ellipse.
The ellipticity can be defined to always be positive. However
negative values can be used to indicate left - handed polarization.
Thus the ellipticity will range from -1 to 0 to 1 as light moves from
LCP to Linear Polarization to RCP.
"""
a, b = ellipse_axes(J)
if _ellipticity_handedness(J) < 0:
return -b / a
return b / a
[docs]
def ellipticity_angle(J):
"""
Return the ellipticity angle of the polarization ellipse.
The tangent of this angle is the ratio of semi - minor:semi - major
radii. It is between -pi / 4 ≤ angle ≤ pi / 4. Positive values
are for right–handed ellipticity. Negative values for left - handed
ellipticity.
"""
a, b = ellipse_axes(J)
if abs(a) >= abs(b):
epsilon = np.arctan2(b, a)
else:
epsilon = np.arctan2(a, b)
if _ellipticity_handedness(J) < 0:
return -epsilon
return epsilon
[docs]
def amplitude_ratio(J):
"""
Return the ratio of electric fields.
This is the amplitude of the vibrations along x measured
relative to the amplitude along y.
"""
Ex0, Ey0 = np.abs(J)
if Ex0 == 0:
return np.inf
return Ey0 / Ex0
[docs]
def amplitude_ratio_angle(J):
"""
Return the ratio of electric fields.
The tangent of this angle is the ratio of electric fields in
the y and x directions.
"""
Ex0, Ey0 = np.abs(J)
psi = np.arctan2(Ey0, Ex0)
return psi
[docs]
def polarization_variable(J):
"""
Return the complex polarization variable, chi.
This reduces the Jones vector to a single complex number and is useful
when the amplitude and absolute - phase are of secondary interest. These
are eliminated and chi is representative of the polarization state
in the complex plane.
"""
return J[1] / J[0]
def poincare_point(J):
"""Return the point on the Poincaré sphere."""
longitude = 2 * ellipse_azimuth(J)
a, b = ellipse_axes(J)
latitude = 2 * np.arctan2(b, a)
return latitude, longitude
[docs]
def jones_op_to_mueller_op(JJ):
"""
Convert a complex 2x2 Jones matrix to a real 4x4 Mueller matrix.
Hauge, Muller, and Smith, "Conventions and Formulas for Using the Mueller-
Stokes Calculus in Ellipsometry," Surface Science, 96, 81 - 107 (1980)
Args:
JJ: Jones matrix
Returns:
equivalent 4x4 Mueller matrix
"""
if alternate_sign_convention:
J = np.conjugate(JJ)
else:
J = JJ
M = np.zeros(shape=[4, 4], dtype=complex)
C = np.conjugate(J)
M[0, 0] = J[0, 0] * C[0, 0] + J[0, 1] * C[0, 1] + J[1, 0] * C[1, 0] + J[1, 1] * C[1, 1]
M[0, 1] = J[0, 0] * C[0, 0] + J[1, 0] * C[1, 0] - J[0, 1] * C[0, 1] - J[1, 1] * C[1, 1]
M[0, 2] = J[0, 1] * C[0, 0] + J[1, 1] * C[1, 0] + J[0, 0] * C[0, 1] + J[1, 0] * C[1, 1]
M[0, 3] = 1j * (J[0, 1] * C[0, 0] + J[1, 1] * C[1, 0] - J[0, 0] * C[0, 1] - J[1, 0] * C[1, 1])
M[1, 0] = J[0, 0] * C[0, 0] + J[0, 1] * C[0, 1] - J[1, 0] * C[1, 0] - J[1, 1] * C[1, 1]
M[1, 1] = J[0, 0] * C[0, 0] - J[1, 0] * C[1, 0] - J[0, 1] * C[0, 1] + J[1, 1] * C[1, 1]
M[1, 2] = J[0, 0] * C[0, 1] + J[0, 1] * C[0, 0] - J[1, 0] * C[1, 1] - J[1, 1] * C[1, 0]
M[1, 3] = 1j * (J[0, 1] * C[0, 0] + J[1, 0] * C[1, 1] - J[1, 1] * C[1, 0] - J[0, 0] * C[0, 1])
M[2, 0] = J[0, 0] * C[1, 0] + J[1, 0] * C[0, 0] + J[0, 1] * C[1, 1] + J[1, 1] * C[0, 1]
M[2, 1] = J[0, 0] * C[1, 0] + J[1, 0] * C[0, 0] - J[0, 1] * C[1, 1] - J[1, 1] * C[0, 1]
M[2, 2] = J[0, 0] * C[1, 1] + J[0, 1] * C[1, 0] + J[1, 0] * C[0, 1] + J[1, 1] * C[0, 0]
M[2, 3] = 1j * (-J[0, 0] * C[1, 1] + J[0, 1] * C[1, 0] - J[1, 0] * C[0, 1] + J[1, 1] * C[0, 0])
M[3, 0] = 1j * (J[0, 0] * C[1, 0] + J[0, 1] * C[1, 1] - J[1, 0] * C[0, 0] - J[1, 1] * C[0, 1])
M[3, 1] = 1j * (J[0, 0] * C[1, 0] - J[0, 1] * C[1, 1] - J[1, 0] * C[0, 0] + J[1, 1] * C[0, 1])
M[3, 2] = 1j * (J[0, 0] * C[1, 1] + J[0, 1] * C[1, 0] - J[1, 0] * C[0, 1] - J[1, 1] * C[0, 0])
M[3, 3] = J[0, 0] * C[1, 1] - J[0, 1] * C[1, 0] - J[1, 0] * C[0, 1] + J[1, 1] * C[0, 0]
MM = M.real / 2
return MM
def _jones_to_stokes(J):
"""
Convert Jones vector to Stokes vector.
Args:
J: Jones vector
Returns:
Stokes vector
"""
Ex = abs(J[0])
Ey = abs(J[1])
phi = phase(J)
S0 = Ex**2 + Ey**2
S1 = Ex**2 - Ey**2
S2 = 2 * Ex * Ey * np.cos(phi)
S3 = 2 * Ex * Ey * np.sin(phi)
return np.array([S0, S1, S2, S3])
[docs]
def jones_to_stokes(J):
"""
Convert Jones vector to Stokes vector.
Args:
J: Jones vector
Returns:
Stokes vector
"""
if J.ndim == 1:
return _jones_to_stokes(J)
n, m = J.shape
if m != 2:
print("Wrong shape ... should be %dx2 not %dx%d" % (m, n, m))
return None
S = np.empty(shape=(n, 4), dtype=np.ndarray)
for i, JJ in enumerate(J):
S[i] = _jones_to_stokes(JJ)
return S