openrocket/core/doc/techdoc/chapter-flight-simulation.tex

\chapter{Flight simulation}
\label{chap-simulation}

In this chapter the actual flight simulation is analyzed.  First in
Section~\ref{sec-atmospheric-properties} methods for simulating
atmospheric conditions and wind are presented.  Then in
Section~\ref{sec-flight-modeling} the actual simulation procedure is
developed.



\section{Atmospheric properties}
\label{sec-atmospheric-properties}

In order to calculate the aerodynamic forces acting on the rocket it
is necessary to know the prevailing atmospheric conditions.  Since the
atmosphere is not constant with altitude, a model must be developed to
account for the changes.  Wind also plays an important role in the
flight of a rocket, and therefore it is important to have a realistic
wind model in use during the simulation.


\subsection{Atmospheric model}

The atmospheric model is responsible to estimating the atmospheric
conditions at varying altitudes.  The properties that are of most
interest are the density of air $\rho$ (which is a scaling parameter
to the aerodynamic coefficients via the dynamic pressure
$\frac{1}{2}\rho v^2$) and the speed of sound $c$ (which affects the
Mach number of the rocket, which in turn affects its aerodynamic
properties).  These may in turn be calculated from the air pressure
$p$ and temperature $T$.

Several models exist that define standard atmospheric conditions as a
function of altitude, including the Internaltional Standard
Atmosphere (ISA)~\cite{international-standard-atmosphere} and the
U.S. Standard Atmosphere~\cite{US-standard-atmosphere}.  These two
models yield identical temperature and pressure profiles for altitudes
up to 32~km.

The models are based on the assumption that air follows the ideal gas
law
%
\begin{equation}
\rho = \frac{Mp}{RT}
\end{equation}
%
where $M$ is the molecular mass of air and $R$ is the ideal gas
constant.  From the equilibrium of hydrostatic forces the differential
equation for pressure as a function of altitude $z$ can be found as
%
\begin{equation}
\dif p = -g_0 \rho \dif z = -g_0 \frac{Mp}{RT} \dif z
\label{eq-pressure-altitude}
\end{equation}
%
where $g_0$ is the gravitational acceleration.  If the temperature of
air were to be assumed to be constant, this would yield an exponential
diminishing of air pressure.

The ISA and U.S. Standard Atmospheres further specity a standard
temperature and pressure at sea level and a temperature profile for
the atmosphere.  The temperature profile is given as eight
temperatures for different altitudes, which are then linearly
interpolated.  The temperature profile and base pressures for the ISA
model are presented in Table~\ref{table-ISA-model}.  These values
along with equation~(\ref{eq-pressure-altitude}) define the
temperature/pressure profile as a function of altitude.

\begin{table}
\caption{Layers defined in the International Standard
  Atmosphere~\cite{wiki-ISA-layers}}
\label{table-ISA-model}
\begin{center}
\begin{tabular}{ccccl}
\hline
Layer & Altitude$^\dagger$ & Temperature & Lapse rate &
                                                \multicolumn{1}{c}{Pressure} \\
      &  m       &  $^\circ$C  & $^\circ$C/km & \multicolumn{1}{c}{Pa} \\
\hline
0     & 0        & $+15.0$     & $-6.5$       & 101\s325 \\
1     & 11\s000  & $-56.5$     & $+0.0$       & \num22\s632 \\
2     & 20\s000  & $-56.5$     & $+1.0$       & \num\num5\s474.9 \\
3     & 32\s000  & $-44.5$     & $+2.8$       & \num\num\num\s868.02 \\
4     & 47\s000  & \num$-2.5$  & $+0.0$       & \num\num\num\s110.91 \\
5     & 51\s000  & \num$-2.5$  & $-2.8$       & \num\num\num\s\num66.939 \\
6     & 71\s000  & $-58.5$     & $-2.0$       & \num\num\num\s\num\num3.9564 \\
7     & 84\s852  & $-86.2$     &              & \num\num\num\s\num\num0.3734 \\
\hline
\end{tabular}
\end{center}
\vspace{-3mm}
{\footnotesize $^\dagger$ Altitude is the geopotential height which
  does not account for the diminution of gravity at high altitudes.}
