Further reorganization of angled sources; fixes to plotting

Author: momchil-flex

requirements.txt

@@ -7,7 +7,7 @@ nlopt
 matplotlib
 shapely==1.7.1
 descartes
-pydantic
+pydantic>=1.9.0
 PyYAML
 boto3
 requests

tidy3d/components/geometry.py

@@ -14,7 +14,7 @@
 from .types import Bound, Size, Coordinate, Axis, Coordinate2D, tidynumpy, Array
 from .types import Vertices, Ax, Shapely
 from .viz import add_ax_if_none, equal_aspect
-from .viz import PLOT_BUFFER, ARROW_LENGTH_FACTOR, ARROW_WIDTH_FACTOR
+from .viz import PLOT_BUFFER, ARROW_LENGTH_FACTOR, ARROW_WIDTH_FACTOR, MAX_ARROW_WIDTH_FACTOR
 from ..log import Tidy3dKeyError, SetupError, ValidationError
 from ..constants import MICROMETER, LARGE_NUMBER
 
@@ -405,15 +405,18 @@ def rotate_points(points: Array[float], axis: Coordinate, angle: float) -> Array
         rot_mat = np.zeros((3, 3))
         cos = np.cos(angle)
         sin = np.sin(angle)
-        rot_mat[0, 0] = cos + ux ** 2 * (1 - cos)
+        rot_mat[0, 0] = cos + ux**2 * (1 - cos)
         rot_mat[0, 1] = ux * uy * (1 - cos) - uz * sin
         rot_mat[0, 2] = ux * uz * (1 - cos) + uy * sin
         rot_mat[1, 0] = uy * ux * (1 - cos) + uz * sin
-        rot_mat[1, 1] = cos + uy ** 2 * (1 - cos)
+        rot_mat[1, 1] = cos + uy**2 * (1 - cos)
         rot_mat[1, 2] = uy * uz * (1 - cos) - ux * sin
         rot_mat[2, 0] = uz * ux * (1 - cos) - uy * sin
         rot_mat[2, 1] = uz * uy * (1 - cos) + ux * sin
-        rot_mat[2, 2] = cos + uz ** 2 * (1 - cos)
+        rot_mat[2, 2] = cos + uz**2 * (1 - cos)
+
+        if len(points.shape) == 1:
+            return rot_mat @ points
 
         return np.einsum("ij,jp...->ip...", rot_mat, points)
 
@@ -423,14 +426,11 @@ def reflect_points(
         polar_axis: Axis,
         angle_theta: float,
         angle_phi: float,
-        plane_point: Coordinate,
     ) -> Array[float]:
-        """Reflect a set of points in 3D at a plane defined by a point on the plane and an axis
-        normal to the plane.
-        all ``points``, array of shape (3, ...) at a plane passing through a given
-        ``plane_point``, normal to an axis defined in polar coordinates (theta, phi) w.r.t. the
+        """Reflect a set of points in 3D at a plane passing through the coordinate origin defined
+        and normal to a given axis defined in polar coordinates (theta, phi) w.r.t. the
         ``polar_axis`` which can be 0, 1, or 2.
-        
+
         Parameters
         ----------
         points : Array[float]
