So in this post I will talk about an interesting application of priority queues. This is actually representative of a whole family of applications and the coming paragraphs can be applied to other fields as well, not only to physics simulation that I will discuss below.

In the first part of this post I will give a general idea about the physics model behind the simulation, and I will outline the algorithm that we will use to run the simulation. In the coming posts I will continue then by showing the implementations behind the fundamental abstractions in the event simulation.

So there are many situations when we want to model a large number of entities which are continuously interacting with each other. Such examples can be gas particles confined to a container, cells of a tissue, stars in a galaxy and many others.

The inherent difficulty in modeling these objects is that their interaction is dynamic and the number of events can grow large. So finding a good model is fundamental to grasp the complexity of the problem.

There are different approaches to handle this complexity. But I will focus on event driven modeling in this chapter (one other major approach is time sliced modeling).

So the basic problem that we are trying to solve is that we have N particles which are constrained to move between the limits of a given rectangular field. And while these particles are moving they can collide with each other and/or the wall and they can bounce off from each other and/or the wall as well.

This is a well defined problem and simulating is easy enough, because by simplifying some of the assumptions about the physical interaction of the particles, we can build a quite sophisticated model using some easy building blocks.

So first of all I will briefly describe the assumptions that we make about the actors in our model.

Our model will assume that the actors are following the hard sphere model:

- N particles are in continuous motion and they are confined to the unit box
- We know the position ($rx_i$, $ry_i$) of each particle
- We know the velocity ($vx_i$, $vy_i$) of each particle
- We know the mass $m_i$ of each particle
- And we know the radius $\sigma_i$ of each particle
- The particles are bouncing off flexibly of each other and their boundaries
- The particles are traveling in straight lines by constant speed

Having these assumptions we can basically use our secondary school physics knowledge to calculate the velocity and position of the given particles.

Then the physics calculation will be basically performed in two steps:

Step #1: Prediction This step is the place where we take the participating particles and we calculate their interactions with all other particles pairwise. The outcome of this prediction is basically the time when the next collision of a given particle will happen.

Step #2: Resolution This step will update the velocity and position of those particles which are participating in the next collision.

But this itself is not enough since we also need a driving force, a controller which will synchronize the predictions and resolutions taking place. And this is where priority queues will help us.

The basic idea is the following: put all the prediction results (the time of next collision) on a priority queue, and keep pulling off events which are based on these time results. Basically we will pull off the next smallest time in each iteration. Then after pulling off the events we can calculate the resolution and run new predictions based on the resolution results, putting these new predictions on the priority queue again.

We keep doing these iterations until we find new events in the queue, or until we reach a given number of iterations. Simple enough.

So let's see how can we implement this. First of all let's take a look at the prediction. We will have to predict three things:

- The time the particles hit each other
- The time the particles hit the vertical wall
- The time the particles hit the horizontal wall

The time the particles hit each other:

```
Particle.prototype.timeToHit = function(that) {
if (this == that) return INFINITY;
var dx = that.rx - this.rx;
var dy = that.ry - this.ry;
var dvx = that.vx - this.vx;
var dvy = that.vy - this.vy;
var dvdr = dx*dvx + dy*dvy;
if (dvdr > 0) return INFINITY;
var dvdv = dvx*dvx + dvy*dvy;
var drdr = dx*dx + dy*dy;
var sigma = this.radius + that.radius;
var d = (dvdr*dvdr) - dvdv * (drdr - sigma*sigma);
// if (drdr < sigma*sigma) console.log("overlapping particles");
if (d < 0) return INFINITY;
return -(dvdr + Math.sqrt(d)) / dvdv;
};
```

The time the particles hit the vertical wall:

```
Particle.prototype.timeToHitVerticalWall = function() {
if (this.vx > 0) return (this.screenW - this.rx - this.radius) / this.vx;
else if (this.vx < 0) return (this.radius - this.rx) / this.vx;
else return INFINITY;
};
```

The time the particles hit the horizontal wall:

```
Particle.prototype.timeToHitHorizontalWall = function() {
if (this.vy > 0) return (this.screenH - this.ry - this.radius) / this.vy;
else if (this.vy < 0) return (this.radius - this.ry) / this.vy;
else return INFINITY;
};
```

Let's break this down. So what's going on? Take a look at the image above first.

The image above contains 2 particles in their initial position (light colored) and the same 2 particles later when their position changed (dark colored). The right particle's initial position is ($rx_i$, $ry_i$) and it's velocity is ($vx_i$, $vy_i$). The other particle's initial position is ($rx_j$, $ry_j$) and it's initial velocity is ($vx_j$, $vy_j$).

