Near2Far

.pylintrc

@@ -1,7 +1,7 @@
 [MASTER]
 extension-pkg-allow-list=pydantic
 ignore=material_library.py, plugins
-good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, f, t, y1, y2, x1, x2, xs, ys, zs, Ax, Nx, Ny, Nz, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c,
+good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, f, t, y1, y2, x1, x2, xs, ys, zs, Ax, Nx, Ny, Nz, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi
 
 [BASIC]
 

tidy3d/convert.py

@@ -359,7 +359,7 @@ def old_json_monitors(sim: Simulation) -> Dict:
                     {
                         "type": "FrequencyMonitor",
                         "store": store,
-                        "interpolate": True,
+                        "interpolate": False,
                     }
                 )
             elif isinstance(monitor, FieldTimeMonitor):

tidy3d/plugins/near2far/near2far.py

@@ -3,63 +3,72 @@
 import numpy as np
 
 from ...constants import C_0, ETA_0
-from ...components.data import FieldData
+from ...components.data import SimulationData
+from ...log import SetupError
+
+# TODO: implement new version with simulation.discretize
 
 
 class Near2Far:
     """Near field to far field transformation tool."""
 
-    def __init__(self, field_data: FieldData):
+    def __init__(self, sim_data: SimulationData, mon_name: str, frequency: float):
         """Constructs near field to far field transformation object from monitor data.
 
         Parameters
         ----------
-        field_data : FieldData
-            Description
+        sim_data : :class:`.SimulationData`
+            Container for simulation data containing a near field monitor.
+        mon_name : str
+            Name of the :class:`.FieldMonitor` to use as source of near field.
+            Must be a :class:`.FieldMonitor` and stored in ``sim_data``.
+        frequency : float
+            Frequency to select from the :class:`.FieldMonitor` to use for projection.
+            Must be a frequency stored in the :class:`FieldMonitor`.
         """
 
-        # get frequency info
-        self.f = float(field_data.data.Ez.f[0].values)
-        self.k0 = 2 * np.pi * self.f / C_0
-
-        # get normal axis (ignore components)
-        xs = field_data.data.Ez.x.values
-        ys = field_data.data.Ez.y.values
-        zs = field_data.data.Ez.z.values
-        mon_size = [xs.shape[-1], ys.shape[-1], zs.shape[-1]]
-        self.axis = mon_size.index(1)
-        assert self.axis == 2, "Currently only works for z normal."
-
-        # get normal and planar coordinates
-        # zs, (xs, ys) = field_data.geometry.pop_axis((xs, ys, zs), axis=self.axis)
-        self.z0 = zs[0]
-        self.xs = xs
-        self.ys = ys
-        self.xx, self.yy = np.meshgrid(self.xs, self.ys, indexing="ij")
-        self.dx = np.mean(np.diff(xs))
-        self.dy = np.mean(np.diff(ys))
-
-        # get tangential near fields
-        # for em_field in "EH":
-        #     for component in "xyz":
-        #         field_name = em_field + component
-        #         assert field_name in field_data.field, f"missing field: {field_name}"
-
-        self.Ex = field_data.data.Ex.values
-        self.Ey = field_data.data.Ey.values
-        self.Ez = field_data.data.Ez.values
-        self.Hx = field_data.data.Hx.values
-        self.Hy = field_data.data.Hy.values
-        self.Hz = field_data.data.Hz.values
-
-        # _, (self.Ex, self.Ey) = field_data.geometry.pop_axis(E, axis=self.axis)
-        # _, (self.Hx, self.Hy) = field_data.geometry.pop_axis(H, axis=self.axis)
+        try:
+            field_data = sim_data[mon_name]
+        except Exception as e:
+            raise SetupError(
+                f"No data for monitor named '{mon_name}' " "found in supplied sim_data."
+            ) from e
+
+        if any(field_name not in field_data for field_name in ("Ex", "Ey", "Hx", "Hy", "Ez")):
+            raise SetupError(f"Monitor named '{mon_name}' doesn't store all field values")
+
+        try:
+            Ex = field_data["Ex"].sel(f=frequency)
+            Ey = field_data["Ey"].sel(f=frequency)
+            # self.Ez = field_data['Ez'].sel(f=frequency)
+            Hx = field_data["Hx"].sel(f=frequency)
+            Hy = field_data["Hy"].sel(f=frequency)
+            Hz = field_data["Hz"].sel(f=frequency)
+        except Exception as e:
+            raise SetupError(
+                f"Frequency {frequency} not found in all fields " f"from monitor '{mon_name}'."
+            ) from e
+
+        self.k0 = 2 * np.pi * frequency / C_0
+
+        # grid sizes
+        self.dx = np.mean(np.diff(Hz.x.values))
+        self.dy = np.mean(np.diff(Hz.y.values))
+
+        # interpolate to centers
+        x_centers = Hz.x.values
+        y_centers = Hz.y.values
+        Ex = Ex.interp(x=x_centers, y=y_centers)
+        Ey = Ey.interp(x=x_centers, y=y_centers)
+        Hx = Hx.interp(x=x_centers, y=y_centers)
+        Hy = Hy.interp(x=x_centers, y=y_centers)
+        self.xx, self.yy = np.meshgrid(x_centers, y_centers, indexing="ij")
 
         # compute equivalent sources
