**1. INTRODUCTION**Solution of the six-degree-of-freedom flight equations for aircraft and missiles continues to represent one of the most important application areas for analog, hybrid, and digital computer systems. Important computer requirements such as accuracy and speed are very much dependent on the choice of axis system for the translation equations of motion. In this connection it is well known that the flight-path axis system (or wind-axis system) makes much lower accuracy and speed demands on the computer than does the body-axis system. /1/, /2/

*Flight-path or wind axis (w), stability axis (s) and body axis (b)*

Despite this, a number of current computer mechanizations continue to use body axes for solving the translational equations of motion. Because of this, unnecessary demands on accuracy or frequency response are placed on the computer, and many mechanizations which could be all-analog or all-digital have shifted to hybrid implementation. Even if the mechanization is hybrid from the outset, there is considerable advantage to be gained by using an efficient axis system.

__The purpose of this paper is to point out again the advantages of flight path axes and to summarize the overall equation requirements for solving the six-degree of freedom flight equations.__

**2. BODY-AXIS TRANSLATIONAL EQUATIONS**

For comparison, we present first the body-axis translational equations. The body axes

**,**

*Xb***, and**

*Yb***are defined as a right-hand set fixed to the vehicle with the**

*Zb***axis along the longitudinal axis and the**

*Xb***axis directed downward for normal level flight. The components of the total vehicle velocity vector**

*Zb***along the**

*Vp***,**

*Xb***, and**

*Yb***axes, respectively, are**

*Zb***,**

*Ub***, and**

*Vb***(see figure 2.1).**

*Wb*The components of the body-axis angular velocity vector

**(and hence the vehicle angular velocity vector) along the**

*Ω***,**

*Xb***, and**

*Yb***axes are**

*Zb***,**

*Pb***, and**

*Qb***, i.e., roll rate, pitch rate, and yaw rate, respectively. Here we assume**

*Rb***and**

*Vp***represent vehicle translational and rotational velocity vectors as viewed from an inertial (nonaccelerating) frame of reference. If we denote the external forces along a set of coordinate axes by**

*Ω***,**

*X***, and**

*Y***, respectively, then Euler’s translational equations of motion, which are obtained by summing forces along the coordinate axes, are the following:**

*Z*where

**is the mass of the vehicle.**

*m*The inefficiency of these equations in body axes is immediately apparent when one considers the approximate size of the various terms. Let the vehicle be, say, a Mach 2 aircraft with

**feet per second. A reasonable upper limit on pitch-rate**

*Vmax = 2000***, might be 2 radians per second. Thus, the term**

*Q***in equation (2.3) might be as large as 4000 ft/s2 or 125 g’s. On the other hand**

*UbQb***, the normal acceleration due to the external force (primarily gravity and aerodynamic lift) may have an upper limit of several g’s. Hence artificial accelerations which are perhaps 20 to 50 times greater than the actual accelerations are introduced because of the high rotation rates which the body-axes experience. This means much less favorable computer scaling and hence much poorer solution accuracy for a given computer precision. Furthermore, the high-speed dynamics of the rotational equations are coupled into the translational equations, thus placing severe dynamic response requirements on the computer.**

*Zb/m*The use of flight-path axes greatly alleviates these problems. As we shall see in the next section, the flight-path axes allow a more efficient calculation of the aerodynamic angle of attack

**(alpha) and the aerodynamic angle of side slip**

*α***(beta) than body axes allow. Using body axes and assuming that the ambient air mass is not moving relative to the inertial frame used to define**

*β***, then the following formulas can be used to obtain**

*Vp***,**

*α***, and velocity magnitude**

*β***, from the body-axis velocity components**

*V***,**

*Ub***, and**

*Vb***:**

*Wb*.

**3. FLIGHT-PATH AXIS TRANSLATIONAL EQUATIONS**

Next consider the flight-path axes

**,**

*Xw***,**

*Yw***, shown in figure 3.1.**

*Zw*These differ from the body axes

**,**

*Xb***,**

*Yb***by the angle of attack**

*Zb***and the angle of side slip**

*α***, as sh own in the figure. To rotate from body axes to flight-path axes one first pitches the body axes about**

*β***through**

*Yb***. This defines an intermediate axis system**

*-α***,**

*Xs***,**

*Ys***, called stability axes. One then yaws about**

*Zs***, through**

*Zs***, which defines the flight-path axes**