The moved particles' position is ($rx_i^\prime$, $ry_i^\prime$) and ($rx_j^\prime$, $ry_j^\prime$) respectively.

Now we can see that when the two particles collide, then the distance of the center of their diameter is $$ \sigma = \sigma_j + \sigma_i $$

Which can be expressed as: $$ \sigma^2 = (rx_i^\prime - rx_j^\prime)^2 + (ry_i^\prime - ry_j^\prime)^2 $$

Note one more thing here though: this last equation basically says that the distance of those particles equals the magnitude of the difference between the two position vectors, that is: $$ |\Delta r| = \sigma $$

This last observation will be important when we derivate the impulse of the particles during collision.

Before the particles collide they move along a straight line getting closer and closer to each other, and their new position is: $$ rx_i^\prime = rx_i + \Delta t \cdot vx_i $$ $$ ry_i^\prime = ry_i + \Delta t \cdot vy_i $$ $$ rx_j^\prime = rx_j + \Delta t \cdot vx_j $$ $$ ry_j^\prime = ry_j + \Delta t \cdot vy_j $$

Substituting these 4 equations into the previous one we get a quadratic equation. Then if we solve the quadratic equation for $\Delta t$, we get the following.

If $\Delta v \cdot \Delta r > 0$, then the quadratic equation has no solution for $\Delta t > 0\:(\Delta t = \infty)$.

If $\Delta d < 0$ then again the quadratic equation has no solution for $\Delta t > 0\:(\Delta t = \infty)$.

But in any other cases $$ \Delta t = -\frac{\Delta v \cdot \Delta r + \sqrt d}{\Delta v \cdot \Delta v} $$ where: $$ d = (\Delta v \cdot \Delta r)^2 - (\Delta v \cdot \Delta v)(\Delta r \cdot \Delta r - \sigma^2) $$ $$ \Delta r = (\Delta rx, \Delta ry) = (rx_j - rx_i, ry_j - ry_i) $$ $$ \Delta v = (\Delta vx, \Delta vy) = (vx_j - vx_i, vy_j - vy_i) $$ $$ \Delta r \cdot \Delta r = (\Delta rx)^2 + (\Delta ry)^2 $$ $$ \Delta v \cdot \Delta v = (\Delta vx)^2 + (\Delta vy)^2 $$ $$ \Delta v \cdot \Delta r = (\Delta vx \cdot \Delta rx) + (\Delta vy \cdot \Delta ry) $$

Note here that $\Delta v$ is basically the relative speed of the observed particles. This will again be important when we derivate the impulse of the particles.

Now if you take a look you can see that these are the same equations that we used to predict the time when the particles will hit each other.

Ok. And what about the collision resolution? Well colliding with the wall is easy. When the particle hits the left or right wall then the new velocity is ($-vx$, $vy$). If it hits the top or bottom walls, then the new velocity is ($vx$, $-vy$).

Colliding with other particles is a bit more complex. But if we reduce our physics model to a simplified version of collision, then it's still doable with simple arithmetics.

