store some progress

src/lumicks/pyoptics/farfield.py

@@ -45,21 +45,20 @@ def transform_to_xyz(self):
 
         return Ex, Ey, Ez
 
+    def get_unit_vectors(
+        self: "FarfieldData",
+    ) -> Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray]]:
+        """Returns the directional (co)sines as unit vector (sx, sy, sz)
 
-def get_unit_vectors(
-    self: "FarfieldData",
-) -> Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray]]:
-    """Returns the directional (co)sines as unit vector (sx, sy, sz)
-
-    Returns
-    -------
-    Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray]]
-        The return values sx, sy, sz
-    """
-    sz = self.cos_theta
-    sx = self.sin_theta * self.cos_phi
-    sy = self.sin_theta * self.sin_phi
-    return sx, sy, sz
+        Returns
+        -------
+        Tuple[Union[float, np.ndarray], Union[float, np.ndarray], Union[float, np.ndarray]]
+            The return values sx, sy, sz
+        """
+        sz = self.cos_theta
+        sx = self.sin_theta * self.cos_phi
+        sy = self.sin_theta * self.sin_phi
+        return sx, sy, sz
 
 
 def get_k_vectors_from_unit_vectors(
@@ -80,7 +79,7 @@ def get_k_vectors_from_unit_vectors(
 
 
 def get_k_vectors_from_cosines(
-    cos_phi, sin_phi, cos_theta, sin_theta, k: float, direction: PropagationDirection
+    cos_theta, sin_theta, cos_phi, sin_phi, k: float, direction: PropagationDirection
 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
     # Calculate properties of the plane waves. As they come from the negative z-direction, a point
     # at infinity with a negative x coordinate leads to a positive value for kx (as the wave is
@@ -89,3 +88,73 @@ def get_k_vectors_from_cosines(
     sign = -1 if direction == PropagationDirection.TO_ORIGIN else 1
     kx, ky, kz = [sign * k * cos_phi * sin_theta, sign * k * sin_phi * sin_theta, k * cos_theta]
     return kx, ky, kz
+
+
+def cosines_from_unit_vectors(sx, sy, sz):
+    """Calculate the sines and cosines involved in creating the set of normalized vectors given by
+    the coordinates (sx, sy, sz), and equivalent to (cos(φ) * sin(θ), sin(φ) * sin(θ), cos(θ)). Here
+    it is assumed that φ is the angle with the x-axis, and θ the angle with the positive z-axis. The
+    domain of θ = [0...π] and the domain of φ = [0...2π].
+
+    Parameters
+    ----------
+    sx : float | np.ndarray
+        Coordinate on the x-axis of a unit vector
+    sy : float | np.ndarray
+        Coordinate on the y-axis of a unit vector
+    sz : float | np.ndarray
+        Coordinate on the z-axis of a unit vector
+
+    Returns
+    -------
+    Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]
+        Returns the values for cos_theta, sin_theta, cos_phi and sin_phi, respectively.
+
+    Notes
+    -----
+    For θ == 0.0, the angle φ is undefined. By definition, this function return cos_phi = 1.0 and
+    sin_phi = 0.0.
+    """
+    cos_theta = sz
+    sin_theta = ((1 + cos_theta) * (1 - cos_theta)) ** 0.5
+    sp = np.hypot(sx, sy)
+    cos_phi = np.ones_like(sp)
+    sin_phi = np.zeros_like(sp)
+    region = sin_theta > 0
+    cos_phi[region] = sx[region] / sin_theta[region]
+    sin_phi[region] = sy[region] / sin_theta[region]
+    return cos_theta, sin_theta, cos_phi, sin_phi
+
+
+def unit_vectors_from_cosines(cos_theta, cos_phi, sin_phi):
+    sp = ((1 + cos_theta)(1 - cos_theta)) ** 0.5
+    sx = sp * cos_phi
+    sy = sp * sin_phi
+    return sx, sy
+
+
+def spherical_to_cartesian(locations, f_theta, f_phi):
+    """Convert farfield spherical field components to cartesian format
+
+    Parameters
+    ----------
+    locations : tuple[np.ndarray, np.ndarray, np.ndarray]
+        List of Numpy arrays describing the (x, y, z) Locations at which E_theta and E_phi are
+        taken.
+    f_theta : np.ndarray
+        Field component in the theta direction
+    f_phi : np.ndarray
+        Field component in the phi direction
+
+    Returns
+    -------
+    tuple[np.ndarray, np.ndarray, np.ndarray]
+        The field components in the x-, y- and z-direction.
+    """
+    locations = np.asarray(locations)
+    n = locations / np.linalg.vector_norm(locations, axis=0)
+    cos_theta, sin_theta, cos_phi, sin_phi = cosines_from_unit_vectors(n[0], n[1], n[2])
+    Ex = f_theta * cos_phi * cos_theta - f_phi * sin_phi
+    Ey = f_theta * sin_phi * cos_theta + f_phi * cos_phi
+    Ez = -f_theta * sin_theta
+    return Ex, Ey, Ez

src/lumicks/pyoptics/farfield_transform/__init__.py

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

src/lumicks/pyoptics/farfield_transform/equivalent_currents.py

@@ -0,0 +1,122 @@
+from dataclasses import dataclass
+
+import numpy as np
+
+
+@dataclass
+class NearfieldData:
+    """A data class to store the near field and everything else that is required to compute the far
+    field by integration over the sampling (integration) points.
+
+    Parameters
+    ----------
+
+    x : np.ndarray
+        The x part of the sampling (integration) locations `x`, `y` and `z`
+    y : np.ndarray
+        Same as `x`, but y component.
+    z : np.ndarray
+        Same as `y`, but z component.
+    normals: np.ndarray
+        The outward normal at every location (x, y, z). First index is x, second index is y, last
+        index is z.
+    weight : np.ndarray
+        The integration weight of every point. Used to calculate the integral by a (weighted)
+        summation of points.
+    E : np.ndarray
+        The electric field at each position, as E[axis, position] and axis = [0, 1, 2] corresponding
+        to x, y & z.
+    H : np.ndarray
+        Like E, but for the magnetic field.
+    k : float
+        The wave number inside the medium, equal to 2 * pi * n_medium / lamnda_vac
+    eta : float
+        Impedance of space inside the medium: (mu_0 / epsilon_0)**0.5 / n_medium
+
+    """
+
+    x: np.ndarray  # x-locations of integration points
+    y: np.ndarray  # y-locations of integration points
+    z: np.ndarray  # z-locations of integration points
+    normals: np.ndarray  # Normals at integration points
+    weight: np.ndarray  # weights of integration points
+    E: np.ndarray  # Ex, Ey, Ez-fields at integration points
+    H: np.ndarray  # Hx, Hy, Hz-fields at integration points
+    k: float
+    eta: float
+
+
+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 near_field_to_far_field(
+    near_field_data: NearfieldData, cos_theta, sin_theta, cos_phi, sin_phi, r, aperture
+):
+    xb, yb, zb, k = [getattr(near_field_data, item) for item in "xyzk"]
+    w = near_field_data.weight * 4 * np.pi
+    eta = near_field_data.eta
+    J_x, J_y, J_z = outer_product(near_field_data.normals, near_field_data.H)
+    M_x, M_y, M_z = outer_product(near_field_data.E, near_field_data.normals)
+    J_x, J_y, J_z, M_x, M_y, M_z = [np.atleast_1d(arr) for arr in (J_x, J_y, J_z, M_x, M_y, M_z)]
+    E_theta, E_phi = _equivalent_currents_to_farfield(
+        [xb, yb, zb],
+        [J_x, J_y, J_z],
+        [M_x, M_y, M_z],
+        w,
+        [cos_theta, sin_theta, cos_phi, sin_phi],
+        r,
+        aperture,
+        k,
+        eta,
+    )
+
+    return E_theta, E_phi
+
+
+def _equivalent_currents_to_farfield(locations, J, M, weights, cosines, r, aperture, k, eta):
+    cos_theta, sin_theta, cos_phi, sin_phi = cosines
+    xb, yb, zb = locations
+    Jx, Jy, Jz = J
+    Mx, My, Mz = M
+    L_phi, L_theta, N_phi, N_theta = [
+        np.zeros_like(cos_theta, dtype="complex128") for _ in range(4)
+    ]
+    for idx in np.flatnonzero(aperture):
+        weighted_phasor = weights * np.exp(
+            (1j * k)
+            * (
+                (xb * cos_phi.flat[idx] + yb * sin_phi.flat[idx]) * sin_theta.flat[idx]
+                + zb * cos_theta.flat[idx]
+            )
+        )
+        L_phi.flat[idx] = (
+            (-Mx * sin_phi.flat[idx] + My * cos_phi.flat[idx]) * weighted_phasor
+        ).sum()
+        L_theta.flat[idx] = (
+            (
+                (Mx * cos_phi.flat[idx] + My * sin_phi.flat[idx]) * cos_theta.flat[idx]
+                - Mz * sin_theta.flat[idx]
+            )
+            * weighted_phasor
+        ).sum()
+        N_phi.flat[idx] = (
+            (-Jx * sin_phi.flat[idx] + Jy * cos_phi.flat[idx]) * weighted_phasor
+        ).sum()
+        N_theta.flat[idx] = (
+            (
+                (Jx * cos_phi.flat[idx] + Jy * sin_phi.flat[idx]) * cos_theta.flat[idx]
+                - Jz * sin_theta.flat[idx]
+            )
+            * weighted_phasor
+        ).sum()
+    G = (1j * k * np.exp(1j * k * r)) / (4 * np.pi * r)
+    E_theta = G * (L_phi + eta * N_theta)
+    E_phi = -G * (L_theta - eta * N_phi)
+    return E_theta, E_phi