Euler Angles

Preface.

One of the central problem in a flight simulator is the correct determination in every instant of all the cinematic parameters of the Aircraft: this includes position (according to a defined reference system), velocity and acceleration both linear and angular.

This implies the resolution of the differential equations of the motion of the aircraft, once some key parameters are known, like all data related to mass, aerodynamics and engine(s), parameters available only with the collaboration of the aircraft manufacturer.

While the resolution of the problem related to linear position, velocities and accelerations is trivial, the same problem related to the angles (aircraft attitude) presents some peculiarities which need some analysis.

A proper reference system.

As a first step let’s define a reference system for the problem under investigation: the universally accepted reference system in aeronautical matters is the following:

1. generally the origin is the center of gravity of the aircraft

2. the x axis is forward, the rotation around the x axis is named “roll” and the angle of roll is named φ or ɸ (phi)

3. the y axis in along the right wing, the rotation around the y axis is named “pitch” and the angle of pitch is named ϑ or θ (theta)

4. the z axis is downward, the rotation around the z axis is named “yaw” and the yaw angle is named psi (ψ)

All those angles (together with the AoA (Angle of Attack or Incidence ) are of paramount importance for the pilot in order to understand where he is going with the aircraft and above all how the aircraft is positioned respect to the ground (attitude), and need to be shown to the pilot in the cockpit in several redundant way depending on the type of aircraft: on the HUD (Head Up Display), and/or on the MHDD (Multifunction Head Down Display), and/or on the EFIS (Electronic Flight Instrument System), and/or on the ADI (Attitude Director Indicator), etc. They need to be carefully calculated in every condition.

The angles defined above are commonly referred as the "Euler angles", from the famous mathematicians of the 18th century who defined them first.

It is very important to add that by convention, changing from one attitude to another is described applying three different rotations: the first rotation by an angle ψ around the Z axis, then the rotation of an angle θ around the Y axis and finally a rotation of an angle φ around the X axis. The exact sequence of rotations is part of the reference system definition and defines the Euler angles, as standardize in the aeronautical world (in fact the same changes in attitude can be described choosing different sequence of rotation with different angles’ values).

Auxiliary reference system

Now let’s pose the following question: is there a standard definition (absolute definition) for an attitude defined as

$\psi =0$   $\theta =0$   $\varphi =0$

The answer is obviously “yes”, and this condition represent an auxiliary reference system in the real world, that we can think always linked to the aircraft, with the x axis oriented to north , the z axis oriented towards the center of the earth and the y axis oriented toward east direction. If we start “measuring” the $\psi$ angle from this reference system, it represents the heading of the aircraft, that is the direction of the fuselage respect to north. In this reference system the plane xy is parallel to the earth surface. When the aircraft changes its attitude, the changes are measured against this auxiliary reference system, according to the Euler angles convention. This auxiliary reference system is parallel to the inertial reference system, a fixed reference system with its origin located at any point on the Earth surface. This frame is commonly referred to as the North-East-Down, or NED frame.

Angular velocities

Attitude changes happen when one or more Euler angles change over time, that is when the aircraft is subject to a rotational velocity. Let’s define the three angular velocities associated to the Euler angles as:

psidot   angular velocity vector around z axis

thetadot  angular velocity vector around y axis after psi rotation

phidot   angular velocity vector around x axis after psi rotation and theta rotation

psidot, thetdot and phidot can be represented as three vectors, whose directions represent the rotational axis and the direction of rotation (according to the right hand rule), and the length represent the value of the rotational velocity.

It’s important to say that phidot, thetadot and psidot are not the three component of the global angular velocity vector responsible for changing the attitude of the aircraft, because they are defined in three different reference systems.

The following video explains better the above concept.

As a side note, I developed the software which displays the animation in C/C++ with the support of the OpenSceneGraph Libraries. Who is interested on the source code can contact me. The model is the Development Aircraft (DA) 7, the seventh prototype of Eurofighter, in charge of Alenia (now Leonardo) where I worked many years in the System Operability and Simulation department, where we were in charge to provide simulation tools (real time Flight Simulators) in support of Aircraft development. The 24th of January 1997 was a very exciting day, when DA7 took off for the maiden flight handled by Maurizio Cheli. I remember the work done with Maurizio in preparing the Display for the 1997 Paris air show on Alenia EF Flight Simulator and all the facilities we developed in order to give him the best to prepare the Display, including an accurate reconstruction of the Le Bourget airport to give him the write references to perform correctly the planned maneuvers respecting the assigned space and time.

Coming back to our problem, the described model (Euler Angles) does not represent the real evolution of the aircraft, that normally changes its attitude according to the instantaneous overall angular velocity represented by a vector as explained later, and not with 3 different sequential rotations. But, as already written in the Euler Parameters post, according to the Euler's rotation theorem, the 3 subsequent rotations can be seen as single overall rotation, and vice versa.

Let’s now define:

$\vec{\omega }$  overall angular velocity vector, that is the vector that with its direction defines the axis of rotation and with its length the angular velocity of rotation, and that is the responsible in the real world, of aircraft attitude changes.

The three components of $\vec{\omega }$ on the aircraft reference systems are universally known as:

P component of $\vec{\omega }$ along x axis

Q component of $\vec{\omega }$ along y axis

R component of $\vec{\omega }$ along z axis

Of course there is a strict relationship among P, Q, R and phidot, thetadot and psidot, as explained later on.

Euler angles rotation matrices

Let’s now consider the aircraft oriented as the auxiliary reference system (defined above) i.e. with $\psi =0\theta =0\varphi =0$ , and consider a generic point (vector) with coordinated P =(x_s, y_s, z_s) in the auxiliary reference system. We want to rotate the aircraft to assume an attitude described by the three angles $\psi ={\psi }_v\theta ={\theta }_v\varphi ={\varphi }_v$ and calculate the coordinates of P respect to this new aircraft reference system.

As a first step let’s apply the first rotation of an angle $\psi ={\psi }_v$ around the z axis.

This defines a new intermediate reference system that we label X1Y1Z1. The situation is described in the following video, where the aircraft is looked at from below.

In the video we highlight the point P, that in the XYZ reference system has coordinates (0.0,yp,0.0).

From the video is quite obvious that the same point P, in the X1Y1Z has coordinates (yp*sin(ψ),yp*cos(ψ),0.0). Adopting a matrix notation, we can write:

$\left(\begin{array}{c}{y}_p\ast \sin (\psi )\\ y_p\ast \cos (\psi )\\ 0\end{array}\right)=\left(\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ c_{31}& c_{32}& c_{33}\end{array}\right)\left(\begin{array}{c}0\\ y_p\\ 0\end{array}\right)$  The effort is to find all c_ij of the matrix

Considering that the rotation is around Z axis, the z coordinates during transformation remain unchanged. This means that c₃₁=c₃₂=0, and c₃₃=1. In addition the matrix is unitary and orthogonal and this ensure that c₂₃=c₁₃=0. So our problem becomes

$\left(\begin{array}{c}{y}_p\ast \sin (\psi )\\ y_p\ast \cos (\psi )\\ 0\end{array}\right)=\left(\begin{array}{ccc}{c}_{11}& c_{12}& 0\\ c_{21}& c_{22}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}0\\ y_p\\ 0\end{array}\right)$

