wip

src/lumicks/pyoptics/trapping/interface.py

@@ -5,6 +5,7 @@
 from scipy.constants import epsilon_0 as EPS0
 from scipy.constants import mu_0 as MU0
 from scipy.constants import speed_of_light as _C
+from scipy.special import roots_legendre
 
 from ..mathutils.lebedev_laikov import get_integration_locations, get_nearest_order
 from ..objective import Objective
@@ -428,9 +429,9 @@ def force_factory(
     objective: Objective,
     bead: Bead,
     bfp_sampling_n: int = 31,
-    num_orders: int = None,
-    integration_orders: int = None,
-    method: Optional[str] = None,
+    num_orders: Optional[int] = None,
+    integration_order: Optional[int] = None,
+    method: str = "lebedev-laikov",
 ):
     """Create and return a function suitable to calculate the force on a bead. Items that can be
     precalculated are stored for rapid subsequent calculations of the force on the bead for
@@ -467,11 +468,20 @@ def force_factory(
         Number of orders that should be included in the calculation the Mie solution. If it is None
         (default), the code will use the Bead.number_of_orders() method to calculate a sufficient
         number.
-    integration_orders : int, optional
-        The order of the integration, following a Lebedev-Laikov integration scheme. If no order is
-        given, the code will determine an order based on the number of orders in the Mie solution.
-        If the integration order is provided, that order or the nearest higher order is used when
-        the provided order does not match one of the available orders.
+    integration_order : int, optional
+        The order of the integration. If no order is given, the code will determine an order based
+        on the number of orders in the Mie solution. If the integration order is provided, that
+        order is used for Gauss-Legendre integration. For Lebedev-Laikov integration, that order or
+        the nearest higher order is used when the provided order does not match one of the available
+        orders. The automatic determination of the integration order sets `integration_order` to
+        `num_orders + 1` for Gauss-Legendre integration, and to `2 * num_orders` for Lebedev-Laikov
+        integration. This typically leads to sufficiently accurate forces, but a convergence check
+        is recommended.
+    method : string, optional
+        Integration method. Choices are "gauss-legendre" and "lebedev-laikov" (default). With
+        automatic determination of the integration order (`integration_order = None`), the
+        Lebedev-Laikov scheme typically has less points to evaulate and therefore could be a bit
+        faster for similar precision. However, the order is limited to 131.
 
     Returns
     -------
@@ -490,49 +500,43 @@ def force_factory(
     ValueError
         Raised if the medium surrounding the bead does not match the immersion medium of the
         objective.
+        Raised when an invalid integration method was specified.
     """
     if bead.n_medium != objective.n_medium:
         raise ValueError("The immersion medium of the bead and the objective have to be the same")
 
     n_orders = bead.number_of_orders if num_orders is None else max(int(num_orders), 1)
-    if method is not None and method not in ["lebedev-laikov", "gauss-legendre", "gauss-chebyshev"]:
+    if method not in ["lebedev-laikov", "gauss-legendre"]:
         raise ValueError(f"Wrong type of integration method specified: {method}")
 
-    if integration_orders is not None:
+    if integration_order is not None:
         # Use user's integration order
-        integration_orders = np.amax((1, int(integration_orders)))
+        integration_order = np.amax((1, int(integration_order)))
         if method == "lebedev-laikov":
             # Find nearest integration order that is equal or greater than the user's
-            integration_orders = get_nearest_order(integration_orders)
+            integration_order = get_nearest_order(integration_order)
     else:
         # Determine reasonable defaults for integrating over a certain bead size
-        if method in ["gauss-chebyshev", "gauss-legendre"]:
-            # Guesstimate based on P_n ** 2 ~ x ** (2 * n + 2) and integration
-            # order m is accurate to x ** (2 * m - 1)
-            integration_orders = n_orders + 2
+        if method == "gauss-legendre":
+            # Guesstimate based on P^1_n ** 2 ~ (1 - x ** 2) * (x ** (n - 1)) ** 2 ~ x ** 2n and
+            # integration order m is accurate to x ** (2 * m - 1)
+            integration_order = n_orders + 1
         # Default integration method is lebedev-laikov
         else:
-            # Get an integration order that is one level higher than the one matching n_orders if no
-            # integration order is specified
-            integration_orders = get_nearest_order(get_nearest_order(get_nearest_order(n_orders) + 1) + 1)
-    if method in ["gauss-legendre", "gauss-chebyshev"]:
-        from scipy.special import roots_legendre, roots_chebyu
-
-        phi = np.arange(1, 2 * integration_orders + 1) * np.pi / integration_orders
-        z, w = (
-            roots_legendre(integration_orders)
-            if method == "gauss-legendre"
-            else roots_chebyu(integration_orders)
-        )
+            # Get an integration order that is one level higher than the one matching 2 * n_orders
+            # if no integration order is specified
+            integration_order = get_nearest_order(get_nearest_order(2 * n_orders))
+    if method == "gauss-legendre":
+        phi = np.arange(1, 2 * integration_order + 1) * np.pi / integration_order
+        z, w = roots_legendre(integration_order)
         x, y = np.cos(phi), np.sin(phi)
         sin_theta = ((1 - z) * (1 + z)) ** 0.5
         x = (sin_theta[:, np.newaxis] * x[np.newaxis, :]).flatten()
         y = (sin_theta[:, np.newaxis] * y[np.newaxis, :]).flatten()
         z = np.repeat(z, phi.size)
