bachelor-project-report/sections/theoretical_background.tex

\section{Theoretical Background}
\label{sec:theory}
The theoretical background is everything related to the physics part of the
project. It covers the calculation of the inertia of different types of polygons;
different algorithms to detect whether there is a collision between
any two polygons; and the resolution of the collision, that is, finding the final
velocity vectors and angular speed of those polygons.

\subsection{Moment of Inertia}
\label{sub:moment}

The inertia of an object refers to the tendency of an object to resist
a change of its state of motion or rest, it describes how the object
behaves when forces are applied to it. In particular, the mass of a
point object (inertial mass) determines the inertia of that objects.
Thus an object with a large mass requires more force to change its
motion, either to make it move from an initial state at rest or to
stop it if it is already moving. On the other hand, an object with a
smaller mass is easier to set in motion or to bring to a halt.

For objects that are not point-like, the motion includes a rotational
component.

The moment of inertia is analogous to the mass in determining the
inertia of an object with respect to its rotational motion. In other
words, the moment of inertia determines the rotational inertia of an
object, that is, the object's resistance to changes in its rotational
motion.  The moment of inertia of an object is therefore computed with
respect to a rotation axis, and it is determined by the way the mass
of the object is distributed with respect to the chosen axis of
rotation.


For the purpose of this project, we are interested in a 2D simulation.
This means that the axis of rotation is always parallel to the
$z$-axis, perpendicular to the plane of the simulation.  Furthermore,
we model the motion of an object using the barycenter of an object as
its reference point.  We therefore compute the moment of inertia of a
polygon with respect to a rotation axis that is perpendicular to the
2D plane and that goes through the center of mass of the polygon.

In general, the moment of inertia about a rotation axis $Z$ of a point
mass $m$ at distance $r$ from $Z$ is $I_m=mr^2$.  Therefore, for a
generic 2D object $Q$, we integrate over the whole 2D object:
\begin{equation}
	\label{eq:moment_general}
	I_Q = \int_{\mathcal{A}}\|\vec r\|^2 \rho(\vec r) \diff \mathcal{A}
\end{equation}
where $\rho$ is the density of object $Q$ (mass per unit area) in the
point $\vec r$ from the across the small pieces of area $\mathcal A$
of the object.

In our case, we use 2D Cartesian coordinates centered in the point
(axis) about which we compute he moment of inertia.  Furthermore, we
consider polygons of uniform density.  We can therefore compute:
\begin{equation}
	\label{eq:moment}
	I_Q = \rho \iint(x^2 + y^2)\diff x\diff y
\end{equation}

The bounds of the integral depend on the shape of the polygon. In the following
sections, we will describe how to compute those bounds, then we will show a
different technique to compute the moment of inertia of arbitrary polygons.

\subsubsection{Rectangle}
The moment of inertia of a rectangle of width $w$ and height $h$ with respect to
the axis of rotation that passes through its barycenter can be visualized in the
\figref{fig:rectangle_inertia}.

\begin{figure}[H]
	\centering
	\hfill
	\begin{subfigure}[]{.4\textwidth}
		\centering
		\inputtikz{rectangle_inertia2d}
		\caption{2d view of rectangle with axis of rotation}
		\label{fig:rectangle_inertia2d}
	\end{subfigure}
	\hfill
	\begin{subfigure}[]{.4\textwidth}
		\centering
		\inputtikz{rectangle_inertia3d}
		\caption{3d view of rectangle with axis of rotation}
		\label{fig:rectangle_inertia3d}
	\end{subfigure}
	\hfill\null
	\caption{Representation of rectangle with respect to axis of rotation $z$}
	\label{fig:rectangle_inertia}
\end{figure}

As figure \figref{fig:rectangle_inertia2d} implies, the bounds of equation
\ref{eq:moment} are trivial to derive:
\begin{equation}
	\label{eq:rect_moment}
	I_{\text{rect}} = \rho \int_{-\frac{h}{2}}^{\frac{h}{2}}
	\int_{-\frac{w}{2}}^{\frac{w}{2}}  x^2 + y^2 \diff x \diff y
	= \frac{\rho wh}{12} \left( w^2 + h^2 \right)
