Skip to content

Commit

Permalink
Merge pull request #14 from sun-data/performance/multilayer-efficienc…
Browse files Browse the repository at this point in the history 
…y-polarization

Performance improvements to `optika.materials.multilayer_efficiency()`
  • Loading branch information
@byrdie
byrdie authored Mar 13, 2024
2 parents a348808 + 615d438 commit fc3e524
Show file tree
Hide file tree
Showing 11 changed files with 200 additions and 348 deletions.
69 changes: 28 additions & 41 deletions optika/materials/_layers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,10 @@ def thickness(self) -> u.Quantity | na.AbstractScalar:
def transfer(
self,
wavelength: u.Quantity | na.AbstractScalar,
direction: na.AbstractCartesian3dVectorArray,
polarization: Literal["s", "p"],
direction: float | na.AbstractScalar,
polarized_s: bool | na.AbstractScalar,
n: float | na.AbstractScalar,
normal: na.AbstractCartesian3dVectorArray,
) -> tuple[na.AbstractScalar, na.Cartesian3dVectorArray, na.Cartesian2dMatrixArray]:
) -> tuple[na.AbstractScalar, na.AbstractScalar, na.Cartesian2dMatrixArray]:
"""
Compute the index of refraction, internal propagation direction,
and the transfer matrix for this layer,
Expand All @@ -54,15 +53,13 @@ def transfer(
wavelength
The wavelength of the incident light in the medium before this layer.
direction
The propagation direction of the incident light in the medium before
this layer.
polarization
Flag controlling whether the incident light is :math:`s` or :math:`p`
polarized.
The component of the incident light's propagation direction in the
medium before this layer antiparallel to the surface normal.
polarized_s
If :obj:`True`, the incident light is :math:`s`-polarized.
If :obj:`False`, the incident light is :math:`p`-polarized.
n
The complex index of refraction of the medium before this layer.
normal
The vector perpendicular to the surface of this layer.
"""

def plot(
Expand Down Expand Up @@ -165,39 +162,35 @@ def _chemical(self) -> optika.chemicals.AbstractChemical:
def n(
self,
wavelength: u.Quantity | na.AbstractScalar,
) -> na.AbstractScalar:
) -> float | na.AbstractScalar:
"""
The complex index of refraction of this layer.
"""
return self._chemical.n(wavelength)
if self.chemical is None:
return 1
else:
return self._chemical.n(wavelength)

def transfer(
self,
wavelength: u.Quantity | na.AbstractScalar,
direction: na.AbstractCartesian3dVectorArray,
polarization: Literal["s", "p"],
direction: float | na.AbstractScalar,
polarized_s: bool | na.AbstractScalar,
n: float | na.AbstractScalar,
normal: na.AbstractCartesian3dVectorArray,
) -> tuple[na.AbstractScalar, na.Cartesian3dVectorArray, na.Cartesian2dMatrixArray]:
) -> tuple[na.AbstractScalar, na.AbstractScalar, na.Cartesian2dMatrixArray]:

n_internal = self.n(wavelength)

direction_internal = snells_law(
wavelength=wavelength / np.real(n),
direction=direction,
index_refraction=np.real(n),
index_refraction_new=np.real(n_internal),
normal=normal,
)
angle_internal = np.arcsin(n * np.sin(np.arccos(direction)) / n_internal)
direction_internal = np.cos(angle_internal)

refraction = matrices.refraction(
wavelength=wavelength,
direction_left=direction,
direction_right=direction_internal,
polarization=polarization,
polarized_s=polarized_s,
n_left=n,
n_right=n_internal,
normal=normal,
interface=self.interface,
)

Expand All @@ -206,7 +199,6 @@ def transfer(
direction=direction_internal,
thickness=self.thickness,
n=n_internal,
normal=normal,
)

transfer = refraction @ propagation
Expand Down Expand Up @@ -376,21 +368,19 @@ def thickness(self) -> u.Quantity | na.AbstractScalar:
def transfer(
self,
wavelength: u.Quantity | na.AbstractScalar,
direction: na.AbstractCartesian3dVectorArray,
polarization: Literal["s", "p"],
direction: float | na.AbstractScalar,
polarized_s: bool | na.AbstractScalar,
n: float | na.AbstractScalar,
normal: na.AbstractCartesian3dVectorArray,
) -> tuple[na.AbstractScalar, na.Cartesian3dVectorArray, na.Cartesian2dMatrixArray]:
) -> tuple[na.AbstractScalar, na.AbstractScalar, na.Cartesian2dMatrixArray]:

result = na.Cartesian2dIdentityMatrixArray()

for layer in self.layers:
n, direction, matrix_transfer = layer.transfer(
wavelength=wavelength,
direction=direction,
polarization=polarization,
polarized_s=polarized_s,
n=n,
normal=normal,
)
result = result @ matrix_transfer

Expand Down Expand Up @@ -485,28 +475,25 @@ def thickness(self) -> u.Quantity | na.AbstractScalar:
def transfer(
self,
wavelength: u.Quantity | na.AbstractScalar,
direction: na.AbstractCartesian3dVectorArray,
polarization: Literal["s", "p"],
direction: float | na.AbstractScalar,
polarized_s: bool | na.AbstractScalar,
n: float | na.AbstractScalar,
normal: na.AbstractCartesian3dVectorArray,
) -> tuple[na.AbstractScalar, na.Cartesian3dVectorArray, na.Cartesian2dMatrixArray]:
) -> tuple[na.AbstractScalar, na.AbstractScalar, na.Cartesian2dMatrixArray]:

period = LayerSequence(self.layers)

n, direction, start = period.transfer(
wavelength=wavelength,
direction=direction,
polarization=polarization,
polarized_s=polarized_s,
n=n,
normal=normal,
)

n, direction, periodic = period.transfer(
wavelength=wavelength,
direction=direction,
polarization=polarization,
polarized_s=polarized_s,
n=n,
normal=normal,
)

periodic = periodic.power(self.num_periods - 1)
Expand Down
161 changes: 61 additions & 100 deletions optika/materials/_multilayers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,13 +30,14 @@

def multilayer_efficiency(
wavelength: u.Quantity | na.AbstractScalar,
direction: None | na.AbstractCartesian3dVectorArray = None,
polarization: None | Literal["s", "p"] = None,
direction: float | na.AbstractScalar = 1,
n: float | na.AbstractScalar = 1,
layers: Sequence[AbstractLayer] | optika.materials.AbstractLayer = None,
substrate: None | Layer = None,
normal: None | na.AbstractCartesian3dVectorArray = None,
) -> tuple[na.AbstractScalar, na.AbstractScalar]:
) -> tuple[
optika.vectors.PolarizationVectorArray,
optika.vectors.PolarizationVectorArray,
]:
r"""
Calculate the reflectivity and transmissivity of a multilayer
film or coating using the method in :cite:t:`Windt1998`.
Expand All @@ -46,12 +47,9 @@ def multilayer_efficiency(
wavelength
The wavelength of the incident light in vacuum.
direction
The propagation direction of incident light in the ambient medium.
If :obj:`None`, the propagation direction is assumed to be
:math:`+\hat{z}`
polarization
The polarization state of the incident light.
If :obj:`None`, the light is assumed to be unpolarized.
The component of the incident light's propagation direction in the
ambient medium antiparallel to the surface normal.
Default is to assume normal incidence.
n
The complex index of refraction of the ambient medium.
layers
Expand All @@ -62,9 +60,6 @@ def multilayer_efficiency(
A layer representing the substrate supporting the multilayer stack.
The thickness of this layer is ignored.
If :obj:`None`, then the substrate is assumed to be a vacuum.
normal
A vector perpendicular to the interface between successive layers.
If :obj:`None`, the normal vector is assumed to be :math:`-\hat{z}`
Examples
--------
Expand Down Expand Up @@ -100,7 +95,7 @@ def multilayer_efficiency(
fig, ax = plt.subplots()
na.plt.plot(
wavelength,
transmissivity,
transmissivity.average,
ax=ax,
axis="wavelength",
label="Zr",
Expand Down Expand Up @@ -155,7 +150,7 @@ def multilayer_efficiency(
fig, ax = plt.subplots()
na.plt.plot(
wavelength,
reflectivity,
reflectivity.average,
ax=ax,
axis="wavelength",
label=rf"Si/Mo $\times$ {layers.num_periods}",
Expand Down Expand Up @@ -211,7 +206,7 @@ def multilayer_efficiency(
fig, ax = plt.subplots()
na.plt.plot(
wavelength,
reflectivity,
reflectivity.average,
ax=ax,
axis="wavelength",
label=r"$\Gamma=" + thickness_ratio.astype(str).astype(object) + "$",
Expand Down Expand Up @@ -353,93 +348,62 @@ def multilayer_efficiency(
The :class:`tuple` :math:`(R, T)` is the quantity returned by this function.
"""
if direction is None:
direction_ambient = na.Cartesian3dVectorArray(0, 0, 1)
else:
direction_ambient = direction

