macromax package

The solver module calculates the solution to the wave equations. More specifically, the work is done in the iteration defined in the Solution.__iter__() method of the Solution class. The convenience function solve() is provided to construct a Solution object and iterate it to convergence using its Solution.solve() method.

Public attributes:

__version__: The MacroMax version number as a str.

solve(): The function to solve the wave problem.

Solution: The class that is used by the solve() function, which can be used for fine-control of the iteration or re-use.

Grid: A class representing uniformly spaced Cartesian grids and their Fourier Transforms.

log: The logging object of the macromax library. This can be used to make the messages more or less verbose.

backend: The sub-package with the back-end specifications.

macromax.solve(grid, vectorial=None, wavenumber=1.0, angular_frequency=None, vacuum_wavelength=None, current_density=None, source_distribution=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, initial_field=0.0, dtype=None, callback=<function <lambda>>)[source]

Function to find a solution for Maxwell’s equations in a media specified by the epsilon, xi, zeta, and mu distributions in the presence of a current source.

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – A Grid object or a Sequence of vectors with uniformly increasing values that indicate the positions in a plaid grid of sample points for the material and solution. In the one-dimensional case, a simple increasing Sequence of uniformly-spaced numbers may be provided as an alternative. The length of the ranges determines the data_shape, to which the source_distribution, epsilon, xi, zeta, mu, and initial_field must broadcast when specified as ndarrays.

  • vectorial (Optional[bool]) – a boolean indicating if the source and solution are 3-vectors-fields (True) or scalar fields (False).

  • wavenumber (Optional[Real]) – the wavenumber in vacuum = 2 pi / vacuum_wavelength. The wavelength in the same units as used for the other inputs/outputs.

  • angular_frequency (Optional[Real]) – alternative argument to the wavenumber = angular_frequency / c

  • vacuum_wavelength (Optional[Real]) – alternative argument to the wavenumber = 2 pi / vacuum_wavelength

  • current_density (Union[Complex, Sequence, ndarray, None]) – (optional, instead of source_distribution) An array or function that returns the free (vectorial) current density input distribution, J. The free current density has units of \(A m^-2\).

  • source_distribution (Union[Complex, Sequence, ndarray, None]) – (optional, instead of current_density) An array or function that returns the (vectorial) source input wave distribution. The source values relate to the current density, J, as 1j * angular_frequency * scipy.constants.mu_0 * J and has units of \(rad s^-1 H m^-1 A m^-2 = rad V m^-3\). More general, non-electro-magnetic wave problems can be solved using the source_distribution, as it does not rely on the vacuum permeability constant, \(mu_0\).

  • epsilon (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (tensor) epsilon that represents the permittivity at the points indicated by the grid specified as its input arguments.

  • xi (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) xi for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • zeta (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) zeta for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • mu (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) permeability at the points indicated by the grid specified as its input arguments.

  • refractive_index (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (complex) (tensor) refractive_index, as the square root of the permittivity, at the points indicated by the grid input argument.

  • bound (Optional[Bound]) – An object representing the boundary of the calculation volume. Default: None, PeriodicBound(grid)

  • initial_field (Union[Complex, Sequence, ndarray]) – optional start value for the E-field distribution (default: all zero E)

  • dtype – optional numpy datatype for the internal operations and results. This must be a complex number type as numpy.complex128 or np.complex64.

  • callback (Callable) – optional function that will be called with as argument this solver. This function can be used to check and display progress. It must return a boolean value of True to indicate that further iterations are required.

Returns:

The Solution object that has the E and H fields, as well as iteration information.

class macromax.Solution(grid, vectorial=None, wavenumber=1.0, angular_frequency=None, vacuum_wavelength=None, current_density=None, source_distribution=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, initial_field=0.0, dtype=None)[source]

Bases: object

__init__(grid, vectorial=None, wavenumber=1.0, angular_frequency=None, vacuum_wavelength=None, current_density=None, source_distribution=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, initial_field=0.0, dtype=None)[source]

Class a solution that can be further iterated towards a solution for Maxwell’s equations in a media specified by the epsilon, xi, zeta, and mu distributions.

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – A Grid object or a Sequence of vectors with uniformly increasing values that indicate the positions in a plaid grid of sample points for the material and solution. In the one-dimensional case, a simple increasing Sequence of uniformly-spaced numbers may be provided as an alternative. The length of the ranges determines the data_shape, to which the source_distribution, epsilon, xi, zeta, mu, and initial_field must broadcast when specified as `numpy.ndarray`s.

  • vectorial (Optional[bool]) – a boolean indicating if the source and solution are 3-vectors-fields (True) or scalar fields (False). Default: True, when vectorial nor the source is specified. Default: vectorial (True), unless the source field is scalar (False if first dimension is a singleton dimension).

  • wavenumber (Optional[Real]) – the wavenumber in vacuum = 2pi / vacuum_wavelength. The wavelength in the same units as used for the other inputs/outputs.

  • angular_frequency (Optional[Real]) – alternative argument to the wavenumber = angular_frequency / c

  • vacuum_wavelength (Optional[Real]) – alternative argument to the wavenumber = 2 pi / vacuum_wavelength

  • current_density (Union[Complex, Sequence, ndarray, None]) – (optional, instead of source_distribution) An array or function that returns the (vectorial) current density input distribution, J. The current density has units of \(A m^-2\).

  • source_distribution (Union[Complex, Sequence, ndarray, None]) – (optional, instead of current_density) An array or function that returns the (vectorial) source input wave distribution. The source values relate to the current density, J, as 1j * angular_frequency * scipy.constants.mu_0 * J and has units of \(rad s^-1 H m^-1 A m^-2 = rad V m^-3\). More general, non-electro-magnetic wave problems can be solved using the source_distribution, as it does not rely on the vacuum permeability constant, \(mu_0\).

  • epsilon (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (tensor) epsilon that represents the permittivity at the points indicated by the grid input argument.

  • xi (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) xi for bi-(an)isotropy at the points indicated by the grid input argument.

  • zeta (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) zeta for bi-(an)isotropy at the points indicated by the grid input argument.

  • mu (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) permeability at the points indicated by the grid input argument.

  • refractive_index (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (complex) (tensor) refractive_index, as the square root of the permittivity, at the points indicated by the grid input argument.

  • bound (Optional[Bound]) – An object representing the boundary of the calculation volume. Default: None, PeriodicBound(grid)

  • initial_field (Union[Complex, Sequence, ndarray]) – optional start value for the E-field distribution (default: all zero E)

  • dtype – optional numpy datatype for the internal operations and results. This must be a complex number type as numpy.complex128 or np.complex64.

property grid: Grid

The sample positions of the plaid sampling grid. This may be useful for displaying result axes.

Returns:

A Grid object representing the sample points of the fields and material.

property vectorial: bool

Boolean to indicates whether calculations happen on vectorial (True) or scalar (False) fields.

property dtype

The numpy equivalent data type used in the calculation. This is either np.complex64 or np.complex128.

property wavenumber: Real

The vacuum wavenumber, \(k_0\), used in the calculation.

Returns:

A scalar indicating the wavenumber used in the calculation.

property angular_frequency: Real

The angular frequency, \(\omega\), used in the calculation.

Returns:

A scalar indicating the angular frequency used in the calculation.

property wavelength: Real

The vacuum wavelength, \(\lambda_0\), used in the calculation.

Returns:

A scalar indicating the vacuum wavelength used in the calculation.

property magnetic: bool

Indicates if this media is considered magnetic.

Returns:

A boolean, True when magnetic, False otherwise.

property bound: Bound

The Bound object that defines the calculation boundaries.

property source_distribution: ndarray

The source distribution, i k0 mu_0 times the current density j.

Returns:

A complex array indicating the amplitude and phase of the source vector field. The dimensions of the array are [1|3, self.grid.shape], where the first dimension is 1 in case of a scalar field, and 3 in case of a vector field.

property j: ndarray

The free current density, j, of the source vector field.

Returns:

A complex array indicating the amplitude and phase of the current density vector field [A m^-2]. The dimensions of the array are [1|3, self.grid.shape], where the first dimension is 1 in case of a scalar field, and 3 in case of a vector field.

property E: ndarray

The electric field for every point in the sample space (SI units).

Returns:

A vector array with the first dimension containing Ex, Ey, and Ez,

while the following dimensions are the spatial dimensions.

property B: ndarray

The magnetic field for every point in the sample space (SI units). This is calculated from H and E.

Returns:

A vector array with the first dimension containing \(B_x, B_y, and B_z\), while the following dimensions are the spatial dimensions.

property D: ndarray

The displacement field for every point in the sample space (SI units). This is calculated from E and H.

Returns:

A vector array with the first dimension containing \(D_x, D_y, and D_z\), while the following dimensions are the spatial dimensions.

property H: ndarray

The magnetizing field for every point in the sample space (SI units). This is calculated from E.

Returns:

A vector array with the first dimension containing \(H_x, H_y, and H_z\), while the following dimensions are the spatial dimensions.

property S: ndarray

The time-averaged Poynting vector for every point in space. :return: A vector array with the first dimension containing \(S_x, S_y, and S_z\), while the following dimensions are the spatial dimensions.

property energy_density: ndarray

Returns the energy density, u.

Returns:

A real array indicating the energy density in space.

property stress_tensor: ndarray

Maxwell’s stress tensor for every point in space.

Returns:

A real and symmetric matrix-array with the stress tensor for every point in space. The units are \(N / m^2\).

property f: ndarray

The electromagnetic force density (force per SI unit volume, not per voxel).

Returns:

A vector array representing the electro-magnetic force exerted per unit volume. The first dimension contains \(f_x, f_y, and f_z\), while the following dimensions are the spatial dimensions. The units are \(N / m^3\).

property torque: ndarray

The electromagnetic force density (force per SI unit volume, not per voxel).

Returns:

A vector array representing the electro-magnetic torque exerted per unit volume. The first dimension contains torque_x, torque_y, and torque_z, while the following dimensions are the spatial dimensions. The units are \(N m / m^3 = N m^{-2}\).

property iteration: int

The current iteration number.

Returns:

An integer indicating how many iterations have been done.

property previous_update_norm: Real

The L2-norm of the last update, the difference between current and previous E-field.

Returns:

A positive scalar indicating the norm of the last update.

property residue: Real

Returns the current relative residue of the inverse problem \(E = H^{-1}S\). The relative residue is return as the l2-norm fraction \(||E - H^{-1}S|| / ||E||\), where H represents the vectorial Helmholtz equation following Maxwell’s equations and S the current density source. The solver searches for the electric field, E, that minimizes the preconditioned inverse problem.

Returns:

A non-negative real scalar that indicates the change in E with the previous iteration normalized to the norm of the current E.

__iter__()[source]

Returns an iterator that on __next__() yields this Solution after updating it with one cycle of the algorithm. Obtaining this iterator resets the iteration counter.

Usage:

for solution in Solution(...):
    if solution.iteration > 100:
        break
print(solution.residue)
solve(callback=<function Solution.<lambda>>)[source]

Runs the algorithm until the convergence criterion is met or until the maximum number of iterations is reached.

Parameters:

callback (Callable) – optional callback function that overrides the one set for the solver. E.g. callback=lambda s: s.iteration < 100

Returns:

This Solution object, which can be used to query e.g. the final field E using Solution.E.

class macromax.ScatteringMatrix(grid, vectorial=True, wavenumber=None, angular_frequency=None, vacuum_wavelength=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, dtype=None, callback=<function ScatteringMatrix.<lambda>>, caching=True, array=None)[source]

Bases: LiteralScatteringMatrix

A class representing scattering matrices.

__init__(grid, vectorial=True, wavenumber=None, angular_frequency=None, vacuum_wavelength=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, dtype=None, callback=<function ScatteringMatrix.<lambda>>, caching=True, array=None)[source]

Construct a scattering matrix object for a medium specified by a refractive index distribution or the corresponding epsilon, xi, zeta, and mu distributions. Each electromagnetic field distribution entering the material is scattered into a certain electromagnetic field distributions propagating away from it from both sides. The complex matrix relates the amplitude and phase of all N propagating input modes to all N propagating output modes. No scattering, as in vacuum, is indicated by the NxN identity matrix. There are N/2 input and output modes on either side of the scattering material. Mode i and i+N/2 correspond to plane wave traveling in opposing directions. The mode directions are taken in raster-scan order and only propagating modes are included. When polarization is considered, the modes come in pairs corresponding to two orthogonal linear polarizations.

The modes are encoded as a vector of length N = 2xMxP for 2 sides, M angles, and P polarizations.

  • First the N/2 modes propagating along the positive x-axis are considered, then those propagating in the reverse direction.

  • In each direction, M different angles (k-vectors) can be considered. We choose propagating modes on a uniformly-spaced plaid grid that includes the origin (corresponding to the k-vector along the x-axis). Modes not propagating along the x-axis, i.e. in the y-z-plane are not considered. The angles are ordered in raster-scan order from negative k_y to positive k_y (slow) and from negative k_z to positive k_z (fast). The grid axes dimensions correspond to x(0), y(1), z(2).

  • When polarization is considered, each angle has a pair of modes, one for each polarization. The first mode has the polarization oriented along the rotated y’-axis and the second mode along the rotated z’-axis. To avoid ambiguity for normal incidence, the Cartesian-coordinate system is rotated along the shortest possible path, i.e. along the axis that is normal to the original x-axis and the mode’s k-vector. All rotations are around the origin of the coordinate system, incurring no phase shift there.

Vectors can be converted to field distributions on the complete grid using the methods:

srcvec2freespace() super-position of free-space plane waves in the whole volume (fast)

srcvec2source() super-position of free-space plane waves at the source planes at the front and back (fast)

source2detfield() calculate the field in the whole volume using the solver.Solution object (slow)

detfield2detvec() vector corresponding to the detected field at the detection planes (fast). The fields

at those planes should only contain the outward propagating waves. Hence, inwards propagating waves should be subtracted before using this method!

srcvec2detfield() calculate the field in the whole volume and convert it to a detection vector (slow)

The latter is used in the matrix multiplication method: matmul, @

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – A Grid object or a Sequence of vectors with uniformly increasing values that indicate the positions in a plaid grid of sample points for the material and solution. In the one-dimensional case, a simple increasing Sequence of uniformly-spaced numbers may be provided as an alternative. The length of the ranges determines the data_shape, to which the source_distribution, epsilon, xi, zeta, mu, and initial_field must broadcast when specified as ndarrays.

  • vectorial (Optional[bool]) – a boolean indicating if the source and solution are 3-vectors-fields (True) or scalar fields (False).

  • wavenumber (Optional[Real]) – the wavenumber in vacuum = 2 pi / vacuum_wavelength. The wavelength in the same units as used for the other inputs/outputs.

  • angular_frequency (Optional[Real]) – alternative argument to the wavenumber = angular_frequency / c

  • vacuum_wavelength (Optional[Real]) – alternative argument to the wavenumber = 2 pi / vacuum_wavelength

  • epsilon (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) epsilon that represents the permittivity at the points indicated by the grid specified as its input arguments.

  • xi (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) xi for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • zeta (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) zeta for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • mu (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) permeability at the points indicated by the grid specified as its input arguments.

  • refractive_index (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) refractive_index = np.sqrt(permittivity) at the points indicated by the grid input argument.

  • bound (Optional[Bound]) – An object representing the boundary of the calculation volume. Default: None, PeriodicBound(grid)

  • dtype – optional numpy datatype for the internal operations and results. This must be a complex number type as numpy.complex128 or np.complex64.

  • callback (Callable) – optional function that will be called with as argument this solver. This function can be used to check and display progress. It must return a boolean value of True to indicate that further iterations are required.

  • caching (bool) – Cache field propagation calculations. By default, the results are cached for multiplications with basis vectors. Numerical errors might accumulate for certain superpositions. Setting this property to False will ensure that field propagations are always used and the constructor argument array is ignored.

  • array (Union[Complex, Sequence, ndarray, LinearOperator, None]) – Optional in case the matrix values have been calculated before and stored. If not specified, the matrix is calculated from the material properties. I specified, this must be a sequence or numpy.ndarray of complex numbers representing the matrix, or a function that returns one.

Returns:

The Solution object that has the E and H fields, as well as iteration information.

property grid: Grid

The calculation grid.

property vectorial: bool

Boolean to indicates whether calculations happen on polarized (True) or scalar (False) fields.

property caching: bool

When set to True, this object uses cached values instead of propagating the field through the scatterer. Otherwise, field propagation is used for all matrix operations. This can help avoid the accumulation of numerical errors.

srcvec2freespace(input_vector)[source]

Convert an input source vector to a superposition of plane waves at the origin and spreading over the whole volume. The input vector specifies the propagating modes in the far-field (the inverse Fourier transform of the fields at the sample origin). Incident waves at an angle will result in higher amplitudes to compensate for the reduction in propagation along the propagation axis through the entrance plane.

Used in ScatteringMatrix.srcvec2source to calculate the source field distribution before entering the scatterer.

Used in ScatteringMatrix.__matmul__ to distinguish incoming from back-scattered light.

Parameters:

input_vector (Union[Complex, Sequence, ndarray, LinearOperator]) – A source vector or array of shape [2, M, P], where the first axis indicates the side (front, back), the second axis indicates the propagation mode (direction, top-bottom-left-right), and the final axis indicates the polarization (1 for scalar, 2 for polarized: V-H).

Return type:

ndarray

Returns:

An nd-array with the field on the calculation grid. Its shape is (1, *self.grid.shape) for scalar

calculations and (3, *self.grid.shape) for vectorial calculations with polarization.

srcvec2source(input_vector, out=None)[source]

Converts a source vector into an (N+1)D-array with the source field at the front and back of the scatterer. The source field is such that it produces E-fields of unit intensity for unit vector inputs.

Used in self.vector2field() and self.__matmul__().

Parameters:
  • input_vector (Union[Complex, Sequence, ndarray, LinearOperator]) – A source vector with self.shape[1] elements. One value per side, per independent polarization (2), and per mode (inwards propagating k-vectors only).

  • out (Optional[ndarray]) – (optional) numpy array to store the result.

Return type:

ndarray

Returns:

The field distribution as an array of shape [nb_pol, *self.grid], where nb_pol = 3 for a vectorial calculation and 1 for a scalar calculation.

source2detfield(source, out=None)[source]

Calculates the (N+1)D-input-field distribution throughout the scatterer for a given source field distribution.

Parameters:
  • source (Union[Complex, Sequence, ndarray, LinearOperator]) – The source field distribution in the whole space.

  • out (Optional[ndarray]) – (optional) numpy array to store the result (shape: self.grid.shape, dtype: self.dtype).

Return type:

ndarray

Returns:

The field distribution as an array of shape [nb_pol, *self.grid], where nb_pol = 3 for a vectorial calculation and 1 for a scalar calculation.

srcvec2detfield(input_vector, out=None)[source]

Calculates the (N+1)D-input-field distribution throughout the scatterer for a given input source vector.

Used in self.__matmul__() and external code.

Parameters:
Return type:

ndarray

Returns:

The field distribution as an array of shape [nb_pol, *self.grid], where nb_pol = 3 for a vectorial calculation and 1 for a scalar calculation.

detfield2detvec(field)[source]

Converts the (N+1)D-output-field defined at the front and back detection planes of the scatterer to a detection vector that describes the resulting far-field distribution (the inverse Fourier transform of the fields at the sample origin). The fields at those planes should only contain the outward propagating waves. Hence, inwards propagating waves should be subtracted before using this method!

This is used in self.__matmul__().

Parameters:

field (Union[Complex, Sequence, ndarray, LinearOperator]) – The detected field in all space (of which only the detection space is used).

Return type:

ndarray

Returns:

Detection vector.

detvec2srcvec(vec)[source]

Convert forward propagating detection vector into a backward propagating (time-reversed) source vector.

Parameters:

vec (Union[Complex, Sequence, ndarray, LinearOperator]) – detection vector obtained from solution or scattering matrix multiplication.

Return type:

ndarray

Returns:

Time-reversed vector, that can be used as a source vector.

__setitem__(key, value)[source]

Updating this matrix is not possible. Use a Matrix object instead.

Parameters:
  • key – Index or slice.

  • value – The new value.

__array__(out=None)[source]

Lazily calculates the scattering matrix as a regular numpy.ndarray

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property back_reflection: BackReflectionMatrix

Select the quarter of the scattering matrix corresponding to the light that is reflected of the back. It indicates how the light coming from positive infinity is back reflected to positive infinity.

Returns:

The back-reflection matrix of shape self.shape // 2.

property backward_transmission: BackwardTransmissionMatrix

Select the backward-transmitted quarter of the scattering matrix. It indicates how the light coming from positive infinity is transmitted to negative infinity.

Returns:

The backward-transmission matrix of shape self.shape // 2.

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

property forward_transmission: ForwardTransmissionMatrix

Select the forward-transmitted quarter of the scattering matrix. It indicates how the light coming from negative infinity is transmitted to positive infinity.

Returns:

The forward-transmission matrix of shape self.shape // 2.

property front_reflection: FrontReflectionMatrix

Select the quarter of the scattering matrix corresponding to the light that is reflected of the front. It indicates how the light coming from negative infinity is back reflected to negative infinity.

Returns:

The front-reflection matrix of shape self.shape // 2.

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transfer(noise_level=0.0)

Calculates the transfer matrix, relating one side of the scatterer to the other side (top, bottom). Each side can have incoming and outgoing waves. This is in contrast to the scattering matrix, self.__array__, which relates incoming waves from both sides to outgoing waves from both sides. One can be calculated from the other using the matrix.convert() function, though this calculation may be ill-conditioned (sensitive to noise). Therefore, the optional argument noise_level should be used to indicate the root-mean-square expectation value of the measurement error. This avoids divisions by near-zero values and obtains a best estimate using Tikhonov regularization.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near) singular backwards transmission matrix.

Return type:

ndarray

Returns:

An nd-array with the transfer matrix relating top-to-bottom instead of in-to-out. This can be converted back into a scattering matrix using the matrix.convert() function.

The first half of the vector inputs and outputs to the scattering and transfer matrices represent fields propagating forward along the positive propagation axis (0) and the second half represents fields propagating backward along the negative direction.

Notation:
p:

positive propagation direction along propagation axis 0

n:

negative propagation direction along propagation axis 0

i:

inwards propagating (from source on either side)

o:

outwards propagating (backscattered or transmitted)

Scattering matrix equation (in -> out):

[po] = [A, B] [pi]

[no] = [C, D] [ni]

Transfer matrix equation (top -> bottom):

[po] = [A - B inv(D) C, B inv(D)] [pi]

[ni] = [ - inv(D) C inv(D)] [no],

where inv(D) is the (regularized) inverse of D.

transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.Grid(shape=None, step=None, extent=None, first=None, center=None, last=None, include_last=False, ndim=None, flat=False, origin_at_center=True, center_at_index=True)[source]

Bases: Sequence

A class representing an immutable uniformly-spaced plaid Cartesian grid and its Fourier Transform.

See also MutableGrid

__init__(shape=None, step=None, extent=None, first=None, center=None, last=None, include_last=False, ndim=None, flat=False, origin_at_center=True, center_at_index=True)[source]

Construct an immutable Grid object.

Parameters:
  • shape – An integer vector array with the shape of the sampling grid.

  • step – A vector array with the spacing of the sampling grid.

  • extent – The extent of the sampling grid as shape * step

  • first – A vector array with the first element for each dimension. The first element is the smallest element if step is positive, and the largest when step is negative.

  • center – A vector array with the center element for each dimension. The center position in the grid is rounded to the next integer index unless center_at_index is set to False for that partical axis.

  • last – A vector array with the last element for each dimension. Unless include_last is set to True for the associated dimension, all but the last element is returned when calling self[axis].

  • include_last – A boolean vector array indicating whether the returned vectors, self[axis], should include the last element (True) or all-but-the-last (False)

  • ndim (Optional[int]) – A scalar integer indicating the number of dimensions of the sampling space.

  • flat (Union[bool, Sequence, ndarray]) – A boolean vector array indicating whether the returned vectors, self[axis], should be flattened (True) or returned as an open grid (False)

  • origin_at_center (Union[bool, Sequence, ndarray]) – A boolean vector array indicating whether the origin should be fft-shifted (True) or be ifftshifted to the front (False) of the returned vectors for self[axis].

  • center_at_index (Union[bool, Sequence, ndarray]) – A boolean vector array indicating whether the center of the grid should be rounded to an integer index for each dimension. If False and the shape has an even number of elements, the next index is used as the center, (self.shape / 2).astype(int).

static from_ranges(*ranges)[source]

Converts one or more ranges of numbers to a single Grid object representation. The ranges can be specified as separate parameters or as a tuple.

Parameters:

ranges (Union[int, float, complex, Sequence, ndarray]) – one or more ranges of uniformly spaced numbers.

Return type:

Grid

Returns:

A Grid object that represents the same ranges.

property ndim: int

The number of dimensions of the space this grid spans.

property shape: array

The number of sample points along each axis of the grid.

property step: ndarray

The sample spacing along each axis of the grid.

property center: ndarray

The central coordinate of the grid.

property center_at_index: array

Boolean vector indicating whether the central coordinate is aligned with a grid point when the number of points is even along the associated axis. This has no effect when the the number of sample points is odd.

property flat: array

Boolean vector indicating whether self[axis] returns flattened (raveled) vectors (True) or not (False).

property origin_at_center: array

Boolean vector indicating whether self[axis] returns ranges that are monotonous (True) or ifftshifted so that the central index is the first element of the sequence (False).

property as_flat: Grid

return: A new Grid object where all the ranges are 1d-vectors (flattened or raveled)

property as_non_flat: Grid

return: A new Grid object where all the ranges are 1d-vectors (flattened or raveled)

property as_origin_at_0: Grid

return: A new Grid object where all the ranges are ifftshifted so that the origin as at index 0.

property as_origin_at_center: Grid

return: A new Grid object where all the ranges have the origin at the center index, even when the number of elements is odd.

swapaxes(axes)[source]

Reverses the order of the specified axes.

Return type:

Grid

transpose(axes=None)[source]

Reverses the order of all axes.

Return type:

Grid

project(axes_to_keep=None, axes_to_remove=None)[source]

Removes all but the specified axes and reduces the dimensions to the number of specified axes.

Parameters:
Return type:

Grid

Returns:

A Grid object with ndim == len(axes) and shape == shape[axes].

property first: ndarray

return: A vector with the first element of each range

property extent: ndarray

The spatial extent of the sampling grid.

property size: int

The total number of sampling points as an integer scalar.

property dtype

The numeric data type for the coordinates.

property f: Grid

The equivalent frequency Grid.

property k: Grid

The equivalent k-space Grid.

__add__(term)[source]

Add a (scalar) offset to the Grid coordinates.

Return type:

Grid

__mul__(factor)[source]

Scales all ranges with a factor.

Parameters:

factor (Union[int, float, complex, Sequence, array]) – A scalar factor for all dimensions, or a vector of factors, one for each dimension.

Return type:

Grid

Returns:

A new scaled Grid object.

__matmul__(other)[source]

Determines the Grid spanning the tensor space, with ndim equal to the sum of both ndims.

Parameters:

other (Grid) – The Grid with the right-hand dimensions.

Return type:

Grid

Returns:

A new Grid with ndim == self.ndim + other.ndim.

__sub__(term)[source]

Subtract a (scalar) value from all Grid coordinates.

Return type:

Grid

__truediv__(denominator)[source]

Divide the grid coordinates by a value.

Parameters:

denominator (Union[int, float, complex, Sequence, ndarray]) – The denominator to divide by.

Return type:

Grid

Returns:

A new Grid with the divided coordinates.

__neg__()[source]

Invert the coordinate values and the direction of the axes.

__len__()[source]

The number of axes in this sampling grid. Or, the number of elements when this object is not multi-dimensional.

Return type:

int

__iter__()[source]
property immutable: Grid

Return a new immutable Grid object.

property mutable: MutableGrid

return: A new MutableGrid object.

__eq__(other)[source]

Compares two Grid objects.

Return type:

bool

property multidimensional: bool

Single-dimensional grids behave as Sequences, multi-dimensional behave as a Sequence of vectors.

classmethod __init_subclass__(*args, **kwargs)

This method is called when a class is subclassed.

The default implementation does nothing. It may be overridden to extend subclasses.

static __new__(cls, *args, **kwds)
count(value) integer -- return number of occurrences of value
index(value[, start[, stop]]) integer -- return first index of value.

Raises ValueError if the value is not present.

Supporting start and stop arguments is optional, but recommended.

Subpackages

Submodules

macromax.bound module

The module provides the abstract Bound to represent the boundary of the simulation, e.g. periodic, or gradually more absorbing. Specific boundaries are implemented as subclasses and can be used directly as the bound argument to macromax.solve() or macromax.Solution. The precludes the inclusion of boundaries in the material description. It is sufficient to leave some space for the boundaries.

class macromax.bound.Electric[source]

Bases: object

Mixin for Bound to indicate that the electric susceptibility is non-zero.

property background_permittivity: Complex

A complex scalar indicating the permittivity of the background.

property electric_susceptibility: ndarray
property permittivity: ndarray

The electric permittivity, epsilon, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

class macromax.bound.Magnetic[source]

Bases: object

Mixin for Bound to indicate that the magnetic susceptibility is non-zero.

property magnetic_susceptibility: ndarray
property permeability: ndarray

The magnetic permeability, mu, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

class macromax.bound.Bound(grid=None, thickness=0.0, background_permittivity=1.0)[source]

Bases: object

A base class to represent calculation-volume-boundaries. Use the subclasses for practical implementations.

__init__(grid=None, thickness=0.0, background_permittivity=1.0)[source]
Parameters:
  • grid (Union[Grid, Sequence, ndarray, None]) – The Grid to which to the boundaries will be applied.

  • thickness (Union[Real, Sequence, ndarray]) – The thickness as a scalar, vector, or 2d-array (axes x side). Broadcasting is used as necessary.

  • background_permittivity (Complex) – The background permittivity of the boundary (default: 1.0 for vacuum). This is only used when the absolute permittivity is requested.

property grid: Grid

The Cartesian grid that indicates the sample positions of this bound and the volume it encompasses.

property thickness: ndarray

The thickness as a 2D-array thickness[axis, front_back] in meters.

property background_permittivity: Complex

A complex scalar indicating the permittivity of the background.

property electric_susceptibility: ndarray

The electric susceptibility, chi_E, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property magnetic_susceptibility: ndarray

The magnetic susceptibility, chi_H, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property between: ndarray

Returns a boolean array indicating True for the voxels between the boundaries where the calculation happens. The inner edge is considered in between the boundaries.

property beyond: ndarray

Returns a boolean array indicating True for the voxels beyond the boundary edge, i.e. inside the boundaries, where the calculation should be ignored.

class macromax.bound.PeriodicBound(grid)[source]

Bases: Bound

__init__(grid)[source]

Constructs an object that represents periodic boundaries.

Parameters:

grid (Union[Grid, Sequence, ndarray]) – The Grid to which to the boundaries will be applied.

property background_permittivity: Complex

A complex scalar indicating the permittivity of the background.

property between: ndarray

Returns a boolean array indicating True for the voxels between the boundaries where the calculation happens. The inner edge is considered in between the boundaries.

property beyond: ndarray

Returns a boolean array indicating True for the voxels beyond the boundary edge, i.e. inside the boundaries, where the calculation should be ignored.

property electric_susceptibility: ndarray

The electric susceptibility, chi_E, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property grid: Grid

The Cartesian grid that indicates the sample positions of this bound and the volume it encompasses.

property magnetic_susceptibility: ndarray

The magnetic susceptibility, chi_H, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property thickness: ndarray

The thickness as a 2D-array thickness[axis, front_back] in meters.

class macromax.bound.AbsorbingBound(grid, thickness=0.0, extinction_coefficient_function=<function AbsorbingBound.<lambda>>, background_permittivity=1.0)[source]

Bases: Bound, Electric

__init__(grid, thickness=0.0, extinction_coefficient_function=<function AbsorbingBound.<lambda>>, background_permittivity=1.0)[source]

Constructs a boundary with depth-dependent extinction coefficient, kappa(rel_depth).

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – The Grid to which to the boundaries will be applied.

  • thickness (Union[Real, Sequence, ndarray]) – The boundary thickness(es) in meters. This can be specified as a 2d-array [axis, side]. Singleton dimensions are broadcast.

  • extinction_coefficient_function (Union[Callable, Sequence, ndarray]) – A function that returns the extinction coefficient as function of the depth in the boundary relative to the total thickness of the boundary.

  • background_permittivity (Complex) – (default: 1.0 for vacuum)

property is_electric: bool
property extinction: ndarray

Determines the extinction coefficient, kappa, of the boundary on a plaid grid. The only non-zero values are found in the boundaries. At the corners, the maximum extinction value of the overlapping dimensions is returned.

Note that the returned array may have singleton dimensions that must be broadcast!

Returns:

An nd-array with the extinction coefficient, kappa.

property electric_susceptibility: ndarray

The electric susceptibility, chi_E, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property background_permittivity: Complex

A complex scalar indicating the permittivity of the background.

property between: ndarray

Returns a boolean array indicating True for the voxels between the boundaries where the calculation happens. The inner edge is considered in between the boundaries.

property beyond: ndarray

Returns a boolean array indicating True for the voxels beyond the boundary edge, i.e. inside the boundaries, where the calculation should be ignored.

property grid: Grid

The Cartesian grid that indicates the sample positions of this bound and the volume it encompasses.

property magnetic_susceptibility: ndarray

The magnetic susceptibility, chi_H, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property permittivity: ndarray

The electric permittivity, epsilon, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property thickness: ndarray

The thickness as a 2D-array thickness[axis, front_back] in meters.

class macromax.bound.LinearBound(grid, thickness=0.0, max_extinction_coefficient=0.25, background_permittivity=1.0)[source]

Bases: AbsorbingBound

__init__(grid, thickness=0.0, max_extinction_coefficient=0.25, background_permittivity=1.0)[source]

Constructs a boundary with linearly increasing extinction coefficient, kappa.

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – The Grid to which to the boundaries will be applied.

  • thickness (Union[Real, Sequence, ndarray]) – The boundary thickness(es) in meters. This can be specified as a 2d-array [axis, side]. Singleton dimensions are broadcast.

  • max_extinction_coefficient (Union[Real, Sequence, ndarray]) – The maximum extinction coefficient, reached at the deepest point of the boundary at the edge of the calculation volume.

  • background_permittivity (Complex) – (default: 1.0 for vacuum)

property background_permittivity: Complex

A complex scalar indicating the permittivity of the background.

property between: ndarray

Returns a boolean array indicating True for the voxels between the boundaries where the calculation happens. The inner edge is considered in between the boundaries.

property beyond: ndarray

Returns a boolean array indicating True for the voxels beyond the boundary edge, i.e. inside the boundaries, where the calculation should be ignored.

property electric_susceptibility: ndarray

The electric susceptibility, chi_E, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property extinction: ndarray

Determines the extinction coefficient, kappa, of the boundary on a plaid grid. The only non-zero values are found in the boundaries. At the corners, the maximum extinction value of the overlapping dimensions is returned.

Note that the returned array may have singleton dimensions that must be broadcast!

Returns:

An nd-array with the extinction coefficient, kappa.

property grid: Grid

The Cartesian grid that indicates the sample positions of this bound and the volume it encompasses.

property is_electric: bool
property magnetic_susceptibility: ndarray

The magnetic susceptibility, chi_H, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property permittivity: ndarray

The electric permittivity, epsilon, at every sample point. Note that the returned array may have singleton dimensions that must be broadcast!

property thickness: ndarray

The thickness as a 2D-array thickness[axis, front_back] in meters.

macromax.matrix module

class macromax.matrix.CachingMatrix(caching=True)[source]

Bases: object

__init__(caching=True)[source]

A mixin for Matrices that can cache the output.

Parameters:

caching (bool) – Cache field propagation calculations. By default, the results are cached for multiplications with basis vectors. Numerical errors might accumulate for certain superpositions. Setting this property to Setting this to False disables the cache so that field propagations are always used and the constructor argument array is ignored.

property caching: bool

When set to True, this object uses cached values instead of propagating the field through the scatterer. Otherwise, field propagation is used for all matrix operations. This can help avoid the accumulation of numerical errors.

class macromax.matrix.Matrix(array=None, shape=None, dtype=<class 'numpy.complex128'>)[source]

Bases: LinearOperator

A class to represent rectangular or square matrices that can be multiplied from the left or right, and pseudo-inverted.

__init__(array=None, shape=None, dtype=<class 'numpy.complex128'>)[source]

Constructs a matrix from a rectangular numpy.ndarray, array-like object, or a function or method that returns one.

Parameters:
  • array (Union[Complex, Sequence, ndarray, LinearOperator, None]) – A sequence or numpy.ndarray of complex numbers representing the matrix, or a function that returns one. If not specified, shape must be specified and _matmul and _rmatmul can be implemented in a subclass.

  • shape (Optional[Sequence[int]]) – The optional shape of the matrix, i.e. the number of rows and columns.

  • dtype – The optional dtype of the matrix elements.

__len__()[source]

The number of rows in the matrix as an integer.

Return type:

int

__setitem__(key, value)[source]

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__array__()[source]

The array values represented by this matrix.

Return type:

ndarray

inv(noise_level=0.0)[source]

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__call__(x)

Call self as a function.

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.SquareMatrix(array=None, side=None, dtype=<class 'numpy.complex128'>)[source]

Bases: Matrix

A class to represent square matrices that can be inverted with or without regularization.

__init__(array=None, side=None, dtype=<class 'numpy.complex128'>)[source]

Constructs a matrix from a square array, array-like object, or a function or method that returns one.

Parameters:
  • array (Union[Complex, Sequence, ndarray, LinearOperator, None]) – A sequence or numpy.ndarray of complex numbers representing the matrix, or a function that returns one. If not specified, the identity matrix is assumed.

  • side (Optional[int]) – The side of the matrix, i.e. the number of rows or columns.

  • dtype – The optional dtype of the matrix elements.

property side: int
property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.LiteralScatteringMatrix(array=None, side=None, dtype=<class 'numpy.complex128'>)[source]

Bases: SquareMatrix

A class to represent scattering matrices constructed from an array of complex numbers.

__init__(array=None, side=None, dtype=<class 'numpy.complex128'>)[source]

Constructs a scattering matrix from a square array, array-like object, or a function or method that returns one.

Parameters:
  • array (Union[Complex, Sequence, ndarray, LinearOperator, None]) – A sequence or numpy.ndarray of complex numbers representing the matrix, or a function that returns one. If not specified, the identity matrix is assumed.

  • side (Optional[int]) – The side of the matrix, i.e. the number of rows or columns.

  • dtype – The optional dtype of the matrix elements.

transfer(noise_level=0.0)[source]

Calculates the transfer matrix, relating one side of the scatterer to the other side (top, bottom). Each side can have incoming and outgoing waves. This is in contrast to the scattering matrix, self.__array__, which relates incoming waves from both sides to outgoing waves from both sides. One can be calculated from the other using the matrix.convert() function, though this calculation may be ill-conditioned (sensitive to noise). Therefore, the optional argument noise_level should be used to indicate the root-mean-square expectation value of the measurement error. This avoids divisions by near-zero values and obtains a best estimate using Tikhonov regularization.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near) singular backwards transmission matrix.

Return type:

ndarray

Returns:

An nd-array with the transfer matrix relating top-to-bottom instead of in-to-out. This can be converted back into a scattering matrix using the matrix.convert() function.

The first half of the vector inputs and outputs to the scattering and transfer matrices represent fields propagating forward along the positive propagation axis (0) and the second half represents fields propagating backward along the negative direction.

Notation:
p:

positive propagation direction along propagation axis 0

n:

negative propagation direction along propagation axis 0

i:

inwards propagating (from source on either side)

o:

outwards propagating (backscattered or transmitted)

Scattering matrix equation (in -> out):

[po] = [A, B] [pi]

[no] = [C, D] [ni]

Transfer matrix equation (top -> bottom):

[po] = [A - B inv(D) C, B inv(D)] [pi]

[ni] = [ - inv(D) C inv(D)] [no],

where inv(D) is the (regularized) inverse of D.

property forward_transmission: ForwardTransmissionMatrix

Select the forward-transmitted quarter of the scattering matrix. It indicates how the light coming from negative infinity is transmitted to positive infinity.

Returns:

The forward-transmission matrix of shape self.shape // 2.

property front_reflection: FrontReflectionMatrix

Select the quarter of the scattering matrix corresponding to the light that is reflected of the front. It indicates how the light coming from negative infinity is back reflected to negative infinity.

Returns:

The front-reflection matrix of shape self.shape // 2.

property back_reflection: BackReflectionMatrix

Select the quarter of the scattering matrix corresponding to the light that is reflected of the back. It indicates how the light coming from positive infinity is back reflected to positive infinity.

Returns:

The back-reflection matrix of shape self.shape // 2.

property backward_transmission: BackwardTransmissionMatrix

Select the backward-transmitted quarter of the scattering matrix. It indicates how the light coming from positive infinity is transmitted to negative infinity.

Returns:

The backward-transmission matrix of shape self.shape // 2.

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.ScatteringMatrix(grid, vectorial=True, wavenumber=None, angular_frequency=None, vacuum_wavelength=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, dtype=None, callback=<function ScatteringMatrix.<lambda>>, caching=True, array=None)[source]

Bases: LiteralScatteringMatrix

A class representing scattering matrices.

__init__(grid, vectorial=True, wavenumber=None, angular_frequency=None, vacuum_wavelength=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, dtype=None, callback=<function ScatteringMatrix.<lambda>>, caching=True, array=None)[source]

Construct a scattering matrix object for a medium specified by a refractive index distribution or the corresponding epsilon, xi, zeta, and mu distributions. Each electromagnetic field distribution entering the material is scattered into a certain electromagnetic field distributions propagating away from it from both sides. The complex matrix relates the amplitude and phase of all N propagating input modes to all N propagating output modes. No scattering, as in vacuum, is indicated by the NxN identity matrix. There are N/2 input and output modes on either side of the scattering material. Mode i and i+N/2 correspond to plane wave traveling in opposing directions. The mode directions are taken in raster-scan order and only propagating modes are included. When polarization is considered, the modes come in pairs corresponding to two orthogonal linear polarizations.

The modes are encoded as a vector of length N = 2xMxP for 2 sides, M angles, and P polarizations.

  • First the N/2 modes propagating along the positive x-axis are considered, then those propagating in the reverse direction.

  • In each direction, M different angles (k-vectors) can be considered. We choose propagating modes on a uniformly-spaced plaid grid that includes the origin (corresponding to the k-vector along the x-axis). Modes not propagating along the x-axis, i.e. in the y-z-plane are not considered. The angles are ordered in raster-scan order from negative k_y to positive k_y (slow) and from negative k_z to positive k_z (fast). The grid axes dimensions correspond to x(0), y(1), z(2).

  • When polarization is considered, each angle has a pair of modes, one for each polarization. The first mode has the polarization oriented along the rotated y’-axis and the second mode along the rotated z’-axis. To avoid ambiguity for normal incidence, the Cartesian-coordinate system is rotated along the shortest possible path, i.e. along the axis that is normal to the original x-axis and the mode’s k-vector. All rotations are around the origin of the coordinate system, incurring no phase shift there.

Vectors can be converted to field distributions on the complete grid using the methods:

srcvec2freespace() super-position of free-space plane waves in the whole volume (fast)

srcvec2source() super-position of free-space plane waves at the source planes at the front and back (fast)

source2detfield() calculate the field in the whole volume using the solver.Solution object (slow)

detfield2detvec() vector corresponding to the detected field at the detection planes (fast). The fields

at those planes should only contain the outward propagating waves. Hence, inwards propagating waves should be subtracted before using this method!

srcvec2detfield() calculate the field in the whole volume and convert it to a detection vector (slow)

The latter is used in the matrix multiplication method: matmul, @

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – A Grid object or a Sequence of vectors with uniformly increasing values that indicate the positions in a plaid grid of sample points for the material and solution. In the one-dimensional case, a simple increasing Sequence of uniformly-spaced numbers may be provided as an alternative. The length of the ranges determines the data_shape, to which the source_distribution, epsilon, xi, zeta, mu, and initial_field must broadcast when specified as ndarrays.

  • vectorial (Optional[bool]) – a boolean indicating if the source and solution are 3-vectors-fields (True) or scalar fields (False).

  • wavenumber (Optional[Real]) – the wavenumber in vacuum = 2 pi / vacuum_wavelength. The wavelength in the same units as used for the other inputs/outputs.

  • angular_frequency (Optional[Real]) – alternative argument to the wavenumber = angular_frequency / c

  • vacuum_wavelength (Optional[Real]) – alternative argument to the wavenumber = 2 pi / vacuum_wavelength

  • epsilon (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) epsilon that represents the permittivity at the points indicated by the grid specified as its input arguments.

  • xi (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) xi for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • zeta (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) zeta for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • mu (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) permeability at the points indicated by the grid specified as its input arguments.

  • refractive_index (Union[Complex, Sequence, ndarray, LinearOperator, None]) – an array or function that returns the (tensor) refractive_index = np.sqrt(permittivity) at the points indicated by the grid input argument.

  • bound (Optional[Bound]) – An object representing the boundary of the calculation volume. Default: None, PeriodicBound(grid)

  • dtype – optional numpy datatype for the internal operations and results. This must be a complex number type as numpy.complex128 or np.complex64.

  • callback (Callable) – optional function that will be called with as argument this solver. This function can be used to check and display progress. It must return a boolean value of True to indicate that further iterations are required.

  • caching (bool) – Cache field propagation calculations. By default, the results are cached for multiplications with basis vectors. Numerical errors might accumulate for certain superpositions. Setting this property to False will ensure that field propagations are always used and the constructor argument array is ignored.

  • array (Union[Complex, Sequence, ndarray, LinearOperator, None]) – Optional in case the matrix values have been calculated before and stored. If not specified, the matrix is calculated from the material properties. I specified, this must be a sequence or numpy.ndarray of complex numbers representing the matrix, or a function that returns one.

Returns:

The Solution object that has the E and H fields, as well as iteration information.

property grid: Grid

The calculation grid.

property vectorial: bool

Boolean to indicates whether calculations happen on polarized (True) or scalar (False) fields.

property caching: bool

When set to True, this object uses cached values instead of propagating the field through the scatterer. Otherwise, field propagation is used for all matrix operations. This can help avoid the accumulation of numerical errors.

srcvec2freespace(input_vector)[source]

Convert an input source vector to a superposition of plane waves at the origin and spreading over the whole volume. The input vector specifies the propagating modes in the far-field (the inverse Fourier transform of the fields at the sample origin). Incident waves at an angle will result in higher amplitudes to compensate for the reduction in propagation along the propagation axis through the entrance plane.

Used in ScatteringMatrix.srcvec2source to calculate the source field distribution before entering the scatterer.

Used in ScatteringMatrix.__matmul__ to distinguish incoming from back-scattered light.

Parameters:

input_vector (Union[Complex, Sequence, ndarray, LinearOperator]) – A source vector or array of shape [2, M, P], where the first axis indicates the side (front, back), the second axis indicates the propagation mode (direction, top-bottom-left-right), and the final axis indicates the polarization (1 for scalar, 2 for polarized: V-H).

Return type:

ndarray

Returns:

An nd-array with the field on the calculation grid. Its shape is (1, *self.grid.shape) for scalar

calculations and (3, *self.grid.shape) for vectorial calculations with polarization.

srcvec2source(input_vector, out=None)[source]

Converts a source vector into an (N+1)D-array with the source field at the front and back of the scatterer. The source field is such that it produces E-fields of unit intensity for unit vector inputs.

Used in self.vector2field() and self.__matmul__().

Parameters:
  • input_vector (Union[Complex, Sequence, ndarray, LinearOperator]) – A source vector with self.shape[1] elements. One value per side, per independent polarization (2), and per mode (inwards propagating k-vectors only).

  • out (Optional[ndarray]) – (optional) numpy array to store the result.

Return type:

ndarray

Returns:

The field distribution as an array of shape [nb_pol, *self.grid], where nb_pol = 3 for a vectorial calculation and 1 for a scalar calculation.

source2detfield(source, out=None)[source]

Calculates the (N+1)D-input-field distribution throughout the scatterer for a given source field distribution.

Parameters:
  • source (Union[Complex, Sequence, ndarray, LinearOperator]) – The source field distribution in the whole space.

  • out (Optional[ndarray]) – (optional) numpy array to store the result (shape: self.grid.shape, dtype: self.dtype).

Return type:

ndarray

Returns:

The field distribution as an array of shape [nb_pol, *self.grid], where nb_pol = 3 for a vectorial calculation and 1 for a scalar calculation.

srcvec2detfield(input_vector, out=None)[source]

Calculates the (N+1)D-input-field distribution throughout the scatterer for a given input source vector.

Used in self.__matmul__() and external code.

Parameters:
Return type:

ndarray

Returns:

The field distribution as an array of shape [nb_pol, *self.grid], where nb_pol = 3 for a vectorial calculation and 1 for a scalar calculation.

detfield2detvec(field)[source]

Converts the (N+1)D-output-field defined at the front and back detection planes of the scatterer to a detection vector that describes the resulting far-field distribution (the inverse Fourier transform of the fields at the sample origin). The fields at those planes should only contain the outward propagating waves. Hence, inwards propagating waves should be subtracted before using this method!

This is used in self.__matmul__().

Parameters:

field (Union[Complex, Sequence, ndarray, LinearOperator]) – The detected field in all space (of which only the detection space is used).

Return type:

ndarray

Returns:

Detection vector.

detvec2srcvec(vec)[source]

Convert forward propagating detection vector into a backward propagating (time-reversed) source vector.

Parameters:

vec (Union[Complex, Sequence, ndarray, LinearOperator]) – detection vector obtained from solution or scattering matrix multiplication.

Return type:

ndarray

Returns:

Time-reversed vector, that can be used as a source vector.

__setitem__(key, value)[source]

Updating this matrix is not possible. Use a Matrix object instead.

Parameters:
  • key – Index or slice.

  • value – The new value.

__array__(out=None)[source]

Lazily calculates the scattering matrix as a regular numpy.ndarray

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property back_reflection: BackReflectionMatrix

Select the quarter of the scattering matrix corresponding to the light that is reflected of the back. It indicates how the light coming from positive infinity is back reflected to positive infinity.

Returns:

The back-reflection matrix of shape self.shape // 2.

property backward_transmission: BackwardTransmissionMatrix

Select the backward-transmitted quarter of the scattering matrix. It indicates how the light coming from positive infinity is transmitted to negative infinity.

Returns:

The backward-transmission matrix of shape self.shape // 2.

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

property forward_transmission: ForwardTransmissionMatrix

Select the forward-transmitted quarter of the scattering matrix. It indicates how the light coming from negative infinity is transmitted to positive infinity.

Returns:

The forward-transmission matrix of shape self.shape // 2.

property front_reflection: FrontReflectionMatrix

Select the quarter of the scattering matrix corresponding to the light that is reflected of the front. It indicates how the light coming from negative infinity is back reflected to negative infinity.

Returns:

The front-reflection matrix of shape self.shape // 2.

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transfer(noise_level=0.0)

Calculates the transfer matrix, relating one side of the scatterer to the other side (top, bottom). Each side can have incoming and outgoing waves. This is in contrast to the scattering matrix, self.__array__, which relates incoming waves from both sides to outgoing waves from both sides. One can be calculated from the other using the matrix.convert() function, though this calculation may be ill-conditioned (sensitive to noise). Therefore, the optional argument noise_level should be used to indicate the root-mean-square expectation value of the measurement error. This avoids divisions by near-zero values and obtains a best estimate using Tikhonov regularization.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near) singular backwards transmission matrix.

Return type:

ndarray

Returns:

An nd-array with the transfer matrix relating top-to-bottom instead of in-to-out. This can be converted back into a scattering matrix using the matrix.convert() function.

The first half of the vector inputs and outputs to the scattering and transfer matrices represent fields propagating forward along the positive propagation axis (0) and the second half represents fields propagating backward along the negative direction.

Notation:
p:

positive propagation direction along propagation axis 0

n:

negative propagation direction along propagation axis 0

i:

inwards propagating (from source on either side)

o:

outwards propagating (backscattered or transmitted)

Scattering matrix equation (in -> out):

[po] = [A, B] [pi]

[no] = [C, D] [ni]

Transfer matrix equation (top -> bottom):

[po] = [A - B inv(D) C, B inv(D)] [pi]

[ni] = [ - inv(D) C inv(D)] [no],

where inv(D) is the (regularized) inverse of D.

transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.QuarterMatrix(matrix, backwards_output, backwards_input)[source]

Bases: SquareMatrix

A base class representing a quarter of a scattering matrix.

__init__(matrix, backwards_output, backwards_input)[source]

Construct an object refering to a quarter of a scattering matrix. Any Scattering matrix should have an even number of rows and columns.

Parameters:
  • matrix (SquareMatrix) – The underlying scattering matrix.

  • backwards_output (bool) – When True, select the bottom half of the matrix, i.e. the quadrants corresponding to the output modes exiting the front side of the scatterer (back-to-front propagation).

  • backwards_input (bool) – When True, select the right half of the matrix, i.e. the quadrants corresponding to the input modes entering from the back side of the scatterer (back-to-front propagation).

property full_matrix: SquareMatrix

The underlying (scattering) matrix of size 2 * self.shape.

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.ForwardTransmissionMatrix(scattering_matrix)[source]

Bases: QuarterMatrix

The forward transmission matrix of the specified scattering matrix. It indicates how the light coming from negative infinity is transmitted to positive infinity.

__init__(scattering_matrix)[source]

Construct an object refering to a quarter of a scattering matrix. Any Scattering matrix should have an even number of rows and columns.

Parameters:
  • matrix – The underlying scattering matrix.

  • backwards_output – When True, select the bottom half of the matrix, i.e. the quadrants corresponding to the output modes exiting the front side of the scatterer (back-to-front propagation).

  • backwards_input – When True, select the right half of the matrix, i.e. the quadrants corresponding to the input modes entering from the back side of the scatterer (back-to-front propagation).

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

property full_matrix: SquareMatrix

The underlying (scattering) matrix of size 2 * self.shape.

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.FrontReflectionMatrix(scattering_matrix)[source]

Bases: QuarterMatrix

The forward reflection matrix of the specified scattering matrix. It indicates how the light coming from positive infinity is back reflected to positive infinity.

__init__(scattering_matrix)[source]

Construct an object refering to a quarter of a scattering matrix. Any Scattering matrix should have an even number of rows and columns.

Parameters:
  • matrix – The underlying scattering matrix.

  • backwards_output – When True, select the bottom half of the matrix, i.e. the quadrants corresponding to the output modes exiting the front side of the scatterer (back-to-front propagation).

  • backwards_input – When True, select the right half of the matrix, i.e. the quadrants corresponding to the input modes entering from the back side of the scatterer (back-to-front propagation).

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

property full_matrix: SquareMatrix

The underlying (scattering) matrix of size 2 * self.shape.

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.BackReflectionMatrix(scattering_matrix)[source]

Bases: QuarterMatrix

The backward reflection matrix of the specified scattering matrix. It indicates how the light coming from negative infinity is back reflected to negative infinity.

__init__(scattering_matrix)[source]

Construct an object refering to a quarter of a scattering matrix. Any Scattering matrix should have an even number of rows and columns.

Parameters:
  • matrix – The underlying scattering matrix.

  • backwards_output – When True, select the bottom half of the matrix, i.e. the quadrants corresponding to the output modes exiting the front side of the scatterer (back-to-front propagation).

  • backwards_input – When True, select the right half of the matrix, i.e. the quadrants corresponding to the input modes entering from the back side of the scatterer (back-to-front propagation).

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

property full_matrix: SquareMatrix

The underlying (scattering) matrix of size 2 * self.shape.

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.BackwardTransmissionMatrix(scattering_matrix)[source]

Bases: QuarterMatrix

The backward transmission matrix of the specified scattering matrix. It indicates how the light coming from positive infinity is transmitted to negative infinity.

__init__(scattering_matrix)[source]

Construct an object refering to a quarter of a scattering matrix. Any Scattering matrix should have an even number of rows and columns.

Parameters:
  • matrix – The underlying scattering matrix.

  • backwards_output – When True, select the bottom half of the matrix, i.e. the quadrants corresponding to the output modes exiting the front side of the scatterer (back-to-front propagation).

  • backwards_input – When True, select the right half of the matrix, i.e. the quadrants corresponding to the input modes entering from the back side of the scatterer (back-to-front propagation).

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

property full_matrix: SquareMatrix

The underlying (scattering) matrix of size 2 * self.shape.

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

property side: int
transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

class macromax.matrix.DepositionMatrix(scattering_matrix, input_operator=None, output_operator=None, caching=True)[source]

Bases: Matrix, CachingMatrix

A rectangular matrix that relates free-space input vectors to an arbitrary subset of fields in the interior.

__init__(scattering_matrix, input_operator=None, output_operator=None, caching=True)[source]

Creates a matrix based on the internally scattered fields of a ScatteringMatrix.

Parameters:
  • scattering_matrix (ScatteringMatrix) – The base ScatteringMatrix

  • input_operator (Optional[LinearOperator]) – The optional field projector, a LinearOperator object that multiplies any input vector to produce a linear combination of source fields. Default: the srcvec2source method of the ScatteringMatrix.

  • output_operator (Union[LinearOperator, Callable[[Union[Complex, Sequence, ndarray, LinearOperator]], ndarray], None]) – The optional detection field projector, a LinearOperator object that multiplies the raveled full-field distribution and returns a projection vector or a function that does the same. Default: the detfield2detvec method of the ScatteringMatrix.

  • caching (bool) – Cache field propagation calculations. By default, the results are cached for multiplications with basis vectors. Numerical errors might accumulate for certain superpositions. Setting this property to False will ensure that field propagations are always used and the constructor argument array is ignored.

property H

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property T

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

__add__(x)
__array__()

The array values represented by this matrix.

Return type:

ndarray

__call__(x)

Call self as a function.

__len__()

The number of rows in the matrix as an integer.

Return type:

int

__matmul__(other)
__mul__(x)
__neg__()
static __new__(cls, *args, **kwargs)
__pow__(p)
__rmatmul__(other)
__rmul__(x)
__setitem__(key, value)

Update (part of) the matrix.

Parameters:
  • key – Index or slice.

  • value – The new value.

__sub__(x)
adjoint()

Hermitian adjoint.

Returns the Hermitian adjoint of self, aka the Hermitian conjugate or Hermitian transpose. For a complex matrix, the Hermitian adjoint is equal to the conjugate transpose.

Can be abbreviated self.H instead of self.adjoint().

Returns:

A_H – Hermitian adjoint of self.

Return type:

LinearOperator

property caching: bool

When set to True, this object uses cached values instead of propagating the field through the scatterer. Otherwise, field propagation is used for all matrix operations. This can help avoid the accumulation of numerical errors.

dot(x)

Matrix-matrix or matrix-vector multiplication.

Parameters:

x (array_like) – 1-d or 2-d array, representing a vector or matrix.

Returns:

Ax – 1-d or 2-d array (depending on the shape of x) that represents the result of applying this linear operator on x.

Return type:

array

inv(noise_level=0.0)

The (Tikhonov regularized) inverse of the scattering matrix.

Parameters:

noise_level (float) – (optional) argument to regularize the inversion of a (near-)singular matrix.

Return type:

ndarray

Returns:

An nd-array with the inverted matrix so that self @ self.inv approximates the identity.

matmat(X)

Matrix-matrix multiplication.

Performs the operation y=A*X where A is an MxN linear operator and X dense N*K matrix or ndarray.

Parameters:

X ({matrix, ndarray}) – An array with shape (N,K).

Returns:

Y – A matrix or ndarray with shape (M,K) depending on the type of the X argument.

Return type:

{matrix, ndarray}

Notes

This matmat wraps any user-specified matmat routine or overridden _matmat method to ensure that y has the correct type.

matvec(x)

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (N,) or (N,1).

Returns:

y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

ndim = 2
rmatmat(X)

Adjoint matrix-matrix multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array, or 2-d array. The default implementation defers to the adjoint.

Parameters:

X ({matrix, ndarray}) – A matrix or 2D array.

Returns:

Y – A matrix or 2D array depending on the type of the input.

Return type:

{matrix, ndarray}

Notes

This rmatmat wraps the user-specified rmatmat routine.

rmatvec(x)

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters:

x ({matrix, ndarray}) – An array with shape (M,) or (M,1).

Returns:

y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.

Return type:

{matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

transpose()

Transpose this linear operator.

Returns a LinearOperator that represents the transpose of this one. Can be abbreviated self.T instead of self.transpose().

macromax.matrix.inv(mat, noise_level=0.0)[source]

ScatteringMatrix inversion with optional Tikhonov regularization.

Return type:

ndarray

macromax.matrix.convert(s, noise_level=0.0)[source]

Converts a scattering matrix into a transfer matrix and vice-versa.

Notation:

p: positive propagation direction along propagation axis 0, n: negative propagation direction along propagation axis 0, i: inwards propagating (from source on either side), o: outwards propagating (backscattered or transmitted).

Scattering matrix equation (in -> out):

[po] = [A, B] [pi], [no] [C, D] [ni].

Transfer matrix equation (top -> bottom):

[po] = [A - B inv(D) C, B inv(D)] [pi], [ni] = [ - inv(D) C inv(D)] [no],

where inv(D) is the (regularized) inverse of D.

The first half of the vector inputs and outputs to the scattering and transfer matrices represent fields propagating forward along the positive propagation axis (0) and the second half represents fields propagating backward along the negative direction.

Parameters:
  • s (Union[Complex, Sequence, ndarray, LinearOperator]) – The scattering (transfer) matrix as a 2D np.ndarray of shape [N, N]. Alternatively, each dimension can be split in two halves for downward and upward propagating fields. In that case the shape would be [2, N//2, 2, N//2].

  • noise_level (float) – The noise level of the measured matrix elements (singular values). If greater than 0, (Tikhonov) regularization will be used.

Return type:

ndarray

Returns:

The transfer (scattering) matrix of the same shape as the scattering matrix.

macromax.solver module

This module calculates the solution to the wave equations. More specifically, the work is done in the iteration defined in the Solution.__iter__() method of the Solution class. The convenience function solve() is provided to construct a Solution object and iterate it to convergence using its Solution.solve() method.

macromax.solver.solve(grid, vectorial=None, wavenumber=1.0, angular_frequency=None, vacuum_wavelength=None, current_density=None, source_distribution=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, initial_field=0.0, dtype=None, callback=<function <lambda>>)[source]

Function to find a solution for Maxwell’s equations in a media specified by the epsilon, xi, zeta, and mu distributions in the presence of a current source.

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – A Grid object or a Sequence of vectors with uniformly increasing values that indicate the positions in a plaid grid of sample points for the material and solution. In the one-dimensional case, a simple increasing Sequence of uniformly-spaced numbers may be provided as an alternative. The length of the ranges determines the data_shape, to which the source_distribution, epsilon, xi, zeta, mu, and initial_field must broadcast when specified as ndarrays.

  • vectorial (Optional[bool]) – a boolean indicating if the source and solution are 3-vectors-fields (True) or scalar fields (False).

  • wavenumber (Optional[Real]) – the wavenumber in vacuum = 2 pi / vacuum_wavelength. The wavelength in the same units as used for the other inputs/outputs.

  • angular_frequency (Optional[Real]) – alternative argument to the wavenumber = angular_frequency / c

  • vacuum_wavelength (Optional[Real]) – alternative argument to the wavenumber = 2 pi / vacuum_wavelength

  • current_density (Union[Complex, Sequence, ndarray, None]) – (optional, instead of source_distribution) An array or function that returns the free (vectorial) current density input distribution, J. The free current density has units of \(A m^-2\).

  • source_distribution (Union[Complex, Sequence, ndarray, None]) – (optional, instead of current_density) An array or function that returns the (vectorial) source input wave distribution. The source values relate to the current density, J, as 1j * angular_frequency * scipy.constants.mu_0 * J and has units of \(rad s^-1 H m^-1 A m^-2 = rad V m^-3\). More general, non-electro-magnetic wave problems can be solved using the source_distribution, as it does not rely on the vacuum permeability constant, \(mu_0\).

  • epsilon (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (tensor) epsilon that represents the permittivity at the points indicated by the grid specified as its input arguments.

  • xi (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) xi for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • zeta (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) zeta for bi-(an)isotropy at the points indicated by the grid specified as its input arguments.

  • mu (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) permeability at the points indicated by the grid specified as its input arguments.

  • refractive_index (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (complex) (tensor) refractive_index, as the square root of the permittivity, at the points indicated by the grid input argument.

  • bound (Optional[Bound]) – An object representing the boundary of the calculation volume. Default: None, PeriodicBound(grid)

  • initial_field (Union[Complex, Sequence, ndarray]) – optional start value for the E-field distribution (default: all zero E)

  • dtype – optional numpy datatype for the internal operations and results. This must be a complex number type as numpy.complex128 or np.complex64.

  • callback (Callable) – optional function that will be called with as argument this solver. This function can be used to check and display progress. It must return a boolean value of True to indicate that further iterations are required.

Returns:

The Solution object that has the E and H fields, as well as iteration information.

class macromax.solver.Solution(grid, vectorial=None, wavenumber=1.0, angular_frequency=None, vacuum_wavelength=None, current_density=None, source_distribution=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, initial_field=0.0, dtype=None)[source]

Bases: object

__init__(grid, vectorial=None, wavenumber=1.0, angular_frequency=None, vacuum_wavelength=None, current_density=None, source_distribution=None, epsilon=None, xi=0.0, zeta=0.0, mu=1.0, refractive_index=None, bound=None, initial_field=0.0, dtype=None)[source]

Class a solution that can be further iterated towards a solution for Maxwell’s equations in a media specified by the epsilon, xi, zeta, and mu distributions.

Parameters:
  • grid (Union[Grid, Sequence, ndarray]) – A Grid object or a Sequence of vectors with uniformly increasing values that indicate the positions in a plaid grid of sample points for the material and solution. In the one-dimensional case, a simple increasing Sequence of uniformly-spaced numbers may be provided as an alternative. The length of the ranges determines the data_shape, to which the source_distribution, epsilon, xi, zeta, mu, and initial_field must broadcast when specified as `numpy.ndarray`s.

  • vectorial (Optional[bool]) – a boolean indicating if the source and solution are 3-vectors-fields (True) or scalar fields (False). Default: True, when vectorial nor the source is specified. Default: vectorial (True), unless the source field is scalar (False if first dimension is a singleton dimension).

  • wavenumber (Optional[Real]) – the wavenumber in vacuum = 2pi / vacuum_wavelength. The wavelength in the same units as used for the other inputs/outputs.

  • angular_frequency (Optional[Real]) – alternative argument to the wavenumber = angular_frequency / c

  • vacuum_wavelength (Optional[Real]) – alternative argument to the wavenumber = 2 pi / vacuum_wavelength

  • current_density (Union[Complex, Sequence, ndarray, None]) – (optional, instead of source_distribution) An array or function that returns the (vectorial) current density input distribution, J. The current density has units of \(A m^-2\).

  • source_distribution (Union[Complex, Sequence, ndarray, None]) – (optional, instead of current_density) An array or function that returns the (vectorial) source input wave distribution. The source values relate to the current density, J, as 1j * angular_frequency * scipy.constants.mu_0 * J and has units of \(rad s^-1 H m^-1 A m^-2 = rad V m^-3\). More general, non-electro-magnetic wave problems can be solved using the source_distribution, as it does not rely on the vacuum permeability constant, \(mu_0\).

  • epsilon (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (tensor) epsilon that represents the permittivity at the points indicated by the grid input argument.

  • xi (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) xi for bi-(an)isotropy at the points indicated by the grid input argument.

  • zeta (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) zeta for bi-(an)isotropy at the points indicated by the grid input argument.

  • mu (Union[Complex, Sequence, ndarray]) – an array or function that returns the (tensor) permeability at the points indicated by the grid input argument.

  • refractive_index (Union[Complex, Sequence, ndarray, None]) – an array or function that returns the (complex) (tensor) refractive_index, as the square root of the permittivity, at the points indicated by the grid input argument.

  • bound (Optional[Bound]) – An object representing the boundary of the calculation volume. Default: None, PeriodicBound(grid)

  • initial_field (Union[Complex, Sequence, ndarray]) – optional start value for the E-field distribution (default: all zero E)

  • dtype – optional numpy datatype for the internal operations and results. This must be a complex number type as numpy.complex128 or np.complex64.

property grid: Grid

The sample positions of the plaid sampling grid. This may be useful for displaying result axes.

Returns:

A Grid object representing the sample points of the fields and material.

property vectorial: bool

Boolean to indicates whether calculations happen on vectorial (True) or scalar (False) fields.

property dtype

The numpy equivalent data type used in the calculation. This is either np.complex64 or np.complex128.

property wavenumber: Real

The vacuum wavenumber, \(k_0\), used in the calculation.

Returns:

A scalar indicating the wavenumber used in the calculation.

property angular_frequency: Real

The angular frequency, \(\omega\), used in the calculation.

Returns:

A scalar indicating the angular frequency used in the calculation.

property wavelength: Real

The vacuum wavelength, \(\lambda_0\), used in the calculation.

Returns:

A scalar indicating the vacuum wavelength used in the calculation.

property magnetic: bool

Indicates if this media is considered magnetic.

Returns:

A boolean, True when magnetic, False otherwise.

property bound: Bound

The Bound object that defines the calculation boundaries.

property source_distribution: ndarray

The source distribution, i k0 mu_0 times the current density j.

Returns:

A complex array indicating the amplitude and phase of the source vector field. The dimensions of the array are [1|3, self.grid.shape], where the first dimension is 1 in case of a scalar field, and 3 in case of a vector field.

property j: ndarray

The free current density, j, of the source vector field.

Returns:

A complex array indicating the amplitude and phase of the current density vector field [A m^-2]. The dimensions of the array are [1|3, self.grid.shape], where the first dimension is 1 in case of a scalar field, and 3 in case of a vector field.

property E: ndarray

The electric field for every point in the sample space (SI units).

Returns:

A vector array with the first dimension containing Ex, Ey, and Ez,

while the following dimensions are the spatial dimensions.

property B: ndarray

The magnetic field for every point in the sample space (SI units). This is calculated from H and E.

Returns:

A vector array with the first dimension containing \(B_x, B_y, and B_z\), while the following dimensions are the spatial dimensions.

property D: ndarray

The displacement field for every point in the sample space (SI units). This is calculated from E and H.

Returns:

A vector array with the first dimension containing \(D_x, D_y, and D_z\), while the following dimensions are the spatial dimensions.

property H: ndarray

The magnetizing field for every point in the sample space (SI units). This is calculated from E.

Returns:

A vector array with the first dimension containing \(H_x, H_y, and H_z\), while the following dimensions are the spatial dimensions.

property S: ndarray

The time-averaged Poynting vector for every point in space. :return: A vector array with the first dimension containing \(S_x, S_y, and S_z\), while the following dimensions are the spatial dimensions.

property energy_density: ndarray

Returns the energy density, u.

Returns:

A real array indicating the energy density in space.

property stress_tensor: ndarray

Maxwell’s stress tensor for every point in space.

Returns:

A real and symmetric matrix-array with the stress tensor for every point in space. The units are \(N / m^2\).

property f: ndarray

The electromagnetic force density (force per SI unit volume, not per voxel).

Returns:

A vector array representing the electro-magnetic force exerted per unit volume. The first dimension contains \(f_x, f_y, and f_z\), while the following dimensions are the spatial dimensions. The units are \(N / m^3\).

property torque: ndarray

The electromagnetic force density (force per SI unit volume, not per voxel).

Returns:

A vector array representing the electro-magnetic torque exerted per unit volume. The first dimension contains torque_x, torque_y, and torque_z, while the following dimensions are the spatial dimensions. The units are \(N m / m^3 = N m^{-2}\).

property iteration: int

The current iteration number.

Returns:

An integer indicating how many iterations have been done.

property previous_update_norm: Real

The L2-norm of the last update, the difference between current and previous E-field.

Returns:

A positive scalar indicating the norm of the last update.

property residue: Real

Returns the current relative residue of the inverse problem \(E = H^{-1}S\). The relative residue is return as the l2-norm fraction \(||E - H^{-1}S|| / ||E||\), where H represents the vectorial Helmholtz equation following Maxwell’s equations and S the current density source. The solver searches for the electric field, E, that minimizes the preconditioned inverse problem.

Returns:

A non-negative real scalar that indicates the change in E with the previous iteration normalized to the norm of the current E.

__iter__()[source]

Returns an iterator that on __next__() yields this Solution after updating it with one cycle of the algorithm. Obtaining this iterator resets the iteration counter.

Usage:

for solution in Solution(...):
    if solution.iteration > 100:
        break
print(solution.residue)
solve(callback=<function Solution.<lambda>>)[source]

Runs the algorithm until the convergence criterion is met or until the maximum number of iterations is reached.

Parameters:

callback (Callable) – optional callback function that overrides the one set for the solver. E.g. callback=lambda s: s.iteration < 100

Returns:

This Solution object, which can be used to query e.g. the final field E using Solution.E.