diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index 5fda683..22248f1 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^{-\imath \beta z} \\ -\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\imath \beta z} \\ +\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} \\ \end{aligned} $$ @@ -40,57 +40,56 @@ Substituting in our expressions for $\vec{E}$, $\vec{H}$ and discretizing: $$ \begin{aligned} --\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_{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_{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 H_y \\ -\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \hat{\partial}_x H_z \\ +\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_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ \end{aligned} $$ Rewrite the last three equations as - $$ \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 \\ +\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 \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 $\imath \beta \tilde{\partial}_x$ to the last equation, -then substitute in for $\imath \beta H_x$ and $\imath \beta H_y$: +Now apply $\gamma \tilde{\partial}_x$ to the last equation, +then substitute in for $\gamma H_x$ and $\gamma H_y$: $$ \begin{aligned} -\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 \\ +\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 \\ &= \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) \\ -\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) \\ +\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) \\ \end{aligned} $$ -With a similar approach (but using $\imath \beta \tilde{\partial}_y$ instead), we can get +With a similar approach (but using $\gamma \tilde{\partial}_y$ instead), we can get $$ \begin{aligned} -\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) \\ +\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) \\ \end{aligned} $$ -We can combine this equation for $\imath \beta \tilde{\partial}_y E_z$ with +We can combine this equation for $\gamma \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} \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 ( +-\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 ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ @@ -101,24 +100,25 @@ and $$ \begin{aligned} --\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 ( +-\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 ( \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 $\imath \beta H_x$ and the so-far unused +However, based on our rewritten equation for $\gamma 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} (\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) \\ +-\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) \\ \end{aligned} $$ @@ -126,7 +126,7 @@ and, similarly, $$ \begin{aligned} --\imath \omega \mu_{yy} (\imath \beta H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +-\imath \omega \mu_{yy} (\gamma 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} -\beta^2 E_x - \tilde{\partial}_x ( +-\gamma^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) \\ --\beta^2 E_y + \tilde{\partial}_y ( +\gamma^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,13 +165,14 @@ $$ E_y \end{bmatrix} $$ -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. +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. 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 $\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 $\gamma$ and $\beta$ will need adjustment +to account for numerical dispersion if the result is introduced into a space with a discretized z-axis. """ @@ -530,37 +531,10 @@ 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` @@ -931,7 +905,7 @@ def inner_product( # TODO documentation prop_phase: float = 0, conj_h: bool = False, trapezoid: bool = False, - ) -> complex: + ) -> tuple[vcfdfield_t, vcfdfield_t]: shape = [s.size for s in dxes[0]] @@ -946,7 +920,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 trapezoid: + if integrate: 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 2c00d02..227008d 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -1,66 +1,14 @@ -r""" +""" Operators and helper functions for cylindrical waveguides with unchanging cross-section. -Waveguide operator is derived according to 10.1364/OL.33.001848. -The curl equations in the complex coordinate system become +WORK IN PROGRESS, CURRENTLY BROKEN -$$ -\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 +As the z-dependence is known, 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, ...]]]`). """ -from typing import Any, cast +# TODO update module docs + +from typing import Any from collections.abc import Sequence import logging @@ -71,7 +19,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__) @@ -82,21 +30,13 @@ def cylindrical_operator( epsilon: vfdfield_t, rmin: float, ) -> sparse.spmatrix: - r""" + """ Cylindrical coordinate waveguide operator of the form - $$ - (\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} - $$ + (NOTE: See 10.1364/OL.33.001848) + TODO: consider 10.1364/OE.20.021583 + + TODO for use with a field vector of the form `[E_r, E_y]`. @@ -106,13 +46,11 @@ 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 (at minimum 'x') + rmin: Radius at the left edge of the simulation domain (minimum 'x') Returns: Sparse matrix representation of the operator @@ -124,9 +62,9 @@ def cylindrical_operator( Ta, Tb = dxes2T(dxes=dxes, rmin=rmin) eps_parts = numpy.split(epsilon, 3) - 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]) + eps_x = sparse.diags(eps_parts[0]) + eps_y = sparse.diags(eps_parts[1]) + eps_z_inv = sparse.diags(1 / eps_parts[2]) omega2 = omega * omega diag = sparse.block_diag @@ -136,6 +74,18 @@ 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 @@ -149,6 +99,7 @@ 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. @@ -191,7 +142,7 @@ def solve_modes( # Wavenumbers assume the mode is at rmin, which is unlikely # Instead, return the wavenumber in inverse radians - angular_wavenumbers = wavenumbers * cast(complex, rmin) + angular_wavenumbers = wavenumbers * rmin order = angular_wavenumbers.argsort()[::-1] e_xys = e_xys[order] @@ -217,7 +168,7 @@ def solve_mode( (e_xy, angular_wavenumber) """ kwargs['mode_numbers'] = [mode_number] - e_xys, angular_wavenumbers = solve_modes(*args, **kwargs) + e_xys, wavenumbers = solve_modes(*args, **kwargs) return e_xys[0], angular_wavenumbers[0] @@ -233,13 +184,11 @@ def linear_wavenumbers( and the mode's energy distribution. Args: - 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` + e_xys: Vectorized mode fields with shape [num_modes, 2 * x *y) + angular_wavenumbers: Angular wavenumbers corresponding to the fields in `e_xys` 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 (at minimum 'x') + rmin: Radius at the left edge of the simulation domain (minimum 'x') Returns: NDArray containing the calculated linear (1/distance) wavenumbers @@ -262,241 +211,3 @@ 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 034d4ba..1b5811d 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 | complexfloating]] | None = None, + dx_e: Sequence[NDArray[floating]] | 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 | complexfloating]] | None = None, + dx_h: Sequence[NDArray[floating]] | 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 | complexfloating]] | None = None, + dx_e: Sequence[NDArray[floating]] | 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 | complexfloating]] | None = None, + dx_h: Sequence[NDArray[floating]] | 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 | complexfloating]] | None = None, + dx_e: Sequence[NDArray[floating]] | 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 | complexfloating]] | None = None, + dx_h: Sequence[NDArray[floating]] | None = None, ) -> Callable: Dx, Dy, Dz = deriv_back(dx_h) diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py index 946eb88..5d50670 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, complexfloating +from numpy import floating 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 | complexfloating]], + dx_e: Sequence[NDArray[floating]], ) -> 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 | complexfloating]], + dx_h: Sequence[NDArray[floating]], ) -> 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 | complexfloating]], + dx_e: Sequence[NDArray[floating]], ) -> 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 | complexfloating]], + dx_h: Sequence[NDArray[floating]], ) -> sparse.spmatrix: """ Curl operator for use with the H field. diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py index 8f5ff39..3871801 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 = 3) -> None: +def unvec(v: None, shape: Sequence[int], nvdim: int) -> None: pass @overload -def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int = 3) -> fdfield_t: +def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int) -> fdfield_t: pass @overload -def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3) -> cfdfield_t: +def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int) -> cfdfield_t: pass def unvec(