After the collision pair is chosen, their post-collision velocities, v i ∗ ∗ and v j ∗ ∗ , are evaluated. Writing the relative velocity in terms of spherical angles, θ θ and ϕ ϕ v r ∗ ∗ = v r [ ( sin ⁡ ⁡ θ θ cos ⁡ ⁡ ϕ ϕ ) x ^ ^ + ( sin ⁡ ⁡ θ θ sin ⁡ ⁡ ϕ ϕ ) y ^ ^ + cos ⁡ ⁡ θ θ z ^ ^ ] these angles are selected by a Monte Carlo process with distributions given by the collision model. For the hard spheres model these angles are uniformly distributed over the unit sphere. The azimuthal angle is uniformly distributed between 0 and 2 π π , so it is selected as ϕ ϕ = 2 π π ℜ ℜ 1 where ℜ ℜ 1 is a uniform deviate in [0, 1). The polar angle is distributed according to the probability density, P θ θ ( θ θ ) d θ θ = 1 2 sin ⁡ ⁡ θ θ d θ θ Using the change of variable q = − − cos ⁡ ⁡ θ θ , we have P q ( q ) d q = ( 1 2 ) d q so − − cos ⁡ ⁡ θ θ = q a n d sin ⁡ ⁡ θ θ = 1 − − q 2 w h e r e q = 2 ℜ ℜ 2 − − 1 The post-collision velocities are set as v i ∗ ∗ = v c m ∗ ∗ + 1 2 v r ∗ ∗ v j ∗ ∗ = v c m ∗ ∗ − − 1 2 v r ∗ ∗ Note that by conservation of linear momentum and energy the center of mass velocity and the relative speed are unchanged in a collision. That is, v c m = 1 2 ( v i + v j ) = 1 2 ( v i ∗ ∗ + v j ∗ ∗ ) = v c m ∗ ∗ and v r = | v i − − v j | = | v i ∗ ∗ − − v j ∗ ∗ | = v r ∗ ∗ This process is repeated for every pair of colliding particles.