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 fprintf('wx = %4.4f[rad/s] wy = %4.4f[rad/s] wz = %4.4f[rad/s] \n\n',... 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; fprintf('wp = %4.4f[rad/s] wn = %4.4f[rad/s] ws = %4.4f[rad/s] \n\n',... 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] wp = -0.4464[rad/s] wn = 2.7704[rad/s] ws = -3.1404[rad/s]
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; fprintf('wx = %6.5f [rad/s] wy = %6.5f [rad/s] wz = %6.5f [rad/s]\n',... wx,wy,wz); fprintf('dwx = %6.5f [rad/s^2] dwy = %6.5d [rad/s^2] dwz = %6.5f [rad/s^2]\n\n',... dwx,dwy,dwz);
wx = -0.09131 [rad/s] wy = 0.09855 [rad/s] wz = 1.19405 [rad/s] dwx = 0.06341 [rad/s^2] dwy = 5.81984e-005 [rad/s^2] dwz = -0.08187 [rad/s^2]
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