diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index 22248f1..5fda683 100644 --- a/meanas/fdfd/waveguide_2d.py +++ b/meanas/fdfd/waveguide_2d.py @@ -18,8 +18,8 @@ $$ \begin{aligned} \nabla \times \vec{E}(x, y, z) &= -\imath \omega \mu \vec{H} \\ \nabla \times \vec{H}(x, y, z) &= \imath \omega \epsilon \vec{E} \\ -\vec{E}(x,y,z) &= (\vec{E}_t(x, y) + E_z(x, y)\vec{z}) e^{-\gamma z} \\ -\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\gamma z} \\ +\vec{E}(x,y,z) &= (\vec{E}_t(x, y) + E_z(x, y)\vec{z}) e^{-\imath \beta z} \\ +\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\imath \beta z} \\ \end{aligned} $$ @@ -40,56 +40,57 @@ Substituting in our expressions for $\vec{E}$, $\vec{H}$ and discretizing: $$ \begin{aligned} --\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \gamma E_y \\ --\imath \omega \mu_{yy} H_y &= -\gamma E_x - \tilde{\partial}_x E_z \\ +-\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta E_y \\ +-\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \tilde{\partial}_x E_z \\ -\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\ -\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \gamma H_y \\ -\imath \omega \epsilon_{yy} E_y &= -\gamma H_x - \hat{\partial}_x H_z \\ +\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta H_y \\ +\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \hat{\partial}_x H_z \\ \imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ \end{aligned} $$ Rewrite the last three equations as + $$ \begin{aligned} -\gamma H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ -\gamma H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ +\imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ +\imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ \imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ \end{aligned} $$ -Now apply $\gamma \tilde{\partial}_x$ to the last equation, -then substitute in for $\gamma H_x$ and $\gamma H_y$: +Now apply $\imath \beta \tilde{\partial}_x$ to the last equation, +then substitute in for $\imath \beta H_x$ and $\imath \beta H_y$: $$ \begin{aligned} -\gamma \tilde{\partial}_x \imath \omega E_z &= \gamma \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - - \gamma \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ +\imath \beta \tilde{\partial}_x \imath \omega E_z &= \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y + - \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z) - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\ &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x) - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y) \\ -\gamma \tilde{\partial}_x E_z &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) - + \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\ +\imath \beta \tilde{\partial}_x E_z &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + + \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\ \end{aligned} $$ -With a similar approach (but using $\gamma \tilde{\partial}_y$ instead), we can get +With a similar approach (but using $\imath \beta \tilde{\partial}_y$ instead), we can get $$ \begin{aligned} -\gamma \tilde{\partial}_y E_z &= \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) - + \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\ +\imath \beta \tilde{\partial}_y E_z &= \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + + \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\ \end{aligned} $$ -We can combine this equation for $\gamma \tilde{\partial}_y E_z$ with +We can combine this equation for $\imath \beta \tilde{\partial}_y E_z$ with the unused $\imath \omega \mu_{xx} H_x$ and $\imath \omega \mu_{yy} H_y$ equations to get $$ \begin{aligned} --\imath \omega \mu_{xx} \gamma H_x &= \gamma^2 E_y + \gamma \tilde{\partial}_y E_z \\ --\imath \omega \mu_{xx} \gamma H_x &= \gamma^2 E_y + \tilde{\partial}_y ( +-\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \imath \beta \tilde{\partial}_y E_z \\ +-\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \tilde{\partial}_y ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ @@ -100,25 +101,24 @@ and $$ \begin{aligned} --\imath \omega \mu_{yy} \gamma H_y &= -\gamma^2 E_x - \gamma \tilde{\partial}_x E_z \\ --\imath \omega \mu_{yy} \gamma H_y &= -\gamma^2 E_x - \tilde{\partial}_x ( +-\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \imath \beta \tilde{\partial}_x E_z \\ +-\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \tilde{\partial}_x ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ \end{aligned} $$ -However, based on our rewritten equation for $\gamma H_x$ and the so-far unused +However, based on our rewritten equation for $\imath \beta H_x$ and the so-far unused equation for $\imath \omega \mu_{zz} H_z$ we can also write $$ \begin{aligned} --\imath \omega \mu_{xx} (\gamma H_x) &= -\imath \omega \mu_{xx} (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\ - &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y - +\imath \omega \mu_{xx} \hat{\partial}_x ( - \frac{1}{-\imath \omega \mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x)) \\ - &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y - -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ +-\imath \omega \mu_{xx} (\imath \beta H_x) &= -\imath \omega \mu_{xx} (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\ + &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y + \imath \omega \mu_{xx} \hat{\partial}_x ( + \frac{1}{-\imath \omega \mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x)) \\ + &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y + -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ \end{aligned} $$ @@ -126,7 +126,7 @@ and, similarly, $$ \begin{aligned} --\imath \omega \mu_{yy} (\gamma H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +-\imath \omega \mu_{yy} (\imath \beta H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ \end{aligned} $$ @@ -135,12 +135,12 @@ By combining both pairs of expressions, we get $$ \begin{aligned} --\gamma^2 E_x - \tilde{\partial}_x ( +\beta^2 E_x - \tilde{\partial}_x ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) ) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ -\gamma^2 E_y + \tilde{\partial}_y ( +-\beta^2 E_y + \tilde{\partial}_y ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) ) &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y @@ -165,14 +165,13 @@ $$ E_y \end{bmatrix} $$ -where $\gamma = \imath\beta$. In the literature, $\beta$ is usually used to denote -the lossless/real part of the propagation constant, but in `meanas` it is allowed to -be complex. +In the literature, $\beta$ is usually used to denote the lossless/real part of the propagation constant, +but in `meanas` it is allowed to be complex. An equivalent eigenvalue problem can be formed using the $H_x$ and $H_y$ fields, if those are more convenient. -Note that $E_z$ was never discretized, so $\gamma$ and $\beta$ will need adjustment -to account for numerical dispersion if the result is introduced into a space with a discretized z-axis. +Note that $E_z$ was never discretized, so $\beta$ will need adjustment to account for numerical dispersion +if the result is introduced into a space with a discretized z-axis. """ @@ -531,10 +530,37 @@ def exy2e( dxes: dx_lists_t, epsilon: vfdfield_t, ) -> sparse.spmatrix: - """ + r""" Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, into a vectorized E containing all three E components + From the operator derivation (see module docs), we have + + $$ + \imath \omega \epsilon_{zz} E_z = \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ + $$ + + as well as the intermediate equations + + $$ + \begin{aligned} + \imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ + \imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ + \end{aligned} + $$ + + Combining these, we get + + $$ + \begin{aligned} + E_z &= \frac{1}{- \omega \beta \epsilon_{zz}} (( + \hat{\partial}_y \hat{\partial}_x H_z + -\hat{\partial}_x \hat{\partial}_y H_z) + + \imath \omega (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y)) + &= \frac{1}{\imath \beta \epsilon_{zz}} (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y) + \end{aligned} + $$ + Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy` @@ -905,7 +931,7 @@ def inner_product( # TODO documentation prop_phase: float = 0, conj_h: bool = False, trapezoid: bool = False, - ) -> tuple[vcfdfield_t, vcfdfield_t]: + ) -> complex: shape = [s.size for s in dxes[0]] @@ -920,7 +946,7 @@ def inner_product( # TODO documentation # Find time-averaged Sz and normalize to it dxes_real = [[numpy.real(dxyz) for dxyz in dxeh] for dxeh in dxes] - if integrate: + if trapezoid: Sz_a = numpy.trapezoid(numpy.trapezoid(E1[0] * H2[1], numpy.cumsum(dxes_real[0][1])), numpy.cumsum(dxes_real[1][0])) Sz_b = numpy.trapezoid(numpy.trapezoid(E1[1] * H2[0], numpy.cumsum(dxes_real[0][0])), numpy.cumsum(dxes_real[1][1])) else: diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index 227008d..2c00d02 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -1,14 +1,66 @@ -""" +r""" Operators and helper functions for cylindrical waveguides with unchanging cross-section. -WORK IN PROGRESS, CURRENTLY BROKEN +Waveguide operator is derived according to 10.1364/OL.33.001848. +The curl equations in the complex coordinate system become -As the z-dependence is known, all the functions in this file assume a 2D grid +$$ +\begin{aligned} +-\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta frac{E_y}{\tilde{t}_x} \\ +-\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \frac{1}{\hat{t}_x} \tilde{\partial}_x \tilde{t}_x E_z \\ +-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\ +\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\ +\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\ +\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ +\end{aligned} +$$ + +where $t_x = 1 + \frac{\Delta_{x, m}}{R_0}$ is the grid spacing adjusted by the nominal radius $R0$. + +Rewrite the last three equations as + +$$ +\begin{aligned} +\imath \beta H_y &= \imath \omega \hat{t}_x \epsilon_{xx} E_x - \hat{t}_x \hat{\partial}_y H_z \\ +\imath \beta H_x &= -\imath \omega \hat{t}_x \epsilon_{yy} E_y - \hat{t}_x \hat{\partial}_x H_z \\ +\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ +\end{aligned} +$$ + +The derivation then follows the same steps as the straight waveguide, leading to the eigenvalue problem + +$$ +\beta^2 \begin{bmatrix} E_x \\ + E_y \end{bmatrix} = + (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\ + 0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} + + \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\ + T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1} + \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + + \begin{bmatrix} \tilde{\partial}_x \\ + \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1} + \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix}) + \begin{bmatrix} E_x \\ + E_y \end{bmatrix} +$$ + +which resembles the straight waveguide eigenproblem with additonal $T_a$ and $T_b$ terms. These +are diagonal matrices containing the $t_x$ values: + +$$ +\begin{aligned} +T_a &= 1 + \frac{\Delta_{x, m }}{R_0} +T_b &= 1 + \frac{\Delta_{x, m + \frac{1}{2} }}{R_0} +\end{aligned} + + +TODO: consider 10.1364/OE.20.021583 for an alternate approach +$$ + +As in the straight waveguide case, all the functions in this file assume a 2D grid (i.e. `dxes = [[[dr_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`). """ -# TODO update module docs - -from typing import Any +from typing import Any, cast from collections.abc import Sequence import logging @@ -19,7 +71,7 @@ from scipy import sparse from ..fdmath import vec, unvec, dx_lists_t, vfdfield_t, vcfdfield_t from ..fdmath.operators import deriv_forward, deriv_back from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration - +from . import waveguide_2d logger = logging.getLogger(__name__) @@ -30,13 +82,21 @@ def cylindrical_operator( epsilon: vfdfield_t, rmin: float, ) -> sparse.spmatrix: - """ + r""" Cylindrical coordinate waveguide operator of the form - (NOTE: See 10.1364/OL.33.001848) - TODO: consider 10.1364/OE.20.021583 - - TODO + $$ + (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\ + 0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} + + \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\ + T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1} + \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + + \begin{bmatrix} \tilde{\partial}_x \\ + \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1} + \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix}) + \begin{bmatrix} E_x \\ + E_y \end{bmatrix} + $$ for use with a field vector of the form `[E_r, E_y]`. @@ -46,11 +106,13 @@ def cylindrical_operator( which can then be solved for the eigenmodes of the system (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields). + (NOTE: See module docs and 10.1364/OL.33.001848) + Args: omega: The angular frequency of the system dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid - rmin: Radius at the left edge of the simulation domain (minimum 'x') + rmin: Radius at the left edge of the simulation domain (at minimum 'x') Returns: Sparse matrix representation of the operator @@ -62,9 +124,9 @@ def cylindrical_operator( Ta, Tb = dxes2T(dxes=dxes, rmin=rmin) eps_parts = numpy.split(epsilon, 3) - eps_x = sparse.diags(eps_parts[0]) - eps_y = sparse.diags(eps_parts[1]) - eps_z_inv = sparse.diags(1 / eps_parts[2]) + eps_x = sparse.diags_array(eps_parts[0]) + eps_y = sparse.diags_array(eps_parts[1]) + eps_z_inv = sparse.diags_array(1 / eps_parts[2]) omega2 = omega * omega diag = sparse.block_diag @@ -74,18 +136,6 @@ def cylindrical_operator( lin0 = sparse.vstack((-Tb @ Dby, Ta @ Dbx)) @ Tb @ sparse.hstack((-Dfy, Dfx)) lin1 = sparse.vstack((Dfx, Dfy)) @ Ta @ eps_z_inv @ sparse.hstack((Dbx @ Tb @ eps_x, Dby @ Ta @ eps_y)) - # op = ( - # # E - # omega * omega * mu_yx @ eps_xy - # + mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) - # + sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy - - # # H - # omega * omega * eps_yx @ mu_xy - # + eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) - # + sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy - # ) - op = sq0 + lin0 + lin1 return op @@ -99,7 +149,6 @@ def solve_modes( mode_margin: int = 2, ) -> tuple[vcfdfield_t, NDArray[numpy.complex128]]: """ - TODO: fixup Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode of the bent waveguide with the specified mode number. @@ -142,7 +191,7 @@ def solve_modes( # Wavenumbers assume the mode is at rmin, which is unlikely # Instead, return the wavenumber in inverse radians - angular_wavenumbers = wavenumbers * rmin + angular_wavenumbers = wavenumbers * cast(complex, rmin) order = angular_wavenumbers.argsort()[::-1] e_xys = e_xys[order] @@ -168,7 +217,7 @@ def solve_mode( (e_xy, angular_wavenumber) """ kwargs['mode_numbers'] = [mode_number] - e_xys, wavenumbers = solve_modes(*args, **kwargs) + e_xys, angular_wavenumbers = solve_modes(*args, **kwargs) return e_xys[0], angular_wavenumbers[0] @@ -184,11 +233,13 @@ def linear_wavenumbers( and the mode's energy distribution. Args: - e_xys: Vectorized mode fields with shape [num_modes, 2 * x *y) - angular_wavenumbers: Angular wavenumbers corresponding to the fields in `e_xys` + e_xys: Vectorized mode fields with shape (num_modes, 2 * x *y) + angular_wavenumbers: Wavenumbers assuming fields have theta-dependence of + `exp(-i * angular_wavenumber * theta)`. They should satisfy + `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` epsilon: Vectorized dielectric constant grid with shape (3, x, y) dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) - rmin: Radius at the left edge of the simulation domain (minimum 'x') + rmin: Radius at the left edge of the simulation domain (at minimum 'x') Returns: NDArray containing the calculated linear (1/distance) wavenumbers @@ -211,3 +262,241 @@ def linear_wavenumbers( return lin_wavenumbers +def exy2h( + angular_wavenumber: complex, + omega: float, + dxes: dx_lists_t, + rmin: float, + epsilon: vfdfield_t, + mu: vfdfield_t | None = None + ) -> sparse.spmatrix: + """ + Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, + into a vectorized H containing all three H components + + Args: + angular_wavenumber: Wavenumber assuming fields have theta-dependence of + `exp(-i * angular_wavenumber * theta)`. It should satisfy + `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` + omega: The angular frequency of the system + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) + rmin: Radius at the left edge of the simulation domain (at minimum 'x') + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representing the operator. + """ + e2hop = e2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, mu=mu) + return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) + + +def exy2e( + angular_wavenumber: complex, + omega: float, + dxes: dx_lists_t, + rmin: float, + epsilon: vfdfield_t, + ) -> sparse.spmatrix: + """ + Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, + into a vectorized E containing all three E components + + Unlike the straight waveguide case, the H_z components do not cancel and must be calculated + from E_x and E_y in order to then calculate E_z. + + Args: + angular_wavenumber: Wavenumber assuming fields have theta-dependence of + `exp(-i * angular_wavenumber * theta)`. It should satisfy + `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` + omega: The angular frequency of the system + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) + rmin: Radius at the left edge of the simulation domain (at minimum 'x') + epsilon: Vectorized dielectric constant grid + + Returns: + Sparse matrix representing the operator. + """ + Dfx, Dfy = deriv_forward(dxes[0]) + Dbx, Dby = deriv_back(dxes[1]) + wavenumber = angular_wavenumber / rmin + + Ta, Tb = dxes2T(dxes=dxes, rmin=rmin) + Tai = sparse.diags_array(1 / Ta.diagonal()) + Tbi = sparse.diags_array(1 / Tb.diagonal()) + + epsilon_parts = numpy.split(epsilon, 3) + epsilon_x, epsilon_y = (sparse.diags_array(epsi) for epsi in epsilon_parts[:2]) + epsilon_z_inv = sparse.diags_array(1 / epsilon_parts[2]) + + n_pts = dxes[0][0].size * dxes[0][1].size + zeros = sparse.coo_array((n_pts, n_pts)) + keep_x = sparse.block_array([[sparse.eye_array(n_pts), None], [None, zeros]]) + keep_y = sparse.block_array([[zeros, None], [None, sparse.eye_array(n_pts)]]) + + mu_z = numpy.ones(n_pts) + mu_z_inv = sparse.diags_array(1 / mu_z) + exy2hz = 1 / (-1j * omega) * mu_z_inv @ sparse.hstack((Dfy, -Dfx)) + hxy2ez = 1 / (1j * omega) * epsilon_z_inv @ sparse.hstack((Dby, -Dbx)) + + exy2hy = Tb / (1j * wavenumber) @ (-1j * omega * sparse.hstack((epsilon_x, zeros)) - Dby @ exy2hz) + exy2hx = Tb / (1j * wavenumber) @ ( 1j * omega * sparse.hstack((zeros, epsilon_y)) - Tai @ Dbx @ Tb @ exy2hz) + + exy2ez = hxy2ez @ sparse.vstack((exy2hx, exy2hy)) + + op = sparse.vstack((sparse.eye_array(2 * n_pts), + exy2ez)) + return op + + +def e2h( + angular_wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + rmin: float, + mu: vfdfield_t | None = None + ) -> sparse.spmatrix: + r""" + Returns an operator which, when applied to a vectorized E eigenfield, produces + the vectorized H eigenfield. + + This operator is created directly from the initial coordinate-transformed equations: + $$ + \begin{aligned} + \imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\ + \imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\ + \imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ + \end{aligned} + $$ + + Args: + angular_wavenumber: Wavenumber assuming fields have theta-dependence of + `exp(-i * angular_wavenumber * theta)`. It should satisfy + `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` + omega: The angular frequency of the system + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) + rmin: Radius at the left edge of the simulation domain (at minimum 'x') + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representation of the operator. + """ + Dfx, Dfy = deriv_forward(dxes[0]) + Ta, Tb = dxes2T(dxes=dxes, rmin=rmin) + Tai = sparse.diags_array(1 / Ta.diagonal()) + Tbi = sparse.diags_array(1 / Tb.diagonal()) + + jB = 1j * angular_wavenumber / rmin + op = sparse.block_array([[ None, -jB * Tai, -Dfy], + [jB * Tbi, None, Tbi @ Dfx @ Ta], + [ Dfy, -Dfx, None]]) / (-1j * omega) + if mu is not None: + op = sparse.diags_array(1 / mu) @ op + return op + + +def dxes2T( + dxes: dx_lists_t, + rmin: float, + ) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]: + r""" + Returns the $T_a$ and $T_b$ diagonal matrices which are used to apply the cylindrical + coordinate transformation in various operators. + + Args: + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) + rmin: Radius at the left edge of the simulation domain (at minimum 'x') + + Returns: + Sparse matrix representations of the operators Ta and Tb + """ + ra = rmin + numpy.cumsum(dxes[0][0]) # Radius at Ey points + rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0]) # Radius at Ex points + ta = ra / rmin + tb = rb / rmin + + Ta = sparse.diags_array(vec(ta[:, None].repeat(dxes[0][1].size, axis=1))) + Tb = sparse.diags_array(vec(tb[:, None].repeat(dxes[1][1].size, axis=1))) + return Ta, Tb + + +def normalized_fields_e( + e_xy: ArrayLike, + angular_wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + rmin: float, + epsilon: vfdfield_t, + mu: vfdfield_t | None = None, + prop_phase: float = 0, + ) -> tuple[vcfdfield_t, vcfdfield_t]: + """ + Given a vector `e_xy` containing the vectorized E_x and E_y fields, + returns normalized, vectorized E and H fields for the system. + + Args: + e_xy: Vector containing E_x and E_y fields + angular_wavenumber: Wavenumber assuming fields have theta-dependence of + `exp(-i * angular_wavenumber * theta)`. It should satisfy + `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` + omega: The angular frequency of the system + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) + rmin: Radius at the left edge of the simulation domain (at minimum 'x') + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction. + Default 0 (continuous propagation direction, i.e. dz->0). + + Returns: + `(e, h)`, where each field is vectorized, normalized, + and contains all three vector components. + """ + e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy + h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy + e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, + mu=mu, prop_phase=prop_phase) + return e_norm, h_norm + + +def _normalized_fields( + e: vcfdfield_t, + h: vcfdfield_t, + omega: complex, + dxes: dx_lists_t, + rmin: float, + epsilon: vfdfield_t, + mu: vfdfield_t | None = None, + prop_phase: float = 0, + ) -> tuple[vcfdfield_t, vcfdfield_t]: + h *= -1 + # TODO documentation for normalized_fields + shape = [s.size for s in dxes[0]] + dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)] + + # Find time-averaged Sz and normalize to it + # H phase is adjusted by a half-cell forward shift for Yee cell, and 1-cell reverse shift for Poynting + phase = numpy.exp(-1j * -prop_phase / 2) + Sz_tavg = waveguide_2d.inner_product(e, h, dxes=dxes, prop_phase=prop_phase, conj_h=True).real # Note, using linear poynting vector + assert Sz_tavg > 0, f'Found a mode propagating in the wrong direction! {Sz_tavg=}' + + energy = numpy.real(epsilon * e.conj() * e) + + norm_amplitude = 1 / numpy.sqrt(Sz_tavg) + norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric + + # Try to break symmetry to assign a consistent sign [experimental] + E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape) + sign = numpy.sign(E_weighted[:, + :max(shape[0] // 2, 1), + :max(shape[1] // 2, 1)].real.sum()) + assert sign != 0 + + norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle) + + print('\nAAA\n', waveguide_2d.inner_product(e, h, dxes, prop_phase=prop_phase)) + e *= norm_factor + h *= norm_factor + print(f'{sign=} {norm_amplitude=} {norm_angle=} {prop_phase=}') + print(waveguide_2d.inner_product(e, h, dxes, prop_phase=prop_phase)) + + return e, h diff --git a/meanas/fdmath/functional.py b/meanas/fdmath/functional.py index 1b5811d..034d4ba 100644 --- a/meanas/fdmath/functional.py +++ b/meanas/fdmath/functional.py @@ -14,7 +14,7 @@ from .types import fdfield_t, fdfield_updater_t def deriv_forward( - dx_e: Sequence[NDArray[floating]] | None = None, + dx_e: Sequence[NDArray[floating | complexfloating]] | None = None, ) -> tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]: """ Utility operators for taking discretized derivatives (backward variant). @@ -38,7 +38,7 @@ def deriv_forward( def deriv_back( - dx_h: Sequence[NDArray[floating]] | None = None, + dx_h: Sequence[NDArray[floating | complexfloating]] | None = None, ) -> tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]: """ Utility operators for taking discretized derivatives (forward variant). @@ -65,7 +65,7 @@ TT = TypeVar('TT', bound='NDArray[floating | complexfloating]') def curl_forward( - dx_e: Sequence[NDArray[floating]] | None = None, + dx_e: Sequence[NDArray[floating | complexfloating]] | None = None, ) -> Callable[[TT], TT]: r""" Curl operator for use with the E field. @@ -94,7 +94,7 @@ def curl_forward( def curl_back( - dx_h: Sequence[NDArray[floating]] | None = None, + dx_h: Sequence[NDArray[floating | complexfloating]] | None = None, ) -> Callable[[TT], TT]: r""" Create a function which takes the backward curl of a field. @@ -123,7 +123,7 @@ def curl_back( def curl_forward_parts( - dx_e: Sequence[NDArray[floating]] | None = None, + dx_e: Sequence[NDArray[floating | complexfloating]] | None = None, ) -> Callable: Dx, Dy, Dz = deriv_forward(dx_e) @@ -136,7 +136,7 @@ def curl_forward_parts( def curl_back_parts( - dx_h: Sequence[NDArray[floating]] | None = None, + dx_h: Sequence[NDArray[floating | complexfloating]] | None = None, ) -> Callable: Dx, Dy, Dz = deriv_back(dx_h) diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py index 5d50670..946eb88 100644 --- a/meanas/fdmath/operators.py +++ b/meanas/fdmath/operators.py @@ -6,7 +6,7 @@ Basic discrete calculus etc. from collections.abc import Sequence import numpy from numpy.typing import NDArray -from numpy import floating +from numpy import floating, complexfloating from scipy import sparse from .types import vfdfield_t @@ -97,7 +97,7 @@ def shift_with_mirror( def deriv_forward( - dx_e: Sequence[NDArray[floating]], + dx_e: Sequence[NDArray[floating | complexfloating]], ) -> list[sparse.spmatrix]: """ Utility operators for taking discretized derivatives (forward variant). @@ -124,7 +124,7 @@ def deriv_forward( def deriv_back( - dx_h: Sequence[NDArray[floating]], + dx_h: Sequence[NDArray[floating | complexfloating]], ) -> list[sparse.spmatrix]: """ Utility operators for taking discretized derivatives (backward variant). @@ -219,7 +219,7 @@ def avg_back(axis: int, shape: Sequence[int]) -> sparse.spmatrix: def curl_forward( - dx_e: Sequence[NDArray[floating]], + dx_e: Sequence[NDArray[floating | complexfloating]], ) -> sparse.spmatrix: """ Curl operator for use with the E field. @@ -235,7 +235,7 @@ def curl_forward( def curl_back( - dx_h: Sequence[NDArray[floating]], + dx_h: Sequence[NDArray[floating | complexfloating]], ) -> sparse.spmatrix: """ Curl operator for use with the H field. diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py index 3871801..8f5ff39 100644 --- a/meanas/fdmath/vectorization.py +++ b/meanas/fdmath/vectorization.py @@ -49,15 +49,15 @@ def vec( @overload -def unvec(v: None, shape: Sequence[int], nvdim: int) -> None: +def unvec(v: None, shape: Sequence[int], nvdim: int = 3) -> None: pass @overload -def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int) -> fdfield_t: +def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int = 3) -> fdfield_t: pass @overload -def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int) -> cfdfield_t: +def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3) -> cfdfield_t: pass def unvec(