\chapter{Physically-based animation} \label{chapter:pba}
\section{Particles} \label{sec:particles}
\subsection{Basic equation}
A particle is a moving point with a mass. Generally the mass is a constant over time while its position and velocity can vary. We therefore consider the mass $ m $ (in $ kg $) as an attribute while its position $ x $ (in $ m $) and its velocity $ v $ (in $m.s^{-1}$) are its state variables. 
Particles obey Newton's law: 
\begin{equation} \label{eq:newton}
ma = \Sigma f
\end{equation} 
where $a$ is the particle's acceleration (in $m.s^{-2}$). This defines the Ordinary Differential Equation (ODE) $a(t) = f(x,v,t)$. Physically-based animation requires us to repeatedly solve it over the interval $[t,t+dt]$ and redisplay.

From equation \ref{eq:newton} we can deduce $a=1/m\;\Sigma f$ and integrate time over a time step dt, and so on. The simplest integration method is Euler's explicit scheme:
\begin{eqnarray}
v(t+dt) &=& v(t)+a(t)dt \nonumber \\
x(t+dt) &=& x(t)+v(t)dt \label{eq:expliciteuler}
\end{eqnarray}

Particles can live in any k-dimensional space $\Re^k$, where state variables and forces are k-dimensional vectors of scalar values. The dynamics equation \ref{eq:newton} applies to each scalar value.
When considering $n$ particles the equations can be conveniently written in matrix form:
\begin{eqnarray*}
\ma M \ve a &=& \ve f\\
\ve a &=& \ma M^{-1} \ve f
\end{eqnarray*}
where \ma M is a diagonal matrix of dimension $kn\times kn$ , \ve a and \ve f are vectors of dimension $kn$ gathering all the scalar components associated with each particle.

Numerical integration can generate instabilities leading to the divergence of the simulated system. To avoid this, one solution is to decrease the time step. The other solution is to use an implicit integration scheme, which takes into account the variation of the forces during the time step. The simplest one is the implicit Euler's method which applies a step based on the forces at the end of the time step instead of the beginning. Equation \ref{eq:expliciteuler} becomes:
\begin{eqnarray}
v(t+dt) &=& v(t)+a(t+dt)dt \nonumber\\
x(t+dt) &=& x(t)+v(t+dt)dt \label{eq:impliciteuler}
\end{eqnarray}
Since $a(t+dt)$ is unknown we have to solve an equation.
If we write $v(t+dt) = v(t)+\Delta v$ the matrix form equation to solve is:
\begin{equation}
\label{eq:matimplicit}
\left( \ma M + dt \ma D + dt^2 \ma K \right) \ve{\Delta v} = dt \left( \ve f(t) + \ma D \ve{v}(t) \right) 
\end{equation}
where the damping matrix $\ma D = \delta \ve f/\delta \ve v$ encodes the variation of force given a variation of velocity, and the stiffness matrix $\ma K = \delta \ve f/\delta \ve p$ encodes the variation of force given a variation of position.

\subsection{Forces}
The forces are responsible for the accelerations of the bodies. Their physical unit is the Newton ($N=kg.m.s^{-2}$).
Here we briefly review the most commonly used forces.

\subsubsection{Weight}
A uniform gravitational field $g$ (in $m.s^{-2}$) applies a force 
\begin{equation} \label{eq:gravity}
f = mg
\end{equation} 
to each particle where $m$ is the mass of the particle.

\subsubsection{Linear damping}
Damping transforms kinetic energy to heat by applying a force opposed to the velocity. It tends to slow down the objects. Linear damping is proportional to the velocity, thus 
\begin{equation}\label{eq:lineardamping}
f=-\nu v
\end{equation}
where $v$ is the velocity of the particle and $\nu$ a positive scalar (in $kg.s^{-1}$).

\subsubsection{Air damping}
Air damping is proportional to the square of the velocity of the body with respect to the air:
\begin{equation}
f = -\rho S_u C_u v^2 u 
\end{equation}
where $\rho$ is the volumic mass the air ($kg.m^{-3}$), $v$ the velocity of the object, $u$ a no-dimensional unit vector in the direction of the velocity, $S$ the area ($m^2$) of the object projected along $u$, and $C_u$ a no-dimensional coefficient associated with the shape of the object and the direction $u$.

