Fix Error in Ball-Ball Frictional Inelastic No-Slip Velocity Calculat…

…ion (Again) (#177)

* clean up ball_ball/frictional_inelastic

* fix frictional inelastic no-slip case (hopefully for the last time)

* add tests for ball-ball frictional inelastic no-slip fix

pooltool/physics/resolve/ball_ball/frictional_inelastic/__init__.py

@@ -22,28 +22,17 @@ def _resolve_ball_ball(rvw1, rvw2, R, u_b, e_b):
     rvw2[1] = ptmath.coordinate_rotation(rvw2[1], -theta)
     rvw2[2] = ptmath.coordinate_rotation(rvw2[2], -theta)
 
-    rvw1_f = rvw1.copy()
-    rvw2_f = rvw2.copy()
-
     # velocity normal component, same for both slip and no-slip after collison cases
     v1_n_f = 0.5 * ((1.0 - e_b) * rvw1[1][0] + (1.0 + e_b) * rvw2[1][0])
     v2_n_f = 0.5 * ((1.0 + e_b) * rvw1[1][0] + (1.0 - e_b) * rvw2[1][0])
     D_v_n_magnitude = abs(v2_n_f - v1_n_f)
 
-    # angular velocity normal component, unchanged
-    w1_n_f = rvw1[2][0]
-    w2_n_f = rvw2[2][0]
-
-    # discard normal components for now
+    # discard velocity normal components for now
     # so that surface velocities are tangent
     rvw1[1][0] = 0.0
-    rvw1[2][0] = 0.0
     rvw2[1][0] = 0.0
-    rvw2[2][0] = 0.0
-    rvw1_f[1][0] = 0.0
-    rvw1_f[2][0] = 0.0
-    rvw2_f[1][0] = 0.0
-    rvw2_f[2][0] = 0.0
+    rvw1_f = rvw1.copy()
+    rvw2_f = rvw2.copy()
 
     v1_c = ptmath.surface_velocity(rvw1, unit_x, R)
     v2_c = ptmath.surface_velocity(rvw2, -unit_x, R)
@@ -71,12 +60,11 @@ def _resolve_ball_ball(rvw1, rvw2, R, u_b, e_b):
     # then slip condition is invalid so we need to calculate no-slip condition
     if not has_relative_velocity or np.dot(v12_c, v12_c_slip) <= 0:  # type: ignore
         # velocity tangent component for no-slip condition
-        D_v1_t = -(1.0 / 9.0) * (
-            2.0 * (rvw1[1] - rvw2[1])
-            + R * ptmath.cross(2.0 * rvw1[2] + 7.0 * rvw2[2], unit_x)
+        D_v1_t = -(1.0 / 7.0) * (
+            rvw1[1] - rvw2[1] + R * ptmath.cross(rvw1[2] + rvw2[2], unit_x)
         )
-        D_w1 = (5.0 / 9.0) * (
-            rvw2[2] - rvw1[2] + ptmath.cross(unit_x, rvw2[1] - rvw1[1]) / R
+        D_w1 = -(5.0 / 14.0) * (
+            ptmath.cross(unit_x, rvw1[1] - rvw2[1]) / R + rvw1[2] + rvw2[2]
         )
         rvw1_f[1] = rvw1[1] + D_v1_t
         rvw1_f[2] = rvw1[2] + D_w1
@@ -86,8 +74,6 @@ def _resolve_ball_ball(rvw1, rvw2, R, u_b, e_b):
     # reintroduce the final normal components
     rvw1_f[1][0] = v1_n_f
     rvw2_f[1][0] = v2_n_f
-    rvw1_f[2][0] = w1_n_f
-    rvw2_f[2][0] = w2_n_f
 
     # rotate everything back to the original frame
     rvw1_f[1] = ptmath.coordinate_rotation(rvw1_f[1], theta)

tests/physics/resolve/ball_ball/test_ball_ball.py

@@ -1,3 +1,4 @@
+import math
 from typing import Tuple
 
 import attrs
@@ -12,27 +13,55 @@
 from pooltool.physics.resolve.ball_ball.frictionless_elastic import FrictionlessElastic
 
 
-def head_on() -> Tuple[Ball, Ball]:
-    cb = Ball.create("cue", xy=(0, 0))
+def velocity_from_speed_and_xy_direction(speed: float, angle_radians: float):
+    """Convert speed and angle to a velocity vector
+
+    Angle is defined CCW from the x-axis in the xy-plane
+    """
+    return speed * np.array([math.cos(angle_radians), math.sin(angle_radians), 0.0])
+
+
+def gearing_z_spin_for_incoming_ball(incoming_ball):
+    """Calculate the amount of sidespin (z-axis spin) required for gearing contact
+    with no relative surface velocity.
+
+    Assumes line of centers from incoming ball to object ball is along the x-axis.
+
+    In order for gearing contact to occur, the sidespin must cancel out any
+    velocity in the tangential (y-axis) direction.
+
+    And from angular velocity equations where
+        'r' is distance from the rotation center
+        'w' is angular velocity
+        'v' is tangential velocity at a distance from the rotation center
+
+        v = w * R
+
+    So, v_y + w_z * R = 0, and therefore w_z = -v_y / R
+    """
+    return -incoming_ball.vel[1] / incoming_ball.params.R
 
-    # Cue ball makes head-on collision with object ball at 1 m/s in +x direction
-    cb.state.rvw[1] = np.array([1, 0, 0])
 
+def ball_collision() -> Tuple[Ball, Ball]:
+    cb = Ball.create("cue", xy=(0, 0))
     ob = Ball.create("cue", xy=(2 * cb.params.R, 0))
     assert cb.params.m == ob.params.m, "Balls expected to be equal mass"
     return cb, ob
 
 
-def translating_head_on() -> Tuple[Ball, Ball]:
-    cb = Ball.create("cue", xy=(0, 0))
-    ob = Ball.create("cue", xy=(2 * cb.params.R, 0))
+def head_on() -> Tuple[Ball, Ball]:
+    cb, ob = ball_collision()
+    # Cue ball makes head-on collision with object ball at 1 m/s in +x direction
+    cb.state.rvw[1] = np.array([1, 0, 0])
+    return cb, ob
+
 
+def translating_head_on() -> Tuple[Ball, Ball]:
+    cb, ob = ball_collision()
     # Cue ball makes head-on collision with object ball at 1 m/s in +x direction
     # while both balls move together at 1 m/s in +y direction
     cb.state.rvw[1] = np.array([1, 1, 0])
     ob.state.rvw[1] = np.array([0, 1, 0])
-
-    assert cb.params.m == ob.params.m, "Balls expected to be equal mass"
     return cb, ob
 
 
@@ -90,7 +119,7 @@ def test_translating_head_on_zero_spin_inelastic(
     cb_f, ob_f = model.resolve(cb_i, ob_i, inplace=False)
 
     # Balls should still be moving together in +y direction
-    assert np.isclose(cb_f.vel[1], ob_f.vel[1], atol=1e-10)
+    assert abs(cb_f.vel[1] - ob_f.vel[1]) < 1e-10
 
 
 @pytest.mark.parametrize("model", [FrictionalInelastic(), FrictionalMathavan()])
@@ -111,3 +140,93 @@ def test_head_on_z_spin(model: BallBallCollisionStrategy, cb_wz_i: float):
     cb_wz_f = cb_f.state.rvw[2][2]
     assert cb_wz_f > 0, "Spin direction shouldn't reverse"
     assert cb_wz_f < cb_wz_i, "Spin should be decay"
+
+
+@pytest.mark.parametrize(
+    "model", [FrictionalInelastic(), FrictionalMathavan(num_iterations=int(1e5))]
+)
+@pytest.mark.parametrize("speed", np.logspace(-1, 1, 5))
+@pytest.mark.parametrize(
+    "cut_angle_radians", np.linspace(0, math.pi / 2.0, 8, endpoint=False)
+)
+def test_gearing_z_spin(
+    model: BallBallCollisionStrategy, speed: float, cut_angle_radians: float
+):
+    """Ensure that a gearing collision causes no throw or induced spin.
+
+    A gearing collision is one where the relative surface speed between the balls is 0.
+    In other words, the velocity of each ball at the contact point is identical, and there is no
+    slip at the contact point.
+    """
+
+    unit_x = np.array([1.0, 0.0, 0.0])
+    cb_i, ob_i = ball_collision()
+
+    cb_i.state.rvw[1] = velocity_from_speed_and_xy_direction(speed, cut_angle_radians)
+    cb_i.state.rvw[2][2] = gearing_z_spin_for_incoming_ball(cb_i)
+
+    # sanity check the initial conditions
+    v_c = (
+        ptmath.surface_velocity(
+            cb_i.state.rvw, np.array([1.0, 0.0, 0.0]), cb_i.params.R
+        )
+        - cb_i.vel[0] * unit_x
+    )
+    assert ptmath.norm3d(v_c) < 1e-10, "Relative surface contact speed should be zero"
+
+    cb_f, ob_f = model.resolve(cb_i, ob_i, inplace=False)
+
+    assert (
+        abs(math.atan2(ob_f.vel[1], ob_f.vel[0])) < 1e-3
+    ), "Gearing english shouldn't cause throw"
+    assert abs(ob_f.avel[2]) < 1e-3, "Gearing english shouldn't cause induced side-spin"
+
+
+@pytest.mark.parametrize("model", [FrictionalInelastic()])
+@pytest.mark.parametrize("speed", np.logspace(0, 1, 5))
+@pytest.mark.parametrize(
+    "cut_angle_radians", np.linspace(0, math.pi / 2.0, 8, endpoint=False)
+)
+@pytest.mark.parametrize("relative_surface_speed", np.linspace(0, 0.05, 5))
+def test_low_relative_surface_velocity(
+    model: BallBallCollisionStrategy,
+    speed: float,
+    cut_angle_radians: float,
+    relative_surface_speed: float,
+):
+    """Ensure that collisions with a "small" relative surface velocity end with 0 relative surface velocity.
+    In other words, that the balls are gearing after the collision.
+
+    Note that how small the initial relative surface velocity needs to be for this condition to be met is dependent
+    on model parameters and initial conditions such as ball-ball friction and the collision speed along the line of centers.
+    """
+
+    unit_x = np.array([1.0, 0.0, 0.0])
+    cb_i, ob_i = ball_collision()
+
+    cb_i.state.rvw[1] = velocity_from_speed_and_xy_direction(speed, cut_angle_radians)
+    cb_i.state.rvw[2][2] = gearing_z_spin_for_incoming_ball(cb_i)
+    cb_i.state.rvw[2][2] += (
+        relative_surface_speed / cb_i.params.R
+    )  # from v = w * R -> w = v / R
+
+    # sanity check the initial conditions
+    v_c = (
+        ptmath.surface_velocity(cb_i.state.rvw, unit_x, cb_i.params.R)
+        - cb_i.vel[0] * unit_x
+    )
+    assert (
+        abs(relative_surface_speed - ptmath.norm3d(v_c)) < 1e-10
+    ), f"Relative surface contact speed should be {relative_surface_speed}"
+
+    cb_f, ob_f = model.resolve(cb_i, ob_i, inplace=False)
+
+    cb_v_c_f = ptmath.surface_velocity(
+        cb_f.state.rvw, unit_x, cb_f.params.R
+    ) - np.array([cb_f.vel[0], 0.0, 0.0])
+    ob_v_c_f = ptmath.surface_velocity(
+        ob_f.state.rvw, -unit_x, ob_f.params.R
+    ) - np.array([ob_f.vel[0], 0.0, 0.0])
+    assert (
+        ptmath.norm3d(cb_v_c_f - ob_v_c_f) < 1e-3
+    ), "Final relative contact velocity should be zero"