Home » 2012 » December (Page 2)

# Monthly Archives: December 2012

## Inertial components of the angular acceleration

Example 10.4, Orbital Mechanics for Engineering Students, 2nd Edition.

The inertial components of the angular momentum of a torque-free rigid body are

```Hg = [320; -375; 450 ];       %[kg*m^2/s] IJK
% the Euler angles [deg] are
fi      =  20;
theta   =  50;
psi     =  75;
% The inertia tensor in the body-fixed principal frame is
Ig = [1000, 0,    0;
0, 2000, 0;
0, 0,    3000]; %[kg*m^2]```

Obtain the inertial components of the (absolute) angular acceleration Matrix of the transformation from body-fixed frame to inertial frame

```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)
]```
```QxX =

0.0309   -0.9646    0.2620
0.6720   -0.1740   -0.7198
0.7399    0.1983    0.6428```

Matrix of the transformation from inertial frame to body-fixed frame

`QXx = QxX'`
```QXx =

0.0309    0.6720    0.7399
-0.9646   -0.1740    0.1983
0.2620   -0.7198    0.6428```

Obtain the components of HG in the body frame

`Hgx = QXx*Hg            % [kg*m^2/s]`
```Hgx =

90.8616
-154.1810
643.0376```

The components of angular velocity in the body frame

```Ig_inv = inv(Ig);
```wx =

0.0909
-0.0771
0.2143```

From Euler ’s equations of motion we calculate angular acceleration in the body frame.

`alfa_x = - Ig_inv*cross(wx,Ig*wx)               % [rad/s^2]`
```alfa_x =

0.0165
0.0195
0.0023```

Angular acceleration in the inertial frame

`alfa_X = QxX*alfa_x                             % [rad/s^2] IJK`
```alfa_X =

-0.0177
0.0060
0.0176```

Published with MATLAB® 7.10

## Spin axis, precession

Example 10.2, Orbital Mechanics for Engineering Students, 2nd Edition.

A cylindrical shell is rotating in torque-free motion about its longitudinal axis.

```r = 1;      % [m]
l = 3;      % [m]
m = 100;    % [kg]
theta = 20; % [deg ]nutation angle
% How long does it take the cylinder to precess through 180 ° if the
% spin rate is 2*pi/radians per minute
C = m*r^2                       % [kg*m^2]
A = 1/2*m*r^2 + 1/12*m*l^2```
```C =

100

A =

125.0000```

A > C, Direct (prograde) precession

`wp = C/(A - C)*ws/cosd(theta)   % [rad/min]`
```wp =

26.7457```

The time for the spin axis to precess through an angle of 180 ° is

`t = pi/wp*60                    % [sec]`
```t =

7.0477```

Published with MATLAB® 7.10

## Moments of inertia about the center of mass of the system of six point masses

Find the moments of inertia about the center of mass of the system of six point masses.

```clc;
clear all;
M = [10,10,8,8,12,12]; % [kg]
X = [1 -1 4 -2 3 -3];  % [m]
Y = [1 -1 -4 2 -3 3];  % [m]
Z = [1 -1 4 -2 -3 3];  % [m]
mt = sum(M);           % The total mass of this system```

Three components of the position vector of the center of mass are

```Xcg =(1/mt)*sum(M.*X);
Ycg =(1/mt)*sum(M.*Y);
Zcg =(1/mt)*sum(M.*Z);
V_cg = [Xcg,Ycg,Zcg]   % [m]```
```V_cg =

0.2667   -0.2667    0.2667```
```Ig = [0.0,0.0,0.0;
0.0,0.0,0.0;
0.0,0.0,0.0];```

The total moment of inertia is the sum of moments of inertia for all point masses in the system

```for i =1:6
x = (X(i) - Xcg);
y = (Y(i) - Ycg);
z = (Z(i) - Zcg);
m = M(i);
Ig = Ig + [  m*(y^2 + z^2),    -m*x*y,         -m*x*z;
-m*y*x,             m*(x^2 + z^2), -m*y*z;
-m*x*z,            -m*y*z,          m*(y^2 + x^2)]; % [kg*m^2]
end
Ig```
```Ig =