*β***,**

*Xw***, and**

*Yw***. Note that**

*Zw***is the**

*Vp***component of vehicle velocity; the**

*Xw***and**

*Yw***components are zero by definition. Let us define the**

*Zw***,**

*Xw***, and**

*Yw***components of flight-path angular velocity relative to inertial space by**

*Zw***,**

*Pw***, and**

*Qw***, and the components of external force along**

*Rw*

*Xw,***and**

*Yw,***by**

*Zw***,**

*Xw***, and**

*Yw***, respectively. Then, since**

*Zw***and**

*Uw = Vp***, the translational equations (2.1), (2.2), and (2.3) referred to the flight-path axes become**

*Vw = Ww = 0*Solution of these three equations results in total velocity

**, flight-path axis yaw rate**

*Vp***, and flight-path axis pitch rate**

*Rw***.**

*Qw*Next consider the formulas for

**and**

*α***. Reference to figure 3.1 shows that**

*β***is directed along**

*α_dot***, with a component**

*Ys***along**

*α_dot cos β***. Thus**

*Yw***is equal to the difference between body-axis and flight-path axis angular rates along**

*α_dot cos β***. Therefore, from equation (3.3) we can write**

*Yw*where

**is the body-axis (not stability axis) angular rate along**

*Pb_s***and is given by (3.5)**

*Xs*Similarly, reference to figure 3.1 shows that

**is directed along**

*β_dot***and is equal to the difference between flight-path axis and body axis angular rates along**

*Zw***. Thus from equation (3.2) we get (3.6) where**

*Zw***is the body-axis angular rate along**

*Rb_s***and is given by (3.7). Note that**

*Zs***since the**

*Qb_s = Qb***body axis and**

*Yb***stability axis are coincident (3.8).**

*Ys*Equations (3.1), (3.4), and (3.6) can be integrated to yield total velocity

**, angle of attack**

*Vp***, and angle of side slip**

*α***. They do not present the scaling difficulty of equations (2.1), (2.2), and (2.3) in the body axes. The body-axis velocity components**

*β***,**

*Ub***, and**

*Vb***can be obtained from**

*Wb***,**

*Vp***, and**

*α***by the formulas (3.9) -3.11). Similarly, the flight-path-axis forces**

*β***,**

*Xw***, and**

*Yw***can be derived from body-axis force components**

*Zw***,**

*Xb***, and**

*Yb***by the formulas (3.12)-(3.14) where**

*Zb***,**

*Xs***, and**

*Ys***are the intermediate stability-axis force components.**

*Zs*Frequently the aerodynamic force components are computed along stability axes, in which case the power plant and gravity forces, computed in body axes, would be resolved into stability axes where the aerodynamic forces are added; then the total forces would be resolved into flight-path axes to allow use of equations (3.1), (3.4), and (3.6).

It should be noted that we have assumed throughout this section that the translational and rotational velocity vectors are velocities relative to an inertial reference frame. If the atmosphere through which the vehicle is flying can be considered to be fixed with respect to this inertial frame, then the velocity magnitude

**is the vehicle velocity relative to the atmosphere, and**

*Vp***and**

*α***represent the aerodynamic angle of attack and side slip, respectively. Using the approximation that the earth is flat, and with constant surface winds, it is possible to define the inertial reference frame as a frame attached to the atmosphere. Then all the formulas, as presented, are correct, and**

*β***,**

*α***and**

*β***can be used for computation of aerodynamic forces and moments.**

*Vp*Unfortunately, for a rotating spherical earth, axes fixed in the atmosphere are not inertial. If we consider such a frame to be inertial, we will make acceleration errors in equations (3.1), (3.2), and (3.3) of the order of

**where**

*V²/r0***is the vehicle velocity relative to an inertial frame with origin at the center of the earth and**

*V***is the radius of the earth. For illustration, consider a vehicle flying eastward in still air with a velocity**

*r0***relative to the atmosphere. Then**

*Va***, where**

*V ~ Va + r0 ωn cos L***is the earth spin rate and**

*ωn***is the latitude. If the vehicle is flying westward,**

*L***. For other headings**

*V ~ Va - r0 ωn cos L***lies between these values.**

*V*For example, if

**= 3000 ft/s and**

*Va***0 degrees,**

*L =***ranges between approximately 1600 and 4400 ft/s. The corresponding acceleration error in equations (3.1), (3.2), and (3.3), given by**