From the above matrix expression is evident that c₁₂ = sin(ψ) and c₂₂=cos(ψ) and thus

$\left(\begin{array}{c}{y}_p\ast \sin (\psi )\\ y_p\ast \cos (\psi )\\ 0\end{array}\right)=\left(\begin{array}{ccc}{c}_{11}& \sin (\psi )& 0\\ c_{21}& \cos (\psi )& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}0\\ y_p\\ 0\end{array}\right)$

But for the orthogonality condition we have

$c_{11}\ast \sin (\psi )+c_{21}\ast \cos (\psi )=0$

and for unitary condition

$c_{11}\ast \cos (\psi )-c_{21}\ast \sin (\psi )=1$

Solving the two condition leads to:

c11=cos(ψ), c21=-sin(ψ)

The transformation matrix becomes:

$\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right)=\left(\begin{array}{ccc}\cos {\psi }_v& \sin {\psi }_v& 0\\ -\sin {\psi }_v& \cos {\psi }_v& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right)=M_{\psi }\left(\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right)$   1)

and x₁, y₁, and z₁ are the coordinates of P in the intermediate reference system X1Y1Z1

Now let’s apply a second rotation of an angle _{$\theta ={\theta }_v$} around the y axis of the X1Y1Z1 intermediate reference system, which identify an additional intermediate reference system that we label X2Y2Z2. Repeating the same reasoning as for psi, the transformation matrix becomes

