diff --git a/CHANGELOG.md b/CHANGELOG.md index 348b91eede..77c4a98f21 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/tests/test_data/test_data_arrays.py b/tests/test_data/test_data_arrays.py index fa975741da..0a32089ef8 100644 --- a/tests/test_data/test_data_arrays.py +++ b/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) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 772cacfd3c..dd8155d8da 100644 --- a/tests/test_data/test_monitor_data.py +++ b/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, ) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index 6ef8875793..b1edb3745a 100644 --- a/tidy3d/__init__.py +++ b/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", diff --git a/tidy3d/components/data/data_array.py b/tidy3d/components/data/data_array.py index 9ab222d94f..34c010d887 100644 --- a/tidy3d/components/data/data_array.py +++ b/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, diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 6775e2afa8..76966f74c9 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/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(