\subsubsection{Linear springs}
A springs applies an elastic force between two points. It is modeled using its rest length $l_0$ (in $m$) and its stiffness $k$ (in $N.m^{-1}$). 
Let $i$ and $j$ be the indices of points linked by a given spring. 
The force applied by a linear spring to point $i$ is given by:
\begin{equation}
\label{eq:spring}
f_i = k( l-l_0 ) u
\end{equation}
where $l=\|x_j - x_i\|$ is the distance between the points and $u$ a unit vector pointing from point $i$ to point $j$. The force $f_j$ applied to point $j$ is the opposite: $f_j = -f_i$. Nonlinear springs can be used to model more complex behaviors.

\subsubsection{Linear damped springs}
Damping forces are commonly associated with springs in order to dissipate energy. They are opposed to the relative velocity of the points. They are typically modeled using a coefficient $\nu$ (in $kg.s^{-1}$).
The force applied by a linear damped spring to point $i$ is given by:
\begin{equation}
\label{eq:dampedspring}
f_i = \left( k( l-l_0 ) + \nu v_{ij} \right) u
\end{equation}
where $v_{ij} = (v_j-v_i).u$ is the relative velocity of the particles along direction $u$.

\subsubsection{Finite elements}
Finite elements is a powerful paradigm for modeling continuous material. At our level, we can see them as springs acting on more than two points simultaneously. For example, a tetrahedral finite element acts on the four vertices of a tetrahedron and allows a more effective control of stiffness and volume than using springs.


%===========================================================================================

\section{Solids}
A solid is a moving reference frame with a mass matrix. In two dimensions it has three degrees of freedom (DOFs), two translations and one rotation, while in three dimensions it has six DOFs, three translations and three independent rotation values. The remainder of this document focuses on three dimensions.

\subsection{Orientation}
There are different ways of modeling orientation ot one frame with respect to another in three dimensions, each of them with advantages and drawbacks:
\begin{itemize}
\item matrices directly define the axes of the solid with respect to a reference frame. They allow fast projections from one frame to another but they contain nine dependent entries and they can not be set up intuitively;
\item Euler angles are compact (three parameters) and intuitive but they have singularities and they can not be easily combined;
 \item (axis, angle) pairs are more intuitive and more compact (four parameters) than matrices but they do not allow projections and combinations;
\item quaternions are compact (four parameters) and allow easy projections and combinations, but they are not easy to set up intuitively.
\end{itemize}

\subsubsection{Orientation matrices}
Orientations matrices are $3\times 3$ matrices, each column gathering the coordinates of one axis of the rotated frame with respect to the reference frame. Each column is thus a unit vector orthogonal to the others. This creates six relations among the nine parameters, leaving three independent DOFs.

\subsubsection{Euler angles}
Euler angles model a sequence of three rotations along three pairwise-independent directions. For example, the following matrix product represents a rotation $\alpha$ along axis $x$ followed by a rotation $\beta$ along rotated axis $y$ followed by a rotation $\gamma$ along the twice rotated axis $z$.
%\begin{equation}\label{eq:angles euler}
$$
\left(\begin{array}{ccc}
1 & 0 & 0 \\
0 & \cos\alpha & -\sin\alpha \\
0 & \sin\alpha &  \cos\alpha
\end{array}\right)
\left(\begin{array}{ccc}
\sin\beta & 0 &  \cos\beta\\
0 & 1 & 0 \\
\cos\beta & 0 & -\sin\beta
\end{array}\right)
\left(\begin{array}{ccc}
\cos\gamma & -\sin\gamma & 0\\
\sin\gamma & \cos\gamma & 0\\
0 & 0 & 1\\
\end{array}\right)
$$
%\end{equation}
Alternatively, this can be seen as a rotation  $\gamma$ along axis $z$ followed by a rotation $\beta$ along the fixed axis $y$ followed by a rotation $\alpha$ along the fixed axis $x$. An example of singularity is the fact that in this system, rotation $(\pi,\pi,0)$ is equivalent with rotation $(0,0,\pi)$. Another example is the fact that rotation $(\alpha,\pi /2, \gamma)$ is equivalent with rotation $(0,\pi /2, \alpha+\gamma)$ for any $\alpha$ and $\gamma$ (this loss of one DOF is called \emph{gimbal lock}).