\vspace{3mm}
\end{table}

These models are totally static and do not take into account any local
flight conditions.  Many rocketeers may be interested in flight
differences during summer and winter and what kind of effect air pressure
has on the flight.  These are also parameters that can easily be
measured on site when launching rockets.  On the other hand, it is
generally hard to know a specific temperature profile for a specific
day.  Therefore the atmospheric model was extended to allow the user
to specify the base conditions either at mean sea level or at the
altitude of the launch site.  These values are simply assigned to the
first layer of the atmospheric model.  Most model rockets do not
exceed altitudes of a few kilometers, and therefore the flight
conditions at the launch site will dominate the flight.

One parameter that also has an effect on air density and the speed of
sound is humidity.  The standard models do not include any definition
of humidity as a function of altitude.  Furthermore, the effect of
humidity on air density and the speed of sound is marginal.  The
difference in air density and the speed of sound between completely dry
air and saturated air at standard conditions are both less than 1\%.
Therefore the effect of humidity has been ignored.




\subsection{Wind modeling}

Wind plays a critical role in the flight of model rockets.  As has
been seen, large angles of attack may cause rockets to lose a
significant amount of stability and even go unstable.  Over-stable
rockets may weathercock and turn into the wind.  In a perfectly static
atmosphere a rocket would, in principle, fly its entire flight
directly upwards at zero angle of attack.  Therefore, the effect of
wind must be taken into account in a full rocket simulation.

Most model rocketeers, however, do not have access to a full wind
profile of the area they are launching in.  Different layers of air
may have different wind velocities and directions.  Modeling such
complex patterns is beyond the scope of this project. Therefore, the
goal is to produce a realistic wind model that can be specified with
only a few parameters understandable to the user and that covers
altitudes of most rocket flights.  Extensions to allow for multiple
air layers may be added in the future.

In addition to a constant average velocity, wind always has some
degree of turbulence in it.  The effect of turbulence can be modeled
by summing the steady flow of air and a random, zero-mean turbulence
velocity.  Two central aspects of the turbulence velocity are the
amplitude of the variation and the frequencies at which they occur.
Therefore a reasonable turbulence model is achieved by a random
process that produces a sequence with a similar distribution and
frequency spectrum as that of real wind.

Several models of the spectrum of wind turbulence at specific
altitudes exist.  Two commonly used such spectra are the {\it Kaimal}
and {\it von K<>rm<72>n} wind turbulence
spectra~\cite[p.~23]{wind-energy-handbook}:
%
\begin{eqnarray}
\mbox{Kaimal:} & & \frac{S_u(f)}{\sigma_u^2} =
    \frac{4 L_{1u} / U}{(1 + 6fL_{1u}/U)^{5/3}} \label{eq-kaimal-wind} \\
%
\mbox{von K<>rm<72>n:} & & \frac{S_u(f)}{\sigma_u^2} =
    \frac{4 L_{2u} / U}{(1 + 70.8(fL_{2u}/U)^2)^{5/6}} \label{eq-karman-wind}
\end{eqnarray}

Here $S_u(f)$ is the spectral density function of the turbulence
velocity and $f$ the turbulence frequency, $\sigma_u$ the standard
deviation of the turbulence velocity, $L_{1u}$ and $L_{2u}$ length
parameters and $U$ the average wind speed.

Both models approach the asymptotic limit 
$S_u(f)/\sigma_u^2 \sim f^{-5/3}$ quite fast.  Above frequencies of
0.5~Hz the difference between equation~(\ref{eq-kaimal-wind}) and the
same equation without the term 1 in the denominator is less than 4\%.
Since the time scale of a model rocket's flight is quite short, the
effect of extremely low frequencies can be ignored.  Therefore
turbulence may reasonably well be modelled by utilizing 
{\it pink noise} that has a spectrum of $1/f^\alpha$ with $\alpha=5/3$.
True pink noise has the additional useful property of being
scale-invariant.  This means that a stream of pink noise samples may
be generated and assumed to be at any sampling rate while maintaining
their spectral properties.

