Introduction to Rotation Matrices

Table of Contents
Getting Started Rotation Matrix to BODY from VEHICLE transform Example 1 Example 2

Getting Started

Let's start by creating an empty script and we'll start with the basic cleaning commands. Add these to your script.
clear;
clear deletes all of the variables in your workspace - this makes sure that every time we run our code, we're always using the data from this run of the script, and not a previous run.
clc;
clc deletes all of the old messages in the console or command window. This makes sure that the only messages we see on our command window come from this run of the script, and not from a previous run.
close all;
close all closes all of the previous plots and graph windows. This way you know that the plots you see generated are new, and not old data.

Rotation Matrix

Let's review our rotation matrices
Given a PSI rotation angle, the corresponding rotation matrix is given by
For a given rotation of
psi = 10;
psi_in_radians = deg2rad(psi);
 
Rz_psi = [ cos(psi_in_radians), sin(psi_in_radians), 0;
-sin(psi_in_radians), cos(psi_in_radians), 0;
0, 0, 1];
disp(Rz_psi);
0.984807753012208 0.17364817766693 0 -0.17364817766693 0.984807753012208 0 0 0 1
We can do the same for the other angles theta and phi
theta = 10;
theta_in_radians = deg2rad(theta);
Ry_the = [cos(theta_in_radians), 0, -sin(theta_in_radians);
0, 1, 0;
sin(theta_in_radians), 0, cos(theta_in_radians)];
disp(Ry_the);
0.984807753012208 0 -0.17364817766693 0 1 0 0.17364817766693 0 0.984807753012208
phi = 10;
phi_in_radians = deg2rad(phi);
Rx_phi = [ 1, 0, 0;
0, cos(phi_in_radians), sin(phi_in_radians);
0, -sin(phi_in_radians), cos(phi_in_radians)];
disp(Rx_phi);
1 0 0 0 0.984807753012208 0.17364817766693 0 -0.17364817766693 0.984807753012208
Now we can rotate a point by and around the x, y and z axes of whatever coordinate reference system we are using.
It is also useful to notate TO and FROM which coordiante reference system we are using.
For example, if we are rotation from coordinate reference system A to coordinate reference system B by rotating along the Z axis, we would notate this rotation as

to BODY from VEHICLE transform

The order of rotation matters significantly so in aerospace, there is a standard order of rotation.
In transforming from the Vehicle reference system to the Body reference system, we use the following order
This reads as the Transform to Body from Vehicle is equal to (from right to left)
  1. Rotation around z to Vehicle-1 from Vehicle
  2. Rotation around y to Vehicle-2 from Vehicle-1
  3. Rotation around z to Body from Vehicle-2
Let's do an example. Let's start by creating our matrix
T_2BFV = Rx_phi*Ry_the*Rz_psi;
disp(T_2BFV);
0.969846310392954 0.171010071662834 -0.17364817766693 -0.141314484355892 0.975082443643152 0.171010071662834 0.198565734023778 -0.141314484355892 0.969846310392954
Next, let's rotation a point by the transform. We'll start by creating a point.
p = [1;0;0];
We set the point as a 3 row, 1 column matrix - each row corresponding to x, y and z coordinates. This point is along the x-axis.
Lets do a rotation for p - we'll rotate is by around the z-axis, then around the y-axis, then around the x-axis.
new_point = T_2BFV*p;
disp(new_point);
0.969846310392954 -0.141314484355892 0.198565734023778
Let's see how we did
While this is an accurate calculation, this probably doesn't mean anything to you. We can rationalize this in our heads a bit though. Imagine a coordinate reference system where the x-axis goes out the right, the y-axis goes up, and the z-axis comes out of the screen. If we wanted to rotate around the z-axis in order to get to the x-axis, our starting point would have to lie below the x-axis. Remember, this is the rotation from the EARTH to BODY. If we wanted to Earth from Body (ie, if the point was at the x-axis and we wanted to see where it ends up after we rotated), then we would have to take the transpose of the transformation matrix.
figure;
quiver3(0,0,0,1,0,0,'Color','#7E2F8E'); hold on;
quiver3(0,0,0,0,1,0,'Color','#7E2F8E');
quiver3(0,0,0,0,0,1,'Color','#7E2F8E');
plot3(1,0,0,'ro','MarkerSize',10);
plot3(0.97,-.141,.199,'bo','MarkerSize',10);
axis([-1.2 1.2 -1.2 1.2]);
title('Points');
xlabel('X axis');
ylabel('Y axis');
zlabel('Z axis');
grid on;
 
view([-0.600 90.000])
And if we rotate with a positive angle around the y-axis (pitching up), then the new point should also have a positive Z values. Go ahead and rotate to see.

Example 1

