`PDielec.GTMcore`#

The pyGTM module.

It has been heavily modified by John Kendrick for inclusion in the PDielec library Quite a lot has been removed as the fields spatial distribution was not needed Thanks to the authors of the original code for all their hard work

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>.

This module implements the generalized 4x4 transfer matrix (GTM) method poposed in Passler, N. C. and Paarmann, A., JOSA B 34, 2128 (2017) and corrected in JOSA B 36, 3246 (2019), as well as the layer-resolved absorption proposed in Passler, Jeannin and Paarman. This code uses inputs from D. Dietze’s FSRStools library https://github.com/ddietze/FSRStools

Please cite the relevant associated publications if you use this code.

Author:

Mathieu Jeannin math.jeannin@free.fr (permanent)

Affiliations:

Laboratoire de Physique de l’Ecole Normale Superieure (2019)
Centre de Nanosciences et Nanotechnologies (2020-2021)

Layers are represented by the Layer class that holds all parameters describing the optical properties of a single layer. The optical system is assembled using the System class.

Change log:

01-10-2024:

Moving back to 63 bit arithmetic as scattering formalism is stable

01-01-2023:

Major changes introduced by John Kendrick to make the package compatible with PDielec and PDGUI.

Allowed only full 3x3 tensors for the permittivity

Changed from 64 to 128 bit arithmetic where possible

Introduced a scattering matrix formalism

Added treatment of incoherence

15-10-2021:

Fixed rounding error bug in lag.eig() causing the program to crash randomly for negligibly small imaginary parts of the wavevectors

Corrected a sign error in gamma32 that lead to field discontinuities

23 June 2021:

integrated the code into pdielec and pdgui

19-03-2020:

