<< Numerical integration of orbital motion Cookbook Simple orbital simulation >>

CelestLab >> - Introduction - > Cookbook > Rotations

Rotations

Use of quaternions and frame transformation matrices

Example: various ways to compute the "qsw" local orbital frame

The following examples show many ways to compute the "qsw" local orbital frame (see Local frames for more details).

The objective is not so much the computation of the local frame as the illustration of the use of functions related to quaternions and frame transformation matrices. The notations used are:

mat_qsw: frame transformation matrix from the inertial frame to the "qsw" local orbital frame. Multiplying by mat_qsw changes coordinates from the inertial frame to the "qsw" local frame.

q_qsw: quaternion of the rotation that transforms the axes of the inertial frame into the axes of the "qsw" local orbital frame. Applying the rotation q_qsw to a vector defined in the inertial frame gives the components of the rotated vector in the same frame. Applying the rotation q_qsw is also equivalent to changing components from the "qsw" local frame to the inertial frame. Note that q_qsw and -q_qsw correspond to the same rotation.

Note: First execute the code snippet below to initialize variables.

// =====================================================
// Initializations of variables
// =====================================================
// position and velocity in inertial frame (m and m/s)
pos = [6000.e3; 3000.e3; 3000.e3]; 
vel = [-4.e3; 3.e3; 6.e3];

// Keplerian orbital elements (using default value for "mu")
kep = CL_oe_car2kep(pos, vel);

inc = kep(3); // inclination
gom = kep(5); // RAAN
pom = kep(4); // arg. of perigee
anv = CL_kp_M2v(kep(2), kep(6)); // true anomaly

1) Local frame from position and velocity (frame transformation matrices)

// =====================================================
// Direct computation (frame transformation matrices)
// =====================================================

// The frame transformation matrix to "qsw" is: 
mat_qsw = CL_fr_qswMat(pos, vel) 

// Another way is to use "CL_fr_locOrbMat": 
mat_qsw = CL_fr_locOrbMat(pos, vel, "qsw") 

// Or "CL_rot_defFrameVec":
mat_qsw = CL_rot_defFrameVec(pos, vel, 1, 2) 
 
 
// -----------------------------------------------------
// Multiplying by mat_qsw yields coordinates relative to "qsw": 

// u: some vector, coordinates in inertial frame 
u = CL_unitVector(CL_cross(pos, vel));   

// u_in_qsw: same vector, coordinates in "qsw"
u_in_qsw = mat_qsw * u 

// The same is obtained by:
u_in_qsw = CL_fr_inertial2qsw(pos, vel, u) 

// Or even by: 
u_in_qsw = CL_rot_pvConvert(u, [], mat_qsw, [])

2) Local frame from orbital elements (frame transformation matrices)

// =====================================================
// LOF from orbital elements (frame transformation matrices)
// => 3 successive rotations (each defined in new axes)
// =====================================================

// The "qsw" local frame can be obtained from the inertial 
// frame by 3 successive rotations:
// - 1st rotation: around Z axis of angle "gom"      => frame1
// - 2nd rotation: around X axis of angle "inc"      => frame2
// - 3rd rotation: around Z axis of angle "pom+anv"  => qsw

// The respective frame transformation matrices are:
mat1 = CL_rot_angles2matrix(3, gom); 
mat2 = CL_rot_angles2matrix(1, inc); 
mat3 = CL_rot_angles2matrix(3, pom+anv); 

// Then:
// (X/qsw = mat3 * X/frame2 
//        = mat3 * mat2 * X/frame1 
//        = mat3 * mat2 * mat1 * X/inertial)
mat_qsw = mat3 * mat2 * mat1 

// Which is the same as: 
mat_qsw = CL_rot_angles2matrix([3,1,3], [gom;inc;pom+anv]) 

// Or even as: 
mat_qsw = CL_rot_compose(mat1,[0;0;0],1,mat2,[0;0;0],1,mat3,[0;0;0],1)

3) Local frame from orbital elements (quaternions)

// =====================================================
// LOF from orbital elements (quaternions)
// => 3 successive rotations (each defined in new axes)
// =====================================================

// The 3 rotations in 2) can be described by 3 quaternions, 
// respectively q1, q2, q3:  

q1 = CL_rot_angles2quat(3, gom); 
q2 = CL_rot_angles2quat(1, inc); 
q3 = CL_rot_angles2quat(3, pom+anv); 

// Or: 
q1 = CL_rot_axAng2quat([0;0;1], gom); 
q2 = CL_rot_axAng2quat([1;0;0], inc); 
q3 = CL_rot_axAng2quat([0;0;1], pom+anv); 
 
 
// The quaternion of the rotation that transforms the inertial 
// frame into "qsw" is then:
q_qsw = q1 * q2 * q3 
 
 
// The quaternion could also be obtained by: 
q_qsw = CL_rot_angles2quat([3,1,3], [gom;inc;pom+anv]) 