\end{equation}
and since $\rho w h$ is the density of the rectangle multiplied by its area, we
can replace this term by its mass $m$, thus
\begin{equation}
	I_{\text{rect}} = \frac{1}{12} m\left(w^2 + h^2\right)
\end{equation}

All the steps to compute equation~\ref{eq:rect_moment} can be found in equation
\ref{eq:rect_moment_long} in Appendix \ref{appendix:calculations}.

\subsubsection{Regular Polygons}
\label{sub:regular_polygons}
A regular polygon is a shape that has sides of equal length and angles between
those sides of equal measure. A polygon of $n$ sides can be subdivided in $n$
congruent (and isosceles since they are all the radius of the circumscribing
circle) triangles that all meet in the polygon's barycenter,
as demonstrated in Figure \ref{fig:pentagon_triangles} with a pentagon.

\begin{figure}[H]
	\centering
	\hfill
	\begin{subfigure}[]{.45\textwidth}
		\centering
		\inputtikz[.6]{pentagon}
		\caption{Regular polygon of 5 sides with its barycenter}
		\label{fig:pentagon}
	\end{subfigure}
	\hfill
	\begin{subfigure}[]{.45\textwidth}
		\centering
		\inputtikz[.6]{pentagon_congruent}
		\caption{Pentagon divided in 5 congruent triangles}
		\label{fig:pentagon_triangles}
	\end{subfigure}
	\hfill\null
	\caption{Subdivision of regular polygons into congruent triangles}
	\label{fig:regular_poly}
\end{figure}

If we define one of the sub-triangle of the regular polygon as $T$, then we can
find the moment of inertia $I_T$ when it is rotating about the barycenter. To
find the bounds of the integral in equation \ref{eq:moment}, we can take the
triangle $T$ and place it along the $x$-axis so that it is symmetric likes shown
in figure. Assuming the side length of the polygon is $l$, the height of the
triangle $T$ is $h$ and the angle of the triangle on the barycenter of the
polygon to be $\theta$, then
\begin{figure}[H]
	\centering
	\inputtikz[.3]{isosceles}
	\caption{Sub-triangle $T$ of regular polygon}
	\label{fig:subtriangle}
\end{figure}
we can see the bounds for the integral
\begin{equation}
	\label{eq:subtriangle_moment}
	I_{T} = \rho \int_0^h\int_{-\frac{lx}{2h}}^{\frac{lx}{2h}}x^2 + y^2 \diff
	y\diff x= \frac{m_Tl^2}{24} \left(1 + 3\cot^2\left(\frac{\theta}{2}\right)\right)
\end{equation}

All the steps to compute equation~\ref{eq:subtriangle_moment} can be found in
equation \ref{eq:subtriangle_moment_long} in Appendix
\ref{appendix:calculations}.

Now that we have the moment of inertia of the sub-triangle, we can make the link
to the overall polygon. Since
$$ \theta = \frac{2\pi}{n} \implies \frac{\theta}{2} = \frac{\pi}{n} $$
and the moment of inertia are additive (as long they are as they are about the
same axis) we can get the moment of inertia with
$$ I_{\text{regular}} = n I_T $$
and since the mass of the regular polygon $m$ is the sum of the masses of the
sub-triangle
$$ m = n m_T $$
we have that
\begin{equation}
	\label{eq:regular_moment}
	I_{\text{regular}} = \frac{ml^2}{24} \left( 1 + 3\cot^2\left(\frac{\pi}{n}\right) \right)
\end{equation}

\subsubsection{Arbitrary Polygons}

For arbitrary polygons, we are taking a slightly different approach. Using the
Cartesian coordinate system to solve the equation \ref{eq:moment} revealed to be
more cumbersome than useful. But similarly to regular polygons (c.f. Section
\ref{sub:regular_polygons}), we can use the additive property of the moment
inertia to divide our arbitrary polygon into sub-triangles. As opposed to
regular polygons, these triangles won't be congruent, so we can't just get the
moment of inertia of one of them and multiply it by the number of sides, but we
need to calculate them individually. So given a polygon of $n$ sides, we can
construct $n$ sub-triangles $T_i$, for $i = 1, \dots, n$. So the moment of
inertia $I$ of the polygon will be
\begin{equation}
	I = \sum_i I_{T_i}