$\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right)=\left(\begin{array}{ccc}\cos {\theta }_v& 0& -\sin {\theta }_v\\ 0& 1& 0\\ \sin {\theta }_v& 0& \cos {\theta }_v\end{array}\right)\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right)=M_{\theta }\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right)$   2)

And now to achieve our result, let’s apply the final rotation of an angle $\varphi ={\varphi }_v$ around the x axis of the X2Y2Z2 reference system. We obtain the final reference system of the aircraft after the three rotation, and we label it as X_vY_vZ_v. The related transformation matrix is:

$\left(\begin{array}{c}{x}_v\\ y_v\\ z_v\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos {\varphi }_v& \sin {\varphi }_v\\ 0& -\sin {\varphi }_v& \cos {\varphi }_v\end{array}\right)\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right)=M_{\varphi }\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right)$   3)

In order to pass from the original coordinates in the original reference systems to the coordinated in the X_vY_vZ_v reference system, we need to substitute relation 1) in 2, then relation 2) in 3), obtaining:

$\left(\begin{array}{c}{x}_v\\ y_v\\ z_v\end{array}\right)=M_{\varphi }\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right)=M_{\varphi }{M}_{\theta }\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right)=M_{\varphi }{M}_{\theta }{M}_{\psi }\left(\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right)=M_{\psi \theta \varphi }\left(\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right)$   4)

where $M_{\psi \theta \varphi }$ is

$M_{\psi \theta \varphi }=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos {\varphi }_v& \sin {\varphi }_v\\ 0& -\sin {\varphi }_v& \cos {\varphi }_v\end{array}\right)\left(\begin{array}{ccc}\cos {\theta }_v& 0& -\sin {\theta }_v\\ 0& 1& 0\\ \sin {\theta }_v& 0& \cos {\theta }_v\end{array}\right)\left(\begin{array}{ccc}\cos {\psi }_v& \sin {\psi }_v& 0\\ -\sin {\psi }_v& \cos {\psi }_v& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right)$  5)

We omit the boring work to perform operation in 5), giving just the result, leaving the reader to verify it (and anyway the result can be found in thousands of paper surfing the internet):

$M_{\psi \theta \varphi }=\left(\begin{array}{ccc}\cos (\theta )\cos (\psi )& \cos (\theta )\sin (\psi )& -\sin (\theta )\\ -\sin (\psi )\cos (\varphi )+\sin (\varphi )\sin (\theta )\cos (\psi )& \cos (\psi )\cos (\varphi )+\sin (\varphi )\sin (\theta )\sin (\psi )& \sin (\varphi )\cos (\theta )\\ \sin (\psi )\sin (\varphi )+\cos (\varphi )\sin (\theta )\cos (\psi )& -\cos (\psi )\sin (\varphi )+\cos (\varphi )\sin (\theta )\sin (\psi )& \cos (\varphi )\cos (\theta )\end{array}\right)$ 6)

Now we have all instrument to put in relation psidot, thetadot and phidot with P,Q and R

Angular velocities relationships

As first approach we try to solve the problem to know the Euler angles, starting from the angular velocities ( psidot, thetadot and phidot). If I know the angular velocities, I should be able, in some way, to find the relative angles, with some sort of integration process. Let’s solve this problem.