Discerete samples of pink noise with spectrum $1/f^\alpha$ can be
generated by applying a suitable digital filter to {\it white noise},
which is simply uncorrelated pseudorandom numbers.  One such filter is
the infinite impulse response (IIR) filter presented by
Kasdin~\cite{pink-filter}:
%
\begin{equation}
x_n = w_n - a_1 x_{n-1} - a_2 x_{n-2} - a_3 x_{n-3} - \ldots
\label{eq-pink-generator}
\end{equation}
%
where $x_i$ are the generated samples, $w_n$ is a generated white
random number and the coefficients are computed using
%
\begin{equation}
\begin{array}{rl}
a_0 & = 1 \\
a_k & = \del{k-1-\frac{\alpha}{2}} \frac{a_{k-1}}{k}.
\end{array}
\label{eq-pink-coefficients}
\end{equation}
%
The infinite sum may be truncated with a suitable number of terms.
In the context of IIR filters these terms are calles {\it poles}.
Experimentation showed that already 1--3 poles provides a reasonably
accurate frequency spectrum in the high frequency range.

One problem in using pink noise as a turbulence velocity
model is that the power spectrum of pure pink noise goes to
infinity at very low frequencies.  This means that a long sequence
of random values may deviate significantly from zero.  However, when
using the truncated IIR filter of equation~(\ref{eq-pink-generator}),
the spectrum density becomes constant below a certain limiting
frequency, dependent on the number of poles used.  By adjusting the
number of poles used, the limiting frequency can be adjusted to a value
suitable for model rocket flight.  Specifically, the number of poles
must be selected such that the limiting frequency is suitable at the
chosen sampling rate.  

It is also desirable that the simulation resolution does not affect
the wind conditions.  For example, a simulation with a time step of
10~ms should experience the same wind conditions as a simulation with
a time step of 5~ms.  This is achieved by selecting a constant
turbulence generation frequency and interpolating between the
generated points when necessary.  The fixed frequency was chosen at
20~Hz, which can still simulate fluctuations at a time scale of 0.1
seconds.

\begin{figure}[p]
\centering
\epsfig{file=figures/wind/pinktime, width=105mm}
\caption{The effect of the number of IIR filter poles on two 20 second
  samples of generated turbulence, normalized so that the two-pole
  sequence has standard deviation one.}
\label{fig-pink-poles}
\end{figure}

\begin{figure}[p]
\centering
\epsfig{file=figures/wind/pinkfreq, width=95mm}
\caption{The average power spectrum of 100 turbulence
  simulations using a two-pole IIR filter (solid) and the Kaimal
  turbulence spectrum (dashed); vertical axis arbitrary.}
\label{fig-pink-spectrum}
\end{figure}



The effect of the number of poles is depicted in
Figure~\ref{fig-pink-poles}, where two pink noise sequences were
generated from the same random number source with two-pole and
ten-pole IIR filters.  A small number of poles generates values strongly
centered on zero, while a larger number of poles introduces more low
frequency variability.  Since the free-flight time of a typical model
rocket is of the order of 5--30 seconds, it is desireable that the
maximum gust length during the flight is substantially shorter than
this.  Therefore the pink noise generator used by the wind model was
chosen to contain only two poles, which has a limiting frequency of
approximately 0.3~Hz when sampled at 20~Hz.  This means that gusts of
wind longer than 3--5 seconds will be rare in the simulted turbulence,
which is a suitable gust length for modeling typical model rocket
flight. Figure~\ref{fig-pink-spectrum} depicts the resulting pink
noise spectrum of the two-pole IIR filter and the Kaimal spectrum of
equation~(\ref{eq-kaimal-wind}) scaled to match each other.