Let's do an example.
An aircraft is moving north at 100 ft per second, with its nose 10 degrees above the horizon. What are the body axis velocities?
First, we'll write the vector for the velocity of the aircraft in Earth reference frame
V_E = [100;0;0];
Next, write the vector for the body axis velocities.
Next we'll calculate the rotation matrix. The angles represent the orientation of the aircraft.
psi = 0;
theta = 10;
phi = 0;
Next, we'll calculate the Transform to Body from World matrix.
psi_in_radians = deg2rad(psi);
Rz_psi = [ cos(psi_in_radians), sin(psi_in_radians), 0;
-sin(psi_in_radians), cos(psi_in_radians), 0;
0, 0, 1];
theta_in_radians = deg2rad(theta);
Ry_the = [cos(theta_in_radians), 0, -sin(theta_in_radians);
0, 1, 0;
sin(theta_in_radians), 0, cos(theta_in_radians)];
phi_in_radians = deg2rad(phi);
Rx_phi = [ 1, 0, 0;
0, cos(phi_in_radians), sin(phi_in_radians);
0, -sin(phi_in_radians), cos(phi_in_radians)];
T_2BFE = Rx_phi*Ry_the*Rz_psi;
disp(T_2BFE);
0.984807753012208 0 -0.17364817766693 0 1 0 0.17364817766693 0 0.984807753012208
Now, we can solve for the aircraft's velocities
V_B = T_2BFE*V_E;
disp(V_B);
98.4807753012208 0 17.364817766693
Let's make some sense of this
figure;
quiver(0,0,100,0,'Color','r'); hold on;
North_vector = [V_B(1)*cos(deg2rad(10));-V_B(1)*sin(deg2rad(10))];
Down_vector = [V_B(3)*sin(deg2rad(10));V_B(3)*cos(deg2rad(10))];
quiver([0,0],[0,0],[North_vector(1),Down_vector(1)],[North_vector(2),Down_vector(2)],'Color','b');
 
 
axis([-20 120 -30 30]);
title('Side Profile');
xlabel('North');
ylabel('Down');
set(gca, 'YDir','reverse')
grid on;
axis equal;
legend('Motion in the World Frame','Body Frame velocities')
When we draw out the vectors, it starts to make some sense. The aircraft is pitched up and moving forward at 98.4 ft per second. But in order to prevent it from gaining altitude, it's also moving in the downward direction at 17.36 ft per second.
What happens if the aircraft was pitched up at instead?
psi = 0;
theta = 18;
phi = 0;
Next, we'll calculate the Transform to Body from World matrix.
psi_in_radians = deg2rad(psi);
Rz_psi = [ cos(psi_in_radians), sin(psi_in_radians), 0;
-sin(psi_in_radians), cos(psi_in_radians), 0;
0, 0, 1];
theta_in_radians = deg2rad(theta);
Ry_the = [cos(theta_in_radians), 0, -sin(theta_in_radians);
0, 1, 0;
sin(theta_in_radians), 0, cos(theta_in_radians)];
phi_in_radians = deg2rad(phi);
Rx_phi = [ 1, 0, 0;
0, cos(phi_in_radians), sin(phi_in_radians);
0, -sin(phi_in_radians), cos(phi_in_radians)];
T_2BFE = Rx_phi*Ry_the*Rz_psi;
disp(T_2BFE);
0.951056516295154 0 -0.309016994374947 0 1 0 0.309016994374947 0 0.951056516295154
Now, we can solve for the aircraft's velocities
V_B = T_2BFE*V_E;
disp(V_B);
95.1056516295154 0 30.9016994374947
figure;
quiver(0,0,100,0,'Color','r'); hold on;
North_vector = [V_B(1)*cos(deg2rad(18));-V_B(1)*sin(deg2rad(18))];
Down_vector = [V_B(3)*sin(deg2rad(18));V_B(3)*cos(deg2rad(18))];
quiver([0,0],[0,0],[North_vector(1),Down_vector(1)],[North_vector(2),Down_vector(2)],'Color','b');
 
 
axis([-20 120 -30 30]);
title('Side Profile');
xlabel('North');
ylabel('Down');
set(gca, 'YDir','reverse')
grid on;
axis equal;
legend('Motion in the World Frame','Body Frame velocities')
As you would expect, if we were pitched up more, we would have a larger downward body velocity to compensate to ensure that we're travelling strictly northward with no change in altitude in the Earth frame.

Example 2

An aircraft is moving forward at 100 ft per second, with a heading of 230 degrees and its nose 5 degrees below the horizon.
What are the velocities of the aircraft in the Earth frame?
Set up what we know:
V_B = [100;0;0];
psi = 230;
theta = -5;
phi = 0;
Calculate the Transform to BODY from EARTH
psi_in_radians = deg2rad(psi);
Rz_psi = [ cos(psi_in_radians), sin(psi_in_radians), 0;
-sin(psi_in_radians), cos(psi_in_radians), 0;
0, 0, 1];
theta_in_radians = deg2rad(theta);
Ry_the = [cos(theta_in_radians), 0, -sin(theta_in_radians);
0, 1, 0;
sin(theta_in_radians), 0, cos(theta_in_radians)];
phi_in_radians = deg2rad(phi);
Rx_phi = [ 1, 0, 0;
0, cos(phi_in_radians), sin(phi_in_radians);
0, -sin(phi_in_radians), cos(phi_in_radians)];
T_2BFE = Rx_phi*Ry_the*Rz_psi;
Next, calculate the Transform to EARTH from BODY
T_2EFB = T_2BFE';
Now we can solve for the aircraft's motion over Earth
V_E = T_2EFB*V_B;
disp(V_E);
-64.0341608768797 -76.312941273777 8.71557427476582
Does this make sense?
figure;
quiver(0,0,100,0,'k--'); hold on;
quiver(0,0,0,100,'k--');
quiver(0,0,V_E(2),V_E(1),'r');
xlabel('East');
ylabel('North');
grid on;
axis equal;
 
figure;
quiver(0,0,100,0,'k--'); hold on;
quiver(0,0,0,100,'k--');
quiver(0,0,99.62,V_E(3),'r');
xlabel('Forward');
ylabel('Down');
set(gca, 'YDir','reverse')
grid on;
axis equal;