Adapted the code to compute the layer-resolved absorption as proposed by Passler et al. (https://arxiv.org/abs/2002.03832), using System.calculate_Poynting_Absorption_vs_z().

Include the correct calculation of intensity transmission coefficients in System.calculate_r_t(). This BREAKS compatibility with the previous definition of the function.

Corrected bugs in System.calculate_Efield() and added magnetic field option

Adapted System.calculate_Efield() to allow hand-defined, irregular grid and a shorthand to compute only at layers interfaces. Regular grid with fixed resolution is left as an option.

20-09-2019:

Added functions in the System class to compute in-plane wavevector of guided modes and dispersion relation for such guided surface modes. This is highly prospective as it depends on the robustness of the minimization procedure (or the lack of thereoff)

Module Contents#

Classes#

`CoherentLayer`	Define a coherent layer inherits from Layer class.
`IncoherentAveragePhaseLayer`	Define an incoherent layer using an average phase in the propagation matrix.
`IncoherentIntensityLayer`	Define an incoherent layer using intensity transfer matrices.
`IncoherentPhaseLayer`	Define an incoherent layer using Arteaga's modification of the phase in the propagation matrix.
`IncoherentThickLayer`	Define an incoherent layer using the thick slab approximation.
`Layer`	Layer class. An instance is a single layer.
`SMatrix`	A class for storing and manipulating scattering matrices.
`ScatteringMatrixSystem`	Define a system of layers which is described by a Scattering Matrix method.
`SemiInfiniteLayer`	Define a semi-infinite layer, inherits from CoherentLayer class.
`System`	System class. An instance is an optical system with substrate, superstrate, and layers.
`TransferMatrixSystem`	Define a system of layers which is described by a Transfer Matrix method.

Functions#

`exact_inv_2x2`(M)	Calculate the inverse of 2x2 complex matrix, M.
`exact_inv_3x3`(M)	Calculate the inverse of a 3x3 complex matrix M.
`exact_inv_4x4`(M)	Compute the 'exact' inverse of a 4x4 matrix using the analytical result.
`vacuum_eps`(f)	Vacuum permittivity function.

Attributes#

`c_const`
`eps0`

class PDielec.GTMcore.CoherentLayer(layer, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Bases: Layer

Define a coherent layer inherits from Layer class.

calculate_scattering_matrix(b)[source]#

Calculate the scattering matrix of this layer with layer b.

This method calculates the scattering matrix based on equations 18-21 in the PyLama paper. However, it’s important to note that the definition of a scattering matrix used here differs from that in the paper.

Parameters#

ba pyGTM layer class: The second layer for which the scattering matrix will be calculated

Returns#

None

Notes#

Explain any important details about the implementation, external references like the PyLama paper, and any differences in definitions or approaches.

class PDielec.GTMcore.IncoherentAveragePhaseLayer(layer, percentage_incoherence=100, number_of_samples=4, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Bases: CoherentLayer

Define an incoherent layer using an average phase in the propagation matrix.

Inherits from CoherentLayer

calculate_propagation_exponents(f)[source]#

Calculate the matrix Ki (or Pi depending on the paper).

A phase shift is included in the calculation. The phase shift can vary between 0 and \(\pi\).

Parameters#

ffloat: Frequency

Returns#

None

class PDielec.GTMcore.IncoherentIntensityLayer(layer, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Bases: CoherentLayer

Define an incoherent layer using intensity transfer matrices.

Inherits from the CoherentLayer class.

class PDielec.GTMcore.IncoherentPhaseLayer(layer, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Bases: CoherentLayer

Define an incoherent layer using Arteaga’s modification of the phase in the propagation matrix.

Inherits from CoherentLayer

calculate_propagation_exponents(f)[source]#

Calculate the matrix Ki (or Pi) depending on the paper.

This routine is taken from Arteaga et al., Thin Solid Films 2014, 571, 701-705. The propagation matrix is suitable for incoherent films. In Passler, the ordering is: Arteaga: 0 p/o -> p -> 1 s/e -> p <- 2 p/o <- s -> 3 s/e <- s <-

Parameters#

ffloat: Frequency

Returns#

None

class PDielec.GTMcore.IncoherentThickLayer(layer, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Bases: CoherentLayer

Define an incoherent layer using the thick slab approximation.

Inherits from CoherentLayer

calculate_propagation_matrix(f)[source]#

Routine to calculate the matrix Ki using the thick slab approximation.

Parameters#

ffloat: Frequency

Returns#

None

calculate_scattering_matrix(b)[source]#

Calculate the scattering matrix of this layer with layer b.

Parameters#

bpyGTM layer: The second layer

Returns#

None

class PDielec.GTMcore.Layer(thickness=1e-06, epsilon=None, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Layer class. An instance is a single layer.

The inherited layer classes wich are used in PDielec are:

CoherentLayer

SemiInfiniteLayer

IncoherentIntensityLayer

IncoherentAveragePhaseLayer

IncoherentPhaseLayer

IncoherentThickLayer

Attributes#

thicknessfloat: Thickness of the layer in meters (m).
epsiloncomplex function: Function epsilon(frequency) representing the full dielectric constant tensor.
thetafloat: Euler angle theta (colatitude) in radians (rad).
phifloat: Euler angle phi in radians (rad).
psifloat: Euler angle psi in radians (rad).

Notes#

If instantiated with default values, it generates a 1 µm thick layer of air. Properties can be checked/changed dynamically using the corresponding get/set methods.

jk_shift#

qsd_thr = 1e-10#

zero_thr = 1e-10#

calculate_ai(zeta)[source]#

Calculate A_i.

Boundary matrix \(A_i\) of the layer.

Parameters#

zetafloat: Reduced in-plane wavevector kx/k0.

calculate_epsilon(f_in)[source]#

Set the value of epsilon in the (rotated) lab frame.

Parameters#

f_infloat: frequency (in Hz)

Returns#

None

Notes#

The values are set according to the epsilon_fi (i=1..3) functions defined using the set_epsilon() method, at the given frequency f. The rotation with respect to the lab frame is computed using the Euler angles.

Use explicitly if you don’t use the Layer.update() function! Modification by JK to allow the use of a full dielectric tensor.

calculate_gamma(zeta)[source]#

Calculate the gamma matrix.

Parameters#

zetacomplex: In-plane reduced wavevector kx/k0.

Returns#

None

calculate_matrices(zeta)[source]#

Calculate the principal matrices necessary for the GTM algorithm.

Parameters#

zetacomplex: In-plane reduced wavevector kx/k0 in the system.

Returns#

None

Notes#

Note that zeta is conserved through the whole system and set externally using the angle of incidence and System.superstrate.epsilon[0,0] value.

Requires prior execution of calculate_epsilon().

calculate_propagation_exponents(f)[source]#

Calculate the propagation exponents.

Parameters#

ffloat: The frequency

calculate_propagation_matrix(f)[source]#

Calculate the matrix Ki (or Pi depending on the paper).

Parameters#

ffloat: Frequency

Returns#

Kindarray: A 4x4-array representing the boundary matrix \(K_i\) of the layer.

calculate_q()[source]#: Calculate the 4 out-of-plane wavevectors for the current layer.

Parameters#

None

Returns#

None

Notes#

From this we also get the Poynting vectors. Wavevectors are sorted according to (trans-p, trans-s, refl-p, refl-s). Birefringence is determined according to a threshold value qsd_thr set at the beginning of the script.

calculate_transfer_matrix(f, zeta)[source]#

Compute the transfer matrix of the whole layer \(T_i=A_iP_iA_i^{-1}\).

Parameters#

ffloat: Frequency (in Hz).
zetacomplex: Reduced in-plane wavevector kx/k0.

Returns#

None

isCoherent()[source]#

Return True if the layer is a coherent layer.

Parameters#

None

Returns#

bool: True if the layer is coherent, False if not.

set_epsilon(epsilon_function)[source]#

Set the dielectric function.

Parameters#

epsilon_functioncomplex function: A function representing the complex dielectric function.

Returns#

None

set_euler(theta, phi, psi)[source]#

Set the values for the Euler rotations angles.

Parameters#

thetafloat: Euler angle theta (colatitude) in rad.
phifloat: Euler angle phi in rad.
psifloat: Euler angle psi in rad.

Returns#

None

set_thickness(thickness)[source]#

Set the layer thickness.

Parameters#

thicknessfloat: The layer thickness (in m).

Returns#

None

update_sm(f, zeta)[source]#

Shortcut to recalculate all layer properties for scattering method.

Appropriate for a scattering matrix method. This avoids the calculation of the exponential; it just calculates the exponents.

Parameters#

ffloat: frequency (in Hz)
zetacomplex: reduced in-plane wavevector kx/k0

Calculates#

Aiarray_like, shape (4, 4): Boundary matrix \(A_i\) of the layer.
Kiarray_like, shape (4, 4): Propagation matrix \(K_i\) of the layer.
Ai_invarray_like, shape (4, 4): Inverse of the \(A_i\) matrix.
Tiarray_like, shape (4, 4): Transfer matrix of the whole layer.

update_tm(f, zeta)[source]#

Shortcut to recalculate all layer properties.

Appropriate for a transfer matrix method

Parameters#

ffloat: frequency (in Hz)
zetacomplex: reduced in-plane wavevector kx/k0

Returns#

Ai4x4-array: Boundary matrix \(A_i\) of the layer
Ki4x4-array: Propagation matrix \(K_i\) of the layer
Ai_inv4x4-array: Inverse of the \(A_i\) matrix
Ti4x4-array: Transfer matrix of the whole layer

class PDielec.GTMcore.SMatrix(S=None)[source]#

A class for storing and manipulating scattering matrices.

S-matrices are square matrices used in physics and engineering to describe the scattering and reflection of waves or particles. They are especially useful in the field of microwave engineering, quantum mechanics, and optical physics.

Attributes#

S11np.ndarray: The top-left 2x2 submatrix of the S-matrix, representing input to input scattering.
S22np.ndarray: The bottom-right 2x2 submatrix of the S-matrix, representing output to output scattering.
S21np.ndarray: The bottom-left 2x2 submatrix of the S-matrix, representing input to output transmission.
S12np.ndarray: The top-right 2x2 submatrix of the S-matrix, representing output to input transmission.
Snp.ndarray: A 4x4 composite S-matrix combining S11, S22, S21, and S12.

Methods#

unitMatrix(): Initializes the S-matrix to a unit matrix with appropriate S11, S22, S21, and S12.
redheffer(b): Performs the Redheffer star product, a specialized matrix multiplication for S-matrices, with another SMatrix instance b.
calculateS(): Reconstructs the composite 4x4 S-matrix from its constituent submatrices (S11, S22, S21, S12).

Examples#

>>> import numpy as np
>>> from smatrix import SMatrix  # Assuming this class is saved in smatrix.py
>>> s_matrix_data = np.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [12, 13, 14, 15]])
>>> smatrix = SMatrix(S=s_matrix_data)
>>> print(smatrix.S)  # This displays the initialized S-matrix

Notes#

The Redheffer star product is particularly useful for cascading two-port networks in microwave engineering and quantum physics. It allows for the efficient calculation of the overall scattering matrix of the combined system.

calculateS()[source]#

Calculate and assign the scattering parameter matrix, S.

Fills a 4 by 4 matrix, S, with the sub-matrices S11, S22, S12, and S21 representing the scattering parameters of a two-port network. The S matrix combines these sub-matrices as follows:

S11 is placed in the top-left quadrant.
S22 is placed in the bottom-right quadrant.
S12 is placed in the top-right quadrant.
S21 is placed in the bottom-left quadrant.

Parameters#

None

Returns#

Sndarray, shape (4, 4), dtype=np.cdouble: The complete scattering parameter matrix of the two-port network.

Notes#

The scattering parameters (S11, S22, S12, S21) must be defined as attributes of the instance before calling this method.

redheffer(b)[source]#

Calculate the product of S-matrices.

Parameters#

bScatteringMatrix: The second scattering matrix to be used the the multiplication

Returns#

SabScatteringMatrix: The product of Sa and Sb

unitMatrix()[source]#

Reset the scattering parameters to form a unit matrix.

Resets the scattering parameters S11, S22, S21, and S12 of the instance to represent a unit matrix, where S11 and S22 are set to zero matrices and S21 and S12 are set to identity matrices. After the reset, calculates the scattering parameters through calculateS.

Parameters#

None

Returns#

None

Notes#

This method is typically used to initialize or reset the scattering parameters in a two-port network representation to a known state before performing further operations or calculations.

class PDielec.GTMcore.ScatteringMatrixSystem(substrate=None, superstrate=None, layers=None)[source]#

Bases: System

Define a system of layers which is described by a Scattering Matrix method.

This implements additional routines needed for the scattering matrix method

calculate_GammaStar(f, zeta_sys)[source]#

Calculate the whole system’s scattering matrix.

This replaces the equivalent routine used in the transfer matrix method.

Parameters#

ffloat: Frequency (Hz)
zeta_syscomplex: In-plane wavevector kx/k0

Returns#

Stotalndarray: System scattering matrix \(S\)

calculate_r_t(zeta_sys)[source]#

Calculate the reflectance and transmittance coefficients using scattering matrix information.

Calculate various field and intensity reflection and transmission coefficients, as well as the 4-valued vector of transmitted field.

Parameters#

zeta_syscomplex: Incident in-plane wavevector

Returns#

r_outarray_like, shape (4,): Complex field reflection coefficients (r_out=([rpp,rps,rss,rsp])).
R_outarray_like, shape (4,): Real intensity reflection coefficients (R_out=([Rpp,Rss,Rsp,Tps])).
t_outarray_like, shape (4,): Complex field transmission coefficients (t=([tpp, tps, tsp, tss])).
T_outarray_like, shape (4,): Real intensity transmission coefficients (T_out=([Tp,Ts])) (mode-inselective).

calculate_scattering_matrices()[source]#

Calculate the scattering matrices in every layer.

Loop through all the layers calculating the scattering matrix of this layer with the next. Care is taken when dealing with the superstrate and the substrate

Parameters#

None

Returns#

None

Note#

The scattering matrix is stored as an attribute If the layer system is changed in anyway all the scattering matrices must be recomputed

class PDielec.GTMcore.SemiInfiniteLayer(layer, theta=0, phi=0, psi=0, exponent_threshold=700)[source]#

Bases: CoherentLayer

Define a semi-infinite layer, inherits from CoherentLayer class.

calculate_propagation_exponents(f)[source]#

Calculate the propagation exponents.

For the semi-infinite case the propagation exponents are set to 0

Parameters#

ffloat: The frequency (ignored for this object)

Returns#

None

calculate_propagation_matrix(f)[source]#

Calculate the matrix Ki (or Pi) depending on the paper.

This routine makes the material transparent. It is therefore suitable for a semi-infinite layer.

Parameters#

ffloat: Frequency

Returns#

Ki : tensor (3x3)

class PDielec.GTMcore.System(substrate=None, superstrate=None, layers=None)[source]#

System class. An instance is an optical system with substrate, superstrate, and layers.

Derived classes used by PDielec are:

TransferMatrixSystem

ScatteringMatrixSystem

Attributes#

thetafloat: Angle of incidence, in radians.
substrateLayer: The substrate layer. Defaults to vacuum (an empty layer instance).
superstrateLayer: The superstrate layer, defaults to vacuum (an empty layer instance).
layerslist of Layer: List of the layers in the system.

Notes#

Layers can be added and removed (not inserted).

The whole system’s transfer matrix is computed using calculate_GammaStar(), which calls Layer.update() for each layer. General reflection and transmission coefficient functions are given; they require the prior execution of calculate_GammaStar(). The electric fields can be visualized in the case of an incident plane wave using calculate_Efield().

add_layer(layer)[source]#

Add a layer instance.

Parameters#

layerLayer: The layer to be added on the stack.

Returns#

None

Notes#

The layers are added in from superstrate to substrate order. Light is incident from the superstrate.

Note that this function adds a reference to layer to the list. If you are adding the same layer several times, be aware that if you change something for one of them, it changes for all of them.

calculate_GammaStar(f, zeta_sys)[source]#

Calculate the whole system’s transfer matrix.

If one layer is incoherent then intensities rather than amplitudes are used

Parameters#

ffloat: Frequency (Hz)
zeta_syscomplex: In-plane wavevector kx/k0

Returns#

GammaStarcomplex ndarray: 4x4 complex matrix representing the system transfer matrix \(\Gamma^{*}\)

calculate_r_t(zeta_sys)[source]#

Calculate various field and intensity reflection and transmission coefficients, as well as the 4-valued vector of transmitted field.

Parameters#

zeta_syscomplex: Incident in-plane wavevector

Returns#

r_outarray_like, shape (4,): Complex field reflection coefficients r_out=([rpp,rps,rss,rsp])
R_outarray_like, shape (4,): Real intensity reflection coefficients R_out=([Rpp,Rss,Rsp,Tps])
t_outarray_like, shape (4,): Complex field transmission coefficients t=([tpp, tps, tsp, tss])
T_outarray_like, shape (4,): Real intensity transmission coefficients T_out=([Tp,Ts]) (mode-inselective)

Notes#

IMPORTANT

As of version 19-03-2020: All intensity coefficients are now well defined. Transmission is defined mode-independently. It could be defined mode-dependently for non-birefringent substrates in future versions. The new definition of this function BREAKS compatibility with the previous one.
As of version 13-09-2019: Note that the field reflectivity and transmission coefficients r and t are well defined. The intensity reflection coefficient is also correct. However, the intensity transmission coefficients T are ill-defined so far. This will be corrected upon future publication of the correct intensity coefficients.

Note also the different ordering of the coefficients, for consistency with Passler’s matlab code.

del_layer(pos)[source]#

Remove a layer at given position. Does nothing for invalid position.

Parameters#

posint: Index of layer to be removed.

Returns#

None

get_all_layers()[source]#

Return the list of all layers in the system.

Returns#

self.layerslist: A list of all layers.

Returns#

None

get_layer(pos)[source]#

Get the layer at a given position.

Parameters#

posint: Position in the stack.

Returns#

self.layers[pos]Layer: The layer at the position pos.

get_substrate()[source]#

Return the System’s substrate.

Parameters#

None

Returns#

self.substrateLayer: The system substrate.

get_superstrate()[source]#

Return the System’s superstrate.

Parameters#

None

Returns#

self.superstrateLayer: The system superstrate.

initialize_sys(f)[source]#

Set the values of epsilon at the given frequency in all the layers.

Parameters#

ffloat: Frequency (Hz)

Returns#

None

Notes#

This function allows to define the in-plane wavevector (\(zeta\)) outside of the class, and thus to explore also guided modes of the system.

overflowErrors()[source]#

Return the total number of overflow errors encountered.

Parameters#

None

Returns#

int: The total number of overflow errors encountered.

set_substrate(sub)[source]#

Set the substrate.

Parameters#

subLayer: Instance of the Layer class, serving as the substrate.

Returns#

None

set_superstrate(sup)[source]#

Set the superstrate.

Parameters#

supLayer: Instance of the layer class, superstrate

Returns#

None

class PDielec.GTMcore.TransferMatrixSystem(substrate=None, superstrate=None, layers=None)[source]#

Bases: System

Define a system of layers which is described by a Transfer Matrix method.

This is basically a place holder for the original pyGTM methods and it is used to distinguish itself from a system which uses the scattering matrix method

PDielec.GTMcore.exact_inv_2x2(M)[source]#

Calculate the inverse of 2x2 complex matrix, M.

Parameters#

Mcomplex 2x2 array: The matrix of which to find the inverse.

Returns#

ndarray: The inverse of M as a 2x2 complex numpy array.

PDielec.GTMcore.exact_inv_3x3(M)[source]#

Calculate the inverse of a 3x3 complex matrix M.

Parameters#

Mcomplex 3x3 array: The matrix to be inverted.

Returns#

numpy.ndarray: The inverse of M as a 3x3 complex numpy array.

PDielec.GTMcore.exact_inv_4x4(M)[source]#

Compute the ‘exact’ inverse of a 4x4 matrix using the analytical result.

Parameters#

Marray-like, shape (4, 4): Matrix to be inverted, consisting of float or complex numbers.

Returns#

array-like, shape (4, 4), complex: Inverse of this matrix or Moore-Penrose approximation if matrix cannot be inverted.

Notes#

This should give a higher precision and speed at a reduced noise, following D.Dietze’s implementation in FSRStools.

Parameters#

ffloat or 1D-array: Frequency (in Hz)

Returns#

epscomplex or 1D-array of complex: Complex value of the vacuum permittivity (1.0 + 0.0j)

PDielec.GTMcore.c_const#

PDielec.GTMcore.eps0#

PDielec.GTMcore#

Module Contents#

Classes#

Functions#

Attributes#

Parameters#

Returns#

Notes#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Returns#

Attributes#

Notes#

Parameters#

Parameters#

Returns#

Notes#

Parameters#

Returns#

Parameters#

Returns#

Notes#

Parameters#

Parameters#

Returns#

Parameters#

Returns#

Notes#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Calculates#

Parameters#

Returns#

Attributes#

Methods#

Examples#

Notes#

Parameters#

Returns#

Notes#

Parameters#

Returns#

Parameters#

Returns#

Notes#

Parameters#

Returns#

Parameters#

Returns#

Parameters#

Returns#

Note#

Parameters#

Returns#

Parameters#

Returns#

Attributes#

Notes#

Parameters#

Returns#

Notes#

Parameters#

Returns#

Parameters#

Returns#

Notes#

`PDielec.GTMcore`#