calibration: add coupling model

lumicks/pylake/force_calibration/detail/drag_models.py

@@ -244,3 +244,39 @@ def calculate_dn(n, k, alpha, beta, delta):
     )
     term3 = (2 * n + 1) * (2 * n + 3) * np.sinh(alpha - beta) * np.sinh(alpha + beta)
     return mul * (term1 + term2 + term3) / delta
+
+
+def coupling_correction_factor(radius1, radius2, distance, summands=5):
+    """Calculate the bead correction factors.
+
+    Parameters
+    ----------
+    radius1, radius2 : float
+        Bead radii
+    distance : float
+        Distance between the bead centers
+    summands : int
+        How many summands to use
+    """
+    a, alpha, beta = to_curvilinear_coordinates(radius1, radius2, distance)
+    # In the paper we have:
+    #   F = 2 sqrt(2) pi nu V / a * summation, whereas the stokes force was given by 6 pi nu r V
+    # We are interested in the factor that translates from the regular flow velocity, to the
+    # effective flow velocity when coupling is present. We divide the force we get by the
+    # one we'd get in the absence of taking it into account.
+    # pre_factor = - (1.0 / 3.0) * np.sqrt(2) / a
+    pre_factor = -(1.0 / 3.0) * np.sqrt(2) / a
+    coupling1, coupling2 = 0, 0
+
+    for n in np.arange(1, summands + 1):
+        k = calculate_k(n, a)
+        delta = calculate_delta(n, alpha, beta)
+        an = calculate_an(n, k, alpha, beta, delta)
+        bn = calculate_bn(n, k, alpha, beta, delta)
+        cn = calculate_cn(n, k, alpha, beta, delta)
+        dn = calculate_dn(n, k, alpha, beta, delta)
+
+        coupling1 += (2 * n + 1) * (an + bn + cn + dn)
+        coupling2 += (2 * n + 1) * (an - bn + cn - dn)
+
+    return pre_factor * coupling1 / radius1, pre_factor * coupling2 / radius2

lumicks/pylake/force_calibration/tests/test_drag.py

@@ -99,3 +99,18 @@ def test_coupling_solution(distance, r1, r2):
     np.testing.assert_allclose(eq2a, eq2b)
     np.testing.assert_allclose(eq3a, eq3b)
     np.testing.assert_allclose(eq4a, eq4b)
+
+
+@pytest.mark.parametrize(
+    "distance, radius1, radius2, ref_factor1, ref_factor2",
+    [
+        [3.0e-6, 0.5e-6, 0.5e-6, 0.8047215852074429, 0.8047215852074429],
+        [6.0e-6, 0.5e-6, 0.8e-6, 0.8224250265105119, 0.8982938795749217],
+        [6.0e-6, 0.9e-6, 0.8e-6, 0.8402922788501829, 0.8147951224475662],
+        [2.0e-6, 1.0e-6, 0.9999e-6, 0.01958146863802238, 0.019580489165418463],
+    ],
+)
+def test_coupling_factors(distance, radius1, radius2, ref_factor1, ref_factor2):
+    f1, f2 = coupling_correction_factor(radius1, radius2, distance, 5)
+    np.testing.assert_allclose(f1, ref_factor1)
+    np.testing.assert_allclose(f2, ref_factor2)