Refractored DirectivitData: radiation characteristics are computed vi…

…a properties

CHANGELOG.md

@@ -31,6 +31,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added an option to simulate plane wave propagation at fixed angles by setting parameter `angular_spec=FixedAngleSpec()` when defining a `PlaneWave` source.
 - The universal `tidy3d.web` can now also be used to handle `ModeSolver` simulations.
 - Support for differentiation with respect to `PolySlab.slab_bounds`.
+- `DirectivityMonitor` can now be used to decompose fields into circular polarization states.
 
 ### Changed
 - Priority is given to `snapping_points` in `GridSpec` when close to structure boundaries, which reduces the chance of them being skipped.
@@ -51,6 +52,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `estimate_cost` is now called at the end of every `web.upload` call.
 - Internal refactor of adjoint shape gradients using `GradientSurfaceMesh`.
 - Enhanced progress bar display in batch operations with better formatting, colors, and status tracking.
+- The behavior of `DirectivityData` was changed, so that only far fields and flux are stored. Radiation characteristics (e.g., directivity, axial ratio) are computed through properties.
 
 ### Fixed
 - Significant speedup for field projection computations.

tests/test_data/test_data_arrays.py

@@ -184,21 +184,16 @@ def make_mode_index_data_array():
     return td.ModeIndexDataArray(values, coords=dict(f=FS, mode_index=MODE_INDICES))
 
 
+def make_far_field_data_array():
+    values = (1 + 1j) * np.random.random((len(PD), len(THETAS), len(PHIS), len(FS)))
+    return td.FieldProjectionAngleDataArray(values, coords=dict(r=PD, theta=THETAS, phi=PHIS, f=FS))
+
+
 def make_flux_data_array():
     values = np.random.random(len(FS))
     return td.FluxDataArray(values, coords=dict(f=FS))
 
 
-def make_directivity_data_array():
-    values = np.random.random((len(PD), len(THETAS), len(PHIS), len(FS)))
-    return td.DirectivityDataArray(values, coords=dict(r=PD, theta=THETAS, phi=PHIS, f=FS))
-
-
-def make_axial_ratio_data_array():
-    values = np.random.random((len(PD), len(THETAS), len(PHIS), len(FS)))
-    return td.AxialRatioDataArray(values, coords=dict(r=PD, theta=THETAS, phi=PHIS, f=FS))
-
-
 def make_flux_time_data_array():
     values = np.random.random(len(TS))
     return td.FluxTimeDataArray(values, coords=dict(t=TS))
@@ -261,13 +256,8 @@ def test_flux_time_data_array():
     data = data.interp(t=1e-13)
 
 
-def test_directivity_data_array():
-    data = make_directivity_data_array()
-    data = data.sel(f=1e14, phi=0)
-
-
-def test_axial_ratio_data_array():
-    data = make_axial_ratio_data_array()
+def test_far_field_data_array():
+    data = make_far_field_data_array()
     data = data.sel(f=1e14, phi=0)
 
 

tests/test_data/test_monitor_data.py

@@ -33,9 +33,8 @@
     PERMITTIVITY_MONITOR,
     SIM,
     SIM_SYM,
-    make_axial_ratio_data_array,
     make_diffraction_data_array,
-    make_directivity_data_array,
+    make_far_field_data_array,
     make_flux_data_array,
     make_flux_time_data_array,
     make_mode_amps_data_array,
@@ -54,9 +53,6 @@
 GRID_CORRECTION = FreqModeDataArray(
     1 + 0.01 * np.random.rand(*N_COMPLEX.shape), coords=N_COMPLEX.coords
 )
-DIRECTIVITY = make_directivity_data_array()
-AXIALRATIO = make_axial_ratio_data_array()
-
 """ Make the montor data """
 
 
@@ -190,8 +186,16 @@ def make_flux_data():
 
 
 def make_directivity_data():
+    data = make_far_field_data_array()
     return DirectivityData(
-        monitor=DIRECTIVITY_MONITOR, directivity=DIRECTIVITY.copy(), axial_ratio=AXIALRATIO.copy()
+        monitor=DIRECTIVITY_MONITOR,
+        flux=FLUX.copy(),
+        Er=data,
+        Etheta=data,
+        Ephi=data,
+        Hr=data,
+        Htheta=data,
+        Hphi=data,
     )
 
 

tidy3d/__init__.py