%, which causes frequency
%components below approximately 0.3~Hz to be subdued.  Therefore, gusts
%of wind longer than 3--5 seconds will be rare in the simulated wind, a
%suitable time scale for the flight of a model rocket.
%Figure~\ref{fig-turbulence}(a) shows a 20 second sample of the
%generated turbulence, normalized to have a standard deviation of one.
%Figure~\ref{fig-turbulence}(b) depicts the actual frequency spectrum
%of the generated turbulence and the Kaimal spectrum of
%equation~(\ref{eq-kaimal-wind}) scaled to match each other.



To simplify the model, the average wind speed is assumed to be
constant with altitude and in a constant direction.  This allows
specifying the model parameters using just the average wind speed and
its standard deviation.  An alternative parameter for specifying the
turbulence amplitude is the {\it turbulence intensity}, which is the
percentage that the standard deviation is of the average wind
velocity,
%
\begin{equation}
I_u = \frac{\sigma_u}{U}.
\end{equation}
%
Wind farm load design standards typically specify turbulence
intensities around 10\ldots20\%~\cite[p.~22]{wind-energy-handbook}.
It is assumed that these intensities are at the top of the range of
conditions in which model rockets are typically flown.

Overall, the process to generate the wind velocity as a function of
time from the average wind velocity $U$ and standard deviation
$\sigma_u$ can be summarized in the following steps:
%
\begin{enumerate}
%\item[Input:]  Average wind velocity $U$ and standard deviation
%  $\sigma_u$.
%
\item Generate a pink noise sample $x_n$ from a Gaussian white noise
  sample $w_n$ using equations~(\ref{eq-pink-generator}) and
  (\ref{eq-pink-coefficients}) with two memory terms included.

\item Scale the sample to a standard deviation one.  This is performed
  by dividing the value by a previously calculated standard deviation
  of a long, unscaled pink noise sequence (2.252 for the two-pole IIR
  filter).

\item The wind velocity at time $n\cdot\Delta t$ ($\Delta t = 0.05\rm~s$)
  is $U_n = U + \sigma_u x_n$.  Velocities in between are interpolated.
\end{enumerate}




\section{Modeling rocket flight}
\label{sec-flight-modeling}


Modeling of rocket flight is based on Newton's laws.  The basic forces
acting upon a rocket are gravity, thrust from the motors and
aerodynamic forces and moments.  These forces and moments are
calculated and integrated numerically to yield a simulation over a
full flight.

Since most model rockets fly at a maximum a few kilometers high, the
curvature of the Earth is not taken into account.  Assuming a flat
Earth allows us to use simple Cartesian coordinates to represent the
position and altitude of the rocket.  As a consequence, the coriolis
effect when flying long distances north or south is not simulated
either.



\subsection{Coordinates and orientation}

During a rocket's flight many quantities, such as the aerodynamical
forces and thrust from the motors, are relative to the rocket itself,
while others, such as the position and gravitational force, are more
naturally described relative to the launch site.  Therefore two sets
of coordinates are defined, the {\it rocket coordinates}, which are
the same as used in Chapter~\ref{chap-aerodynamics}, and 
{\it world coordinates}, which is a fixed coordinate system with the
origin at the position of launch.

The position and velocity of a rocket are most naturally maintained as
Cartesian world coordinates.  Following normal convensions, the
$xy$-plane is selected to be parallel to the ground and the $z$-axis
is chosen to point upwards.  In flight dynamics of aircraft the
$z$-axis often points towards the earth, but in the case of rockets it
is natural to have the rocket's altitude as the $z$-coordinate.

Since the wind is assumed to be unidirectional and the Coriolis effect
is ignored, it may be assumed that the wind is directed along the
$x$-axis.  The angle of the launch rod may then be positioned relative
to the direction of the wind without any loss of generality.