*V***, ranges between approximately 0.1 and 1 ft/s2. For many flight vehicles this is a negligible error. On the other hand, for a supersonic transport cruising eastward at 3000 ft,/s this lowers the required steady-state lift by about 3per cent, which could lower the drag significantly and hence make a noticeable difference in maximum range.**

*V²/r0*There appears to the authors to be no simple way to take these accelerations into account and still use a flight-path axis system referenced to the atmosphere. One could add an approximate correction acceleration in the vertical direction given by

**to the translatory forces in equations (3.1), (3.2), and (3.3). In fact, one could further simplify the computation by adding the term only to equation (3.3), based on the argument that most of the time a supersonic aircraft will be in near-level flight at cruise and that a moderate acceleration error during transient conditions can be accepted.**

*V²/r0*In this case equation (3.3) becomes

Note that this equation will exhibit acceleration errors of the order of

**for conditions far from level flight.**

*V²/r0***4. ROTATIONAL EQUATIONS**

The only reasonable axes to use for the rotational equations of motion are the body axes. If we let the external moment components along

**,**

*Xb***, and**

*Yb***, be**

*Zb***,**

*Lb***, and**

*Mb***, respectively, then summation of moments about the three body axes of a body symmetrical about the**

*Nb***plane leads to the equations:**

*XbZb*Here

**,**

*Ixx***, and**

*Iyy***are the moments of inertia about**

*Izz***,**

*Xb***, and**

*Yb***respectively, and**

*Zb***is the product of inertia of the symmetrical body. Note in equations (4.1), (4.2), and (4.3) that the second term in each equation represents a nonlinear inertial coupling term. For flight vehicles such as large transport aircraft which do not generate relatively high angular rates these terms often can be neglected. For many flight vehicles the roll rate**

*Ixz***has a maximum value which is considerably higher than pitch-rate**

*Pb***or yaw-rate**

*Qb***. Hence the**

*Rb***term in equation (4.1 ) often can be neglected in comparison with the**

*QbRb***and**

*RbPb***terms in equations (4.2) and (4.3), respectively.**

*PbQb*The third term in each of equations (4.1), (4.2), and (4.3) represents the effect of the product of inertia

**. If the**

*Ixz***,**

*x***,**

*y***body axes have been chosen to be almost coincident with the principal axes, this term may be negligible in all three equations, since**

*z***will be very small compared with the principal moments of inertia. For relatively low angular rates the nonlinear terms (**

*Izz***,**

*PbQb***² -**

*Pb***², and**

*Rb***) usually can be neglected and in any event**

*QbRb***² frequently can be neglected as small compared with**

*Rb***² in equation (4.2).**

*Pb***5. COMPUTATION OF EULER ANGLES FOR A FLAT EARTH**

Solution of equations (4.1), (4.2), and (4.3) results in computation of the body-axis angular velocity components

**, and**

*Pb Qb***. These must be integrated a second time to obtain the orientation of the vehicle body axes with respect to the desired references axes, typically Euler axes which point north, east, and toward the center of the earth (**

*Rb***,**

*Xe***, and**

*Ye***in figure 5.1 ).**

*Ze*This orientation usually is expressed in terms of the conventional aircraft Euler angles, i.e., heading-angle

**(psi), pitch angle**

*Ψ***(theta), and bank angle**

*Θ***(phi). These angles usually are computed from**

*Φ***,**

*Pb***and**

*Qb***by the following well-known equations:**

*Rb*Note that, since

**,**

*Pb***,**

*Qb***are the body-axis components of the vehicle angular velocity relative to an inertial reference system, there is a small error introduced by the angular velocity of axes which point north, east, and down relative to a spherical earth. This small error can be corrected, if necessary, using equations (7.1)-(7.3) in the next section.**

*Rb*The well known singularity of the Euler angle system at

**=**

*Θ***±π/2**can be avoided, if necessary, by computing direction cosines or quaternions /4/ instead of Euler angles, or by introducing a fourth angle. /3/

The use of quaternions rather than direction cosines should be considered if a system free of singularities is needed, since there are only four quaternions with a single redundancy as compared with nine direction cosines with six redundancies.