@@ -37,11 +37,9 @@
 
 # data
 from .components.data.data_array import (
-    AxialRatioDataArray,
     CellDataArray,
     ChargeDataArray,
     DiffractionDataArray,
-    DirectivityDataArray,
     EMECoefficientDataArray,
     EMEModeIndexDataArray,
     EMEScalarFieldDataArray,
@@ -417,8 +415,6 @@ def set_logging_level(level: str) -> None:
     "FieldProjectionCartesianDataArray",
     "FieldProjectionKSpaceDataArray",
     "DiffractionDataArray",
-    "DirectivityDataArray",
-    "AxialRatioDataArray",
     "HeatDataArray",
     "ChargeDataArray",
     "FieldDataset",

tidy3d/components/data/data_array.py

@@ -921,55 +921,6 @@ class FieldProjectionAngleDataArray(DataArray):
     _data_attrs = {"long_name": "radiation vectors"}
 
 
-class DirectivityDataArray(DataArray):
-    """Directivity in the frequency domain as a function of angles theta and phi.
-    Directivity is a dimensionless quantity defined as the ratio of the radiation
-    intensity in a given direction to the average radiation intensity over all directions.
-
-    Example
-    -------
-    >>> f = np.linspace(1e14, 2e14, 10)
-    >>> r = np.atleast_1d(5)
-    >>> theta = np.linspace(0, np.pi, 10)
-    >>> phi = np.linspace(0, 2*np.pi, 20)
-    >>> coords = dict(r=r, theta=theta, phi=phi, f=f)
-    >>> values = np.random.random((len(r), len(theta), len(phi), len(f)))
-    >>> data = DirectivityDataArray(values, coords=coords)
-    """
-
-    __slots__ = ()
-    _dims = ("r", "theta", "phi", "f")
-    _data_attrs = {"long_name": "radiation intensity"}
-
-
-class AxialRatioDataArray(DataArray):
-    """Axial Ratio (AR) in the frequency domain as a function of angles theta and phi.
-    AR is a dimensionless quantity defined as the ratio of the major axis to the minor
-    axis of the polarization ellipse.
-
-    Note
-    ----
-    The axial ratio computation is based on:
-
-    Balanis, Constantine A., "Antenna Theory: Analysis and Design,"
-    John Wiley & Sons, Chapter 2.12 (2016).
-
-    Example
-    -------
-    >>> f = np.linspace(1e14, 2e14, 10)
-    >>> r = np.atleast_1d(5)
-    >>> theta = np.linspace(0, np.pi, 10)
-    >>> phi = np.linspace(0, 2*np.pi, 20)
-    >>> coords = dict(r=r, theta=theta, phi=phi, f=f)
-    >>> values = np.random.random((len(r), len(theta), len(phi), len(f)))
-    >>> data = AxialRatioDataArray(values, coords=coords)
-    """
-
-    __slots__ = ()
-    _dims = ("r", "theta", "phi", "f")
-    _data_attrs = {"long_name": "axial ratio"}
-
-
 class FieldProjectionCartesianDataArray(DataArray):
     """Far fields in frequency domain as a function of local x and y coordinates.
 
@@ -1267,8 +1218,6 @@ class IndexedDataArray(DataArray):
     FieldProjectionCartesianDataArray,
     FieldProjectionKSpaceDataArray,
     DiffractionDataArray,
-    DirectivityDataArray,
-    AxialRatioDataArray,
     FreqModeDataArray,
     FreqDataArray,
     TimeDataArray,

tidy3d/components/data/monitor_data.py

@@ -63,10 +63,8 @@
 )
 from ..validators import enforce_monitor_fields_present, required_if_symmetry_present
 from .data_array import (
-    AxialRatioDataArray,
     DataArray,
     DiffractionDataArray,
-    DirectivityDataArray,
     EMEFreqModeDataArray,
     FieldProjectionAngleDataArray,
     FieldProjectionCartesianDataArray,
@@ -97,6 +95,7 @@
 
 # how much to shift the adjoint field source for 0-D axes dimensions
 SHIFT_VALUE_ADJ_FLD_SRC = 1e-5
+AXIAL_RATIO_CAP = 100
 
 
 class MonitorData(AbstractMonitorData, ABC):
@@ -2028,13 +2027,14 @@ class DirectivityData(MonitorData):
 
     Example
     -------
-    >>> from tidy3d import DirectivityDataArray, AxialRatioDataArray
+    >>> from tidy3d import DirectivityDataArray, AxialRatioDataArray, FieldProjectionAngleDataArray
     >>> f = np.linspace(1e14, 2e14, 10)
     >>> r = np.atleast_1d(1e6)
     >>> theta = np.linspace(0, np.pi, 10)
     >>> phi = np.linspace(0, 2*np.pi, 20)
     >>> coords = dict(r=r, theta=theta, phi=phi, f=f)
     >>> values = np.random.random((len(r), len(theta), len(phi), len(f)))
+    >>> complex_values = (1+1j) * np.random.random((len(r), len(theta), len(phi), len(f)))
     >>> scalar_directivity_field = DirectivityDataArray(values, coords=coords)
     >>> scalar_axial_ratio_field = AxialRatioDataArray(values, coords=coords)
     >>> monitor = DirectivityMonitor(center=(1,2,3), size=(2,2,2), freqs=f, name='n2f_monitor', phi=phi, theta=theta)
@@ -2047,14 +2047,104 @@ class DirectivityData(MonitorData):
         description="Monitor describing the angle-based projection grid on which to measure directivity data.",
     )
 
-    directivity: DirectivityDataArray = pd.Field(
-        ..., title="Directivity", description="Directivity with an angle-based projection grid."
+    flux: FluxDataArray = pd.Field(..., title="Flux", description="Flux values.")
+    Er: FieldProjectionAngleDataArray = pd.Field(
+        ...,
+        title="Er",
+        description="Spatial distribution of r-component of the electric field.",
     )
-
-    axial_ratio: AxialRatioDataArray = pd.Field(
-        ..., title="Axial Ratio", description="Axial ratio with an angle-based projection grid."
+    Etheta: FieldProjectionAngleDataArray = pd.Field(
+        ...,
+        title="Etheta",
+        description="Spatial distribution of the theta-component of the electric field.",
+    )
+    Ephi: FieldProjectionAngleDataArray = pd.Field(
+        ...,
+        title="Ephi",
+        description="Spatial distribution of phi-component of the electric field.",
+    )
+    Hr: FieldProjectionAngleDataArray = pd.Field(
+        ...,
+        title="Hr",
+        description="Spatial distribution of r-component of the magnetic field.",
+    )
+    Htheta: FieldProjectionAngleDataArray = pd.Field(
+        ...,
+        title="Htheta",
+        description="Spatial distribution of theta-component of the magnetic field.",
+    )
+    Hphi: FieldProjectionAngleDataArray = pd.Field(
+        ...,
+        title="Hphi",
+        description="Spatial distribution of phi-component of the magnetic field.",
     )
 
+    @property
+    def directivity(self):
+        """Directivity in the frequency domain as a function of angles theta and phi.
+        Directivity is a dimensionless quantity defined as the ratio of the radiation
+        intensity in a given direction to the average radiation intensity over all directions."""
+
+        power_theta = 0.5 * np.real(self.Etheta * np.conj(self.Hphi))
+        power_phi = 0.5 * np.real(-self.Ephi * np.conj(self.Htheta))
+        power = power_theta + power_phi
+
+        # Normalize the aligned flux by dividing by (4 * pi * r^2) to adjust the flux for
+        # spherical surface area normalization
+        flux_normed = self.flux / (4 * np.pi * self.monitor.proj_distance**2)
+
+        return power / flux_normed
+
+    @property
+    def axial_ratio(self):
+        """Axial Ratio (AR) in the frequency domain as a function of angles theta and phi.
+        AR is a dimensionless quantity defined as the ratio of the major axis to the minor
+        axis of the polarization ellipse.
+        Note
+        ----
+        The axial ratio computation is based on:
+
+        Balanis, Constantine A., "Antenna Theory: Analysis and Design,"
+        John Wiley & Sons, Chapter 2.12 (2016).
+        """
+
+        # Calculate the terms of the equation
+        E1_abs_squared = np.abs(self.Etheta) ** 2
+        E2_abs_squared = np.abs(self.Ephi) ** 2
+        E1_squared = self.Etheta**2
+        E2_squared = self.Ephi**2
+
+        # Axial ratio calculations based on equations (2-65) to (2-67)
+        # from Balanis, Constantine A., "Antenna Theory: Analysis and Design,"
+        # John Wiley & Sons, 2016. These calculations use complex numbers
+        # directly and are equivalent to the referenced equations.
+        AR_numerator = E1_abs_squared + E2_abs_squared + np.abs(E1_squared + E2_squared)
+        AR_denominator = E1_abs_squared + E2_abs_squared - np.abs(E1_squared + E2_squared)
+
+        inds_zero = AR_numerator == 0
+        axial_ratio_inverse = xr.zeros_like(AR_numerator)
+        # Perform the axial ratio inverse calculation where the numerator is non-zero
+        axial_ratio_inverse = axial_ratio_inverse.where(
+            inds_zero, np.sqrt(np.abs(AR_denominator / AR_numerator))
+        )
+
+        # Cap the axial ratio values at 1 / AXIAL_RATIO_CAP
+        axial_ratio_inverse = axial_ratio_inverse.where(
+            axial_ratio_inverse >= 1 / AXIAL_RATIO_CAP, 1 / AXIAL_RATIO_CAP
+        )
+
+        return 1 / axial_ratio_inverse
+
+    @property
+    def left_polarization(self):
+        "Left polarization (counter clockwise component of circular polarization) electric far field with an angle-based projection grid."
+        return (self.Etheta - 1j * self.Ephi) / np.sqrt(2)
+
+    @property
+    def right_polarization(self):
+        "Right polarization (clockwise component of circular polarization) electric far field with an angle-based projection grid."
+        return (self.Etheta + 1j * self.Ephi) / np.sqrt(2)
+
 
 ProjFieldType = Union[
     FieldProjectionAngleDataArray,
@@ -2083,7 +2173,7 @@ class AbstractFieldProjectionData(MonitorData):
 
     Er: ProjFieldType = pd.Field(
         ...,
-        title="Ephi",
+        title="Er",
         description="Spatial distribution of r-component of the electric field.",
     )
     Etheta: ProjFieldType = pd.Field(
@@ -2098,7 +2188,7 @@ class AbstractFieldProjectionData(MonitorData):
     )
     Hr: ProjFieldType = pd.Field(
         ...,
-        title="Hphi",
+        title="Hr",
         description="Spatial distribution of r-component of the magnetic field.",
     )
     Htheta: ProjFieldType = pd.Field(