Home » Space Flight/Orbital Mechanics Excercise » Euler Angles, Euler’s Equations, Yaw, Pitch and Roll Angles

# Euler Angles, Euler’s Equations, Yaw, Pitch and Roll Angles

Contents

## Euler Angles

The direction cosine matrix of an orthogonal transformation from XYZ to xyz is Q. Find the Euler angles fi , theta and psi for this transformation.

```clc;
clear all;
% Direction cosine matrix
Q = [-0.32175 0.89930 -0.29620;
0.57791 -0.061275 -0.81380;
-0.7500 -0.43301, -0.5000];
fi = atan2d(Q(3,1),-Q(3,2)); % Precession angle [0..360]
fprintf('Precession angle   fi    = %4.2f [deg]\n',fi);
theta = acosd(Q(3,3));        % Nutation angle [0..180]
fprintf('Precession angle   theta = %4.2f [deg]\n',theta);
psi = atan2d(Q(1,3),Q(2,3)); % Spin angle [0..360]
fprintf('Precession angle   fi    = %4.2f [deg]\n\n',psi);```
```Precession angle   fi    = 300.00 [deg]
Precession angle   theta = 120.00 [deg]
Precession angle   fi    = 200.00 [deg]```

At a given instant, the unit vectors of a body frame are

```i = [ 0.40825 -0.40825   0.81649];  % IJK inertial frame
j = [-0.10102 -0.90914  -0.40405];
k = [ 0.90726  0.082479 -0.41240];
% The angular velocity is
wX = [-3.1; 2.5; 1.7];                 % [rad/s] IJK inertial frame
% Calculate wp , wn and ws (the precession, nutation and spin rates) at
% this instant.
QXx = [i;j;k];
wx = QXx*wX;
w_x = wx(1);w_y = wx(2);w_z = wx(3);
% The components of the angular velocity in the body frame are
wx(1),wx(2),wx(3));
% Calculating the Euler angles fi,theta and psi from direction cosine matrix.
fi = atan2d(QXx(3,1),-QXx(3,2)); % Precession angle [0..360]
fprintf('Precession angle   fi    = %4.2f [deg]\n',fi);
theta = acosd(QXx(3,3));        % Nutation angle [0..180]
fprintf('Precession angle   theta = %4.2f [deg]\n',theta);
psi = atan2d(QXx(1,3),QXx(2,3)); % Spin angle [0..360]
fprintf('Precession angle   fi    = %4.2f [deg]\n\n',psi);
% Now we can obtain Euler angle rates from the angular velocitis wx,wy,wz
wp = 1/sind(theta)*(w_x*sind(psi) + w_x*cosd(psi));
wn = w_x*cosd(psi) - w_y*sind(psi);
ws = -1/tand(theta)*(w_x*sind(psi) + w_y*cosd(psi)) + w_z;
wp,wn,ws);```
```wx = -0.8982[rad/s]  wy = -2.6466[rad/s]  wz = -3.3074[rad/s]

Precession angle   fi    = 95.19 [deg]
Precession angle   theta = 114.36 [deg]
Precession angle   fi    = 116.33 [deg]

The mass moments of inertia of a body about the principal body frame axes with origin at the center of mass G are

```A = 1000;           % [kg*m^2]
B = 2000;           % [kg*m^2]
C = 3000;           % [kg*m^2]
% The Euler angles in radians are given as functions of time in seconds
% as follows:
dt = 0.01;
t = 0:dt:(10+2*dt);
fi    = rem(2*t.*exp(-0.05*t),2*pi);
theta = rem(0.02 + 0.3*sin(0.25*t),2*pi);
psi   = rem(0.6*t,2*pi);
% Compute derivatives for Euler angles
wp = diff(fi)/dt;      % d(fi)/d(t)
wn = diff(theta)/dt;   % d(theta)/d(t)
ws = diff(psi)/dt;     % d(psi)/d(t)

dwp = diff(wp)/dt;     % d(wp)/d(t)
dwn = diff(wn)/dt;     % d(wn)/d(t)
dws = diff(ws)/dt;     % d(ws)/d(t)
% Evaluating all of these quantities  at t =10 [s] yields
n = size(t);
fprintf('fi = %4.2f [deg]       wp = %6.5f [rad/s]     dwp = %6.5f [rad/s^2]\n',...
fi(n(2)-2)*180/pi,wp(n(2)-2),dwp(n(2)-2));
fprintf('theta = %4.2f [deg]     wn = %6.5f [rad/s]     dwn = %6.5f [rad/s^2]\n',...
theta(n(2)-2)*180/pi,wn(n(2)-2),dwn(n(2)-2));
fprintf('psi = %4.2f [deg]       ws = %6.5f [rad/s]     dws = %6.5f [rad/s^2]\n\n',...
psi(n(2)-2)*180/pi,ws(n(2)-2),dws(n(2)-2));```
```fi = 335.03 [deg]       wp = 0.60608 [rad/s]     dwp = -0.09090 [rad/s^2]
theta = 11.43 [deg]     wn = -0.06014 [rad/s]     dwn = -0.01118 [rad/s^2]
psi = 343.77 [deg]       ws = 0.60000 [rad/s]     dws = 0.00000 [rad/s^2]```

At t = 10 [s], find the net moment about G Calculating angular velocity components using the Euler angle rates at t = 10

```fi    =  fi(n(2)-2)*180/pi;     wp = wp(n(2)-2);    dwp = dwp(n(2)-2);
theta =  theta(n(2)-2)*180/pi;  wn = wn(n(2)-2);    dwn = dwn(n(2)-2);
psi   =  psi(n(2)-2)*180/pi;    ws = ws(n(2)-2);    dws = dws(n(2)-2);
% Angular velocity components in body fixed xyx frame
wx = wp*sind(theta)*sind(psi) + wn*cosd(psi);
wy = wp*sind(theta)*cosd(psi) - wn*sind(psi);
wz = ws + wp*cosd(theta);
% Taking the time derivative of each of these equations we will obtain
% angular acceleration of the space craft in the body fixed frame
dwx = wp*wn*cosd(theta)*sind(psi) + wp*ws*sind(theta)*cosd(psi)-...
wn*ws*sind(psi) + dwp*sind(theta)*sind(psi) + dwn*cosd(psi);
dwy = wp*wn*cosd(theta)*cosd(psi) - wp*ws*sind(theta)*sind(psi)-...
wn*ws*cosd(psi) + dwp*sind(theta)*cosd(psi) - dwn*sind(psi);
dwz = -wp*wn*sind(theta) + dwp*cosd(theta) + dws;
wx,wy,wz);
dwx,dwy,dwz);```
```wx = -0.09131  [rad/s]     wy = 0.09855 [rad/s]      wz = 1.19405  [rad/s]