\end{equation}
\begin{figure}[H]
	\centering
	\begin{subfigure}[]{.5\textwidth}
		\centering
		\inputtikz[.7]{arbitrary}
		\caption{An arbitrary 6-sided polygon}
		\label{fig:arbitrary}
	\end{subfigure}
	\begin{subfigure}[]{.49\textwidth}
		\centering
		\inputtikz[.7]{arbitrary_divided}
		\caption{Arbitrary polygon divided into 6 sub-triangles}
		\label{fig:abitrary_divded}
	\end{subfigure}
\end{figure}

To calculate the moment of inertia $I_{T_i}$, instead of using the classical
$x$- and $y$-axis as we did before, we decided to use the edges of the triangle
as axis and therefore express what we need to integrate in function of those as
can be seen in Figure \ref{fig:abitrary_subtriangle}.
\begin{figure}[H]
	\centering
	\inputtikz[.7]{arbitrary_subtriangle}
	\caption{Sub-triangle of arbitrary polygon}
	\label{fig:abitrary_subtriangle}
\end{figure}

In Figure \ref{fig:abitrary_subtriangle}, $C$ represent the barycenter of the
polygon (as is shown in Figure \ref{fig:abitrary_divded}). The axis we are going
to integrate on are $\vv{CA}$ and $\vv{AB}$.
We can now define
\begin{equation} \label{eq:alpha}
	\vv{CP_1} = \alpha \vv{CA}, \qquad \vv{CP_2} = \alpha
	\vv{CB}, \qquad \forall \alpha \in [0, 1]
\end{equation}
and
$$ \vv{P_1Q} = \beta \vv{P_1P_2}, \qquad \forall \beta \in [0, 1] $$

From \ref{eq:alpha}, it quickly follows that
\[
	\vv{P_1P_2} = \alpha \vv{AB}
\]
therefore
\begin{equation} \label{eq:beta_alpha}
	\vv{P_1Q} = \beta \alpha \vv{AB}
\end{equation}
Finally, if we put together equations \ref{eq:alpha} and \ref{eq:beta_alpha}, we
have that
\begin{equation}\label{eq:r}
	\vec r = \vv{CP_1} + \vv{P_1Q} = \alpha \vv{CA} + \beta \alpha \vv{AB}
\end{equation}

Now we got the first part equation \ref{eq:moment_general}. To find the $\diff
	\mathcal A$, we
just need to get the area of the square that contains $Q$ in Figure
\ref{fig:abitrary_subtriangle}. Since $\|\vv{AB}\|$ represents the base of the
triangle $T_i$, we can define
$$ b = \| \vv{AB}\| $$
we consequently have that
\begin{equation} \label{eq:dA}
	\diff \mathcal{A} = b \alpha \diff \beta h\diff \alpha
\end{equation}
where $h = \| \vv{CH} \|$ is the height of triangle.
We can now assemble \ref{eq:r} and \ref{eq:dA}
\begin{equation}
	\label{eq:subtriangle_arbitrary_moment}
	I_{T_i} = \rho \int_0^1 \int_0^1 \vec r^2 hb \alpha  \diff \alpha
	\diff \beta = \frac{\rho h b}{4} \left(\frac{1}{3} \vv{AB}^2 +
	\vv{AB} \cdot \vv{CA} + \vv{CA}^2\right)
\end{equation}

Since $\frac{\rho h b}{2}$ is the mass of the triangle we can write the result
as
\begin{equation}
	I_{T_i} = \frac{m_{T_i}}{2} \left(\frac{1}{3} \vv{AB}^2 + \vv{AB} \cdot \vv{CA}
	+ \vv{CA}^2\right)
\end{equation}


All the steps to compute equation~\ref{eq:subtriangle_arbitrary_moment} can be
found in equation \ref{eq:subtriangle_arbitrary_moment_long} in Appendix
\ref{appendix:calculations}.

Now that we have the moment of inertia of the sub-triangle, we can make the link
to the overall polygon.
\begin{equation}
	I_{\text{arbitrary}} = \sum_i I_{T_i} = \sum_{i=1}^n \frac{m_{T_i}}{2}
	\left(\frac{1}{3} \vv{P_iP_{i+1}}^2 + \vv{CP_i} \cdot \vv{P_iP_{i+1}} +
	\vv{CP_i}^2\right)
