Home » Space Flight/Orbital Mechanics Excercise (Page 2)

# Category Archives: Space Flight/Orbital Mechanics Excercise

## 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

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

# Contents

## Free-Space Path Loss vs. Cable Loss

Click for more…! Consider a space link with 100 million kilometer distance and a transmit frequency of 2 GHz (S-Band). 1. Calculate the Free-Space Path Loss [dB].

```clear all;
clc;
close all;
F = 2;                                %transmit frequency (S-Band)[GHz]
d = 1e8;                              %[km]
L = 92.4 + 20*log10(F) + 20*log10(d); %[dB]
fprintf('Free-Space Path Loss L = %4.2f [dB] \n',L);
% Assume a typical coaxial cable with a loss of 0.3[dB/m] at S-Band.
% How long may this cable be for a total loss equal to the Free-Space
% Path Loss from above [m]?
c_l = 0.3;          %[dB/m]
l_cabel = L/c_l;    %[m]
fprintf('Equivalent Cabel Length = %4.2f [m] \n\n',l_cabel);```
```Free-Space Path Loss L = 258.42 [dB]
Equivalent Cabel Length = 861.40 [m]```

## Propagation Delays

1.Calculate the two-way propagation delays[min] between Earth and spacecrafts at diferent planets (from Mercury to Saturn; consider the following average distances between Sun-Planet: 0.3871AU for Mercury, 0.723AU for Venus, 1.524AU for Mars, 5.203AU for Jupiter, 9.582AU for Saturn). Assume conjunction Sun-Planet-Earth (Mercury, Venus), or Sun-Earth-Planet

```%(Mars,Jupiter, Saturn) for an easy calculation of the minimum distances
%(Note that this is not the worst case in terms of maximum distances to be
% considered for the actual link design).
v_light   = 300000;% Speed of light in vacuum [km/s]
R_mercury = 0.3871;    %[AU]
R_venus   = 0.7230;    %[AU]
R_earth   = 1.0000;    %[AU]
R_mars    = 1.5240;    %[AU]
R_jupiter = 5.2030;    %[AU]
R_saturn  = 9.5820;    %[AU]
AU        = 149597871; %[km]
mu_sun = 132712440018;          %(km^3s?2)```

Earth – Mercury

```p_d = 2*(R_earth - R_mercury)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Mercury = %4.2f [min]\n',p_d);
figure(1);
T_mercury  = 2*pi*((R_mercury*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:10*T_mercury;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_mercury = rem(time,T_mercury)/(T_mercury)*360;
phase = fi_earth - fi_mercury;
dist = (R_earth^2 + R_mercury^2 - 2*R_earth*R_mercury*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Mercury Comunication');
ylabel('Distance Earth - Mercury [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Mercury = 10.19 [min]` Earth – Venus

```p_d = 2*(R_earth - R_venus)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Venus   = %4.2f  [min]\n',p_d);
figure(2);
T_venus  = 2*pi*((R_venus*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:7*T_venus;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_venus = rem(time,T_venus)/(T_venus)*360;
phase = fi_earth - fi_venus;
dist = (R_earth^2 + R_venus^2 - 2*R_earth*R_venus*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Venus Comunication');
ylabel('Distance Earth - Venus [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Venus   = 4.60  [min]` Earth – Mars

```p_d = 2*(R_mars - R_earth)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Mars    = %4.2f  [min]\n',p_d);
figure(3);
T_mars  = 2*pi*((R_mars*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:6*T_mars;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_mars  = rem(time,T_mars)/(T_mars)*360;
phase = fi_earth - fi_mars;
dist = (R_earth^2 + R_mars^2 - 2*R_earth*R_mars*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Mars Comunication');
ylabel('Distance Earth - Mars [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Mars    = 8.71  [min]` Earth – Jupiter

```p_d = 2*(R_jupiter- R_earth)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Jupiter    = %4.2f  [min]\n',p_d);
figure(4);
T_jupiter  = 2*pi*((R_mars*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:3*T_mars;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_jupiter = rem(time,T_jupiter)/(T_jupiter)*360;
phase = fi_earth - fi_jupiter;
dist = (R_earth^2 + R_jupiter^2 - 2*R_earth*R_jupiter*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title('Earth - Jupiter Comunication');
ylabel('Distance Earth - Jupiter [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Jupiter    = 69.86  [min]` Earth – Saturn

```p_d = 2*(R_saturn- R_earth)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Saturn    = %4.2f  [min]\n',p_d);
figure(5);
T_saturn  = 2*pi*((R_saturn*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:3*T_mars;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_saturn = rem(time,T_saturn)/(T_saturn)*360;
phase = fi_earth - fi_saturn;
dist = (R_earth^2 + R_saturn^2 - 2*R_earth*R_saturn*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Saturn Comunication');
ylabel('Distance Earth - Saturn [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Saturn    = 142.65  [min]` ```%2. For the various cases calculate the Free-Space Path Loss[dB], assuming
% an RF frequency of 6 GHz.
F = 6; %transmit frequency [GHz]
% Earth - Mercury
L = 92.4 + 20*log10(F) + 20*log10((R_earth - R_mercury)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Mercury L = %4.2f [dB] \n',L);
% Earth - Venus
L = 92.4 + 20*log10(F) + 20*log10((R_earth - R_venus)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Venus   L = %4.2f  [dB] \n',L);
% Earth - Mars
L = 92.4 + 20*log10(F) + 20*log10((R_mars - R_earth)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Mars    L = %4.2f [dB] \n',L);
% Earth - Jupiter
L = 92.4 + 20*log10(F) + 20*log10((R_jupiter - R_earth)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Jupiter L = %4.2f [dB] \n',L);
% Earth - Saturn
L = 92.4 + 20*log10(F) + 20*log10((R_saturn - R_earth)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Saturn  L = %4.2f [dB] \n\n',L);```
```Free-Space Path Loss, Earth - Mercury L = 267.21 [dB]
Free-Space Path Loss, Earth - Venus   L = 260.31  [dB]
Free-Space Path Loss, Earth - Mars    L = 265.85 [dB]
Free-Space Path Loss, Earth - Jupiter L = 283.93 [dB]
Free-Space Path Loss, Earth - Saturn  L = 290.13 [dB]```

## Satellite Design

Your satellite designer wants to reduce the satellite transmitter output power from 50 W to 25 W to save weight. How much is this reduction expressed in dB scale?

```P1 = 50;                % [W]
P1_db = 10*log10(P1);   % 17 dBW
P2 = 25;                % [W]
P2_db = 10*log10(P2);   % 14 dBW
red = (P2_db - P1_db);  % - 3  dBW

fprintf('Reduction expressed in DB scale %4.0f\n',red);

% If you want to maintain the satellite-to-ground
% station data link at the same data rate, you could achieve this by
% modifying the antenna on ground: By which factor do you have to increase
% the diameter of a your dish antenna then?
% Pt1*A1r = Pt2*A2r  A2r/A1r = Pt1/Pt2 = 2
% Pt - Transmited Power ;
% Ar - Effective Area of receiving antena;

R_factor = sqrt(P1/P2);
fprintf('Diameter increase your dish antenna %4.2f \n',R_factor);```
```Reduction expressed in DB scale   -3
Diameter increase your dish antenna 1.41```

Published with MATLAB® 7.10