\subsubsection{Axis, angle}
The rotation $\theta$ along an axis defined by a unit vector $n$ has the following matrix:
$$
\rot{\theta}{u} = \mat{I} + \sin\theta\oppvec{n} + (1-\cos\theta)\oppvec{n}^2
$$
where matrix $\oppvec{n}$  is the vector product matrix operator: 
%\begin{equation}\label{eq:oppvec}
$$
\oppvec{n} = \left(\begin{array}{ccc}
0 & -n_z & n_y \\
n_z & 0 & -n_x \\
-n_y & n_x & 0
\end{array}\right)
$$
%\end{equation}
It is possible to convert a matrix back to (axis, angle) by noticing that  $\trace{\mat R}=1+2\cos\theta$ et que $\mat R -\transp{\mat R} = 2\sin\theta\oppvec{n}$

\subsubsection{Quaternions}
Quaternions are an extension of complex numbers: $ \bm q = w + x\bm i + y \bm j + z \bm k = (w,\vect v)$.\\
w is the real part, \vect v the imagianry part.
Properties of \bm i, \bm j, \bm k:
$$
\begin{array}{l}
  \bm i^2 = \bm j^2 = \bm k^2 = -1\\
  \bm{ij}=\bm k,\;\bm{ji} = -\bm k\\
  \bm{jk}=\bm i,\;\bm{kj} = -\bm i\\
  \bm{ki}=\bm j,\;\bm{ik} = -\bm j\\
\end{array}
$$
A 3d vector is a pure imaginary quaternion:
$$
\bm p = (0,x,y,z)
$$
Product of quaternions (not commutative):
$$
\bm{q_1q_2} = (w_1w_2 - \bm v_1.\bm v_2, \; w_1\bm v_2 + w_2\bm v_1 + \bm v_1 \wedge \bm v_2)
$$
Conjugate quaternion:
$$
\begin{array}{l}
  \bm{\bar q} = w - x\bm i - y \bm j - z \bm k\\
  \bm{q\bar q} = w^2 + x^2 + y^2 + z^2
\end{array}
$$
Unit quaternions used to model rotations:
$$
\bm{q\bar q} = 1
$$
Rotation $(\theta, \bm u)$: {\em ($\bm u^2=1$)}
$$
\bm q_{(\theta,u)} = ( \cos{\frac{\theta}{2}}, u_x\sin{\frac{\theta}{2}}, u_y\sin{\frac{\theta}{2}}, u_z\sin{\frac{\theta}{2}})
$$
Rotation of a vector \bm p:$\;\;\;\bm{ qp\bar{q} }$\\
Rotation matrix associated with a unit quaternion:
$$
\left( \begin{array}{ccc}
  1 - 2y^2 - 2z^2 & 2xy - 2wz & 2xz + 2wy \\
  2xy + 2wz & 1-2x^2-2z^2 & 2yz - 2wx \\
  2xz - 2wy & 2yz + 2wx & 1 - 2x^2 - 2y^2
\end{array} \right)
$$
Combination of rotations: $\rot{\alpha}{u}\rot{\beta}{v} \longrightarrow q_{(\alpha,u)}q_{(\beta,v)}$\\
Inverse rotation: $\inv{ q_{(\theta,u)}} = q_{(-\theta,u)} = q_{(\theta,-u)} = (-w,\vect v) = (w,\vect -v)$ \\
Conversion $(w, \vect v) \longrightarrow ( \theta, \vect u)$ :
\begin{eqnarray*}
\cos(\theta/2) &=& w\\
\sin(\theta/2) &=& \|\vect v\|\\
\vect u &=& \vect v/\|\vect v\|
\end{eqnarray*}

