\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á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á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}