So the important thing is that we assume our particles are hard and perfectly elastic, and there is neither friction nor spin. This will simplify a few things and we can make assumptions, for example that the linear and angular momentum will be conserved (so we can apply Newton's second law).

To formalize this we can write: $$ m_i v_i + m_j v_j = m_i v_i^\prime + m_j v_j^\prime \:\:(\star) $$

Just like kinetic energy (which is not true in the general case) because an elastic collision between two objects is one in which the total kinetic energy (as well as total momentum) of the system is the same before and after the collision.

So let's formalize this with an equation. The conservation of kinetic energy means that: $$ m_i(v_i^2 - {v_i^\prime}^2)\hat{n} = m_j(v_j^2 - {v_j^\prime}^2)\hat{n} $$ If we factor both sides then we get: $$ m_i(v_i - v_i^\prime)(v_i + v_i^\prime)\hat{n} = m_j(v_j - v_j^\prime)(v_j + v_j^\prime)\hat{n} $$ Next if we separate the terms containing $m_i$ and $m_j$ in ($\star$) we get: $$ m_i(v_i - v_i^\prime) = m_j(v_j - v_j^\prime) $$ Now if we divide the factored kinetic energy equation with this last equation we get: $$ (v_i + v_i^\prime)\hat{n} = (v_j + v_j^\prime)\hat{n} $$ Or the equivalent: $$ (v_i - v_j)\hat{n} = -(v_i^\prime - v_j^\prime)\hat{n} \:\:(\star\star) $$

And one more thing finally: assuming no friction or spin, the normal force acts along the line connecting their centers. And the impulse as well.

So if we want to express the velocity after the collision we can write respectively: $$ v_i^\prime = v_i - \frac{J}{m_i}\cdot \hat{n} $$ $$ v_j^\prime = v_j + \frac{J}{m_j}\cdot \hat{n} $$

Now if we substitute these equations into (*\star\star*):
$$
(v_j + \frac{J}{m_j}\hat{n} - v_i + \frac{J}{m_i}\hat{n})\hat{n} = -(v_j - v_i)\hat{n}
$$

Now if we divide by $\hat{n}$ and factoring out the inner $\hat{n}$ we get: $$ (v_j + \frac{J}{m_j} - v_i + \frac{J}{m_i})\hat{n} = -(v_j - v_i) $$

Which is: $$ v_j + \frac{J}{m_j}\hat{n} - v_i + \frac{J}{m_i}\hat{n} = v_i - v_j $$ And then: $$ 2v_j - 2v_i = -\frac{J}{m_j}\hat{n} - \frac{J}{m_i}\hat{n} $$ $$ 2v_j - 2v_i = -(\frac{J}{m_j} + \frac{J}{m_i})\hat{n} $$ $$ 2v_i - 2v_j = \frac{J}{m_j}\hat{n} + \frac{J}{m_i}\hat{n} $$ $$ 2 m_j m_i(v_i - v_j) = m_i J \hat{n} + m_j J \hat{n} $$ $$ \frac{2 m_j m_i(v_i - v_j)}{m_i + m_j} = J\hat{n} $$ $$ \frac{2 m_j m_i(v_i - v_j)}{m_i + m_j}\cdot\frac{1}{\hat{n}} = J $$

Now notice that: $$ \hat{n}^-1 = \frac{\hat{n}}{|\hat{n}|^2} = \hat{n} $$

So we can write: $$ \frac{2 m_j m_i\Delta v\hat{n}}{m_i + m_j} = J $$

And we also know that: $$ \frac{\Delta r}{\sigma} = \hat{n} \:\:because\:\: |\Delta r| = \sigma $$

So we get: $$ \frac{2 m_j m_i\Delta v\Delta r}{\sigma(m_i + m_j)} = J $$

Now we can decompose the impulse into the following $x$ and $y$ components: $$ J_x = J \cdot \Delta rx / \sigma $$ because: $$ J_x = J \cdot \cos(\theta) $$ And $$ J_y = J \cdot \Delta ry / \sigma $$ because: $$ J_y = J \cdot \sin(\theta) $$

where J is the earlier computed impulse and ($J_x$, $J_y$) is the impulse, $m_i$ and $m_j$ are masses of the particles, and the rest is the same as above.

Once we know the impulse we can apply Newton's 2nd law to compute the velocities of the particles $$ vx_i^\prime = vx_i + \frac{J_x}{m_i} $$ $$ vx_j^\prime = vx_j - \frac{J_x}{m_j} $$ $$ vy_i^\prime = vy_i + \frac{J_y}{m_i} $$ $$ vy_j^\prime = vy_j - \frac{J_y}{m_j} $$

The code basically reflects these equations in a condensed form. We will have to compute three things then:

- The speed after the particles hit each other
- The speed after the particles hit the vertical wall
- The speed after the particles hit the horizontal wall

The speed after the particles hit each other:

```
Particle.prototype.bounceOff = function(that) {
var dx = that.rx - this.rx;
var dy = that.ry - this.ry;
var dvx = that.vx - this.vx;
var dvy = that.vy - this.vy;
var dvdr = dx*dvx + dy*dvy; // dv dot dr
var dist = this.radius + that.radius; // distance between particle centers at collison
// normal force F, and in x and y directions
var F = 2 * this.mass * that.mass * dvdr / ((this.mass + that.mass) * dist);
var fx = F * dx / dist;
var fy = F * dy / dist;
// update velocities according to normal force
this.vx += fx / this.mass;
this.vy += fy / this.mass;
that.vx -= fx / that.mass;
that.vy -= fy / that.mass;
// update collision counts
this.cnt++;
that.cnt++;
};
```

The speed after the particles hit the vertical wall:

```
Particle.prototype.bounceOffVerticalWall = function() {
this.vx = -this.vx;
this.cnt++;
};
```

The speed after the particles hit the horizontal wall:

```
Particle.prototype.bounceOffHorizontalWall = function() {
this.vy = -this.vy;
this.cnt++;
};
```

So now that we know how to model the physical interaction, we can start building our collision system. This is where we will continue in the coming post.