if normal is None:
normal = na.Cartesian3dVectorArray(0, 0, -1)

direction_ambient = direction
n_ambient = n

if substrate is not None:
n_substrate = substrate.n(wavelength)
direction_substrate = snells_law(
wavelength=wavelength / np.real(n_ambient),
direction=direction_ambient,
index_refraction=np.real(n_ambient),
index_refraction_new=np.real(n_substrate),
normal=normal,
)
interface_substrate = substrate.interface
substrate = dataclasses.replace(substrate, thickness=0 * u.nm)
else:
n_substrate = 1
direction_substrate = direction_ambient
interface_substrate = None

cos_theta_ambient = direction_ambient @ normal
cos_theta_substrate = direction_substrate @ normal

if polarization is None:
polarization = ["s", "p"]
elif polarization == "s" or polarization == "p":
polarization = [polarization]
else: # pragma: nocover
raise ValueError(
f"Unrecognized polarization state '{polarization}'."
f"Only 's', 'p', or `None` is supported."
)
substrate = optika.materials.Layer(thickness=0 * u.nm)

polarized_s = na.ScalarArray(
ndarray=np.array([True, False]),
axes="_polarization",
)

if layers is None:
layers = []
if not isinstance(layers, optika.materials.AbstractLayer):
layers = optika.materials.LayerSequence(layers)

reflectivity = []
transmissivity = []

for pol in polarization:
n, direction, m = layers.transfer(
wavelength=wavelength,
direction=direction_ambient,
polarization=pol,
n=n_ambient,
normal=normal,
)
m_substrate = matrices.refraction(
wavelength=wavelength,
direction_left=direction,
direction_right=direction_substrate,
polarization=pol,
n_left=n,
n_right=n_substrate,
normal=normal,
interface=interface_substrate,
)
m = m @ m_substrate

r = m.y.x / m.x.x
t = 1 / m.x.x

if pol == "s":
q_ambient = cos_theta_ambient * n_ambient
q_substrate = cos_theta_substrate * n_substrate
else:
q_ambient = cos_theta_ambient / n_ambient
q_substrate = cos_theta_substrate / n_substrate

R = np.square(np.abs(r))
T = np.square(np.abs(t)) * np.real(q_substrate / q_ambient)

reflectivity.append(R)
transmissivity.append(T)

reflectivity = na.stack(reflectivity, axis="polarization")
transmissivity = na.stack(transmissivity, axis="polarization")

reflectivity = reflectivity.mean("polarization")
transmissivity = transmissivity.mean("polarization")
n, direction, m = layers.transfer(
wavelength=wavelength,
direction=direction_ambient,
polarized_s=polarized_s,
n=n_ambient,
)
n_substrate, direction_substrate, m_substrate = substrate.transfer(
wavelength=wavelength,
direction=direction,
polarized_s=polarized_s,
n=n,
)
m = m @ m_substrate

r = m.y.x / m.x.x
t = 1 / m.x.x

impedance_ambient = np.where(polarized_s, n_ambient, 1 / n_ambient)
impedance_substrate = np.where(polarized_s, n_substrate, 1 / n_substrate)

q_ambient = direction_ambient * impedance_ambient
q_substrate = direction_substrate * impedance_substrate

reflectivity = np.square(np.abs(r))
transmissivity = np.square(np.abs(t)) * np.real(q_substrate / q_ambient)

index_s = dict(_polarization=0)
index_p = dict(_polarization=1)

reflectivity = optika.vectors.PolarizationVectorArray(
s=reflectivity[index_s],
p=reflectivity[index_p],
)

transmissivity = optika.vectors.PolarizationVectorArray(
s=transmissivity[index_s],
p=transmissivity[index_p],
)

return reflectivity, transmissivity

Expand Down Expand Up @@ -530,14 +494,12 @@ def efficiency(

reflectivity, transmissivity = multilayer_efficiency(
wavelength=wavelength,
direction=rays.direction,
polarization=None,
direction=-rays.direction @ normal,
n=n,
layers=self.layers,
substrate=None,
normal=normal,
)
return transmissivity
return transmissivity.average

def plot_layers(
self,
Expand Down Expand Up @@ -588,13 +550,12 @@ def efficiency(

reflectivity, transmissivity = multilayer_efficiency(
wavelength=wavelength,
direction=rays.direction,
polarization=None,
direction=-rays.direction @ normal,
n=n,
layers=self.layers,
substrate=self.substrate,
)
return reflectivity
return reflectivity.average

def plot_layers(
self,
Expand Down
Loading

0 comments on commit fc3e524

Please sign in to comment.