**5.1 Additional notes about Euler angles /6/**

The most popular application is to describe aircraft attitudes, normally using a Tait–Bryan convention so that zero degrees elevation represents the horizontal attitude. Tait–Bryan angles represent the orientation of the aircraft respect a reference axis system (world frame) with three angles which in the context of an aircraft are normally called

*Heading*

**(psi),**

*Ψ**Elevation*

**(theta)**

*Θ**and*

*Bank*

**(phi). When dealing with vehicles, different axes conventions are possible.**

*Φ**Tait–Bryan Euler angles in aircraft application (all zero degrees yield normal flight attitude north)*

When studying rigid bodies in general, one calls the xyz system space coordinates, and the XYZ system body coordinates. The space coordinates are treated as unmoving, while the body coordinates are considered embedded in the moving body. Calculations involving acceleration, angular acceleration, angular velocity, angular momentum, and kinetic energy are often easiest in body coordinates, because then the moment of inertia tensor does not change in time. If one also diagonalizes the rigid body's moment of inertia tensor (with nine components, six of which are independent), then one has a set of coordinates (called the principal axes) in which the moment of inertia tensor has only three components.

The angular velocity of a rigid body takes a simple form using Euler angles in the moving frame. Also the Euler's rigid body equations are simpler because the inertia tensor is constant in that frame.

"It is important to appreciate the range of the Euler angles:

−π ≤ ψ ≤ +π

−π/2 ≤ θ ≤ +π/2

−π ≤ φ ≤ +πIt is not obvious that pitch should be constrained to the range ±90◦. Consider a loop; as the pilot pulls up into the loop, the pitch increases to 80◦, then85◦, 86◦, 87◦, 88◦ and 89◦. At this point, the aircraft is still pointing forwards, but almost vertically upwards. In the next two degrees of motion, the aircraft passes through 90◦ and subsequently is pointing backwards, but at 89◦ to the horizon. Notice also that during this maneuverer, as the pitch angle passes through 90◦, the aircraft heading changes from forwards to backwards (an 180◦ instantaneous change in yaw) and also, the aircraft changes from upright to inverted (an 180◦ instantaneous change in roll)." /8/

German words: "Die Lage des körperfesten Achsensystems wird mit Hilfe der drei Eulerwinkel

**(**

*Ψ**Azimut, Gierwinkel*),

**(**

*Θ**Längsneigung, Nickwinkel*) und

**(**

*Φ**Hängewinkel, Rollwinkel*) beschrieben." /7/

**6. COMPUTATION OF VEHICLE POSITION FOR A FLAT EARTH**

Once the orientation of the flight vehicle with respect to the Euler axes has been established e.g., by means of the Euler angles, then it is possible to compute the velocity north

**, the velocity east**

*Ue***, and the velocity downward**

*Ve***, or its negative, the rate of climb**

*We***. Direct integration then yields the vehicle position.**

*h_dot*Determination of

**,**

*Ue***,**

*Ve***is complicated by the fact that the vehicle velocity vector**

*We***lies along the**

*Vp***wind axis and therefore must be resolved from wind to earth axes. Unfortunately, the complete orientation of the wind axes is known only relative to body axes; hence it is necessary to perform the resolution of**

*x***into earth axes by first resolving it into body axes,. then from body axes to earth axes. The resolution of**

*Vp***from wind axes to body axes is accomplished by the transformation given in equations (3.9), (3.10), and (3.11 ).**

*Vr*The resolution of vehicle velocity from body axes to earth axes can be accomplished by using direction cosines. Thus let

**,**

*l1***, and**

*l2***be the projections of a unit vector along the**

*l3***body axis onto the**

*x***,**

*Xe***, and**

*Ye***earth axes, respectively. Similarly, let**

*Ze***,**

*m1***, and**

*m2***be the projections of a unit vector along the**

*m3***body axis onto the**

*y***,**

*Xe***, and**

*Ye***earth axes, respectively. In the same way let**

*Ze***,**

*n1***, and**

*n2***be the projections of a unit vector the**

*n3***body axis along the**

*z***,**

*Xe***, and**

*Ye***earth axes, respectively. Then by definition**

*Ze*It is easy to show that the direction cosines are related to Euler angles by the following formulas:

(Equivalent formulas for direction cosines in terms of quaternions are given in reference 4.)

An alternative mechanization of equations (6.1 ) through (6.6) avoids computation of the direction cosines by instead performing successive resolutions of the velocity components

**,**

*Ub***, and**

*Vb***through the angles -**

*Wb***, -**

*Φ***, and -**

*Θ***. Consider the intermediate axis system**

*Ψ***’,**

*x***’,**

*y***’ in figure 5.1. Clearly the velocity components**

*z***’,**

*U***’, and**

*V***’ along**

*W***’,**

*x***’, and**

*y***’ are given by**

*z*Defining

**,**

*U"***, and**

*V"***as the velocity components along the intermediate axes**

*W"***,**

*x"***, and**

*y"***, we have (6.8). Finally, from figure 4.1 we see (6.9).**

*z"*Successive application of equations (6.7), (6.8), and (6.9) for

**,**

*U***, and**

*V***to obtain**

*W***,**

*Ue***, and**

*Ye***requires fewer mathematical operations than using equations (6.1) through (6.6). It therefore has computational advantages using either an analog or digital mechanization. In computing ground coordinates from**

*We***and**

*Ue***in equations (6.1 ) and (6.2), or equation (6.9), it is important to note that**

*Ve***and**

*Ue***represent airspeed components north and east, respectively. To convert them to groundspeed components, the north component of wind**

*Ve***must be subtracted from**

*Wz***, and the east component of wind**

*Ue***must be subtracted from**

*Wy***. Thus if**

*Ve***and**

*Sz***represent distance traveled north and east, respectively, then**

*Sy*Equations (6.10) and (6.11) are valid only for steady winds, since the implicit assumption was made that axes stationary relative to the atmosphere are inertial. Equations (3.1), (3.2), and (3.3) are referred to inertial space and correction terms must be added if the reference axes are not inertial.

**7. COMPUTATION OF VEHICLE EULER ANGLES AND POSITION FOR A ROTATING SPHERICAL EARTH**

In the previous section we presented the formulas for computing vehicle position over a flat earth with steady winds. We can use the same position formulas to obtain velocity north

**and velocity east**

*Ue***over a rotating spherical earth with radius**

*Ve***and angular rate**

*r0***. However,**

*ωN***and**

*Ue***will represent airspeed components and must be corrected to yield groundspeed components**

*Ve***and**

*Sx_dot***respectively. Furthermore, when a spherical earth is considered, it may be necessary to correct the vehicle angular rates used to compute Euler angles in order to take into account the rotating reference frame, as pointed out earlier in section 6. It can be shown that the body-axis components of the vehicle angular velocity relative to the Euler reference frame,**

*Sy_dot***,**

*Pbe***,**

*Qbe***are given by the following formulas: /5/**

*Rbe*Here

**is the radial distance of the vehicle from the center of the earth and is given by (7.4) where**

*r***is the radius of the earth and h is the vehicle altitude. In many cases we can substitute**

*r0***for**

*r0***in equations (7.1), (7.2), and (7.3) and still obtain sufficient accuracy.**

*r*The values of

**,**

*Pbe***, and**

*Qbe***in these equations are then used to compute the Euler angle rates. Thus by analogy with equations (5.1), (5.2), and (5.4)**

*Rbe*In many cases the computations involved in equations (7.1), (7.2), and (7.3) can be neglected i.e., we can assume that

**,**

*Pbe ~ Pb***,**

*Qbe ~ Qb***. This is particularly true if the overall six-degree-of-freedom computation involves a control system (automatic or human) which attempts to maintain**

*Rbe ~ Rb***,**

*Ψ***, and**

*Θ***at specified values. In any case the correction rates are the order of**

*Φ***. For example, if**

*Vp/r0***= 2000 ft/s, the correction rate is equal to approximately 0.005 degree per second. On the other hand, if the flight-vehicle problem includes a stable platform, the rate corrections given by equations (7.1), (7.2), and (7.3) may be important.**

*Vp*It should be noted that

**and**

*Sz_dot***in equations (7.1 ), (7.2), and (7.3) represent vehicle velocity components north and east, respectively, over the surface of a nonrotating earth with steady winds. On the other hand equations (6.1 ) and (6.2), or alternatively, equation (6.9), gives us**

*Sy_dot***and**

*Ue***, i.e., vehicle velocity components north and east, respectively, relative to the inertial reference frame for the translational equations of motion. We made the approximation in section 3 that this reference frame is fixed relative to the ambient atmosphere. We can account for the linear velocity of the atmospheric reference frame (but not the angular velocity) by noting that relative to the surface of a nonrotating earth it is moving northward with the northerly component of wind and eastward with the sum of the rate due to earth spin and that due to eastward component of wind. Thus we can write the following equations:**

*Ve*It can also be shown /5/ that the time rate of change of latitude and longitude are given by the formulas (7.10) and (7.11).

It should be noted that motion of a flight-vehicle over a rotating earth can be treated exactly, /5/ but that the exact translational equations referred to axes fixed relative to the ambient atmosphere are very complicated. Thus many of the computer mechanization advantages for flight-path axes are lost, and the computation of

**,**

*α***, and**

*β***is much less elegant. We have attempted in this section to describe how one can utilize the flight-path axis system and still correct approximately for the fact that the vehicle is flying over a rotating earth with surface winds. This approach should be adequate for all but the most exacting requirements, unless the vehicle reaches hypersonic speeds. For subsonic vehicles or supersonic vehicles traveling over relatively short distances the flat-earth equations in sections 3 through 6 should be adequate.**

*Vp***8. COMPUTATION OF AERODYNAMIC FORCES AND MOMENTS**

Aerodynamic forces and moments for flight vehicles normally are computed in stability axes. Thus let

**,**

*Xa***, and**

*Ya***be the aerodynamic forces along the**

*Za***.,**

*Xs***, and**

*Ys***stability axes shown in figure 3.1. Here**

*Zs***corresponds to drag**

*-Xa***, and**

*D***corresponds to lift**

*-Za***.**

*L*In addition to aerodynamic forces the flight vehicle experiences propulsion and gravity forces. These are most conveniently defined in body axes. Let us denote the propulsion force components along the

**and**

*x***body axes by**

*z***and**

*Xp***, respectively. The gravity force components along**

*Zp***,**

*x***, and**

*y***will by definition be**

*z***,**

*mgl3***, and**

*mgm3***respectively. The sum of the propulsive forces and gravity forces along body axes can then be resolved to stability axes, where the aerodynamic forces are added to obtain the total force component**

*mgn3***,**

*Xs***and**

*Ys***along the stability axes. Using equations (6.4), (6.5), and (6.6) to express the direction cosines**

*Zs***,**

*l3***,**

*m3***in terms of Euler angles, we obtain:**

*n3*Finally, reference to figure 3.1 shows that the force components

**,**

*Xw***, and**

*Yw***along the flight-path axes can be computed using the formulas (8.4) - (8.6).**

*Zw*These force components are used to mechanize the translational equations given in section 3. If the aerodynamic forces are given in body axes rather than stability axes, the modification of equations (8.1 ), (8.2), and (8.3) is apparent.

The external moments acting on the flight vehicle consist of aerodynamic moments

**,**

*La***, and**

*Ma***, normally given in stability axes, and power plant moments**

*Na***,**

*Lp***, and**

*Mp***, given in body axes. The total moments**

*Np***,**

*L***, and**

*M***in body axes become the following:**

*N*Again, if the aerodynamic moments are given in body axes, the simplification of equations (8.7), (8.8), and (8.9) is obvious (set

**= 0). Figure 8.1 shows a block diagram of the over-all six-degree-of-freedom equations for the case of a flat earth.**

*α**Figure 8.1*

*L. E. FOGARTY*

**,***R. M. HOWE***University of Michigan**

(REPRINT of the original with minor fixes and additions.)

**REFERENCES**

/1/ Howe R M

Coordinate systems for solving the three-dimensional flight equations

WADC Technical Note 55-747 June 1956

/2/ Howe R M

Coordinate systems and methods of coordinate transformation for three-dimensional flight equations

Proceedings of the First Flight Simulation Symposium WSPG

Special Report 9 September 1957

/3/ Greenwood D T

An extended Euler angle coordinate system for use with all-altitude aircraft simulators

WADC Technical Report 60-372 August 1960

/4/ Robinson A C

On the use of quaternions in simulation of rigid-body motion

WADC Technical Report 58-17 December 1958

/5/ Fogarty L E, Howe R M

Analog computer solution of the orbital flight equations

SIMULATION vol 1 no 1 Fall 1963 pp R-41 through R-47

/6/ http://en.wikipedia.org/wiki/Euler_angles

/7/ http://tumb1.biblio.tu-muenchen.de/publ/diss/mw/2004/holzapfel.pdf

/8/ David Allerton,

PRINCIPLES OF FLIGHT SIMULATION, Department of Automatic Control and Systems Engineering, The University of Shefﬁeld