The time derivative of the unit quaternion $q$ defining the orientation of a solid with angular velocity $\omega$ (vector of $\RRR$, see section \ref{sec:omega}) is: $\dot q = \frac{1}{2}\omega q$.


\subsection{Kinematics}
\todo{choose omega or Omega. Simplify notations where possible.}
\todo{choose n or u.}
\subsubsection{Derivative in \Rep{0} of a vector fixed in \Rep{1}. Angular velocity.} \label{sec:omega}
Consider vector \fixedans{\vect u}{1}, fixed in frame \rep{1}. Frame \rep{1} rotates with respect to frame \rep{0}. Consider the projection \vecin{\fixedans{\vect u}{1}}{0} of this vector to \rep{0}, which we sometimes call \vect u for clarity, and its derivative in \rep{0} which we write $\derivedans{\fixedans{\bm u}{1}}{0}$.


Let \mat{R(dt)} be the rotation of \rep{1} between time $t$ and $t+dt$. We write:
\begin{eqnarray}
 \vect u(t+dt) &=& \mat{R(dt)} \vect u(t)\\
 \vect u(t+dt) - \vect u(t)&=& (\mat{R(dt)}-\ident{}) \vect u(t) \label{eq ri}
\end{eqnarray}
Let $\dot{\theta}$ be the angular velocity along the rotation axis, which we set to \vect z for clarity. The first-order Taylor series is: 
$$
 \mat{R(dt)}-\ident{} = 
 \left(\begin{array}{ccc}
  cos(\dot{\theta}dt)-1 & -sin(\dot{\theta}dt) & 0\\   
   sin(\dot{\theta}dt)& cos(\dot{\theta}dt)-1& 0\\ 
  0 & 0 & 0 
 \end{array}\right)
 \longrightarrow
 \left(\begin{array}{ccc}
  0 & -\dot{\theta}dt & 0\\
  \dot{\theta}dt & 0 & 0\\
  0 & 0 & 0
 \end{array}\right) 
 =
 \dot{\theta}dt \oppvec{z}
$$
which can be easily extended to any rotation axis. Let $\vecrot{1}{0}=\dot{\theta}\vect n$. Dividing expression \ref{eq ri} by $dt$ and decreasing $dt$ to $0$ gives $\dot{\mat R} = \oppvec{ \vecrot{1}{0} }$.

We can write the time derivative in \rep{0}: $\derivedans{\fixedans{\bm u}{1}}{0}  = \vecrot{1}{0} \wedge  \fixedans{\vect u}{1}$, or more simply:

\begin{equation}\label{vrot}
\begin{array}{rcl}
 \dot{\vect u} &=& \dot{\mat R} \vect u \\
               &=& \oppvec{ \vecrot{1}{0} } \vect u \\
               &=& \vecrot{1}{0} \wedge \vect u
\end{array}
\end{equation}
 

\subsection{Velocity in \Rep{0} of a point fixed in \Rep{1}. Velocity field.}
We consider the velocity $\vfdans{A}{1}{0}$ in \rep{0} of a point $A$ fixed in \rep{1} while \rep{1} moves with respect to \rep{0}. Let $O_0$ be the origin of \rep{0} and $O_1$ the origin of \rep{1}. The following relation holds:
$$ \vfdans{A}{1}{0} = \vfdans{O_1}{1}{0} + \vecrot{1}{0} \wedge \vecf{O_1A} \label{eq vit} \label{eq vit solide}
$$

\subsection{Acceleration in \Rep{0} of  point fixed in \Rep{1}. Acceleration field. }
By deriving equation \ref{eq vit}, and based on the fact that $\vecf{O_1A}$ is fixed in \rep{1}, we get the acceleration of A, fixed in \rep{1}, with respect to \rep{0}:
\begin{equation}\label{eq acc}
 \afdans{A}{1}{0} = \afdans{O_1}{1}{0} + \accrot{1}{0}\wedge \vecf{O_1A} + \vecrot{1}{0} \wedge \left( \vecrot{1}{0} \wedge \vecf{O_1A} \right)
\end{equation}