Determining the orientation of a rocket is more complicated.  A
natural choise for defining the orientation would be to use the
spherical coordinate zenith and azimuth angles $(\theta, \phi)$ and an
additional roll angle parameter.  Another choise common in aviation is
to use {\it Euler angles}~\cite{wiki-euler-angles}.  However, both of
these systems have notable shortcomings.  Both systems have
singularity points, in which the value of some parameter is
ambiguous.  With spherical coordinates, this is the direction of the
$z$-axis, in which case the azimuth angle $\phi$ has no effect on the
position.  Rotations that occur near these points must often be
handled as special cases.  Furthermore, rotations in spherical
coordinate systems contain complex trigonometric formulae which are
prone to programming errors.

The solution to the singularity problem is to introduce an extra
parameter and an additional constraint to the system.  For example,
the direction of a rocket could be defined by a three-dimensional unit
vector $(x,y,z)$ instead of just the zenith and azimuth angles.  The
additional constraint is that the vector must be of unit length.  This
kind of representation has no singularity points which would require
special consideration.

Furthermore, Euler's rotation theorem states that a rigid body can be
rotated from any orientation to any other orientation by a single
rotation around a specific axis~\cite{wiki-euler-rotation-theorem}.
Therefore instead of defining quantities that define the orientation
of the rocket we can define a three-dimensional rotation that rotates
the rocket from a known reference orientation to the current
orientation. This has the additional advantage that the same rotation
and its inverse can be used to transform any vector between world
coordinates and rocket coordinates.

A simple and efficient way of descibing the 3D rotation is by using
{\it unit quaternions}.  Each unit quaternion corresponds to a unique
3D rotation, and they are remarkably simple to combine and use.  The
following section will present a brief overview of the properties of
quaternions.

The fixed reference orientation of the rocket defines the rocket
pointing towards the positive $z$-axis in world coordinates and an
arbitrary but fixed roll angle.  The orientation of the rocket is then
stored as a unit quaternion that rotates the rocket from this
reference orientation to its current orientation.  
This rotation can also be used to transform vectors from world
coordinates to rocket coordinates and its inverse from rocket
coordinates to world coordinates.  (Note that the rocket's initial
orientation on the launch pad may already be different than its
reference orientation if the launch rod is not completely vertical.)




\subsection{Quaternions}

{\it Quaternions} are an extension of complex numbers into four
dimensions.  The usefulness of quaternions arises from their use in
spatial rotations.  Similar to the way multiplication with a complex
number of unit length $e^{i\phi}$ corresponds to a rotation of angle
$\phi$ around the origin on the complex plane, multiplication with
unit quaternions correspond to specific 3D rotations around an axis.
A more thorough review of quaternions and their use in spatial
rotations is available in Wikipedia~\cite{wiki-quaternion-rotations}.

The typical notation of quaternions resembles the addition of a scalar
and a vector:
%
\begin{equation}
q = w + x\vi + y\vj + z\vk = w + \vect v
\end{equation}
%
Addition of quaternions and multiplication with a scalar operate as
expected.  However, the multiplication of two quaternions is
non-commutative (in general $ab \neq ba$) and follows the rules
%
\begin{equation}
\vi^2 = \vj^2 = \vk^2 = \vi\vj\vk = -1.
\end{equation}
%
As a corollary, the following equations hold:
%
\begin{equation}
\begin{array}{rl}
\vi\vj = \vk  \hspace{15mm}& \vj\vi = -\vk \\
\vj\vk = \vi  \hspace{15mm}& \vk\vj = -\vi \\
\vk\vi = \vj  \hspace{15mm}& \vi\vk = -\vj 
\end{array}
\end{equation}
%
The general multiplication of two quaternions becomes
%
\begin{equation}
\begin{array}{rl}
(a + b\vi + c\vj + d\vk)(w + x\vi + y\vj + z\vk)\;\; =
 &   (aw-bx-cy-dz) \\
 & + (ax+bw+cz-dy)\;\vi \\
 & + (ay-bz+cw+dx)\;\vj \\
 & + (az+by-cx+dw)\;\vk
