# Encoding: utf-8
from dataclasses import dataclass
import chardet
import numpy as np
import numpy.typing as npt
import pandas as pd
import scipy.constants as sc
from numpy.lib.scimath import sqrt
from scipy.linalg import expm as scipy_expm
[docs]
@dataclass(frozen=True)
class DeltaRange:
lower: int
upper: int
def __post_init__(self):
if not (isinstance(self.lower, int) or isinstance(self.upper, int)):
raise TypeError("delta range values must be an integer")
if not -360 <= self.lower <= 0:
raise ValueError(
f"Lower bound ({self.lower}) is not in the valid range (-360, 0)."
)
if not 0 <= self.upper <= 360:
raise ValueError(
f"Upper bound ({self.upper}) is not in the valid range (0, 360)."
)
if not (
((self.lower, self.upper) == (0, 180))
or abs(self.upper - self.lower) == 360
):
raise ValueError(f"Invalid delta range ({self.lower}, {self.upper})")
[docs]
def convert_delta_range(
delta_values: npt.ArrayLike, lower: int, upper: int
) -> npt.ArrayLike:
"""Shifts the delta values to a given range.
Args:
delta_values (npt.ArrayLike): Delta values to be converted.
lower (int): The lower delta range. Should be in the range of -360 and 0.
upper (int): The upper delta range. Should be in the range of 0 and 360.
Returns:
npt.ArrayLike: Delta values with shifted range.
"""
delta_range = DeltaRange(lower, upper)
return np.mod((delta_values + delta_range.lower), 360) + delta_range.lower
[docs]
def calc_pseudo_diel(rho, angle: float, output: str = "eps") -> pd.DataFrame:
"""Calculates the pseudo dielectric function of a measurement from rho.
Args:
rho (pandas.DataFrame):
Measurement DataFrame containing rho as complex number as column and wavelength as index
angle (float): Angle of measurement in degree
output (str, optional): Output format for dielectric function.
'n': refractive index,
'eps': Dielectric function as two-column pandas.DataFrame,
'epsi': Dielectric function as imaginary number.
Defaults to 'eps'.
Returns:
pandas.DataFrame: Frame containing the pseudo dielectric function or refractive index.
"""
theta = angle * np.pi / 180
eps = np.sin(theta) ** 2 * (1 + np.tan(theta) ** 2 * ((1 - rho) / (1 + rho)) ** 2)
if output == "n":
n = sqrt(eps)
return pd.DataFrame({"n": n.real, "k": n.imag}, index=eps.index)
if output == "epsi":
return eps
return pd.concat(
{"ϵ1": eps.apply(lambda x: x.real), "ϵ2": eps.apply(lambda x: x.imag)}, axis=1
)
[docs]
def calc_rho(psi_delta: pd.DataFrame) -> pd.DataFrame:
"""Calculate rho from a Psi-Delta DataFrame.
The Psi-Delta DataFrame should be structured as follows:
:index: Wavelength
:column 'Ψ': Psi from measurement
:column 'Δ': Delta from measurement
This format is as returned from SpectraRay.read_psi_delta_file(...).
Args:
psi_delta (pandas.DataFrame): DataFrame containing Psi+Delta Measurement data
Returns:
pandas.DataFrame: Frame containing rho as an imaginary number.
"""
return psi_delta.apply(
lambda x: np.tan(np.deg2rad(x["Ψ"])) * np.exp(-1j * np.deg2rad(x["Δ"])), axis=1
)
[docs]
def convert_psi_delta_to_isotropic_mueller_matrix(psi_delta_df) -> pd.DataFrame:
r"""Extract aois and wavelengths values from psi_delta pandas.DataFrame and convert it to a Muellermatrix;
only works for isotropic media.
Args:
psi_delta_df (pd.DataFrame): dataFrame returned from data importers'
Returns:
pd.DataFrame: pyElli-compatible DataFrame of the Mueller matrix
"""
aois = psi_delta_df.index.get_level_values("Angle of Incidence").unique().to_numpy()
wavelengths = psi_delta_df.index.get_level_values("Wavelength").unique().to_numpy()
# create empty Mueller matrix
MM = pd.DataFrame(
{},
columns=[
"M11",
"M12",
"M13",
"M14",
"M21",
"M22",
"M23",
"M24",
"M31",
"M32",
"M33",
"M34",
"M41",
"M42",
"M43",
"M44",
],
index=pd.MultiIndex.from_product(
[aois, wavelengths], names=["Angle of Incidence", "Wavelength"]
),
dtype=float,
)
# fill in the Mueller matrix coefficients for an isotropic material based equation in https://www.jawoollam.com/resources/ellipsometry-faq#toggle-id-15
for aoi in aois:
for λ in wavelengths:
Δ = psi_delta_df.loc[aoi, λ]["Δ"]
Ψ = psi_delta_df.loc[aoi, λ]["Ψ"]
N = np.cos(2 * Ψ)
C = np.sin(2 * Ψ) * np.cos(2 * Δ)
S = np.sin(2 * Ψ) * np.cos(2 * Δ)
MM.loc[(aoi, λ), "M11"] = 1
MM.loc[(aoi, λ), "M12"] = -N
MM.loc[(aoi, λ), "M13"] = 0
MM.loc[(aoi, λ), "M14"] = 0
MM.loc[(aoi, λ), "M21"] = -N
MM.loc[(aoi, λ), "M22"] = 1
MM.loc[(aoi, λ), "M23"] = 0
MM.loc[(aoi, λ), "M24"] = 0
MM.loc[(aoi, λ), "M31"] = 0
MM.loc[(aoi, λ), "M32"] = 0
MM.loc[(aoi, λ), "M33"] = C
MM.loc[(aoi, λ), "M34"] = S
MM.loc[(aoi, λ), "M41"] = 0
MM.loc[(aoi, λ), "M42"] = 0
MM.loc[(aoi, λ), "M43"] = -S
MM.loc[(aoi, λ), "M44"] = C
return MM
[docs]
def get_qwp_thickness(material: "Material", lbda: float) -> float:
"""Return the thickness of a material in nm for a quarter wave plate at wavelength 'lbda'.
Args:
material (Material): Material object of the quarter wave plate
lbda (float): Wavelength (in nm) at which the quarter wave plate is calculated
Returns:
float: Thickness (in nm) of quarter wave plate
"""
nr = np.real(material.get_refractive_index(lbda)[0, 0, 0])
return lbda / (4.0 * nr)
# base vectors
E_X = np.array([1, 0, 0]).reshape((3,))
E_Y = np.array([0, 1, 0]).reshape((3,))
E_Z = np.array([0, 0, 1]).reshape((3,))
[docs]
def conversion_wavelength_energy(value: npt.ArrayLike) -> npt.ArrayLike:
r"""Converts wavelength values to energy values and vice versa.
.. math::
value_{\text{target}} = c \cdot \hbar / \boldsymbol{value}
Args:
value (npt.ArrayLike): Single value or array of wavelengths in nm or energy in eV.
Returns:
npt.ArrayLike: Energy in eV or wavelength in nm.
"""
return sc.speed_of_light * sc.value("Planck constant in eV/Hz") / (value * 1e-9)
[docs]
def conversion_frequency2energy(f: npt.ArrayLike) -> npt.ArrayLike:
r"""Converts frequency values to energy values.
.. math::
E = \boldsymbol{f} \cdot \hbar
Args:
f (npt.ArrayLike): Single value or array of frequencies in Hz.
Returns:
npt.ArrayLike: Energy in eV.
"""
return f * sc.value("Planck constant in eV/Hz")
[docs]
def conversion_energy2frequency(E: npt.ArrayLike) -> npt.ArrayLike:
r"""Converts energy values to frequency values.
.. math::
f = \boldsymbol{E} / \hbar
Args:
E (npt.ArrayLike): Single value or array of energies in eV.
Returns:
npt.ArrayLike: Frequency in Hz.
"""
return E / sc.value("Planck constant in eV/Hz")
[docs]
def conversion_wavelength_frequency(value: npt.ArrayLike) -> npt.ArrayLike:
r"""Converts wavelength values to frequency values and vice versa.
.. math::
value_{\text{target}} = c / \boldsymbol{value}
Args:
value (npt.ArrayLike): Single value or array of wavelengths in nm or frequencies in Hz.
Returns:
npt.ArrayLike: Frequencies in Hz or wavelengths in nm.
"""
return sc.speed_of_light / (value * 1e-9)
[docs]
def conversion_wavelength_wavenumber(value: npt.ArrayLike) -> npt.ArrayLike:
r"""Converts wavelength values to wavenumber values and vice versa.
.. math::
value_{\text{target}} = 1 / \boldsymbol{value}
Args:
value (npt.ArrayLike): Single value or array of wavelengths in nm
or wavenumbers in :math:`\text{cm}^{-1}`.
Returns:
npt.ArrayLike: Wavenumbers in :math:`\text{cm}^{-1}` or wavelengths in nm.
"""
return 1e7 / value
#########################################################
# Rotations
[docs]
def rotation_euler(p: float, n: float, r: float) -> npt.NDArray:
"""Returns rotation matrix defined by Euler angles p, n, r.
Successive rotations : z,x',z'
Note : The inverse rotation is -r, -n, -p
Args:
p (float): precession angle, 1st rotation, around z (0..360°).
n (float): nutation angle, 2nd rotation, around x' (0..180°).
r (float): 3rd rotation, around z' (0..360°).
Returns:
npt.NDArray: rotation matrix :math:`M_R`
"""
p = np.deg2rad(p)
n = np.deg2rad(n)
r = np.deg2rad(r)
c1 = np.cos(p)
s1 = np.sin(p)
c2 = np.cos(n)
s2 = np.sin(n)
c3 = np.cos(r)
s3 = np.sin(r)
return np.array(
[
[c1 * c3 - s1 * c2 * s3, -c1 * s3 - s1 * c2 * c3, s1 * s2],
[s1 * c3 + c1 * c2 * s3, -s1 * s3 + c1 * c2 * c3, -c1 * s2],
[s2 * s3, s2 * c3, c2],
]
)
[docs]
def rotation_v(v: npt.ArrayLike) -> npt.NDArray:
r"""Returns rotation matrix defined by a rotation vector v.
The calculation is made with the matrix exponential
:math:`M_R = \exp(W)`, with :math:`W_{ij} = - ε_{ijk} V_{k}`,
where :math:`ε_{ijk}` is the Levi-Civita antisymmetric tensor.
If V is separated in a unit vector v and a magnitude θ, V = θ·v, with
θ = ∥V∥, the calculation of the matrix exponential is avoided, and only
sin(θ) and cos(θ) are needed instead.
Note : The inverse rotation is -v
Args:
v (npt.ArrayLike): rotation vector (list or array)
Returns:
npt.NDArray: rotation matrix :math:`M_R`
"""
# fmt: off
m_w = np.array([[0, -v[2], v[1]],
[v[2], 0, -v[0]],
[-v[1], v[0], 0]])
# fmt: on
return scipy_expm(m_w)
[docs]
def rotation_v_theta(v: npt.ArrayLike, theta: float) -> npt.NDArray:
"""Returns rotation matrix defined by a unit rotation vector and an angle.
Notes : The inverse rotation is (v,-theta)
Args:
v (npt.ArrayLike): unit vector orienting the rotation (list or array)
theta (float): rotation angle around v in degrees
Returns:
npt.NDArray: rotation matrix :math:`M_R`
"""
# fmt: off
m_w = np.array([[0, -v[2], v[1]],
[v[2], 0, -v[0]],
[-v[1], v[0], 0]])
# fmt: on
return (
np.identity(3)
+ m_w * np.sin(np.deg2rad(theta))
+ np.linalg.matrix_power(m_w, 2) * (1 - np.cos(np.deg2rad(theta)))
)
[docs]
def detect_encoding(fname: str) -> str:
r"""Detects the encoding of file fname.
Args:
fname (str): Filename
Returns:
str: Encoding identifier string.
"""
with open(fname, "rb") as f:
raw_data = f.read()
result = chardet.detect(raw_data)
return result["encoding"]