\end{equation}
where, $P_{n+1} = P_1$ in the case of $i = n$.

\subsection{Collision detection}
\label{sub:collision-detection}

Collision detection, as the name suggests, are the algorithms used to detect
whether two polygons are colliding. The result of this procedure must be an
impact point and a normal vector, that will then be used for the  collision
resolution \ref{sub:resolution}.

\subsubsection{Separating Axis Theorem}

This algorithm was the first one studied for this project and was inspired by
the works of David Eberly \cite{convexcollisionsSAT}. The separating axis
theorem (SAT) states that if you can draw a line between two convex objects,
they do not overlap. We will call this line a \textit{separating line}. More
technically, two convex shapes do not overlap if there exists an axis onto which
the two objects' projections do not overlap. We'll call this axis a \textit{separating
	axis}. This concept can be visualized in Figure \ref{fig:SAT-intro}.

\begin{figure}[H]
	\centering
	\inputtikz[.7]{SAT_intro}
	\caption{SAT: Separating axis ($A$) vs non-separating axis ($B$), with
		separating line ($C$)}
	\label{fig:SAT-intro}
\end{figure}

As we can see in Figure \ref{fig:SAT-intro}, the axis $B$ show that the
projections of the both polygons overlap, but we were able to find an axis $A$
where this is not the case. As soon as we find an axis for which the projections
do not overlap, it means that the polygons are not colliding. For 2D objects, we
only need to consider the axes that are orthogonal to each edge. In Figure
\ref{fig:SAT-intro}, only two of those axes are shown for better readability,
but they would be 7, one for each edge.

To move (or push) one polygon away from the other, we also need to find a vector
that, when added to the polygons position, will make the shapes not overlap. We
want the minimum displacement possible, We'll call this vector the minimum push
vector (MPV). For 2-dimensional polygons, this vector will lie in some of the
orthogonal axes.

\begin{figure}[H]
	\centering
	\inputtikz[.7]{SAT_mvp}
	\caption{SAT: Minimum push vector $\vec v_i$ on axis defined by $\vec o_i$,
		orthogonal to edge $e_i$}
	\label{fig:SAT-mpv}
\end{figure}

The candidate MPVs $\vec v_i$ are the vectors that define the axis $\vec o_i$
(orthogonal to edge $e_i$), with $\| \vec o_i \| = 1$, multiplied by the minimum
overlap between the two polygons, as shown in \ref{fig:SAT-mpv}. The final MPV
is simply the $\vec v_i$ with the smallest norm.

\paragraph{Pitfalls of SAT} The issue with the SAT algorithm is that although it
is good to find whether two polygons are colliding and the MPV, it isn't trivial
to gather the point of impact, i.e. the vertex that is penetrating the other
polygon. It is doable, but during the implementation, it came with some caveats
that introduced some bugs, so we decided to switch strategy and go with an
algorithm of our own. Moreover, SAT only supports convex polygons, which limits
the original objective of the project, which was to have any arbitrary polygon.

\subsubsection{Vertex collisions}
\label{sub:vertex-collision}

The solution that was adopted for the project, after trying SAT, was a more
intuitive one, developed by Prof. Carzaniga. The idea is simple: check if a
vertex of a polygon is colliding with an edge of another polygon.

\begin{figure}[H]
	\centering
	\inputtikz[.5]{vertex_intro}
	\caption{Vertex-edge collision between polygons $A$ and $B$}
	\label{fig:vertex-edge}
\end{figure}

If we have a polygon defined as a set of points $P \subset \mathbb{R}^2$, we
define a vertex as a pair of segments $\left(\overline{V_{i-1}V_i},
	\overline{V_{i+1}V_i}\right)$. To check if vertex $V_i$ of polygon $A$ is inside
polygon $B$, we just check if the both segments $\overline{V_{i-1}V_i}$ and
$\overline{V_{i+1}V_i}$ intersect edge $e_j$. If such is the case, we know that
$V_i$ is inside, and we can use it as impact point. We can now take the vector
perpendicular to the edge $e_j$ and normalize it, which is the normal we need
for the collision resolution.