As a first step, we need to “transport” all the three velocities in the same reference system. In fact, as explained previously and shown in a previous movie, the three velocities are defined in three different reference system: in order to operate on them all three must be ported in the same reference system, like the last after all three rotations. Fortunately we have the instrument to do that by means of the transformation matrices just defined.

Let’s start considering that for phidot we don’t need to do anything, being phidot already defined in the (X_v, Y_v, Z_v). Remember that phidot is responsible for the final rotation around the X axis, from (X₂, Y₂, Z₂) to (X_v, Y_v, Z_v), but having phidot only the X component, it remains unchanged during rotation. The vector phidot , in the reference system (X_v, Y_v, Z_v), has thus the expression

$\dot{\overrightarrow{{\varphi }_v}}=\left(\begin{array}{c}\dot{\varphi }\\ 0\\ 0\end{array}\right)$   7)

As far as thetadot is concerned, we need to apply the theta transformation to transport it in the (X_v, Y_v, Z_v) reference system, i.e:

$\dot{\overrightarrow{{\theta }_v}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \varphi & \sin \varphi \\ 0& -\sin \varphi & \cos \varphi \end{array}\right)\left(\begin{array}{c}0\\ \dot{\theta }\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ \dot{\theta }\cos (\varphi )\\ -\dot{\theta }\sin (\varphi )\end{array}\right)$   8)

To transport psidot we need to apply the complete transformation matrix. Anyway, let’s consider that psidot has only the third component different from 0 whose value is $\dot{\psi }$ . This implies that the result is $\dot{\psi }$ multiplied by the last column of matrix 6) and thus

$\dot{\overrightarrow{{\psi }_v}}=\left(\begin{array}{c}-\dot{\psi }\sin (\theta )\\ \dot{\psi }\sin (\varphi )\cos (\theta )\\ \dot{\psi }\cos (\varphi )\cos (\theta )\end{array}\right)$   9)

Now the three angular velocities are defined in the same reference system. Summing up (vectorially) these three velocities results in the total angular velocity, whose components are, by definition, P,Q and R, as defined at the beginning of this post.

We can thus write:

$\vec{\omega }=\left(\begin{array}{c}P\\ Q\\ R\end{array}\right)=\left(\begin{array}{c}\dot{\varphi }\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}0\\ \dot{\theta }\cos (\varphi )\\ -\dot{\theta }\sin (\varphi )\end{array}\right)+\left(\begin{array}{c}-\dot{\psi }\sin (\theta )\\ \dot{\psi }\sin (\varphi )\cos (\theta )\\ \dot{\psi }\cos (\varphi )\cos (\theta )\end{array}\right)$   10)

From 10) the following equation can be derived:

$P=\dot{\varphi }-\dot{\psi }\sin (\theta )$      11a)

$Q=\dot{\theta }\cos (\varphi )+\dot{\psi }\sin (\varphi )\cos (\varphi )$    11b)

$R=-\dot{\theta }\sin (\varphi )+\dot{\psi }\cos (\varphi )\cos (\theta )$   11c)

Fantastic, we reach some results! We have now a relationship which links (P,Q,R) and $(\dot{\psi },\dot{\theta },\dot{\psi })$ . Unfortunately we are looking for a relationship which links $(\dot{\psi },\dot{\theta },\dot{\psi })$ to (P,Q,R), because the last ones, as discussed in other posts, are the known entities. To solve the problem let’s write 11a), 11b) and 11c) in a matricial form:

$\vec{\omega }=\left(\begin{array}{c}P\\ Q\\ R\end{array}\right)=\left(\begin{array}{ccc}1& 0& -\sin (\theta )\\ 0& \cos (\varphi )& \sin (\varphi )\cos (\theta )\\ 0& -\sin (\varphi )& \cos (\varphi )\cos (\theta )\end{array}\right)\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)$  12)

Relation 12) is in form

$\left(\begin{array}{c}P\\ Q\\ R\end{array}\right)=M\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)$ where $M=\left(\begin{array}{ccc}1& 0& -\sin (\theta )\\ 0& \cos (\varphi )& \sin (\varphi )\cos (\theta )\\ 0& -\sin (\varphi )& \cos (\varphi )\cos (\theta )\end{array}\right)$   13)