783.4667  351.7333   40.2667
351.7333  783.4667  -80.2667
40.2667  -80.2667  783.4667```

Published with MATLAB® 7.10

## The Quaternion, Euler principal rotation angle & Euler axis of rotation.

Contents

```% For the yaw-pitch-roll sequence calculate
clc;
clear all;
yaw   = 50;  % [deg]
pitch = 90;  % [deg]
roll  = 120; % [deg]```

## The quaternion

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

K =1/3*[ QXx(1,1)-QXx(2,2)-QXx(3,3), QXx(2,1)+QXx(1,2), QXx(3,1)+QXx(1,3),...
QXx(2,3)-QXx(3,2);
QXx(2,1)+QXx(1,2), -QXx(1,1)+QXx(2,2)-QXx(3,3), QXx(3,2)+QXx(2,3),...
QXx(3,1) - QXx(1,3);
QXx(3,1)+QXx(1,3), QXx(3,2)+QXx(2,3), -QXx(1,1)-QXx(2,2)-QXx(3,3),...
QXx(1,2)-QXx(2,1);
QXx(2,3) - QXx(3,2), QXx(3,1)- QXx(1,3), QXx(1,2)- QXx(2,1),...
QXx(1,1)+QXx(2,2)+QXx(3,3)
]
% The eigenvalues and eigenvectors of this matrix are
[V,D] = eig(K);
% Finding the largest eigenvalue
eigv = [D(1,1),D(2,2),D(3,3),D(4,4)]
[max_eig,i] = max(eigv);
% The quaternion is the eigenvector coresponding to the maximum eigenvalue
q = V(:,i)
% Checking the norm of quaternion.Should be equal to 1
q_norm = norm(q)```
```R1_roll =

1.0000         0         0
0   -0.5000    0.8660
0   -0.8660   -0.5000

R2_pitch =

0     0    -1
0     1     0
1     0     0

R3_yaw =

0.6428    0.7660         0
-0.7660    0.6428         0
0         0    1.0000

QXx =

0         0   -1.0000
0.9397    0.3420         0
0.3420   -0.9397         0

K =

-0.1140    0.3132   -0.2193    0.3132
0.3132    0.1140   -0.3132    0.4473
-0.2193   -0.3132   -0.1140   -0.3132
0.3132    0.4473   -0.3132    0.1140

eigv =

-0.3333   -0.3333   -0.3333    1.0000

q =

-0.4056
-0.5792
0.4056
-0.5792

q_norm =

1.0000```

Lets now calculate quaternion using another approach.This procedure is easy but fails when q4 = 0.

```q4 = 0.5*sqrt(1 + QXx(1,1) + QXx(2,2) + QXx(3,3));
q1 = (QXx(2,3) - QXx(3,2))/(4*q4);
q2 = (QXx(3,1) - QXx(1,3))/(4*q4);
q3 = (QXx(1,2) - QXx(2,1))/(4*q4);

q = [q1;q2;q3;q4]
q_norm = norm(q)```
```q =

0.4056
0.5792
-0.4056
0.5792

q_norm =

1```

## Euler principal rotation angle and Euler axis of rotation.

Euler principal rotation angle [deg]

```theta = 2*acos(q(4))*180/pi
% The unit vector along the Euler axis around which the inertial reference
% frame is rotated into the body fixed frame

u = [q1 q2 q3]/(sind(theta/2))```
```theta =

109.2075

u =

0.4975    0.7106   -0.4975```

Published with MATLAB® 7.10

Contents

Consider the following example of a spacecraft telemetry link for which we wish to find the maximum range. Suppose we have a 1 W telemetry transmitter onboard a spacecraft and that the data rate requires a channel capacity corresponding to a signal-to-noise ratio of at least unity in a 1 Hz bandwidth. This link uses a frequency of 3 GHz (10 cm wavelength). The transmitting antenna on the spacecraft is a 2 m diameter dish. The ground station antenna is a 10 m diameter dish. Assume that both these dish antennas have an eective area equal to 60% of their physical apertures. Assume also that there is no pointing error, that is, the antennas always point directly at each other and that the system temperature of the ground station receiver is 25 K. (The system temperature is the sum of the equivalent receiver noise temperature, the antenna noise temperature and the sky noise temperature.) Boltzmann constant k = 1.38*10-23 W/(Hz*K)

```clc;
clear all;
k    = 1.38e-23;
c    = 3e8;         % Speed of Light [m/s]
F    = 3;           % Operating Frequency[GHz]
lm   = c/(F*10^9);  % Wavelength at the operating frequency  [m/s]
E_gs = 10;          % Ground Station antenn diameter [m]
E_sc = 2 ;          % Spacecraft antenn diameter [m]
eff  = 0.6;         % 60[%] Aperture Efficiency```