\end{array}
\end{equation}
%
while the norm of a quaternion is defined in the normal manner
%
\begin{equation}
|q| = \sqrt{w^2+x^2+y^2+z^2}.
\end{equation}

The usefulness of quaternions becomes evident when we consider a
rotation around a vector $\vect u$, $|\vect u|=1$ by an angle $\phi$.
Let
%
\begin{equation}
q = \cos\frac{\phi}{2} + \vect u \sin\frac{\phi}{2}.
\label{eq-rotation-quaternion}
\end{equation}
%
Now the previously mentioned rotation of a three-dimensional vector
$\vect v$ defined by $\vi$, $\vj$ and $\vk$ is equivalent to the
quaternion product
%
\begin{equation}
\vect v \mapsto q\vect v q^{-1}.
\end{equation}
%
Similarly, the inverse rotation is equivalent to the transformation
%
\begin{equation}
\vect v \mapsto q^{-1} \vect v q.
\end{equation}
%
The problem simplifies even further, since for unit quaternions
%
\begin{equation}
q^{-1} = (w + x\vi + y\vj + z\vk)^{-1} = w - x\vi - y\vj - z\vk.
\end{equation}
%
Vectors can therefore be considered quaternions with no scalar
component and their rotation is equivalent to the left- and right-sided
multiplication with unit quaternions, requiring a total of 24
floating-point multiplications.  Even if this does not make the
rotations more efficient, it simplifies the trigonometry considerably
and therefore helps reduce programming errors.


\subsection{Mass and moment of inertia calculations}
\label{sec-mass-inertia}

Converting the forces and moments into linear and angluar acceleration
requires knowledge of the rocket's mass and moments of inertia.  The
mass of a component can be easily calculated from its volume and
density.  Due to the highly symmetrical nature of rockets, the rocket
centerline is commonly a principal axis for the moments of inertia.
Furthermore, the moments of inertia around the in the $y$- and
$z$-axes are very close to one another.  Therefore as a simplification
only two moments of inertia are calculated, the longitudal and
rotational moment of inertia.  These can be easily calculated for each
component using standard formulae~\cite{wiki-moments-of-inertia} and
combined to yield the moments of the entire rocket.

This is a good way of calculating the mass, CG and inertia of a rocket
during the design phase.  However, actual rocket components often have
a slightly different density or additional sources of mass such as
glue attached to them.  These cannot be effectively modeled by the
simulator, since it would be extremely tedious to define all these
properties.  Instead, some properties of the components can be
overridden to utilize measured values.

Two properties that can very easily be measured are the mass and
CG position of a component.  Measuring the moments of inertia is a
much harder task.  Therefore the moments of inertia are still computed
automatically, but are scaled by the overridden measurement values.

If the mass of a component is overridden by a measured value, the
moments of inertia are scaled linearly according to the mass.  This
assumes that the extra weight is distributed evenly along the
component.  If the CG position is overridden, there is no knowledge
where the extra weight is at.  Therefore as a best guess the moments
of inertia are updated by shifting the moment axis according to the
parallel axis theorem.

As the components are computed individually and then combined, the
overriding can take place either for individual components or larger
combinations.  It is especially useful is to override the mass and/or CG
position of the entire rocket.  This allows constructing a rocket from
components whose masses are not precisely known and afterwards scaling
the moments of inertia to closely match true values.



\subsection{Flight simulation}

The process of simulating rocket flight can be broken down into the
following steps:

\begin{enumerate}
\setcounter{enumi}{-1}
\item Initialize the rocket in a known position and orientation at
  time $t=0$.
\item Compute the local wind velocity and other atmospheric conditions.
\item Compute the current airspeed, angle of attack, lateral wind
  direction and other flight parameters.
\item Compute the aerodynamic forces and moments affecting the rocket.
\item Compute the effect of motor thrust and gravity.
\item Compute the mass and moments of inertia of the rocket and from
  these the linear and rotational acceleration of the rocket.
