From 9069d7f24917a044f37d77135eb0f46671d89882 Mon Sep 17 00:00:00 2001 From: Derek McBlane Date: Tue, 7 Jan 2025 11:43:26 -0500 Subject: [PATCH] Fix Error in Ball-Ball Frictional Inelastic No-Slip Velocity Calculation (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 --- .../frictional_inelastic/__init__.py | 28 +--- .../resolve/ball_ball/test_ball_ball.py | 139 ++++++++++++++++-- 2 files changed, 136 insertions(+), 31 deletions(-) diff --git a/pooltool/physics/resolve/ball_ball/frictional_inelastic/__init__.py b/pooltool/physics/resolve/ball_ball/frictional_inelastic/__init__.py index 0f14f6fe..342e35c7 100644 --- a/pooltool/physics/resolve/ball_ball/frictional_inelastic/__init__.py +++ b/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) diff --git a/tests/physics/resolve/ball_ball/test_ball_ball.py b/tests/physics/resolve/ball_ball/test_ball_ball.py index 2e843916..2623d291 100644 --- a/tests/physics/resolve/ball_ball/test_ball_ball.py +++ b/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"