wip

src/lumicks/pyoptics/farfield_transform/__init__.py

@@ -1 +1,2 @@
 from .nf2ff_czt import nearfield_to_farfield_czt
+from .nf2ff_currents import farfield_factory

src/lumicks/pyoptics/farfield_transform/nf2ff_currents.py

@@ -1,7 +1,29 @@
+from typing import Optional, Tuple
+
 import numpy as np
+from scipy.constants import epsilon_0, mu_0
+import matplotlib.pyplot as plt
+
+from ..mathutils.integration import (
+    determine_integration_order,
+    get_integration_locations,
+    get_nearest_order,
+)
 from ..objective import Objective
 from ..trapping.bead import Bead
 from ..trapping.interface import focus_field_factory
+from ..trapping.local_coordinates import LocalBeadCoordinates
+from ..field_distributions.dipole import field_dipole_y, field_dipole_z, field_dipole_x
+
+
+def outer_product(x1: np.ndarray, x2: np.ndarray) -> np.ndarray:
+    x, y, z = x1
+    vx, vy, vz = x2
+
+    px = y * vz - z * vy
+    py = z * vx - x * vz
+    pz = x * vy - y * vx
+    return [px, py, pz]
 
 
 def farfield_factory(
@@ -43,54 +65,85 @@ def farfield_factory(
         local_coordinates,
         False,
     )
-    _eps = EPS0 * bead.n_medium**2
-    _mu = MU0
 
-    n = np.vstack((x, y, z))
+    eta = (mu_0 / epsilon_0) ** 0.5 / bead.n_medium
+
     bfp = condenser.get_back_focal_plane_coordinates(condenser_bfp_sampling_n)
     cos_theta, sin_theta, cos_phi, sin_phi, aperture = condenser.get_farfield_cosines(bfp)
+    # fig, ax = plt.subplots(1, 5)
+    # for idx, field in enumerate([cos_theta, sin_theta, cos_phi, sin_phi, aperture]):
+    #     p = ax[idx].imshow(field)
+    #     fig.colorbar(p)
+    # fig.show()
+
     k = 2 * np.pi * bead.n_medium / bead.lambda_vac
 
     def farfield(bead_center: Tuple[float, float, float], num_threads: Optional[int] = None):
         bead_center = np.atleast_2d(bead_center)
         # shape = (numbers of bead positions, number of field coordinates)
-        Ex, Ey, Ez, Hx, Hy, Hz = external_fields_func(bead_center, True, True, True, num_threads)
-        E = np.stack((Ex, Ey, Ez), axis=2)
-        H = np.stack((Hx, Hy, Hz), axis=2)
-        Ms, Js = [field - np.sum(field * n.reshape(bead_center.shape[0], x.size, 3), axis=2) * field for field in (E, H)]
-        Ms *= -1.0
+        Ex, Ey, Ez, Hx, Hy, Hz = field_dipole_z(1, bead.n_medium, bead.lambda_vac, xb, yb, zb)
+        # external_fields_func(bead_center, True, True, True, num_threads)
+        J_x, J_y, J_z = outer_product([x, y, z], [Hx, Hy, Hz])
+        M_x, M_y, M_z = outer_product([Ex, Ey, Ez], [x, y, z])
         L_phi, L_theta, N_phi, N_theta = [
             np.zeros_like(cos_theta, dtype="complex128") for _ in range(4)
         ]
-        for idx in range(len(cos_phi)):
-            weighted_phasor = w * np.exp(
-                (-1j * k)
-                * (
-                    ((xb + bead_center[0]) * cos_phi[idx] + (yb + bead_center[1]) * sin_phi[idx])
-                    * sin_theta[idx]
-                    + (zb + bead_center[2]) * cos_theta[idx]
+        E_theta, E_phi = [], []
+        for bead_idx, (bead_pos_x, bead_pos_y, bead_pos_z) in enumerate(bead_center):
+            for idx in np.flatnonzero(aperture):
+                weighted_phasor = w * np.exp(
+                    (1j * k)
+                    * (
+                        (
+                            (xb + bead_pos_x) * cos_phi.flat[idx]
+                            + (yb + bead_pos_y) * sin_phi.flat[idx]
+                        )
+                        * sin_theta.flat[idx]
+                        + (zb + bead_pos_z) * cos_theta.flat[idx]
+                    )
                 )
+                L_phi.flat[idx] = (
+                    (-M_x[bead_idx] * sin_phi.flat[idx] + M_y[bead_idx] * cos_phi.flat[idx])
+                    * weighted_phasor
+                ).sum()
+                L_theta.flat[idx] = (
+                    (
+                        (M_x[bead_idx] * cos_phi.flat[idx] + M_y[bead_idx] * sin_phi.flat[idx])
+                        * cos_theta.flat[idx]
+                        - M_z[bead_idx] * sin_theta.flat[idx]
+                    )
+                    * weighted_phasor
+                ).sum()
+                N_phi.flat[idx] = (
+                    (-J_x[bead_idx] * sin_phi.flat[idx] + J_y[bead_idx] * cos_phi.flat[idx])
+                    * weighted_phasor
+                ).sum()
+                N_theta.flat[idx] = (
+                    (
+                        (J_x[bead_idx] * cos_phi.flat[idx] + J_y[bead_idx] * sin_phi.flat[idx])
+                        * cos_theta.flat[idx]
+                        - J_z[bead_idx] * sin_theta.flat[idx]
+                    )
+                    * weighted_phasor
+                ).sum()
+
+            fig, ax = plt.subplots(1, 4)
+            for idx, field in enumerate([L_phi, L_theta, N_phi, N_theta]):
+                fig.colorbar(ax[idx].imshow(np.abs(field)), ax=ax[idx])
+            fig.show()
+            E_theta.append(
+                (1j * k * np.exp(1j * k * condenser.focal_length))
+                / (4 * np.pi * condenser.focal_length)
+                * (4 * np.pi)
+                * (L_phi + eta * N_theta)
             )
-            L_phi[idx] = ((-Ms[0] * sin_phi[idx] + Ms[1] * cos_phi[idx]) * weighted_phasor).sum()
-            L_theta[idx] = (
-                (
-                    (Ms[0] * cos_phi[idx] + Ms[1] * sin_phi[idx]) * cos_theta[idx]
-                    - Ms[2] * sin_theta[idx]
-                )
-                * weighted_phasor
-            ).sum()
-            N_phi[idx] = ((-Js[0] * sin_phi + Js[1] * cos_phi) * weighted_phasor).sum()
-            N_theta[idx] = (
-                ((Js[0] * cos_phi + Js[1] * sin_phi) * cos_theta - Js[2] * sin_theta)
-                * weighted_phasor
-            ).sum()
-
-        eta = (_mu / _eps) ** 0.5
-        E_theta = (1j * k * np.exp(1j * k * condenser.focal_length)) / (
-            (4 * np.pi * condenser.focal_length) * (4 * np.pi) * (L_phi + eta * N_theta)
-        )
-        E_phi = (-1j * k * np.exp(1j * k * condenser.focal_length)) / (
-            (4 * np.pi * condenser.focal_length) * (4 * np.pi) * (L_theta - eta * N_phi)
-        )
+            E_phi.append(
+                (-1j * k * np.exp(1j * k * condenser.focal_length))
+                / (4 * np.pi * condenser.focal_length)
+                * (4 * np.pi)
+                * (L_theta - eta * N_phi)
+            )
+
+        return np.squeeze(np.asarray(E_theta)), np.squeeze(np.asarray(E_phi))
 
     return farfield

tests/mathutils/test_cross.py

@@ -0,0 +1,44 @@
+import numpy as np
+import pytest
+
+from lumicks.pyoptics.mathutils.integration import get_integration_locations
+from lumicks.pyoptics.farfield_transform.nf2ff_currents import outer_product
+
+
+
+@pytest.mark.parametrize(
+    "e1,e2,result",
+    [
+        ((1, 0, 0), (0, 1, 0), (0, 0, 1)),
+        ((0, 1, 0), (1, 0, 0), (0, 0, -1)),
+        ((0, 1, 0), (0, 0, 1), (1, 0, 0)),
+        ((0, 0, 1), (0, 1, 0), (-1, 0, 0)),
+        ((1, 0, 0), (0, 0, 1), (0, -1, 0)),
+        ((0, 0, 1), (1, 0, 0), (0, 1, 0)),
+        ((1, 0, 0), (1, 0, 0), (0, 0, 0)),
+        ((0, 1, 0), (0, 1, 0), (0, 0, 0)),
+        ((0, 0, 1), (0, 0, 1), (0, 0, 0)),
+    ],
+)
+def test_cross(e1, e2, result):
+    np.testing.assert_allclose(outer_product(e1, e2), result)
+
+
+@pytest.mark.parametrize("order", [3, 15, 31])
+def test_cross2(order: int):
+    x, y, z, _ = get_integration_locations(order, "lebedev-laikov")
+
+    rho = np.hypot(x, y)
+    mask = rho > 0
+    cos_phi = np.ones_like(rho)
+    sin_phi = np.zeros_like(rho)
+    cos_phi[mask] = x[mask] / rho[mask]
+    sin_phi[mask] = y[mask] / rho[mask]
+    cos_theta = z
+    sin_theta = ((1 - z) * (1 + z)) ** 0.5
+    er = [x, y, z]
+    et = [cos_phi * cos_theta, sin_phi * cos_theta, -sin_theta]
+    ep = [-sin_phi, cos_phi, np.zeros(sin_phi.shape)]
+    np.testing.assert_allclose(outer_product(ep, er), et, atol=1e-7)
+    np.testing.assert_allclose(outer_product(er, et), ep, atol=1e-7)
+    np.testing.assert_allclose(outer_product(et, ep), er, atol=1e-7)