If we left multiply the first and second member of 13) by M^-1 we have:

$M^{-1}\left(\begin{array}{c}P\\ Q\\ R\end{array}\right)=M^{-1}M\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)=\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)$  and thus

$\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)=M^{-1}\left(\begin{array}{c}P\\ Q\\ R\end{array}\right)$  14)

We need to find the inverse of matrix M. It is not a big problem: we need first to obtain the transpose matrix of M, then to make its adjoint matrix a then divide all the members by the determinant of M. Applying these rules leads to:

$\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)=\left(\begin{array}{ccc}1& \tan (\theta )\sin (\varphi )& \tan (\theta )\cos (\varphi )\\ 0& \cos (\theta )& -\sin (\varphi )\\ 0& \frac{\sin (\varphi )}{\cos (\theta )}& \frac{\cos (\varphi )}{\cos (\theta )}\end{array}\right)\left(\begin{array}{c}P\\ Q\\ R\end{array}\right)$  15)

which leads to the following system of differential equations:

$\dot{\varphi }=P+\tan (\theta )\sin (\varphi )Q+\tan (\theta )\cos (\varphi )R$  16a)

$\dot{\theta }=\cos (\varphi )Q-\sin (\varphi )R$   16b)

$\dot{\psi }=\frac{\sin (\varphi )}{\cos (\theta )}Q+\frac{\cos (\varphi )}{\cos (\theta )}R$  16c)

From equation 16a), 16b) and 16c) it is evident the difficulty to solve the problem (remember: find $\varphi ,\theta ,\psi$ : even if we find a smart algorithm to solve numerically the system of differential equations, we are not able to face the condition where $\theta =\frac{\pi }{2}$ . In fact in this case we have equations 16a) with $\tan (\theta )$ going to $\mathrm{\infty }$ , and in equation 16c) two denominators going to 0. In simple words the system of differential equations becomes undefined: note that this not true only for the condition $\theta =\frac{\pi }{2}$ , but also in an interval around this value: being the word length of a computer limited, there is the risk on overflow during computation. This problem reflects the fact that for $\theta =\frac{\pi }{2}$ the angle $\varphi$ is indeterminate, condition known as “gimbal lock”. The condition $\theta =\frac{\pi }{2}$ seems to be not quite frequent, but impose the pilot not to reach that condition on a flight simulator could be reasonable for a transport aircraft, but think about a flight simulator for advanced fighter: this limitation would be simply unacceptable. This means that we need to find some other solution for calculation of $\varphi ,\theta ,\psi$

From Matrix transformation to Euler angles

The work done in this post has not been a waste of time, and it will be useful as soon as we develop a new strategy to solve the problem in the following posts. Let’s assume now to manage to calculate the matrix 6) ( $M_{\psi \theta \varphi }$ ) in an efficient and reliable way (this will be done in other posts), and let $M_{\psi \theta \varphi }$ be;

$M_{\psi \theta \varphi }=\left(\begin{array}{ccc}{m}_{11}& m_{12}& m_{13}\\ m_{21}& m_{22}& m_{23}\\ m_{31}& m_{32}& m_{33}\end{array}\right)$   7)

we would have immediately that

$\theta =-\arcsin (m_{13})\varphi =\arctan (\frac{m_{23}}{m_{33}})\psi =\arctan (\frac{m_{12}}{m_{11}})$

Remember also that we are looking for solution able to be integrated in a numerical computer: in all modern HLL (High Level Language) there is the built-in mathematical instruction atant2(x,y) to calculate the $\arctan (\frac{x}{y})$ with no risk to have an overflow condition in case of y=0. In the next post we will introduce the Euler parameters, another way to see the $\varphi ,\theta ,\psi$ matrix transformation, then the quaternions, a smart interpretation of the Euler parameters, than a smart algorithm to find the quaternion components knowing P,Q and R, and final an alternative way to find the matrix $M_{\psi \theta \varphi }$ and thus, for what written above, to solve our problem to find $\varphi ,\theta ,\psi$