# Applying impulses simultaneously at each contact point or sequentially?

I'm coding a 2d physics engine in python, and I'm struggling to understand the right way to implement collision resolution when there are multiple contact points.

Consider a very simple collision case between 2 axis-aligned bounding boxes with the same mass. The collision will end up having two contact points, each corresponding to an endpoint of the collision manifold.

As far as I can tell, there are two ways of doing collision resolution:

1. Apply an impulse to the first contact point, updating the velocities, and then apply the second impulse.
2. Compute both impulses independently and store them in an accumulator, and apply them simultaneously.

Here's code for both options: option 1 corresponds to:

for arbiter in collisions:
# arbiter is (body1, body2, normal, penetration, contact_points)
for contact in arbiter[-1]:



and option 2 is roughly:

for arbiter in collisions:
# arbiter is (body1, body2, normal, penetration, contact_points)
accum_vx = np.array([0.0, 0.0])
accum_vy = np.array([0.0, 0.0])
accum_ax = 0
accum_ay = 0
for contact in arbiter[-1]:

arbiter[0].velocity += accum_vx
arbiter[1].velocity += accum_vy

arbiter[0].angular_velocity += accum_ax
arbiter[1].angular_velocity += accum_ay


I'm unsure of which one is correct / better. Implementing the first option causes the angular velocities to not cancel out correctly, meaning that a full on collision causes the boxes to spin away - which definitely should not be happening.

For the second option, I'm not sure how to handle friction impulses. Compute / apply them after both impulses are applied? Or compute them along with the collision impulses and apply them also at the same time?

Here's the collision response code:

def resolve_velocities(x: Body, y: Body, normal, penetration: float, contact):
ra = contact - x.pos
rb = contact - y.pos

v_xy = y.velocity + cross_scalar(y.angular_velocity, rb) - x.velocity - cross_scalar(x.angular_velocity, ra)
vn = np.dot(v_xy, normal)

if vn > 0:
return (0, 0), (0, 0)

e = min(x.restitution, y.restitution)

ra_dot_n = np.dot(ra, normal)
rb_dot_n = np.dot(rb, normal)

effective_mass = (x.inv_mass + y.inv_mass + ra_dot_n * ra_dot_n * x.inv_inertia + rb_dot_n * rb_dot_n * y.inv_inertia)

j = -(1 + e) * vn / effective_mass
impulse = j * normal

return (-impulse * x.inv_mass, impulse * y.inv_mass), (-x.inv_inertia * cross(ra, impulse), y.inv_inertia * cross(rb, impulse))
$$$$
`