## Euler’s equations

```Mx_net  = A*dwx + (C - B)*wy*wz;
My_net  = B*dwy + (A - C)*wx*wz;
Mz_net  = C*dwz + (B - A)*wx*wy;
fprintf('Mx_net = %4.2f [N*m], My_net = %4.2f [N*m], Mz_net = %4.2f [N*m]\n\n',...
Mx_net,My_net,Mz_net);
% Now we will project the components of relative acceleration onto the
% axes of the inertial frame.
% Required orthogonal transformation matrix
QxX = [-sind(fi)*cosd(theta)*sind(psi) + cosd(fi)*cosd(psi), ...
-sind(fi)*cosd(theta)*cosd(psi) - cosd(fi)*sind(psi),sind(fi)*sind(theta);
cosd(fi)*cosd(theta)*sind(psi) + sind(fi)*cosd(psi),...
cosd(fi)*cosd(theta)*cosd(psi) - sind(fi)*sind(psi),-cosd(fi)*sind(theta);
sind(theta)*sind(psi),   sind(theta)*cosd(psi),    cosd(theta)
]
% Angular velocity rates in inertial frame;
% Transition from Body Fiexed to Inertial Frame
w_xyz = [dwx;dwy;dwz];
w_XYZ = QxX*w_xyz;
fprintf('w_xyz = [ %6.5f  %6.5f  %6.5f ] rad/s^2\n',w_xyz);
fprintf('w_XYZ = [ %6.5f  %6.5f  %6.5f ] rad/s^2\n\n',w_XYZ);
% To check our results
% Transition from Inertial Frame to Body Fiexed
w_xyz = QxX'*w_XYZ;
fprintf('w_XYZ = [ %6.5f  %6.5f  %6.5f ] rad/s^2\n',w_XYZ);
fprintf('w_xyz = [ %6.5f  %6.5f  %6.5f ] rad/s^2\n\n',w_xyz);
% We get the desired values.```
```Mx_net = 181.08 [N*m], My_net = 218.18 [N*m], Mz_net = -254.62 [N*m]

QxX =

0.7548    0.6505   -0.0837
-0.6536    0.7352   -0.1797
-0.0554    0.1903    0.9802

w_xyz = [ 0.06341  0.00006  -0.08187 ] rad/s^2
w_XYZ = [ 0.05475  -0.02669  -0.08375 ] rad/s^2

w_XYZ = [ 0.05475  -0.02669  -0.08375 ] rad/s^2
w_xyz = [ 0.06341  0.00006  -0.08187 ] rad/s^2```

## Yaw, Pitch and Roll Angles

For the given direction cosine matrix obtain YAW, PITCH AND ROLL ANGLES

```QXx =[-0.32175   0.89930   -0.29620;
0.57791  -0.061275  -0.81380;
-0.75     -0.43301   -0.5  ];
yaw   = atan2d(QXx(1,2),QXx(1,1));
pitch = asin(-QXx(1,3))*180/pi;
roll  = atan2d(QXx(2,3),QXx(3,3));
fprintf('Yaw = %6.5f [deg]  Pitch = %6.5f [deg]   Roll = %6.5f [deg]\n'...
,yaw,pitch,roll);
% Now lets optain direction cosine matrix from Yaw,Pitch and Roll angles
R1_roll  =  [1            0           0;
0            cosd(roll)  sind(roll);
0           -sind(roll)  cosd(roll);
]
R2_pitch = [cosd(pitch)  0     -sind(pitch);
0            1       0;
sind(pitch)  0      cosd(pitch)]

R3_yaw   =   [ cosd(yaw)   sind(yaw)   0;
-sind(yaw)   cosd(yaw)   0;
0            0          1]
QXx      = R1_roll*R2_pitch*R3_yaw```
```Yaw = 109.68611 [deg]  Pitch = 17.22951 [deg]   Roll = 238.43337 [deg]

R1_roll =

1.0000         0         0
0   -0.5235   -0.8520
0    0.8520   -0.5235

R2_pitch =

0.9551         0   -0.2962
0    1.0000         0
0.2962         0    0.9551

R3_yaw =

-0.3369    0.9416         0
-0.9416   -0.3369         0
0         0    1.0000

QXx =

-0.3218    0.8993   -0.2962
0.5779   -0.0613   -0.8138
-0.7500   -0.4330   -0.5000```

Published with MATLAB® 7.10