@@ -441,20 +441,15 @@ def reflect_points(
             Polar angle w.r.t. the polar axis.
         angle_phi : float
             Azimuth angle around the polar axis.
-        plane_point : Coordinate
-            A point lying on the reflection plane.
         """
 
-        # Offset coordinates such that ``plane_point`` is at the origin
-        points_new = points - plane_point
-
         # Rotate such that the plane normal is along the polar_axis
         axis_theta, axis_phi = [0, 0, 0], [0, 0, 0]
         axis_phi[polar_axis] = 1
         plane_axes = [0, 1, 2]
         plane_axes.pop(polar_axis)
         axis_theta[plane_axes[1]] = 1
-        points_new = self.rotate_points(points_new, axis_phi, -angle_phi)
+        points_new = self.rotate_points(points, axis_phi, -angle_phi)
         points_new = self.rotate_points(points_new, axis_theta, -angle_theta)
 
         # Flip the ``polar_axis`` coordinate of the points, which is now normal to the plane
@@ -617,7 +612,7 @@ def _intersect_dist(self, position, z0) -> float:
         dz = np.abs(z0 - position)
         if dz > self.radius:
             return None
-        return 2 * np.sqrt(self.radius ** 2 - dz ** 2)
+        return 2 * np.sqrt(self.radius**2 - dz**2)
 
 
 """ importable geometries """
@@ -786,7 +781,7 @@ def _plot_arrow(  # pylint:disable=too-many-arguments, too-many-locals
 
         plot_axis, _ = self.parse_xyz_kwargs(x=x, y=y, z=z)
         arrow_axis = [component == 0 for component in direction]
-        arrow_length = self._arrow_length(ax, length_factor)
+        arrow_length, arrow_width = self._arrow_dims(ax, length_factor, width_factor)
 
         # only add arrow if the plotting plane is perpendicular to the source
         if arrow_axis.count(0.0) > 1 or arrow_axis.index(0.0) != plot_axis:
@@ -800,7 +795,7 @@ def add_arrow(sign=1.0):
                     y=y0,
                     dx=sign * arrow_length * dx,
                     dy=sign * arrow_length * dy,
-                    width=width_factor * arrow_length,
+                    width=arrow_width,
                     color=color,
                     alpha=alpha,
                     zorder=np.inf,
@@ -812,17 +807,28 @@ def add_arrow(sign=1.0):
 
         return ax
 
-    def _arrow_length(self, ax: Ax, length_factor: float = ARROW_LENGTH_FACTOR) -> float:
-        """Length of arrow is the minimum size of the axes times the length factor."""
+    def _arrow_dims(
+        self,
+        ax: Ax,
+        length_factor: float = ARROW_LENGTH_FACTOR,
+        width_factor: float = ARROW_WIDTH_FACTOR,
+    ) -> Tuple[float, float]:
+        """Length and width of arrow based on axes size and length and width factors."""
 
         # get the sizes of the matplotlib axes
         xmin, xmax = ax.get_xlim()
         ymin, ymax = ax.get_ylim()
         ax_width = xmax - xmin
         ax_height = ymax - ymin
 
-        # apply length factor to the minimum size to get arrow width
-        return length_factor * min(ax_width, ax_height)
+        # apply length factor to the minimum size to get arrow length
+        arrow_length = length_factor * min(ax_width, ax_height)
+
+        # constrain arrow width by the maximum size and the max arrow width factor
+        arrow_width = width_factor * arrow_length
+        arrow_width = min(arrow_width, MAX_ARROW_WIDTH_FACTOR * max(ax_width, ax_height))
+
+        return arrow_length, arrow_width
 
 
 class Sphere(Circular):
@@ -854,7 +860,7 @@ def inside(self, x, y, z) -> bool:
         dist_x = np.abs(x - x0)
         dist_y = np.abs(y - y0)
         dist_z = np.abs(z - z0)
-        return (dist_x ** 2 + dist_y ** 2 + dist_z ** 2) <= (self.radius ** 2)
+        return (dist_x**2 + dist_y**2 + dist_z**2) <= (self.radius**2)
 
     def intersections(self, x: float = None, y: float = None, z: float = None):
         """Returns shapely geometry at plane specified by one non None value of x,y,z.
@@ -976,7 +982,7 @@ def inside(self, x, y, z) -> bool:
         dist_x = np.abs(x - x0)
         dist_y = np.abs(y - y0)
         dist_z = np.abs(z - z0)
-        inside_radius = (dist_x ** 2 + dist_y ** 2) <= (self.radius ** 2)
+        inside_radius = (dist_x**2 + dist_y**2) <= (self.radius**2)
         inside_height = dist_z < (self.length / 2)
         return inside_radius * inside_height
 

tidy3d/components/mode.py

@@ -3,8 +3,9 @@
 from typing import Tuple
 
 import pydantic as pd
+import numpy as np
 
-from ..constants import MICROMETER, RADIAN
+from ..constants import MICROMETER, RADIAN, GLANCING_CUTOFF
 from .base import Tidy3dBaseModel
 from .types import Axis2D, Literal
 from ..log import SetupError
@@ -83,3 +84,13 @@ def bend_axis_given(cls, val, values):
         if val is None and values.get("bend_radius") is not None:
             raise SetupError("bend_axis must also be defined if bend_radius is defined.")
         return val
+
+    @pd.validator("angle_theta", allow_reuse=True, always=True)
+    def glancing_incidence(cls, val):
+        """Warn if close to glancing incidence."""
+        if np.abs(np.pi / 2 - val) < GLANCING_CUTOFF:
+            raise SetupError(
+                "Mode propagation axis too close to glancing angle for accurate injection. "
+                "For best results, switch the injection axis."
+            )
+        return val

tidy3d/components/monitor.py

@@ -183,6 +183,34 @@ class AbstractModeMonitor(PlanarMonitor, FreqMonitor):
         description="Parameters to feed to mode solver which determine modes measured by monitor.",
     )
 