\subsection{Derivative in \rep{0} of a vector defined in \rep{1}}
Let $(\vect e_1, \vect e_e, \vect e_3)$ be a base of \rep{1}. We have:
\begin{eqnarray*}
 \vecin{u}{1} &=& \sum_i x_i \vect e_i\\
 \dot{\vect u} &=& \sum_i \dot x_i \vect e_i + \sum_i x_i \dot{\vect e}_i
\end{eqnarray*}
and thus:
\begin{equation}\label{eq vec mob}
 \derivedans{u}{0} = \derivedans{u}{1} + \vecrot{1}{0} \wedge \vect u
\end{equation}


\subsection{Velocity in \rep{0} of a point moving in \rep{1}.}
Let \vmdans{A}{1} be the velocity of point $A$ with respect to \rep{1}. We have:
\begin{equation}\label{eq vit mob}
\vmdans{A}{0} = \vmdans{A}{1} + \vfdans{O_1}{1}{0} + \vecrot{1}{0} \wedge \vecf{O_1A}
\end{equation}
Note that $O_1$ being the origin of frame \rep{1}, we have $\vmdans{O_1}{0} = \vfdans{O_1}{1}{0}$.

\subsection{Acceleration in \rep{0} of a point moving in \rep{1}. Coriolis acceleration.}
EBy deriving equation \ref{eq vit mob} nous we get:
$$
 \amdans{A}{0} = \underbrace{\amdans{A}{1} + \vecrot{1}{0}\wedge \vmdans{A}{1}}_{\overset{\circ}{\vmdans{A}{1}}} + \amdans{O_1}{0} + \underbrace{\accrot{1}{0}\wedge \vect{O_1}{A} + \vecrot{1}{0}\wedge \vmdans{A}{1} + \vecrot{1}{0} \wedge (\vecrot{1}{0}\wedge \vecf{O_1A})}_{\overset{\circ}{\vecrot{1}{0} \wedge \vecf{O_1A}}}
$$
or:
\begin{equation}\label{eq acc mob}
 \amdans{A}{0} = \amdans{A}{1} +  \amdans{O_1}{0} + \vecrot{1}{0} \wedge (\vecrot{1}{0}\wedge \vecf{O_1A}) + 2\vecrot{1}{0}\wedge \vmdans{A}{1}
\end{equation}
with:
\begin{itemize}
\item $\amdans{A}{1} = \sum_i \ddot x_i \vect e_i$ relative acceleration
\item $\amdans{O_1}{0}$ frame acceleration (?)
\item $\vecrot{1}{0} \wedge (\vecrot{1}{0}\wedge \vecf{O_1A})$ centripetal acceleration
\item $2\vecrot{1}{0}\wedge \vmdans{A}{1}$ Coriolis acceleration
\end{itemize}


\subsection{Dynamics}
Solids accelerate linearly due to forces, and accelerate angularly due to torques. A given torque has the same value everywhere in the solid. However, a given force generates different torques at different points. The torque $\tau$ applied at point $c$ generated by a force $f$ applied at point $b$ is: $\tau=cb\wedge f$.

The acceleration $\ddot c$ of the mass center of a solid and its angular acceleration $\dot \omega$ with respect to the world are given by the relations:

    \begin{eqnarray*}
        m \bf{ \ddot c } = \sum\bf{f_{ext}} \\%\label{eq:PFD1}\\
        \bf{ I_M \dot \omega } + \bm{ \omega \times I_M \omega} = \sum\bf{\tau_{ext}} %\label{eq:PFD2}
    \end{eqnarray*}
where $m$ is the mass of the solid, $\sum {\bf f_{ext}}$ is the sum of the forces applied to the solid, ${\bf I_M}$ is the inertia matrix et $\sum {\bf \tau_{ext}}$ the sum of the torques applied to the solid and expressed at its mass center. The inertia matrix is given by:
$${\bf I_M} = \int_x \int_y \int_z \rho (x,y,z)\begin{bmatrix}  y^2+z^2 & -xy & -xz \\ -xy & x^2+z^2 & -yz \\ -xz & -yz & x^2+y^2 \end{bmatrix} dx dy dz$$
where $\rho(x,y,z)$ is the volumic mass of the naterial (in $kg.m^{-3}$).