\paragraph{Edge cases} With this approach (and many other collision detection
algorithms), it is easy to see that the case described in the Figure
\ref{fig:vertex-edge} is only a general one. It will ultimately happen most
often, but there are especially two edge cases that are note-worthy and occurred
during the simulation a greater number of times than expected.

\paragraph{Parallel collision} Parallel collision occur when two edge collide
with each other, an example can be seen in Figure \ref{fig:vertex-parallel}.
This collision does not get detected by the vertex-edge detection method because
the edge of $A$ parallel to the edge of $B$ do not collide, since they are
parallel.

\begin{figure}[H]
	\centering
	\inputtikz[.5]{vertex_parallel}
	\caption{Parallel edges collision between polygons $A$ and $B$}
	\label{fig:vertex-parallel}
\end{figure}

To determine if polygon $A$ is having a parallel collision with polygon $B$ we
check if only the segment $\overline{V_{i-1}V_i}$ intersects with edge $e_j$ and
if the segment $\overline{V_{i+1}V_i}$ is parallel to edge $e_j$.

To find the impact point, we first need to find the minimum overlap between the
parallel edges. We can calculate it by projecting the points that make the
parallel edges of $A$ and $B$ on the axis generated by the edge of $A$ (c.f.
Figure \ref{fig:parallel-impact}).Finally, the impact point, is the midpoint of
the said overlap.

\begin{figure}[H]
	\centering
	\begin{subfigure}[]{.5\textwidth}
		\centering
		\inputtikz[.7]{vertex_parallel_full}
		\caption{Edge fully contained by other edge}
	\end{subfigure}
	\begin{subfigure}[]{.49\textwidth}
		\centering
		\inputtikz[.7]{vertex_parallel_partial}
		\caption{Edge partly contained by other edge}
	\end{subfigure}
	\caption{Parallel collision, finding the impact point}
	\label{fig:parallel-impact}
\end{figure}

The normal vector is given in the same way as the vertex-edge collision, i.e.
taking the perpendicular vector to~$e_j$ and normalizing it.

% TODO: ask Carza about the possibility of resolving vertex-vertex normal with SAT
\paragraph{Vertex on vertex collision} These collision happen when, at the
moment of the frame, two polygons are "inside each other", i.e. they both have
have a vertex present inside the area of the other polygon, as shown in Figure
\ref{fig:vertex-vertex}.

\begin{figure}[H]
	\centering
	\inputtikz[.5]{vertex_vertex}
	\caption{Vertex on vertex collision between polygons $A$ and $B$}
	\label{fig:vertex-vertex}
\end{figure}

\noindent
Two problems arise when trying to deal with this edge case:
\begin{enumerate}
	\item determining which of the two vertices we chose as the impact point;
	      \label{itm:impact}
	\item determining what normal vector.
	      \label{itm:normal}
\end{enumerate}

% \noindent
For point \ref{itm:impact}, it actually doesn't really matter. What is
represented in Figures~\ref{fig:vertex-edge}-\ref{fig:vertex-vertex} are an
over exaggeration of what happens in the engine. Since the time delta between
two frames is so small, the collisions look more like what is in
Figure~\ref{fig:real-vertex-vertex}, where the concerned vertices are located
practically in the same spot, and the difference is negligible. So for convenience
with the general vertex-edge case, we will take the vertex of $A$ being inside
$B$ as the impact point.

For point \ref{itm:normal}, that's where the real problem of the edge case
occur. We are unable to determine what vectors is supposed to be the normal
vector. One might say "just take the vector that goes from the vertex of $A$
(that is inside $B$) to the one of $B$ (that is inside $A$)", but it doesn't
work out in general. It works (kind of) when we have non-rotating objects going
straight at each other, but as soon as you have rotational motion, this
reasoning collapses. The solution that was decided (together with advisor Prof.
Carzaniga) was to treat the collision as a vertex-edge collision, and choosing
whatever the first edge of $B$ comes up first in the calculations as the edge to
find the normal. The results look realistic enough to be accepted.