\item Numerically integrate the acceleration to the rocket's position
  and orientation during a time step $\Delta t$ and update the current
  time $t \mapsto t+\Delta t$.
\end{enumerate}

Steps 1--6 are repeated until an end criteria is met, typically until
the rocket has landed.

The computation of the atmospheric properties and instantaneous wind
velocity were discussed in Section~\ref{sec-atmospheric-properties}.
The local wind velocity is added to the rocket velocity to get the
airspeed velocity of the rocket.  By inverse rotation this quantity is
obtained in rocket coordinates, from which the angle of attack and
other flight parameters can be computed.

After the instantaneous flight parameters are known, the aerodynamic
forces can be computed as discussed in
Chapter~\ref{chap-aerodynamics}.  The computed forces are in the
rocket coordinates, and can be converted to world coordinates by
applying the orientation rotation.  The thrust from the motors is
similarly calculated from the thrust curves and converted to world
coordinates, while the direction of gravity is already in world
coordinates.  When all of the the forces and moments acting upon the
rocket are known, the linear and rotational accelerations can be
calculated using the mass and moments of inertia discussed in
Section~\ref{sec-mass-inertia}.

The numerical integration is performed using the Runge-Kutta~4 (RK4)
integration method.  In order to simulate the differential equations
%
\begin{equation}
\begin{split}
x''(t) &= a(t) \\
\phi''(t) &= \alpha(t)
\end{split}
\end{equation}
%
the equation is first divided into first-order equations using the
substitutions $v(t)=x'(t)$ and $\omega(t)=\phi'(t)$:
%
\begin{equation}
\begin{split}
v'(t) &= a(t) \\
x'(t) &= v(t) \\
\omega'(t) &= \alpha(t) \\
\phi'(t)   &= \omega(t)
\end{split}
\end{equation}
%
For brevity, this is presented in the first order representation
%
\begin{equation}
y' = f(y,\; t)
\end{equation}
%
where $y$ is a vector function containing the position and orientation
of the rocket.

Next the right-hand side is evaluated at four positions, dependent on
the previous evaluations:
%
\begin{equation}
\begin{split}
k_1 &= f(y_0,\; t_0) \\
k_2 &= f(y_0 + k_1\:\mbox{$\frac{\Delta t}{2}$},\; 
       t_0 + \mbox{$\frac{\Delta t}{2}$}) \\
k_3 &= f(y_0 + k_2\:\mbox{$\frac{\Delta t}{2}$},\; 
       t_0 + \mbox{$\frac{\Delta t}{2}$}) \\
k_4 &= f(y_0 + k_3\:\Delta t,\; t_0 + \Delta t)
\end{split}
\end{equation}
%
Finally, the result is a weighted sum of these values:
%
\begin{align}
y_1 &= y_0 + \frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\,\Delta t \\
t_1 &= t_0 + \Delta t
\end{align}

Computing the values $k_1\ldots k_4$ involves performing steps~1--5
four times per simulation iteration, but results in significantly
better simulation precision.  The method is a fourth-order integration
method, meaning that the error incurred during one simulation step is
of the order  $O(\Delta t^5)$ and of the total simulation 
$O(\Delta t^4)$.  This is a considerable improvement
over, for example, simple Euler integration, which has a total error
of the order $O(\Delta t)$.  Halving the time step in an Euler
integration only halves the total error, but reduces the error of a
RK4 simulation 16-fold.

The example above used a total rotation vector $\phi$ to contain the
orientation of the rocket.  Instead, this is replaced by the rotation
quaternion, which can be utilized directly as a transformation between
world and rocket coordinates.  Instead of updating the total rotation
vector,
%
\begin{equation}
\phi_1 = \phi_0 + \omega\,\Delta t,
\end{equation}
%
the orientation quaternion $o$ is updated by the same amount by
%
\begin{equation}
o_1 = \del{\cos\del{|\omega|\,\Delta t} +
\hat\omega\sin\del{|\omega|\,\Delta t}} \cdot o_0.
\end{equation}
%
The first term is simply the unit quaternion corresponding to the
3D rotation $\omega\,\Delta t$ as in
equation~(\ref{eq-rotation-quaternion}).  It is applied to the
previous value $o_0$ by multiplying the quaternion from the left.
This update is performed both during the calculation of 
$k_2\ldots k_4$ and when computing the final step result.  Finally, in
order to improve numerical stability, the quaternion is normalized to
unit length.