// Or by: 
q_qsw = CL_rot_defFrameVec(pos, vel, 1, 2, res="q") 

// Or:
q_qsw = CL_rot_defRotVec([1;0;0], [0;1;0], pos, vel) 
 
 
// -----------------------------------------------------
// The X axis of the "qsw" frame is the image of [1;0;0] by the rotation: 
x_axis_of_qsw = CL_rot_rotVect(q_qsw, [1;0;0]) // coord in inertial frame

// The Euler angles can be derived from the quaternion by: 
// (only one solution considered here!)
angles = CL_rot_quat2angles(q_qsw, [3,1,3]) 
 
 
// -----------------------------------------------------
// The same could be done by using the corresponding 
// transformation matrix: 
mat_qsw = CL_rot_quat2matrix(q_qsw) 

// mat_qsw': frame transformation matrix from "qsw" to inertial frame 
x_axis_of_qsw = mat_qsw' * [1;0;0] 

angles = CL_rot_matrix2angles(mat_qsw, [3,1,3])

4) Local frame from orbital elements (quaternions in inertial frame)

// =====================================================
// LOF from orbital elements (quaternions)
// => 3 successive rotations, all defined in the inertial frame
// =====================================================

// The 3 rotations in 2) and 3) can be defined by quaternions relative 
// to the inertial frame: 

q1 = CL_rot_axAng2quat([0;0;1], gom); 
x1 = CL_rot_rotVect(q1, [1;0;0]); // new x-axis in inertial frame
z1 = CL_rot_rotVect(q1, [0;0;1]); // new z-axis in inertial frame

q2 = CL_rot_axAng2quat(x1, inc); 
z2 = CL_rot_rotVect(q2, z1); // new z-axis in inertial frame

q3 = CL_rot_axAng2quat(z2, pom+anv); 

// The combined rotation quaternion is obtained by multiplying
// the quaternions in reverse order: 
q_qsw = q3 * q2 * q1

5) Local frame from right ascension, declination, "bearing" (quaternions)

// =====================================================
// LOF from orbital elements (quaternions)
// => 3 successive rotations (right ascension, declination, "bearing")
// =====================================================

// We first compute the rotation angles in "right ascension", "declination" 
// and "bearing" (see function for more details):  
[rra, decl, bear] = CL_gm_orbSphTriangle(inc, "aol", pom+anv, ..
                    output=["rra", "decl", "bear"]);

// The "inertial to qsw" rotation quaternion can be computed by 3 successive 
// rotations around Z, then Y, then X: 
q_qsw = CL_rot_angles2quat([3,2,1], [gom+rra;-decl;%pi/2-bear])

6) Local frame in 2 successive rotations (quaternions)

// =====================================================
// LOF (quaternions) 
// => Determination of 2 successive rotations 
// =====================================================

// The rotation that transforms the inertial frame into "qsw" can be considered as 
// the combination of 2 rotations R1 followed by R2: 
// R1: rotation such that X=[1;0;0] becomes "pos" 
// (R1 is itself a combination of 2 rotations: around Z then Y)
// R2: rotation around "new X" = "pos" so that "new Y" becomes "vel"

q1 = CL_rot_defRot2Ax([1;0;0], pos, 3, 2); 
x1 = CL_rot_rotVect(q1, [1;0;0]); // new X axis after R1 (components in inertial frame)
y1 = CL_rot_rotVect(q1, [0;1;0]); // new Y axis after R1 (components in inertial frame)

q2 = CL_rot_defRot1Ax(y1, x1, vel, 0); 

// The quaternion of the combined rotation is:
q_qsw = q2 * q1
 
 
// -----------------------------------------------------
// The second rotation could be defined in the new axes (obtained after R1):

q1 = CL_rot_defRot2Ax([1;0;0], pos, 3, 2); 

// Compute velocity components in the new axes: 
vel1 = CL_rot_rotVect(q1', vel);  
// Same as: 
vel1 =  CL_rot_quat2matrix(q1) * vel; 

q2 = CL_rot_defRot1Ax([0;1;0], [1;0;0], vel1, 0); 

// The quaternion of the combined rotation is:
q_qsw = q1 * q2
 
 
// -----------------------------------------------------
// The rotations R1 and R2 are not unique. 
// Below, R1 transforms the X axis into "pos", in a different way:
q1 = CL_rot_defFrameVec(pos, [0;0;1], 1, 2, res="q");  
q2 = CL_rot_defRot1Ax([0;1;0], [1;0;0], CL_rot_rotVect(q1', vel), 0); 

// The quaternion of the combined rotation is:
q_qsw = q1 * q2


Report an issue
<< Numerical integration of orbital motion Cookbook Simple orbital simulation >>