\begin{figure}[H]
	\centering
	\begin{subfigure}[]{.33\textwidth}
		\centering
		\inputtikz[.7]{vertex_edge_real}
		\caption{Realistic vertex-edge collision}
	\end{subfigure}
	\begin{subfigure}[]{.32\textwidth}
		\centering
		\inputtikz[.7]{vertex_parallel_real}
		\caption{Realistic edge-edge collision}
	\end{subfigure}
	\begin{subfigure}[]{.34\textwidth}
		\centering
		\inputtikz[.7]{vertex_vertex_real}
		\caption{Realistic vertex-vertex collision}
		\label{fig:real-vertex-vertex}
	\end{subfigure}
	\caption{Realistic collisions}
	\label{fig:real-collisions}
\end{figure}

\subsection{Collision resolution}
\label{sub:resolution}

The collision resolution is the last step in the processing of the collision.
Algorithmically, it is much less heavy than collision detection, since, once the
simulation has two colliding polygons, a point of impact and a normal vector,
it's just a case of applying the rigid body physics formulas to the polygons
that are colliding. This part has been helped a lot by the works of Erik Neumann
\cite{collision:resolution-site} and Chris Hecker
\cite{collision:resolution-paper}.

\begin{figure}[H]
	\centering
	\inputtikz[.5]{collision_resolution}
	\caption{Collision resolution between polygons $A$ and $B$}
	\label{fig:collision_resolution}
\end{figure}

\paragraph{Variable definition} Before getting into any maths, let's define some
variables that we are going to use

\begin{multicols}{2}
	\begin{itemize}
		\item $m_a, m_b = $ mass of the bodies $A$ and $B$
		\item $\vec r_{ap} = $ distance vector from center of mass of body $A$ to
		      point $P$
		\item $\vec r_{bp} = $ distance vector from center of mass of body $A$ to
		      point $P$
		\item $\omega_{a1}, \omega_{b1} = $ initial angular velocity of bodies
		      $A, B$
		\item $\omega_{a2}, \omega_{b2} = $ final angular velocity of bodies
		      $A, B$
		\item $\vec v_{a1}, \vec v_{b1} =$ initial velocities of center of mass
		      bodies $A, B$
		\item $\vec v_{a2}, \vec v_{b2} =$ final velocities of center of mass
		      bodies $A, B$
		\item $\vec v_{ap1}=$ initial velocity of impact point $P$ on body $A$
		\item $\vec v_{bp1}=$ initial velocity of impact point $P$ on body $B$
		\item $\vec v_{p1}=$ initial relative velocity of impact points on body $A, B$
		\item $\vec v_{p2}=$ final relative velocity of impact points on body $A, B$
		\item $\vec n=$ normal vector
		\item $ e=$ elastic coefficient (0 = inelastic, 1 = perfectly elastic)
	\end{itemize}
\end{multicols}

Some of those variables, like the ones in the left column, are already given by
the simulation, some of them have been computed thanks to the algorithms in
section \ref{sub:collision-detection} (normal vector and position of point $P$),
and some have still to be defined mathematically, such as $\vec v_{ap1},\vec
	v_{bp1},\vec v_{p1}$ and $\vec v_{p2}$. So the velocity vectors of impact
point $P$ on both bodies, before the collision, are
\begin{equation} \label{eq:vabp1}
	\begin{split}
		\vec v_{ap1} = \vec v_{a1} + \omega_{a1} \times \vec r_{ap} \\
		\vec v_{bp1} = \vec v_{b1} + \omega_{b1} \times \vec r_{bp}
	\end{split}
\end{equation}

Similarly, the final velocities are
\begin{equation} \label{eq:vabp2}
	\begin{split}
		\vec v_{ap2} = \vec v_{a2} + \omega_{a2} \times \vec r_{ap}\\
		\vec v_{bp2} = \vec v_{b2} + \omega_{b2} \times \vec r_{bp}
	\end{split}
\end{equation}
Here we are regarding the angular velocity as a 3-dimensional vector
perpendicular to the plane, so that the cross product is calculated as

$$ \omega \times \vec r = \begin{pmatrix} 0\\0\\\omega \end{pmatrix} \times
	\begin{pmatrix} r_x\\r_y\\0 \end{pmatrix} = \begin{pmatrix} -\omega r_y \\
		\omega r_x  \\0\end{pmatrix}  $$

We these variables, we can finally define the relative velocities $\vec
	v_{p1}$ and $\vec v_{p2}$

