finesse.knm.integrators module

A collection of methods to compute overlap integrals for modal scattering matrices. Essentially this involves solving the following integral

\[K_{abnm} = \iint^{\infty}_{-\infty} U_{nm}(x,y,q_i) M(x,y) U^*_{ab}(x,y,q_o) \, dy \, dx \]

\(U_{nm}\) is the initial modes in the basis \(q_i\) we are converting from and \(U_{ab}\) are the target modes in a basis \(q_o\) we are projecting into.

TODO

Should explore if decomposing compute_map_knm_matrix_riemann_optimised into real and imaginary integrals might be faster. In cases where q_in == q_out then integrals are real, apart from the map component which can be complex.

Explore use of zgemm3m which is 25% faster than zgemm

Probably look into using CUDA if necessary for more speed.

finesse.knm.integrators.composite_newton_cotes_weights(N, order)[source]

Constructs the weights for a composite Newton-Cotes rule for integration along 1-dimensional line with equally spaced steps. Newton-Cotes are a generalisation of a various discrete integration methods that are approximating an integrated by some polynomial. Common methods are:

N = 0 : Riemann sum

N = 1 : Trapezoid rule

N = 2 : Simpsons rule

N = 4 : Simpsons 3/8 rule

Approximating a large bound with a high order polynomial can be numerically problematic. Thus a composite rule is generated by subdividing a larger area into multiple chunks and applying each rule along it.

If the order Newton-Cotes order specified does not produce a rule that fits into N, the order is decremented and that is used to fill gaps that do no fit.

See https://mathworld.wolfram.com/Newton-CotesFormulas.html

Parameters

Ninteger >= 0: 1D size of array being ingrated over
orderinteger >= 0: Order of Newton-Cotes rule to use

Returns

weightsarray: Array of weights to multiply data with before summing

finesse.knm.integrators.compute_map_scattering_coeffs_riemann_optimised(double dA, double complex[, : :1] Z, double complex[, , : :1] Unn_, double complex[, , : :1] Umm_, double complex[, , : :1] tmp, double complex[, , , : :1] result)[source]

Calculates a mode scattering matrix using a Riemann sum. This method uses an computationally optimised approach making use of fast BLAS functions. This requires the input modes to be specified in specific formats and memory layouts. What this functions computes is the following via a Riemann sum:

\[K_{abnm} = \int^{\infty}_{-\infty} u_{m}(y,q^y_i) u^*_{b}(y,q^y_o) \Bigg[ \int^{\infty}_{-\infty} Z(x,y) u_{n}(x,q^x_i) u^*_{a}(x,q^x_o) \, dx \Bigg] \, dy \]

This integral is not actually performed to infinity, it is bound by the dimensions of the discretised map \(Z\). The map bound and uniform discretisation must be chosen to efficiently sample the size of any of the beams and maximum mode order being used.

Due to the optimised calculation method, the result indexing is not [a,b,n,m] but [m,a,n,b]. To convert it back to a more usable indexing use:

>>> result = np.transpose(result, (2,1,0,3))

Notes

Nx - number of x samples
Ny - number of y samples
Nm - number of modes (n, m) being calculated

Parameters

dAdouble: Area of discrete integral, dx * dy
Zarray[complex]: 2D map of size [Ny, Nx]
Unnarray[complex]: 3D array of size [Nm, Nm, Nx]. This should contain the \(U_n(x) * U_n'(x)^*\) products
Ummarray[complex]: 3D array of size [Nm, Nm, Ny]. This should contain the \(U_m(y) * U_m'(y)^*\) products
tmparray[complex]: Temporary storage that can be used to compute the dot products between Z and Unn. Should be of size (Nm, Nm, Ny)
resultarray[complex]: Resulting Knmnm output of size [Nm, Nm, Nm, Nm]. IMPORTANT: Note output indexing of result K[m,n,a,b], where (n,m) are the input mode indices and (a,b) are the output mode indices (m,b sagittal; n,a tangential).

finesse.knm.integrators.map_coupling_matrix_riemann(double complex[, : :1] Y, double dx, double dy, double complex[, : :1] Un, double complex[, : :1] Um, long[, : :1] index_map)[source]

finesse.knm.integrators.outer_conj_product(double complex[, : :1] U, double complex[, , : :1] result) → void: Computes U * U^C and returns the output into the result array. Result array must be of shape (N,N,M) where U is shape (N,M).

finesse.knm.integrators.outer_conj_product_2(double complex[, : :1] U1, double complex[, : :1] U2, double complex[, , : :1] result) → void: Computes U1 * U2**C and returns the output into the result array. Result array must be of shape (N,N,M) where U is shape (N,M).