Since most of a rocket's flight occurs in a straight line, rather
large time steps can be utilized.  However, the rocket may encounter
occasional oscillation, which may affect its flight notably.
Therefore the time step utilized is dynamically reduced in cases where
the angular velocity or angular acceleration exceeds a predefined
limit.  This allows utilizing reasonably large time steps for most of
the flight, while maintaining the accuracy during oscillation.


\subsection{Recovery simulation}

All model rockets must have some recovery system for safe landing.
This is typically done either using a parachute or a streamer.  When a
parachute is deployed the rocket typically splits in half, and it is
no longer practical to compute the orientation of the rocket.
Therefore at this point the simulation changes to a simpler, three
degree of freedom simulation, where only the position of the rocket is
computed.

The entire drag coefficient of the rocket is assumed to come from the
deployed recovery devices.  For parachutes the drag coefficient is
by default 0.8~\cite[p.~13-23]{hoerner} with the reference area being the
area of the parachute.  The user can also define their own drag
coefficient.

The drag coefficient of streamers depend on the material, width and
length of the streamer.  The drag coefficient and optimization of
streamers has been an item of much intrest within the rocketry
community, with competitions being held on streamer descent time
durations~\cite{streamer-optimization}.  In order to estimate the drag
coefficient of streamers, a series of experiments were perfomed using
the $40\times40\times120$~cm wind tunnel of
Pollux~\cite{pollux-wind-tunnel}.  The experiments were performed
using various materials, widths and lengths of streamers and at
different wind speeds.  From these results an empirical formula was
devised that estimates the drag coefficient of streamers.  The
experimental results and the derivation of the empirical formula are
presented in Appendix~\ref{app-streamers}.  Validation performed with
an independent set of measurements indicates that the drag coefficient
is estimated with an accuracy of about 20\%, which translates to a
descent velocity accuracy within 15\% of the true value.




\subsection{Simulation events}

Numerous different events may cause actions to be taken during a
rocket's flight.  For example in high-power rockets the burnout or
ignition charge of the first stage's motor may trigger the ignition of
a second stage motor.  Similarly a flight computer may deploy a small
drogue parachute when apogee is detected and the main parachute is
deployed later at a predefined lower altitude.  To accomodate
different configurations a simulation event system is used, where
events may cause other events to be triggered.

Table~\ref{tab-simulation-events} lists the available simulation
events and which of them can be used to trigger motor ignition or recovery
device deployment.  Each trigger event may additionally include a
delay time.  For example, one motor may be configured to ignite at
launch and a second motor to ignite using a timer at 5 seconds after
launch.  Alternatively, a short delay of 0.5--1 seconds may be used to
simulate the delay of an ejection charge igniting the upper stage
motors.

The flight events are also stored along with the simulated flight data
for later analysis.  They are also available to the simulation
listeners, described in Section~\ref{sec-listeners}, to act upon
specific conditions.

\begin{table}
\caption{Simulation events and the actions they may trigger (motor
  ignition or recovery device deployment).}
\label{tab-simulation-events}
%
\begin{center}
\begin{tabular}{ll}
Event description & Triggers \\
\hline
Rocket launch at $t=0$          & Ignition, recovery \\
Motor ignition                  & None \\
Motor burnout                   & Ignition \\
Motor ejection charge           & Ignition, recovery \\
Launch rod cleared              & None \\
Apogee detected                 & Recovery \\
Change in altitude              & Recovery \\
Touchdown after flight          & None \\
Deployment of a recovery device & None \\
End of simulation               & None \\
\hline
\end{tabular}
\end{center}
\end{table}