\begin{equation}
	\label{eq:vp1}
	\begin{split}
		\vec v_{p1} = \vec v_{ap1} - \vec v_{bp2}\\
		\vec v_{p2} = \vec v_{ap2} - \vec v_{bp2}
	\end{split}
\end{equation}
If we expand by using \ref{eq:vabp1} and \ref{eq:vabp2}, we get
\begin{equation}
	\begin{split}
		\vec v_{p1} = \vec v_{a1} + \omega_{a1} \times \vec r_{ap} - \vec v_{b1}
		+ \omega_{b1} \times \vec r_{bp}\\
		\vec v_{p2} = \vec v_{a2} + \omega_{a2} \times \vec r_{ap} - \vec v_{b2} + \omega_{b2} \times \vec r_{bp}
	\end{split}
\end{equation}

The relative velocity of point $P$ along the normal vector $\vec n$ is
$$ \vec v_{p1} \cdot \vec n $$
Note that for a collision to occur this relative normal velocity must be
negative (that is, the objects must be approaching each other). Let e be the
elasticity of the collision, having a value between 0 (inelastic) and 1
(perfectly elastic). We now make an important assumption in the form of the
following relation
\begin{equation}
	\label{eq:vp2n}
	\vec v_{p2} \cdot \vec n = - e \vec v_{p1} \cdot \vec n
\end{equation}
This says that the velocity at which the objects fly apart is proportional to
the velocity with which they were coming together. The proportionality factor is
the elasticity $e$.

\paragraph{Collision Impulse} In simple terms, collision impulse refers to the
change in momentum experienced by an object during a collision. It is a measure
of the force applied to an object over a short period of time. We imagine that
during the collision there is a very large force acting for a very brief period
of time. If you integrate (sum) that force over that brief time, you get the
impulse.

In our simulation, we are assuming that no friction is happening during the
collision, so that the only force we have to consider is the one of along the
normal vector $\vec n$. The friction would create a force perpendicular to the
normal, and it would make things a little too complicated for the scope of this
project.

Since the only force we consider is the normal one, we can consider the net
impulse to be $j \vec n$, where $j$ is the impulse parameter. The body $A$ will
experience the net impulse $j \vec n$ and body $B$ will experience it's negation
$- j \vec n$ since the force that $B$ experience is equal an opposite to the one
experienced by $A$. The impulse is a change in momentum. Momentum has units of
velocity times mass, so if we divide the impulse by the mass we get the change
in velocity. We can relate initial and final velocities as

\begin{equation}
	\label{eq:va2}
	\vec v_{a2} = \vec v_{a1} +  \frac{j\vec n}{m_a}
\end{equation}
\begin{equation}
	\label{eq:vb2}
	\vec v_{b2} = \vec v_{b1} -  \frac{j\vec n}{m_b}
\end{equation}


For the final angular speed, it's change from the impulse $j\vec n$ is given by
$\vec r_{ap} \times j \vec n$. We can divide the result by the moment of
inertia, which was calculated in section \ref{sub:moment}, to get convert the
change in angular momentum into change in angular speed. Similarly as above, we
can relate the initial and final angular velocity as

\begin{equation}
	\label{eq:omega_a2}
	\omega_{a2} = \omega_{a1} + \frac{\vec r_{ap} \times j\vec n}{I_a}
\end{equation}
\begin{equation}
	\label{eq:omega_b2}
	\omega_{b2} = \omega_{b1} - \frac{\vec r_{bp} \times j\vec n}{I_b}
\end{equation}

\paragraph{Solving for the impulse parameter}
Now we have everything we need to solve for the impulse parameter $j$. The bulk
of the calculations can be found in section \ref{app:impulse_long} of the
appendix \ref{appendix:calculations}. But basically we start with equation
\ref{eq:vp2n} and we end up with
\begin{equation}
	\label{eq:j}
	j = \frac{ - (1+e) \cdot \vec v_{ap1} \cdot \vec n }{\frac{1}{m_a} + \frac{1}{m_b} +
		\frac{\left( \vec r_{ap} \times \vec n \right)^2}{I_a} + \frac{\left( \vec
			r_{bp} \times \vec n \right)^2}{I_b}}
\end{equation}