1. What is the equivalent input Noise Power of the receiver[dBm]?

```%N = 10*log10()
Te = 25;  % System (receiver) noise temperature in [K];
B = 1;    % System noise bandwidth in              [Hz]
N = 10*log10(k*Te*B); % or
N = -228.6 + 10*log10(Te) + 10*log10(B);         % [dBW]
N = N + 30;                                      % [dBm]
fprintf('Equivalent input Noise Power of the receiver %4.2f [dBm] \n\n',N);```
`Equivalent input Noise Power of the receiver -184.62 [dBm]`

2. What is the Effective Area of the receiving antenna [m2]?

```A_gs = pi*E_gs^2/4;        % Area of antenna aperture
G_gs = 20*log10(E_gs) + 20*log10(F) + 20*log10(10*pi/3*sqrt(eff));
% or
G_gs = 10*log10(4*pi*A_gs*eff/lm^2); % [dBi]
Aeff_gs = 10^(G_gs/10)*lm^2/(4*pi);          % [m^2]
fprintf('Gain of the receiving antenna = %4.2f [dBi] \n',G_gs);
fprintf('Effective Area of the receiving antenna = %4.2f [m^2]\n\n',Aeff_gs);```
```Gain of the receiving antenna = 47.72 [dBi]
Effective Area of the receiving antenna = 47.12 [m^2]```

3. What is the gain of the transmitting antenna [dBi]?

```A_sc = pi*E_sc^2/4;                  % Area of SC antenna aperture
G_sc = 10*log10(4*pi*A_sc*eff/lm^2); % [dBi]
Aeff_sc = 10^(G_sc/10)*lm^2/(4*pi);          % [m^2]
fprintf('Gain of the transmitting antenna = %4.2f [dBi] \n',G_sc);
fprintf('Effective Area of the transmitting antenna = %4.2f [m^2]\n',Aeff_sc);```
```Gain of the transmitting antenna = 33.75 [dBi]
Effective Area of the transmitting antenna = 1.88 [m^2]```

4. What is the maximum range, R [km] for the spacecraft to maintain the required signal-to-noise ratio?

```L = -(N-30); % dBw maximum path Loss
Power = G_sc + G_gs; %dBi
R = 10^((Power + L - 92.4 - 20*log10(F))/20);            %[km]
fprintf('Maximum range, R = %4.2d [km]\n\n',R);```
`Maximum range, R = 5.10e+009 [km]`

## Coherent-Frequency Turnaround

Consider a spacecraft transponder which operates in frequency-coherent mode. The transmitted uplink frequency is 2 GHz (S-Band), and the downlink frequency is in X-Band. The turn-around ratio r inside the spacecraft transponder is 880/221 (downlink-freq/uplink-freq). The assumed velocity of the spacecraft relative to the groundstation is 10km/sec. Assume further that the downlink telemetry (TM) rate of 10 kbit/sec is modulated onto a sub-carrier with frequency 150 kHz.

1. Calculate the Doppler shift[kHz] on the uplink and downlink carriers, which are received by the transponder and the groundstation,respectively.

```f_uplink   =  2;                            %[GHz]
f_sub = 150;                                %[kHz]
v_sat = 10000;                              %[m/s]
dr = 10;                                    %[kbit/sec]
```Doppler shift Uplink Carrier  =  66.67 [kHz]
Doppler shift Downlink Carrier = 265.46 [kHz]```

2. Calculate the Doppler shifts on the received downlink sub-carrier[Hz], and the TM data rate[bit/s]

```df_sub= v_sat/c*f_sub*10^3;                %[Hz]
fprintf('Doppler shift Sub-Carrier = %4.2f [Hz] \n',df_sub);
% Assuming that the spacecraft was opserved during the recession
% Shifted frequency
f_shif = f_sub - df_sub;
%TM data shift rate[bit/s]
dr_sub= v_sat/c* dr*1000;
fprintf('Downlink telemetry (TM) rate shift = %4.2f [bit/sec] \n',dr_sub);```
```Doppler shift Sub-Carrier = 5.00 [Hz]
Downlink telemetry (TM) rate shift = 0.33 [bit/sec]```

Published with MATLAB® 7.10