-        w = np.repeat(w, phi.size) * 0.25 / integration_orders
+        w = np.repeat(w, phi.size) * 0.25 / integration_order
     else:
-        x, y, z, w = [np.asarray(c) for c in get_integration_locations(integration_orders)]
-
+        x, y, z, w = [np.asarray(c) for c in get_integration_locations(integration_order)]
     xb, yb, zb = [c * bead.bead_diameter * 0.51 for c in (x, y, z)]
 
     local_coordinates = LocalBeadCoordinates(
@@ -573,6 +577,9 @@ def force_on_bead(bead_center: Tuple[float, float, float], num_threads: Optional
         Th13 = np.atleast_2d(_mu * np.real(Hx * np.conj(Hz)))
         Th23 = np.atleast_2d(_mu * np.real(Hy * np.conj(Hz)))
 
+        # F = 1/2 Int Re[T dot n] da ~ 4 pi r^2 Sum w * T dot n, with r the diameter of the
+        # imaginary sphere surrounding the bead on which the force (density) is calculated
+
         T = np.empty((len(bead_center), x.size, 3, 3))
         T[..., 0, 0] = Te11 + Th11
         T[..., 0, 1] = Te12 + Th12
@@ -585,10 +592,11 @@ def force_on_bead(bead_center: Tuple[float, float, float], num_threads: Optional
         T[..., 2, 2] = Te33 + Th33
 
         # T is [num_bead_positions, num_coords_around_bead, 3, 3] large, nw is [num_coords, 3] large
-        # F is the 3D force on the bead at each position `num_bead_positions`.
+        # F is the 3D force on the bead at each position `num_bead_positions`, except for the factor
+        # 1/2 * 4 pi r^2
         F = np.matmul(T, nw[np.newaxis, ..., np.newaxis])[..., 0].sum(axis=1)
 
-        # Note: factor 1/2 incorporated as 2 pi instead of 4 pi
+        # Note: the factor 1/2 is incorporated as 2 pi instead of 4 pi
         return np.squeeze(F) * (bead.bead_diameter * 0.51) ** 2 * 2 * np.pi
 
     return force_on_bead
@@ -601,7 +609,7 @@ def forces_focus(
     bead_center=(0, 0, 0),
     bfp_sampling_n=31,
     num_orders=None,
-    integration_orders=None,
+    integration_order=None,
     method=None,
 ):
     """
@@ -657,7 +665,7 @@ def forces_focus(
         bead=bead,
         bfp_sampling_n=bfp_sampling_n,
         num_orders=num_orders,
-        integration_orders=integration_orders,
+        integration_order=integration_order,
         method=method,
     )
     return force_fun(bead_center)
@@ -670,7 +678,7 @@ def absorbed_power_focus(
     bead_center=(0, 0, 0),
     bfp_sampling_n=31,
     num_orders=None,
-    integration_orders=None,
+    integration_order=None,
     verbose=False,
 ):
     """
@@ -726,14 +734,14 @@ def absorbed_power_focus(
 
     n_orders = bead.number_of_orders if num_orders is None else max(int(num_orders), 1)
 
-    if integration_orders is None:
+    if integration_order is None:
         # Go to the next higher available order of integration than should
         # be strictly necessary
         order = get_nearest_order(get_nearest_order(n_orders) + 1)
         x, y, z, w = get_integration_locations(order)
     else:
         x, y, z, w = get_integration_locations(
-            get_nearest_order(np.amax((1, int(integration_orders))))
+            get_nearest_order(np.amax((1, int(integration_order))))
         )
 
     # Enforce floats to ensure Numba has all floats when doing @ / np.matmul
@@ -778,7 +786,7 @@ def scattered_power_focus(
     bead_center=(0, 0, 0),
     bfp_sampling_n=31,
     num_orders=None,
-    integration_orders=None,
+    integration_order=None,
     verbose=False,
 ):
     """
@@ -831,14 +839,14 @@ def scattered_power_focus(
 
     n_orders = bead.number_of_orders if num_orders is None else max(int(num_orders), 1)
 
-    if integration_orders is None:
+    if integration_order is None:
         # Go to the next higher available order of integration than should
         # be strictly necessary
         order = get_nearest_order(get_nearest_order(n_orders) + 1)
         x, y, z, w = get_integration_locations(order)
     else:
         x, y, z, w = get_integration_locations(
-            get_nearest_order(np.amax((1, int(integration_orders))))
+            get_nearest_order(np.amax((1, int(integration_order))))
         )
 
     # Enforce floats to ensure Numba has all floats when doing @ / np.matmul