+    def plot(
+        self, x: float = None, y: float = None, z: float = None, ax: Ax = None, **kwargs
+    ) -> Ax:
+
+        # call the monitor.plot() function first
+        ax = super().plot(x=x, y=y, z=z, ax=ax, **kwargs)
+
+        # and then add an arrow using the direction comuputed from `_dir_arrow`.
+        ax = self._plot_arrow(
+            x=x,
+            y=y,
+            z=z,
+            ax=ax,
+            direction=self._dir_arrow,
+            color=ARROW_COLOR_MONITOR,
+            alpha=ARROW_ALPHA,
+            both_dirs=True,
+        )
+        return ax
+
+    @property
+    def _dir_arrow(self) -> Tuple[float, float, float]:
+        """Source direction normal vector in cartesian coordinates."""
+        dx = np.cos(self.mode_spec.angle_phi) * np.sin(self.mode_spec.angle_theta)
+        dy = np.sin(self.mode_spec.angle_phi) * np.sin(self.mode_spec.angle_theta)
+        dz = np.cos(self.mode_spec.angle_theta)
+        return self.unpop_axis(dz, (dx, dy), axis=self.size.index(0.0))
+
 
 class FieldMonitor(AbstractFieldMonitor, FreqMonitor):
     """:class:`Monitor` that records electromagnetic fields in the frequency domain.
@@ -357,34 +385,6 @@ def storage_size(self, num_cells: int, tmesh: int) -> int:
         # stores 3 complex numbers per frequency, per mode.
         return 3 * BYTES_COMPLEX * len(self.freqs) * self.mode_spec.num_modes
 
-    def plot(
-        self, x: float = None, y: float = None, z: float = None, ax: Ax = None, **kwargs
-    ) -> Ax:
-
-        # call the monitor.plot() function first
-        ax = super().plot(x=x, y=y, z=z, ax=ax, **kwargs)
-
-        # and then add an arrow using the direction comuputed from `_dir_arrow`.
-        ax = self._plot_arrow(
-            x=x,
-            y=y,
-            z=z,
-            ax=ax,
-            direction=self._dir_arrow,
-            color=ARROW_COLOR_MONITOR,
-            alpha=ARROW_ALPHA,
-            both_dirs=True,
-        )
-        return ax
-
-    @property
-    def _dir_arrow(self) -> Tuple[float, float, float]:
-        """Source direction normal vector in cartesian coordinates."""
-        normal = [0.0, 0.0, 0.0]
-        normal_axis = self.size.index(0.0)
-        normal[normal_axis] = 1.0
-        return tuple(normal)
-
 
 class ModeSolverMonitor(AbstractModeMonitor):
     """:class:`Monitor` that stores the mode data (field profiles and effective index)