-        self.Jx = -self.Hy
-        self.Jy = self.Hx
-        self.Mx = self.Ey
-        self.My = -self.Ex
+        self.Jx = -np.squeeze(Hy.values)
+        self.Jy = np.squeeze(Hx.values)
+        self.Mx = np.squeeze(Ey.values)
+        self.My = -np.squeeze(Ex.values)
 
     def _radiation_vectors(self, theta, phi):
         """Compute radiation vectors at an angle in spherical coordinates
@@ -78,49 +87,32 @@ def _radiation_vectors(self, theta, phi):
         """
 
         # precompute trig functions and add extra dimensions
-        theta = np.array(theta)
-        phi = np.array(phi)
-        sin_theta = np.sin(theta).reshape((-1, 1, 1))
-        cos_theta = np.cos(theta).reshape((-1, 1, 1))
-        sin_phi = np.sin(phi).reshape((-1, 1, 1))
-        cos_phi = np.cos(phi).reshape((-1, 1, 1))
+        sin_theta = np.sin(theta)
+        cos_theta = np.cos(theta)
+        sin_phi = np.sin(phi)
+        cos_phi = np.cos(phi)
 
         # precompute fourier transform phase term {dx dy e^(ikrcos(psi))}
-        FT_phase_x = np.exp(1j * self.k0 * self.xx * sin_theta * cos_phi)
-        FT_phase_y = np.exp(1j * self.k0 * self.yy * sin_theta * sin_phi)
-        FT_phase = self.dx * self.dy * FT_phase_x * FT_phase_y
-
-        # multiply the phase terms with the current sources
-        Jx_phased = np.sum(self.Jx * FT_phase, axis=(-2, -1))
-        Jy_phased = np.sum(self.Jy * FT_phase, axis=(-2, -1))
-        Mx_phased = np.sum(self.Mx * FT_phase, axis=(-2, -1))
-        My_phased = np.sum(self.My * FT_phase, axis=(-2, -1))
-
-        # get rid of extra dimensions
-        cos_phi = np.squeeze(cos_phi)
-        sin_phi = np.squeeze(sin_phi)
-        cos_theta = np.squeeze(cos_theta)
-        sin_theta = np.squeeze(sin_theta)
+        phase_x = np.exp(1j * self.k0 * self.xx * sin_theta * cos_phi)
+        phase_y = np.exp(1j * self.k0 * self.yy * sin_theta * cos_phi)
+        phase = self.dx * self.dy * phase_x * phase_y
+
+        Jx_k = np.sum(self.Jx * phase)
+        Jy_k = np.sum(self.Jy * phase)
+        Mx_k = np.sum(self.Mx * phase)
+        My_k = np.sum(self.My * phase)
 
         # N_theta (8.33a)
-        integrand_x = +Jx_phased * cos_theta * cos_phi
-        integrand_y = +Jy_phased * cos_theta * sin_phi
-        N_theta = integrand_x + integrand_y
+        N_theta = Jx_k * cos_theta * cos_phi + Jy_k * cos_theta * sin_phi
 
         # N_phi (8.33b)
-        integrand_x = -Jx_phased * sin_phi
-        integrand_y = +Jy_phased * cos_phi
-        N_phi = integrand_x + integrand_y
+        N_phi = -Jx_k * sin_phi + Jy_k * cos_phi
 
         # L_theta  (8.34a)
-        integrand_x = +Mx_phased * cos_theta * cos_phi
-        integrand_y = +My_phased * cos_theta * sin_phi
-        L_theta = integrand_x + integrand_y
+        L_theta = Mx_k * cos_theta * cos_phi + My_k * cos_theta * sin_phi
 
         # L_phi  (8.34b)
-        integrand_x = -Mx_phased * sin_phi
-        integrand_y = +My_phased * cos_phi
-        L_phi = integrand_x + integrand_y
+        L_phi = -Mx_k * sin_phi + My_k * cos_phi
 
         return N_theta, N_phi, L_theta, L_phi
 
@@ -143,16 +135,23 @@ def fields_spherical(self, r, theta, phi):
             (Er, Etheta, Ephi), (Hr, Htheta, Hphi), fields in polar
             coordinates.
         """
+
+        # project radiation vectors to distance r away for given angles
         N_theta, N_phi, L_theta, L_phi = self._radiation_vectors(theta, phi)
         scalar_proj_r = 1j * self.k0 * np.exp(-1j * self.k0 * r) / (4 * np.pi * r)
+
+        # assemble E felds
         E_theta = -scalar_proj_r * (L_phi + ETA_0 * N_theta)
         E_phi = scalar_proj_r * (L_theta - ETA_0 * N_phi)
         E_r = np.zeros_like(E_phi)
         E = np.stack((E_r, E_theta, E_phi))
+
+        # assemble H fields
         H_theta = -E_phi / ETA_0
         H_phi = E_theta / ETA_0
         H_r = np.zeros_like(H_phi)
         H = np.stack((H_r, H_theta, H_phi))
+
         return E, H
 
     def fields_cartesian(self, x, y, z):